C++ TFile::Get方法代码示例

void PP_MC_normResponseMatrix(int radius = 3){

  gStyle->SetOptStat(0);
  
  bool printDebug = true;
  
  TFile * fin = TFile::Open(Form("pp_MC_new_ak%d_20_eta_20.root",radius));
  
  TH2F * mPP_Matrix;
  TH2F * mPP_Response;
  TH2F * mPP_ResponseNorm;
  TH2F * mPP_ResponseNorm_recoTrunc;

  // TH2F * mPbPb_Matrix_recoTrunc[nbins_cent];
  // TH2F * mPbPb_Response_recoTrunc[nbins_cent];
  // TH2F * mPbPb_ResponseNorm_recoTrunc[nbins_cent];

  TFile fout(Form("PP_kineticEfficiency_correctionFactors_R%d.root",radius),"RECREATE");
  fout.cd();
  
  TH1F * hGen_Projection;
  TH1F * hGen_Projection_recoTrunc;
  TH1F * hKineticEfficiency;

  TCanvas * cGenProj = new TCanvas("cGenProj","",800,600);
  TCanvas * cKineticEfficiency = new TCanvas("cKineticEfficiency","",800,600);  

  TCanvas * cMatrix = new TCanvas("cMatrix","",800,600);
  TCanvas * cMatrix_recoTrunc = new TCanvas("cMatrix_recoTrunc","",800,600);
  
  mPP_Matrix = (TH2F*)fin->Get(Form("hpp_matrix_HLT_R%d_20_eta_20",radius));
  mPP_Matrix->Rebin2D(10, 10);

  for (int y=1;y<=mPP_Matrix->GetNbinsY();y++) {
    double sum=0;
    for (int x=1;x<=mPP_Matrix->GetNbinsX();x++) {
      if (mPP_Matrix->GetBinContent(x,y)<=1*mPP_Matrix->GetBinError(x,y)) {
	mPP_Matrix->SetBinContent(x,y,0);
	mPP_Matrix->SetBinError(x,y,0);
      }
      sum+=mPP_Matrix->GetBinContent(x,y);
    }
    
    for (int x=1;x<=mPP_Matrix->GetNbinsX();x++) {	   
      double ratio = 1;
      // if (hGenSpectraCorrPP->GetBinContent(x)!=0) ratio = 1e5/hGenSpectraCorrPP->GetBinContent(x);
      mPP_Matrix->SetBinContent(x,y,mPP_Matrix->GetBinContent(x,y)*ratio);
      mPP_Matrix->SetBinError(x,y,mPP_Matrix->GetBinError(x,y)*ratio);
    }
  }
  // mPbPb_Matrix[i]->Smooth(0);
  
  // Ok major differences here between my code and Kurt in b-jet Tools under Unfold - lines 469 and above.  

  if(printDebug) cout<<"getting the response matrix"<<endl;

  mPP_Response = (TH2F*)mPP_Matrix->Clone("mPP_Response");
  TH1F *hProjPP = (TH1F*)mPP_Response->ProjectionY()->Clone("hProjPP");
  
  
  for (int y=1;y<=mPP_Response->GetNbinsY();y++) {
    double sum=0;
    for (int x=1;x<=mPP_Response->GetNbinsX();x++) {
      if (mPP_Response->GetBinContent(x,y)<=1*mPP_Response->GetBinError(x,y)) {
	// in the above if statement, kurt has 1*error and my old has 0*error
	mPP_Response->SetBinContent(x,y,0);
	mPP_Response->SetBinError(x,y,0);
      }
      sum+=mPP_Response->GetBinContent(x,y);
    }
    
    for (int x=1;x<=mPP_Response->GetNbinsX();x++) {  	
      if (sum==0) continue;
      double ratio = 1;
      //if(dPbPb_TrgComb[i]->GetBinContent(y)==0) ratio = 1e-100/sum;
      // else ratio = dPbPb_TrgComb[i]->GetBinContent(y)/sum
      ratio = 1./sum;
      if (hProjPP->GetBinContent(y)==0) ratio = 1e-100/sum;
      else ratio = hProjPP->GetBinContent(y)/sum;
      mPP_Response->SetBinContent(x,y,mPP_Response->GetBinContent(x,y)*ratio);
      mPP_Response->SetBinError(x,y,mPP_Response->GetBinError(x,y)*ratio);
    }
  }
  if(printDebug) cout<<"getting the normalized response matrix"<<endl;
  mPP_ResponseNorm = (TH2F*)mPP_Matrix->Clone("mPP_ResponseNorm");
  for (int x=1;x<=mPP_ResponseNorm->GetNbinsX();x++) {
    double sum=0;
    for (int y=1;y<=mPP_ResponseNorm->GetNbinsY();y++) {
      if (mPP_ResponseNorm->GetBinContent(x,y)<=1*mPP_ResponseNorm->GetBinError(x,y)) {
	mPP_ResponseNorm->SetBinContent(x,y,0);
	mPP_ResponseNorm->SetBinError(x,y,0);
      }
      sum+=mPP_ResponseNorm->GetBinContent(x,y);
    }
    
    for (int y=1;y<=mPP_ResponseNorm->GetNbinsY();y++) {  	
      if (sum==0) continue;
      double ratio = 1./sum;
      mPP_ResponseNorm->SetBinContent(x,y,mPP_ResponseNorm->GetBinContent(x,y)*ratio);
      mPP_ResponseNorm->SetBinError(x,y,mPP_ResponseNorm->GetBinError(x,y)*ratio);
//.........这里部分代码省略.........

int main() {
  float min_logL1 = 5986.94;
  float min_logL0 = 5987.16;
  
  string filepath="FINAL_RESULT_AB.root_RESULT__RESULT";
  filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_175622__RNSG_46692.4.root__RESULT__RESULT"; // final MCwS 
  
  filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_175622__RNSG_46692.4.root__RESULT__RESULT";


  filepath="/shome/buchmann/KillerKoala/CBAF/Development/exchange/RooFit__WorkSpace__Exchange_201417_141126__RNSG_97048.1.root__RESULT__RESULT";
//  *************************************************************************************

  setlumi(PlottingSetup::luminosity);
  setessentialcut(PlottingSetup::essential);  // this sets the essential cut; this one is used in the draw command so it is AUTOMATICALLY applied everywhere. IMPORTANT: Do NOT store weights here!
  stringstream resultsummary;

  // write_analysis_type(PlottingSetup::RestrictToMassPeak,PlottingSetup::DoBTag);
  do_png(true);
  do_pdf(true);
  do_eps(false);
  do_C(true);
  do_root(false);

  PlottingSetup::directoryname = "pValuePlot";
  gROOT->SetStyle("Plain");
  bool do_fat_line = false;  // if you want to have HistLineWidth=1 and FuncWidth=1 as it was before instead of 2
  setTDRStyle(do_fat_line);
  gStyle->SetTextFont(42);
  bool showList = true;
  set_directory(PlottingSetup::directoryname);  // Indicate the directory name where you'd like to save the output files in Setup.C
  set_treename("events");        // you can set the treename here to be used; options are "events" (for reco) for "PFevents" (for particle flow)


  
  TFile *f = new TFile(filepath.c_str());

  if(f->IsZombie()) {
    cout << "Seems to be a zombie. goodbye." << endl;
    return -1;
  }

  RooWorkspace *wa = (RooWorkspace*)f->Get("transferSpace");
  RooPlot *plot = (RooPlot*) wa->obj("frame_mlledge_109fde50");
  
  
//  cout << plot << endl;
wa->Print("v");
  


   TCanvas *can = new TCanvas("can","can");
   
cout << "Address of plot : " << plot << endl;
//   plot->Draw();
   
   float pVal_mllmin=35;
   float pVal_mllmax=90;
   int is_data=PlottingSetup::data;
   
   vector < std::pair < float, float> > loglikelihoods;
   
   string function="";
   for(int i=0; i< plot->numItems();i++){
     string name = plot->getObject(i)->GetName();
     if (plot->getObject(i)->IsA()->InheritsFrom( "RooCurve" ))function=name;
   }
   
   RooCurve* curve = (RooCurve*) plot->findObject(function.c_str(),RooCurve::Class()) ;
   if (!curve) { 
     dout << "RooPlot::residHist(" << plot->GetName() << ") cannot find curve" << endl ;
     return 0 ;
   }
   
   int iMinimum=0;
   float min=1e7;

   for(int i=0;i<curve->GetN();i++) {
     double x,y;
     curve->GetPoint(i,x,y);
     if(y<min & y>=0) {
       min=y;
       iMinimum=i;
     }
   }
   
   double x,y;
   curve->GetPoint(iMinimum,x,y);
   cout << "Minimum is at " << x << " : " << y << endl;
   loglikelihoods.push_back(make_pair(x,y+min_logL1));
   
   //move right starting from the minimum
   for(int i=iMinimum+1;i<curve->GetN();i++) {
     float yold=y;
     curve->GetPoint(i,x,y);
     //if(abs((y-yold)/yold)>0.5) continue;
     loglikelihoods.push_back(make_pair(x,y+min_logL1));
   }
   
   /*  
//.........这里部分代码省略.........

void PhotonIsolation()

{
  const float chargeHadronIsoCut_LWP_EB = 1.3;
  const float chargeHadronIsoCut_MWP_EB = 0.44;
  const float chargeHadronIsoCut_TWP_EB = 0.2;

  const float chargeHadronIsoCut_LWP_EE = 1.36;
  const float chargeHadronIsoCut_MWP_EE = 0.82;
  const float chargeHadronIsoCut_TWP_EE = 0.27;

  const float HoverECut_LWP_EB = 0.06;
  const float HoverECut_MWP_EB = 0.04;
  const float HoverECut_TWP_EB = 0.027;

  const float HoverECut_LWP_EE = 0.027;
  const float HoverECut_MWP_EE = 0.021;
  const float HoverECut_TWP_EE = 0.020;

  const float SigIEtaIEtaCut_LWP_EB = 0.01031;
  const float SigIEtaIEtaCut_MWP_EB = 0.01023;
  const float SigIEtaIEtaCut_TWP_EB = 0.00994;

  const float SigIEtaIEtaCut_LWP_EE = 0.0269;
  const float SigIEtaIEtaCut_MWP_EE = 0.0259;
  const float SigIEtaIEtaCut_TWP_EE = 0.0258;

  const float beamSpotZ = 0.282329;
  const float beamSpotSigmaZ = 43.131299;
  const float NPU = 140.0;

  const int n_density = 5;

  float IsoEff_NoTiming[n_density] = {0.0};// for n_density pu density bins, ranging from 0 to 2.0, bin width 0.1
  float Err_IsoEff_NoTiming[n_density] = {0.0};// for n_density pu density bins, ranging from 0 to 2.0, bin width 0.1
  float IsoEff_Timing50_TrkVtx[n_density] = {0.0}; 
  float IsoEff_Timing80_TrkVtx[n_density] = {0.0}; 
  float Err_IsoEff_Timing80_TrkVtx[n_density] = {0.0}; 
  float IsoEff_Timing120_TrkVtx[n_density] = {0.0}; 

  float IsoEff_Timing50_TrkPho[n_density] = {0.0}; 
  float IsoEff_Timing80_TrkPho[n_density] = {0.0}; 
  float IsoEff_Timing120_TrkPho[n_density] = {0.0}; 
  float Err_IsoEff_Timing120_TrkPho[n_density] = {0.0}; 


  int N_Pho_Total_NoBin = 0;
  int N_Pho_Total_NoBin200 = 0;
  int N_Pho_Total[n_density] = {0};
  int N_Pho_Total200[n_density] = {0};
  int N_Pho_PassIso_NoTiming[n_density] = {0}; 
  int N_Pho_PassIso_Timing50_TrkVtx[n_density] = {0}; 
  int N_Pho_PassIso_Timing80_TrkVtx[n_density] = {0}; 
  int N_Pho_PassIso_Timing120_TrkVtx[n_density] = {0}; 
 
  int N_Pho_PassIso_Timing50_TrkPho[n_density] = {0}; 
  int N_Pho_PassIso_Timing80_TrkPho[n_density] = {0}; 
  int N_Pho_PassIso_Timing120_TrkPho[n_density] = {0}; 
 

  TTree * tree = 0;
  TFile *f  =new TFile("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU140_Timing_Iso.root");
  TDirectory * dir = (TDirectory*)f->Get("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU140_Timing_Iso.root:/ntuples");
  dir->GetObject("RazorEvents",tree);
  //cout<<tree->GetEntries()<<endl;  
  int NEntries = tree->GetEntries();

  TTree * tree200 = 0;
  TFile *f200  =new TFile("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU200_Timing_Iso.root");
  TDirectory * dir200 = (TDirectory*)f200->Get("/afs/cern.ch/work/z/zhicaiz/public/release/CMSSW_8_1_0_pre15/src/SUSYBSMAnalysis/RazorTuplizer/python/razorNtuple_PU200_Timing_Iso.root:/ntuples");
  dir200->GetObject("RazorEvents",tree200);
  //cout<<tree->GetEntries()<<endl;  
  int NEntries200 = tree200->GetEntries();


  Int_t           nPhotons;

  Float_t         pvZ;
  Float_t         phoSigmaIetaIeta[700];
  Float_t         pho_superClusterEta[700];
  Float_t         pho_superClusterEnergy[700];
  Float_t         pho_HoverE[700];
  Bool_t          pho_passEleVeto[700];
  
  Float_t         pho_sumChargedHadronPt_NewPV_NoTiming[700];
  Float_t         pho_sumChargedHadronPt_NewPV_Timing50_TrkVtx[700];
  Float_t         pho_sumChargedHadronPt_NewPV_Timing80_TrkVtx[700];
  Float_t         pho_sumChargedHadronPt_NewPV_Timing120_TrkVtx[700];
 
  Float_t         pho_sumChargedHadronPt_NewPV_Timing50_TrkPho[700];
  Float_t         pho_sumChargedHadronPt_NewPV_Timing80_TrkPho[700];
  Float_t         pho_sumChargedHadronPt_NewPV_Timing120_TrkPho[700];
  
  tree->SetBranchAddress("nPhotons", &nPhotons);
  tree->SetBranchAddress("pvZ_New", &pvZ);
  tree->SetBranchAddress("phoSigmaIetaIeta", phoSigmaIetaIeta);
  tree->SetBranchAddress("pho_superClusterEta", pho_superClusterEta);
  tree->SetBranchAddress("pho_superClusterEnergy", pho_superClusterEnergy);
  tree->SetBranchAddress("pho_HoverE", pho_HoverE);
  tree->SetBranchAddress("pho_passEleVeto", pho_passEleVeto);
//.........这里部分代码省略.........

int mumucross()
{
  TCanvas *c1 = new TCanvas("c1","di-mu cross section");
  TF1 *f1 = new TF1("f1",mumuxsdang,2,3.1,0);
  f1->GetXaxis()->SetTitle("Ecm (GeV)");
  f1->GetYaxis()->SetTitle("#sigma (nbar)");
  f1->Draw();
  
  TCanvas *c2 = new TCanvas("c2","BhaBha cross section");
  TF1 *f2 = new TF1("f2",bhabhaxsdang,2,3.1,0);
  f2->GetXaxis()->SetTitle("Ecm (GeV)");
  f2->GetYaxis()->SetTitle("#sigma (nbar)");
  f2->Draw();
  //return 0;

  double energys[22] = {
  2.0000, 2.0500, 2.1000, 2.1500, 2.1750,
  2.2000, 2.2324, 2.3094, 2.3864, 2.3960,
  2.5000, 2.6444, 2.6464, 2.7000, 2.8000,
  2.9000, 2.9500, 2.9810, 3.0000, 3.0200,
  3.0800, 2.125
  };
  double luminosity[22] = {
  10074,   3343,   12167,  2841,   10625,
  13699,  11856,   21089, 22549,   66869,
   1098,  33722,   34003,  1034,    1008,
 105253,  15942,   16071, 15881,   17290,
  126188, 108490
  };
  double mccross[22] = {
  22.3763923, 21.2939032, 20.4018357, 19.4157221, 18.9842334,
  18.5074456, 18.0584889, 16.9121088, 15.8285505, 15.6795155,
  14.3957219, 12.9117968, 12.9257851, 12.3962790, 11.5661025,
  10.7698318, 10.4156161, 10.2122681, 10.0576468, 9.9391192,
  9.5904453,  21.1564523
  };
  double nevt[22] = {
  5.0e5, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
  5.0e5, 5.0e5, 5.0e5, 5.0e5, 1.07e6,
  5.0e5, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
  1.13e6, 5.0e5, 5.0e5, 5.0e5, 5.0e5,
  1.16e6, 1.93e6
  };

  TFile *file = new TFile("output/output_bck.root");

  TH1D *hdata;
  TH1D *hmcmu;
  TH1D *hmcKK;
  gStyle->SetOptStat(0);
  for (int i=0 ;i<22;i++) {
    //double nmumu = luminosity[i]*f1->Eval(energys[i]);
    //cout<<nmumu<<"\t events at "<< energys[i] <<" GeV."<<endl;
    //double nee = luminosity[i]*f2->Eval(energys[i]);
    //cout<<nee<<"\t events at "<< energys[i] <<" GeV."<<endl;
    cout<<"ene "<<energys[i]<<", scale is "<< luminosity[i]/(nevt[i]/mccross[i])<<endl;
    //continue;
    char hname[100];
    // get data spectrum
    //sprintf(hname,"hp_KK_%d",(int)(energys[i]*10000+0.5));
    sprintf(hname,"hp_mcKK_%.4f",energys[i]);
    hdata = (TH1D*)file->Get(hname);
    hdata->SetFillColor(0);
    //hdata->SetTitle("momentum spectrum of data and di-#mu MC");
    hdata->SetTitle("");
    hdata->GetXaxis()->SetTitle("p (GeV/c)");
    hdata->GetYaxis()->SetTitle("Counts");
    // get dimu spectrum
    sprintf(hname,"hp_mcmumu_%.4f",energys[i]);
    hmcmu = (TH1D*)file->Get(hname);

    // get mc KK spectrum
    sprintf(hname,"hp_mcKK2_%.4f",energys[i]);
    hmcKK = (TH1D*)file->Get(hname);
    int binid;
    double maxdata = getMax(hdata,1,50,binid);
    double maxmcKK = getMax(hmcKK,1,50,binid);

    cout<< hdata <<"\t"<<hmcmu<<" ndata,nmcKK "<< maxdata<<" "<<maxmcKK<<endl;
    sprintf(hname,"Canvas_%02d",i+1);
    
    TCanvas *chp = new TCanvas(hname,"Background analysis");
    chp->SetMargin(0.15,0.1,0.15,0.1);
    hdata->SetLineColor(1);
    hdata->SetTitle("");
    hdata->GetXaxis()->SetNdivisions(505);
    hdata->GetXaxis()->SetTitleSize(0.05);
    hdata->GetXaxis()->SetLabelSize(0.05);
    hdata->GetYaxis()->SetNdivisions(505);
    hdata->GetYaxis()->SetTitleSize(0.05);
    hdata->GetYaxis()->SetLabelSize(0.05);
    hdata->GetYaxis()->SetTitleOffset(1.1);

    hdata->Draw("E");
    // scale dimu to data according luminosity
    hmcmu->Scale(luminosity[i]/(nevt[i]/mccross[i]));
    hmcmu->SetLineColor(2);
    hmcmu->SetMarkerColor(2);
    hmcmu->SetFillColor(0);
    hmcmu->Draw("same");
//.........这里部分代码省略.........

//注意，主函数名称和文件名称要一样
int rootelectron(void)
{
    const double PI = 3.1415926535897932384626433832795;
    ///////////////////////%%%%%%输入参数%%%%%%%%%%%/////////////////////////////////////////////////////////////
    string filename="spectrum.root";                                                        //数据文件名称
    //root文件中tree的名称，需要和Geant4中的一致
    string treename="selectree";
    string branchname = "particledata";                                                 //root文件中branch的名称，需要和Geant4中的一致
    //    string treename[4]={"pelectree","pgamatree","selectree","sgamatree"};
    double Ebins=200;   //统计能谱的直方图中bin的个数，
    double Elow=-75.0;      //能谱的低能截止,单位keV
    double Eup=95.5555673; //能谱的高能能截止,单位keV
    double cbins[3] = {180,180,180};//统计角分布的三个直方图中bin的个数【x,y,z】
    double chlow[3] = {0,-90,0};//统计角分布的三个直方图的左边界【x,y,z】
    double chup[3] = {180,90,90};//统计角分布的三个直方图的左边界【x,y,z】
    double eneconv=1;//1000->MeV,1->keV将输入的能量转换单位

    ///////////////////////%%%%%%读入ROOT文件和Tree%%%%%%%%%%%/////////////////////////////////////////////////////////////
    gROOT->Reset();
    cout<<endl;
    cout << "\nStart\n" <<endl;
    TFile *myfile = new TFile(&filename[0]);    //新建对象，指向数据文件
if (!myfile)    {
        cout <<"no input file!" <<endl;
        return;
    }
    cout<<"Open file: "<<filename<<endl;
    TTree *mytree = (TTree *)myfile->Get(&treename[0]);//新建tree,用于处理数据
    if (!mytree)    {
        cout <<"no  tree named:  "<<treename<<endl;
        return;
    }
    cout << "input tree: "<<treename <<endl;

    //////////////////////////////%%设置新tree，用于数据处理%%%%%%/////////////////////////////////////%%%%%%%%%%%%%%%

    double data[4]={0};//保存从root文件中取出的数据。
    mytree->SetBranchAddress(&branchname[0],data);//设置新建的tree中branch的指向。
    //create two histograms
    double EEup=Eup/eneconv;
    double EElow=Elow/eneconv;
    double sintheta=0;
    TH1D * h1     = new TH1D("E","Energy Spectrum of source",Ebins,EElow,EEup);             //统计能谱
   // TH2D *h2Ex = new TH2D("h2Ex","h2 Energy vs x",cbins[0],chlow[0],chup[0],Ebins,EElow,EEup);//统计角分布，出射方向与X方向的夹角
   // TH2D *h2Ey = new TH2D("h2Ey","h2 Energy vs y",cbins[1],chlow[1],chup[1],Ebins,EElow,EEup);//统计角分布，出射方向与Y方向的夹角
    TH2D *h2Ez = new TH2D("h2Ez","h2 Energy vs z",cbins[2],chlow[2],chup[2],Ebins,EElow,EEup);//统计角分布，出射方向与Z方向的夹角（横坐标名，纵坐标名，横坐标bins，横坐标低截止，横坐标高截止，纵坐标bins，纵坐标低截止，纵坐标高截止）

    //read all entries and fill the histograms
    Long64_t nentries = mytree->GetEntries();   //注意：因为事例entry很多，所以需要Long64_t
    for (Long64_t i=0;i<nentries;i++) {
        mytree->GetEntry(i);
        //将处理数据的tree指向数据文件中的第i个事例，由于上面设置了的tree中branch的指针为变量data，所以处理data变量就处理了数据。
        h1->Fill(data[3]);                      //统计能谱
      //  sintheta= -sqrt(1-data[3]*data[3]);//data[3]=cos(theta)
       // h2Ex->Fill(acos(data[1]/sintheta)/PI*180,data[0]);
       // h2Ey->Fill(asin(data[2]/sintheta)/PI*180,data[0]);
        h2Ez->Fill(acos(data[1])/PI*180,data[0]);
    }

    ///////////////////////%%%%%%%%%Profile Histogram%%%%%%%/////////////////////////////////////
    TH1D* az;
  //  ax=h2Ex->ProjectionX();   //将二维散点图投影到1维就得到了一维的角分布，
  //  ay=h2Ey->ProjectionX();   //将二维散点图投影到1维就得到了一维的角分布，
    az=h2Ez->ProjectionX();   //将二维散点图投影到1维就得到了一维的角分布，

    ///////////////////////%%%%%%%%%Plot Histogram%%%%%%%/////////////////////////////////////
    //gROOT->SetStyle("Plain");   // uncomment to set a different style
    gStyle->SetPalette(1);        // use precomputed color palette 1
    gStyle->SetOptStat("ner");

 /*   TCanvas *c1 = new TCanvas();
    c1->cd();
    h2Ex->Draw("C");  //energy
    h2Ex->SetTitle("Energy Spectrum");     //散点图
    h2Ex->GetXaxis()->SetTitle("alpha(x)");
    h2Ex->GetYaxis()->SetTitle("Energy (keV)");
    h2Ex->GetXaxis()->CenterTitle();
    h2Ex->GetYaxis()->CenterTitle();
    gPad->Update();
    //    eneh->GetXaxis()->SetRangeUser(0, 1200);
    //    eneh->GetYaxis()->SetRangeUser(1e-8,1e-1);
    //        gPad->SetLogx();
    //            gPad->SetLogy();

    TCanvas *c2 = new TCanvas();
    c2->cd();
    h2Ey->Draw("C");  //energy
    h2Ey->SetTitle("Energy Spectrum");
    h2Ey->GetXaxis()->SetTitle("beta(y)");
    h2Ey->GetYaxis()->SetTitle("Energy (MeV)");
    h2Ey->GetXaxis()->CenterTitle();
    h2Ey->GetYaxis()->CenterTitle();
    gPad->Update();            */
   /*
    TCanvas *c3 = new TCanvas();           
    c3->cd();
    h2Ez->Draw("C");  //energy
    h2Ez->SetTitle("Energy Spectrum");
    h2Ez->GetXaxis()->SetTitle("gamma(z)");
//.........这里部分代码省略.........

void hackedbob(bool eff = false)
{
    using namespace std;
    
    double integral, preInt = 1.0;
    
    vector<pair<string,string> > fnames;
    vector<pair<string,string> > hnames;

    //TFile * tfile = new TFile("/home/ugrad/pastika/cms/HeavyNu/CMSSW_4_2_7/src/HeavyNu/AnalysisModules/analysis.root");
    TFile * tfile = new TFile("/local/cms/user/pastika/heavynu/heavynu_2011Bg_summer11_DYToLL_M-50_7TeV-sherpa_nov3.root");
    hnames.push_back(pair<string,string>("none","hNu%s/cut0_none/mWR"));
    hnames.push_back(pair<string,string>("LLJJ Pt","hNu%s/cut1_LLJJpt/mWR"));
    hnames.push_back(pair<string,string>("trig","hNu%s/cut2_TrigMatches/mWR"));
    hnames.push_back(pair<string,string>("vertex","hNu%s/cut3_Vertex/mWR"));
    hnames.push_back(pair<string,string>("mu1 pt","hNu%s/cut4_Mu1HighPt/mWR"));
    hnames.push_back(pair<string,string>("Mll","hNu%s/cut5_diLmass/mWR"));
    hnames.push_back(pair<string,string>("MWR","hNu%s/cut6_mWRmass/mWR"));

    fnames.push_back(pair<string,string>("nom Mu24", "Mu24"));
    fnames.push_back(pair<string,string>("nom Mu40", "Mu40"));
    //fnames.push_back(pair<string,string>("JES Hi Mu24", "Mu24jesHi"));
    //fnames.push_back(pair<string,string>("JES Lo Mu24", "Mu24jesLo"));
    //fnames.push_back(pair<string,string>("JES Hi Mu40", "Mu40jesHi"));
    //fnames.push_back(pair<string,string>("JES Lo Mu40", "Mu40jesLo"));
    //fnames.push_back(pair<string,string>("MES Hi Mu24", "Mu24mesHi"));
    //fnames.push_back(pair<string,string>("MES Lo Mu24", "Mu24mesLo"));
    //fnames.push_back(pair<string,string>("MES Hi Mu40", "Mu40mesHi"));
    //fnames.push_back(pair<string,string>("MES Lo Mu40", "Mu40mesLo"));
    //fnames.push_back(pair<string,string>("MuID Hi Mu24", "Mu24midHi"));
    //fnames.push_back(pair<string,string>("MuID Lo Mu24", "Mu24midLo"));
    //fnames.push_back(pair<string,string>("MuID Hi Mu40", "Mu40midHi"));
    //fnames.push_back(pair<string,string>("MuID Lo Mu40", "Mu40midLo"));
    //fnames.push_back(pair<string,string>("Trig Hi Mu24", "Mu24trigHi"));
    //fnames.push_back(pair<string,string>("Trig Lo Mu24", "Mu24trigLo"));
    //fnames.push_back(pair<string,string>("Trig Hi Mu40", "Mu40trigHi"));
    //fnames.push_back(pair<string,string>("Trig Lo Mu40", "Mu40trigLo"));
    fnames.push_back(pair<string,string>("Pileup Hi Mu24", "Mu24puHi"));
    fnames.push_back(pair<string,string>("Pileup Lo Mu24", "Mu24puLo"));
    fnames.push_back(pair<string,string>("Pileup Hi Mu40", "Mu40puHi"));
    fnames.push_back(pair<string,string>("Pileup Lo Mu40", "Mu40puLo"));
    
	printf("%8s", "cut");
	for(vector<pair<string,string> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
	{
		printf(eff?" & %15s        ":" & %15s", i->first.c_str());
	}
	printf(" \\\\ \\hline\n");
	for(vector<pair<string,string> >::const_iterator j = hnames.begin(); j != hnames.end(); j++)
	{
		printf("%8s", j->first.c_str());
		for(vector<pair<string,string> >::const_iterator i = fnames.begin(); i != fnames.end(); i++)
		{
		    char histname[256];
		    sprintf(histname, j->second.c_str(), i->second.c_str());
			TH1* h = (TH1*)tfile->Get(histname);
			integral = h->Integral(0, h->GetNbinsX()+1);
			
			if(integral > 10) printf((j == hnames.begin() && eff)?" & %15.0f        ":" & %15.0f", integral);
			else if(integral > 1) printf((j == hnames.begin() && eff)?" & %15.2f        ":" & %15.2f", integral);
			else printf((j == hnames.begin() && eff)?" & %15.3f        ":" & %15.3f", integral);
			//if(j != hnames.begin() && eff) 
			//{
			//	TH1* hpre = (TH1*)tfile->Get((j-1)->second.c_str());
			//	preInt = hpre->Integral(0, h->GetNbinsX()+1);
			//	printf("  (%4.2f)", integral/preInt);
			//}
		} 
		printf(" \\\\ \\hline\n");
	}
}

void plotMVAOutput( bool printgif = false ){

  //gROOT->ProcessLine(".x selection.h");

  //char* path    = "Trainings/H130_WWTo2L2Nu/output";
  //char* path    = "Trainings/H130_WWTo2L2Nu_WJetsToLNu/output";
  char* path    = "Trainings/H130_allbkg_4vars/output";

  //char* mvaname = "MVA_PDERS";
  //char* mvaname = "MVA_MLPBNN";
  //char* mvaname = "MVA_BDT";
  //char* mvaname = "LikelihoodPCA";

  vector<char*> mvanames;
  mvanames.push_back("BDT");
  mvanames.push_back("MLPBNN");
  const unsigned int nmva = mvanames.size();

  int rebin     = 10;
  int colors[]  = { 5 , 2 , 4 , 3 , 7 , 8 , 6 , 9 , 10};

  vector<char*> samples;
  samples.push_back("WWTo2L2Nu");
  samples.push_back("GluGluToWWTo4L");
  samples.push_back("WZ");
  samples.push_back("ZZ");
  samples.push_back("TTJets");
  samples.push_back("tW");
  samples.push_back("WJetsToLNu");
  samples.push_back("DY");
  const unsigned int nsamples = samples.size();

  char* higgssample = "Higgs130";
  //char* higgssample = "Higgs160";
  //char* higgssample = "Higgs200";

  TCanvas *can[nmva];


  for( unsigned int imva = 0 ; imva < nmva ; ++imva ){


    TFile*   file     = new TFile();
    TH1F*    hist     = new TH1F();
    TH1F*    bkghist  = new TH1F();
    THStack* bkgstack = new THStack("bkgstack","bkgstack");
    TLegend *leg      = new TLegend(0.3,0.7,0.5,0.9);
    leg->SetBorderSize(1);
    leg->SetFillColor(0);

    //loop over backgrounds
    for( unsigned int i = 0 ; i < nsamples ; ++i ){

      file = TFile::Open(Form("%s/%s.root",path,samples.at(i)));
      hist = cloneHist( (TH1F*) file->Get( Form("MVA_%s",mvanames.at(imva) ) ) );
   
      hist->SetFillColor(colors[i]);
    
      leg->AddEntry(hist,samples.at(i),"f");
  
      if( i == 0 ) bkghist = (TH1F*) hist->Clone();
      else         bkghist -> Add(hist);

      hist->Rebin( rebin );
      bkgstack->Add(hist);

    }

    //higgs sample
    file = TFile::Open(Form("%s/%s.root",path,higgssample));
    TH1F* higgshist = cloneHist( (TH1F*) file->Get( Form("MVA_%s",mvanames.at(imva) ) ) );
    higgshist->SetLineWidth(2);
    leg->AddEntry(higgshist,higgssample,"l");



    float bkg    = 0;
    float sig    = 0;
    float minbkg = 1.48;
    //float minbkg = 1.10;
    float cut    = 0.;

    for( int ibin = 1 ; ibin < bkghist->GetNbinsX() ; ibin++ ){

      bkg = bkghist->Integral( ibin , 10000 );
      sig = higgshist->Integral( ibin , 10000 );

      if( bkg < minbkg ){
        cut = bkghist->GetBinCenter(ibin);
        cout << endl;
        cout << "S/B       " << sig/bkg    << endl;
        cout << "Sig       " << sig        << endl;
        cout << "Bkg       " << bkg        << endl;
        cout << "cut value " << cut        << endl;
        break;
      }

    }

    float cutsig = sig;
//.........这里部分代码省略.........

float ComputeTestStat(TString wsfile, double mu_susy_sig_val) {

  gROOT->Reset();

  TFile* wstf = new TFile( wsfile ) ;

  RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
  ws->Print() ;
  
  ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
  
  modelConfig->Print() ;

  RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;
  
  rds->Print() ;
  rds->printMultiline(cout, 1, kTRUE, "") ;
  
  RooAbsPdf* likelihood = modelConfig->GetPdf() ;
  
  RooRealVar* rrv_mu_susy_sig = ws->var("mu_susy_all0lep") ;
  if ( rrv_mu_susy_sig == 0x0 ) {
    printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
    return ;
  } else {
    printf(" current value is : %8.3f\n", rrv_mu_susy_sig->getVal() ) ; cout << flush ;
    rrv_mu_susy_sig->setConstant(kFALSE) ;
  }

  /*
  // check the impact of varying the qcd normalization:

  RooRealVar *rrv_qcd_0lepLDP_ratioH1 = ws->var("qcd_0lepLDP_ratio_H1");
  RooRealVar *rrv_qcd_0lepLDP_ratioH2 = ws->var("qcd_0lepLDP_ratio_H2");
  RooRealVar *rrv_qcd_0lepLDP_ratioH3 = ws->var("qcd_0lepLDP_ratio_H3");
  
  rrv_qcd_0lepLDP_ratioH1->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH2->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH3->setVal(0.3);
  
  rrv_qcd_0lepLDP_ratioH1->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH2->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH3->setConstant(kTRUE);
  */
  
  printf("\n\n\n  ===== Doing a fit with SUSY component floating ====================\n\n") ;

  RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
  double logLikelihoodSusyFloat = fitResult->minNll() ;
  
  double logLikelihoodSusyFixed(0.) ;
  double testStatVal(-1.) ;
  if ( mu_susy_sig_val >= 0. ) {
    printf("\n\n\n  ===== Doing a fit with SUSY fixed ====================\n\n") ;
    printf(" fixing mu_susy_sig to %8.2f.\n", mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setVal( mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setConstant(kTRUE) ;
    
    fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
    logLikelihoodSusyFixed = fitResult->minNll() ;
    testStatVal = 2.*(logLikelihoodSusyFixed - logLikelihoodSusyFloat) ;
    printf("\n\n\n ======= test statistic : -2 * ln (L_fixed / ln L_max) = %8.3f\n\n\n", testStatVal ) ;
  }


  return testStatVal ;

}

void makeStackedPlot_Jet1pT(){

  gROOT->SetStyle("Plain");
  gStyle->SetErrorX(0);

  TCanvas c1("c1","Stacked Histogram",10,10,700,800);
  TPad *p_2=new TPad("p_2", "p_2", 0, 0, 1, 0.35);
  TPad *p_1=new TPad("p_1", "p_1", 0, 0.35, 1, 1);
  p_1->SetBottomMargin(0.005);
  p_1->SetFillStyle(4000);
  p_1->SetFrameFillColor(0);
  p_2->SetFillStyle(4000);
  p_2->SetFrameFillColor(0);
  p_1->SetLogy();
  p_1->Draw();
  p_2->Draw();
  p_1->cd();

  double qcd_100_200 = 27540000.0*209.0/80093092.0;
  double qcd_200_300 = 1735000.0*209.0/18717349.0;
  double qcd_300_500 = 366800.0*209.0/20086103.0;
  double qcd_500_700 = 29370.0*209.0/19542847.0;
  double qcd_700_1000 = 6524.0*209.0/15011016.0;
  double qcd_1000_1500 = 1064.0*209.0/4963895.0;
  double qcd_1500_2000 = 121.5*209.0/3848411.0;
  double qcd_2000_inf =  25.4*209.0/1961774.0;

  THStack hs("hs","Jet pT > 50 GeV");

  TFile* file0 = TFile::Open("output_QCD_HT_100_200.root");
  TH1F *h0 = (TH1F*)file0->Get("general/jet_pT50");
  h0->Rebin(1);
  h0->GetXaxis()->SetRangeUser(0, 300);
  h0->Scale(qcd_100_200); 
  h0->SetFillColor(kYellow);
  h0->Draw();
  hs.Add(h0);

  TFile* file1 = TFile::Open("output_QCD_HT_200_300.root");
  TH1F *h1 = (TH1F*)file1->Get("general/jet_pT50");
  h1->Rebin(1);
  h1->GetXaxis()->SetRangeUser(0, 300);
  h1->Scale(qcd_200_300);
  h1->SetFillColor(kYellow);
  h1->Draw();
  hs.Add(h1);

  TFile* file2 = TFile::Open("output_QCD_HT_300_500.root");
  TH1F *h2 = (TH1F*)file2->Get("general/jet_pT50");
  h2->Rebin(1);
  h2->GetXaxis()->SetRangeUser(0, 300);
  h2->Scale(qcd_300_500);
  h2->SetFillColor(kYellow);
  h2->Draw();
  hs.Add(h2);

  TFile* file3 = TFile::Open("output_QCD_HT_500_700.root");
  TH1F *h3 = (TH1F*)file3->Get("general/jet_pT50");
  h3->Rebin(1);
  h3->GetXaxis()->SetRangeUser(0, 300);
  h3->Scale(qcd_500_700);
  h3->SetFillColor(kYellow);
  h3->Draw();
  hs.Add(h3);

  TFile* file4 = TFile::Open("output_QCD_HT_700_1000.root");
  TH1F *h4 = (TH1F*)file4->Get("general/jet_pT50");
  h4->Rebin(1);
  h4->GetXaxis()->SetRangeUser(0, 300);
  h4->Scale(qcd_700_1000);
  h4->SetFillColor(kYellow);
  h4->Draw();
  hs.Add(h4);

  TFile* file5 = TFile::Open("output_QCD_HT_1000_1500.root");
  TH1F *h5 = (TH1F*)file5->Get("general/jet_pT50");
  h5->Rebin(1);
  h5->GetXaxis()->SetRangeUser(0, 300);
  h5->Scale(qcd_1000_1500);
  h5->SetFillColor(kYellow);
  h5->Draw();
  hs.Add(h5);

  TFile* file6 = TFile::Open("output_QCD_HT_1500_2000.root");
  TH1F *h6 = (TH1F*)file6->Get("general/jet_pT50");
  h6->Rebin(1);
  h6->GetXaxis()->SetRangeUser(0, 300);
  h6->Scale(qcd_1500_2000);
  h6->SetFillColor(kYellow);
  h6->Draw();
  hs.Add(h6);

  TFile* file7 = TFile::Open("output_QCD_HT_2000_inf.root");
  TH1F *h7 = (TH1F*)file7->Get("general/jet_pT50");
  h7->Rebin(1);
  h7->GetXaxis()->SetRangeUser(0, 300);
  h7->Scale(qcd_2000_inf);
  h7->SetFillColor(kYellow);
  h7->Draw();
  hs.Add(h7);
//.........这里部分代码省略.........

void Plot_AM_events_07Sep_Susy_1_auto() {
 TString cutNameBefore = Form("Data/");

 //                            cut_variable 
 TString cutNameAfter  = Form("_1_6_Tot_temp");

 gROOT->LoadMacro("PlotVHqqHggH.C+");
 gInterpreter->ExecuteMacro("LatinoStyle2.C");

 TCanvas* c1 = new TCanvas("events","events",500,600);
 TFile* f = new TFile("~/Cern/Code/VBF/qqHWW/AnalysisPackage_qqHWWlnulnu/out_test_Latinos_07Sep2012_2000_Run2012AB_8TeV_SUSY.root");


 PlotVHqqHggH* hs = new PlotVHqqHggH();

 hs->setLumi(5.063);
 hs->setLabel("event");
 hs->addLabel("    #sqrt{s} = 8 TeV");

 TString name;

 std::vector<int> vectColourBkg;        
 std::vector<double> vectSystBkg;       
 std::vector<std::string> vectNameBkg;  
 std::vector<TH1F*> vectTHBkg;          

 std::vector<int> vectColourSig;      
 std::vector<double> vectSystSig;       
 std::vector<std::string> vectNameSig;  
 std::vector<TH1F*> vectTHSig;          

 ///==== signal (begin) ====

 name = Form("%sT2tt-350-70%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop350-lsp70");
 vectColourSig.push_back(6);

 name = Form("%sT2tt-500-70%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop500-lsp70");
 vectColourSig.push_back(97);

 name = Form("%sT2tt-350-100%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop350-lsp100");
 vectColourSig.push_back(70);

 name = Form("%sT2tt-500-100%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop500-lsp100");
 vectColourSig.push_back(65);

 name = Form("%sT2tt-500-200%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop500-lsp200");
 vectColourSig.push_back(5);

 name = Form("%sT2tt-200-0%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHSig.push_back ( (TH1F*) f->Get(name) );
 vectNameSig.push_back ("T2tt-stop200-lsp0");
 vectColourSig.push_back(7);

 ///==== signal (end)  ====

 name = Form("%sDATA%s",cutNameBefore.Data(),cutNameAfter.Data());
 hs->setDataHist ((TH1F*)f->Get(name));



 hs->setBlindBinSx(0);
 hs->setBlindBinDx(0);

 hs->setCutSx(-999,">");
 hs->setCutDx(-999,"<");


 ///==== background (begin)  ====

 name = Form("%sH-125%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHBkg.push_back ( (TH1F*) f->Get(name) );
 vectNameBkg.push_back ("H-125");
 vectColourBkg.push_back(633);
 vectSystBkg.push_back(0.00);


 name = Form("%sVV%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHBkg.push_back ( (TH1F*) f->Get(name) );
 vectNameBkg.push_back ("WZ/ZZ");
 vectColourBkg.push_back(858);
 vectSystBkg.push_back(0.00);

 name = Form("%sWJets%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHBkg.push_back ( (TH1F*) f->Get(name) );
 vectNameBkg.push_back ("W+jets");
 vectColourBkg.push_back(921);
 vectSystBkg.push_back(0.36);

 name = Form("%sTop%s",cutNameBefore.Data(),cutNameAfter.Data());
 vectTHBkg.push_back ( (TH1F*) f->Get(name) );
//.........这里部分代码省略.........

void DrawStatUncertainty()
{
  TFile *f;
  TString ss("ExpStat_PYTHIA_canvasOutFile.root");
  f = new TFile(ss);
  TCanvas *c
  c = (TCanvas*)f->Get("_R92777");
  TString algoTitle;
  algoTitle = "k_{T} D = 0.6";
  TList *li = (TList*)gPad->GetListOfPrimitives();
  TObject *obj;
  
  TIter next(li);
  TH1D *histList[100], *htmp;
  int i,j,N;
  

  while ((obj = (TObject*)next()))
    {
      TString cname = obj->ClassName();
      TString name = obj->GetName();
      
      cout << cname <<" "<<name<<endl;    
      if (cname=="TH1D")
        {
          
          histList[N] = (TH1D*)gPad->FindObject(obj);
          histList[N]->SetFillColor(0);
          histList[N]->SetFillStyle(0);
          N++;
        } 
    }
  
  TCanvas *c = new TCanvas("c","c");
  gPad->SetLogx();
  gPad->SetLogy();
  //histList[0]->SetMaximum(1e+10);
  histList[0]->SetTitle("");
  histList[0]->GetXaxis()->SetTitle("jet p_{T} (GeV)");
  histList[0]->GetXaxis()->SetTitleFont(42);
  histList[0]->GetXaxis()->SetLabelFont(42);
  histList[0]->GetXaxis()->SetTitleSize(0.05);
  histList[0]->GetXaxis()->SetLabelSize(0.05);
  histList[0]->GetYaxis()->SetTitle("Relative Statistical Uncertainty");
  histList[0]->GetYaxis()->SetTitleFont(42);
  histList[0]->GetYaxis()->SetLabelFont(42);
  histList[0]->GetYaxis()->SetTitleSize(0.05);
  histList[0]->GetYaxis()->SetLabelSize(0.05);
  histList[0]->GetYaxis()->SetNdivisions(505);

  histList[0]->SetLineWidth(1);
  histList[0]->SetLineStyle(1);
  histList[0]->SetLineColor(1);
  histList[0]->SetMarkerColor(1);
  histList[0]->SetMarkerStyle(20);
  histList[0]->SetMarkerSize(1.2); 

  histList[1]->SetLineColor(2);
  histList[1]->SetMarkerColor(2);
  histList[1]->SetLineWidth(1);
  histList[1]->SetMarkerStyle(25);
  histList[1]->SetLineStyle(1);
  histList[1]->SetMarkerSize(1.2); 

  histList[2]->SetLineWidth(1);
  histList[2]->SetLineStyle(1);
  histList[2]->SetLineColor(4);
  histList[2]->SetMarkerColor(4);
  histList[2]->SetMarkerStyle(22);
  histList[2]->SetMarkerSize(1.2); 

  histList[0]->Draw("P");
  histList[1]->Draw("sameP");
  histList[2]->Draw("sameP");

  TLegend *leg = new TLegend(0.47,0.2,0.92,0.35);
  leg->SetTextFont(42);
  leg->AddEntry(histList[0],"|y| < 0.55","P");
  leg->AddEntry(histList[1],"1.10 < |y| < 1.70","P");
  leg->AddEntry(histList[2],"1.70 < |y| < 2.50","P");
  leg->SetFillColor(0);
  leg->SetBorderSize(0);
  leg->Draw();

  TPaveText *pave3 = new TPaveText(0.2,0.65,0.45,0.9,"NDC");
  pave3->SetTextFont(42);
  pave3->AddText("CMS preliminary");
  pave3->AddText(algoTitle);
  pave3->AddText("#sqrt{s} = 10 TeV");
  pave3->AddText("");
  pave3->AddText("#int L dt = 10 pb^{-1}");
  pave3->SetFillColor(0);
  pave3->SetBorderSize(0);
  pave3->Draw();

  c->Print("KT6_StatUncertainty.eps");
  
}

int Ana(const char *filename, const char* outdir, TFile *fileout)
{
  //TFile *file = new TFile("mc_KPI_22324.root");
  TFile *file;
  if (filename==0) file= new TFile("KK_22324.root");
  else file = new TFile(filename);
  std::cout<<"File name is "<<file->GetName()<<std::endl;
  if (file==0) return -1;
  TTree *tree = (TTree*)file->Get("TwoProng");
  if (tree==0) return -2;

  double ene[22];
  double pcut[22];
  double epcut[22];
  double thecut[22];
  double m_pcut;
  double m_epcut;
  double m_thecut;
  ifstream cutin("cutpar");
  char line[1000];
  if (cutin.is_open()) cutin.getline(line,1000);
  int iene=0;
  // set E/p and p cut, mark: set cut
  while (!cutin.eof()){
    cutin.getline(line,1000);
    istringstream iss;
    iss.str(line);
    iss>>ene[iene]>>pcut[iene]>>epcut[iene];
    iene++;
    if (iene==22) break;
  }
  for (int i=0;i<22;i++) epcut[i] = 1.64295 - 0.629622*ene[i] + 0.104755 *pow(ene[i],2);
  //for (int i=0;i<21;i++) thecut[i] = 173.946 + 1.74736*ene[i];
  //for (int i=0;i<22;i++) thecut[i] = 179;

  double kappx,kappy,kappz,kampx,kampy,kampz;
  int nneu;
  int run;
  int idxmc;
  int pdgid[100];
  int motheridx[100];
  int emcstatusInt;
  short emcstatusShort;
  double emctrk1;
  double emctrk2;
  double epratio1;
  double epratio2;
  int ntof1;
  int ntof2;
  int tofl1[5];
  int tofl2[5];
  double tof1[5];
  double tof2[5];
  tree->SetBranchAddress("run",&run);
  tree->SetBranchAddress("indexmc",&idxmc);
  tree->SetBranchAddress("pdgid",pdgid);
  tree->SetBranchAddress("motheridx",motheridx);
  tree->SetBranchAddress("kappx",&kappx);
  tree->SetBranchAddress("kappy",&kappy);
  tree->SetBranchAddress("kappz",&kappz);
  tree->SetBranchAddress("kampx",&kampx);
  tree->SetBranchAddress("kampy",&kampy);
  tree->SetBranchAddress("kampz",&kampz);
  //tree->SetBranchAddress("costheta",&costheta);
  tree->SetBranchAddress("nneu",&nneu);
  if (strncmp(tree->GetBranch("emcstatus")->GetLeaf("emcstatus")->GetTypeName(),"Short_t",7)==0)
  	tree->SetBranchAddress("emcstatus",&emcstatusShort);
  else 
  	tree->SetBranchAddress("emcstatus",&emcstatusInt);
  tree->SetBranchAddress("epratio1",&epratio1);
  tree->SetBranchAddress("epratio2",&epratio2);
  tree->SetBranchAddress("emctrk1",&emctrk1);
  tree->SetBranchAddress("emctrk2",&emctrk2);
  tree->SetBranchAddress("ntof1",&ntof1);
  tree->SetBranchAddress("toflayer1",tofl1);
  tree->SetBranchAddress("tof1",tof1);
  tree->SetBranchAddress("ntof2",&ntof2);
  tree->SetBranchAddress("toflayer2",tofl2);
  tree->SetBranchAddress("tof2",tof2);
  double pi = TMath::Pi();
  double mka = 0.493677;
  //double mka = 0.13957;
  //double mka = 0.1057;
  

  tree->GetEntry(0);
  double Ebeam = GetEnergy(run);
  if (Ebeam<1.0) Ebeam = getEne(filename);
  //Ebeam = 3.08;
  for (int iene=0;iene<22;iene++){
    if (fabs(Ebeam-ene[iene])<1e-4) {
      m_pcut = pcut[iene];
      m_epcut = epcut[iene];
      m_thecut = thecut[iene];
      break;
    }
  }
  //m_epcut = m_epcut +0.05; // change E/p to determine uncertainty from this cut
  //m_thecut = 179+0.2; // uncertainty from theta cut
  //m_pcut = 0.956409;
//.........这里部分代码省略.........

void Draweff(){
	int sth=0, Gth=0;
	TFile *f = TFile::Open(outG);
	if(sth==0){TString dirname = "std";}
        else if(sth==1){TString dirname ="Gri055";}
        else {TString dirname ="Gri101";}
	gStyle->SetErrorX(0);
	TString name;
        TObjString* dataname = (TObjString*)f->Get(Form("dataname"));
        TObjString* histoname = (TObjString*)f->Get(Form("histoname"));
        if(Gth==0)
                name = "G0";
        else if(Gth<nGlau)
                name = Form("Glau_%d",Gth);
        else
                name = Form("bin_%d",Gth-nGlau+1);
        TObjString* Glaubername = (TObjString*)f->Get(Form("%s/%s/Glaubername",dirname.Data(),name.Data()));
        TVectorD* xmin = (TVectorD*)f->Get(Form("%s/%s/xmin",dirname.Data(),name.Data()));
        TVectorD* xmax = (TVectorD*)f->Get(Form("%s/%s/xmax",dirname.Data(),name.Data()));
        TVectorD* mubest = (TVectorD*)f->Get(Form("%s/%s/mubest",dirname.Data(),name.Data()));
        TVectorD* kbest = (TVectorD*)f->Get(Form("%s/%s/kbest",dirname.Data(),name.Data()));
        TVectorD* Ndf = (TVectorD*)f->Get(Form("%s/%s/Ndf",dirname.Data(),name.Data()));
        TVectorD* chis = (TVectorD*)f->Get(Form("%s/%s/chis",dirname.Data(),name.Data()));
	TVectorD *kpoint = (TVectorD*)f->Get(Form("%s/%s/kpoint",dirname.Data(),name.Data()));
	TFile *fdata = TFile::Open(dataname->GetString());
	TH1D *histo_obs = (TH1D*)fdata->Get(histoname->GetString());
	histo_obs->Sumw2();
	TFile *fGlauber = TFile::Open(Glaubername->GetString());
        int binnum = histo_obs->GetNbinsX();
        double Minx = histo_obs->GetXaxis()->GetXmin();
        double Maxx = histo_obs->GetXaxis()->GetXmax();
	double binsize = (Double_t)(Maxx-Minx)/binnum;
        int xbinmin=(int)(((*xmin)[0]-Minx)/binsize);
        int xbinmax=(int)(((*xmax)[0]-Minx)/binsize);
	TH1D *histo_exp = new TH1D("histo_exp","Simulated distribution;Multiplicity;Event Fraction",binnum,Minx,Maxx);
	histo_exp->Sumw2();
	Int_t ibin;
		
	TH1D *histo_obs_norm = (TH1D*)histo_obs->Clone();
        histo_obs_norm->Scale(1/histo_obs->Integral(xbinmin,xbinmax));

	TF1 *NBD_fun = new 
	TF1("NBD_fun","[0]*TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1]))*TMath::Power([2]/[1],x)/TMath::Power([2]/[1]+1,x+[1])",0,100);
	NBD_fun->SetParameter(0,1);	//[0]: Normalized constant
	NBD_fun->SetParameter(1,(*kbest)[0]);	//[1]: k value
	NBD_fun->SetParameter(2,(*mubest)[0]);	//[2]: mu value
		
	TTree *t = (TTree*) fGlauber->Get("nt_Pb_Pb");
	Long_t Nevent;

	Nevent = (Long_t) t->GetEntries();

	Long_t Ev;	Int_t Bino;	Double_t Para, Bi_Para, Mult;
	Float_t Ncoll;
	t->SetBranchAddress("Ncoll",&Ncoll);

	for(Ev=0; Ev<Nevent; Ev++){
		if(Ev%100000==0)	 cout<<"Have run "<<Ev<<" events"<<endl;
		t->GetEntry(Ev);
		Para = 0; //make sure that Para doesn't accumulate through loops
		for(Bino=0; Bino<Ncoll; Bino++){
                         Bi_Para = NBD_fun->GetRandom();
                         Para += Bi_Para;
		}	
		histo_exp->Fill(Para);
	}
	Double_t SumEvent, scale;
	SumEvent = histo_exp->Integral(xbinmin,xbinmax);
	scale = 1/SumEvent;
	TH1D *histo_exp_norm = (TH1D*) histo_exp->Clone();
	histo_exp_norm->Scale(scale);

	TCanvas *c1 = new TCanvas();
        gStyle->SetOptStat(kFALSE);

	double hfbin[]={0,1,2,3,4,6,8,10,13,16,20,25,30,40,55,70,90};
	int nhfbin = 16;
	rehisto_obs_norm = (TH1D*)histo_obs_norm->Rebin(nhfbin,"rehisto_obs_norm",hfbin);
	normalizeByBinWidth(rehisto_obs_norm);
	rehisto_exp_norm = (TH1D*)histo_exp_norm->Rebin(nhfbin,"rehisto_exp_norm",hfbin);
	normalizeByBinWidth(rehisto_exp_norm);
	TH1D* ratio = (TH1D*)rehisto_obs_norm->Clone("ratio");
	ratio->Divide(rehisto_exp_norm);
        ratio->SetMaximum(1.2);
        ratio->SetMinimum(0);
        ratio->GetXaxis()->SetTitle("HF #Sigma E_{T}");
       	ratio->GetYaxis()->SetTitle("ratio");
/*

void make_rValues(std::string model, int m1, int m2){

  //setup file to hold r-values if it doesn't already exist
  TFile *f_temp = new TFile(Form("r-values_%s.root", model.c_str()), "NEW");
  if(f_temp){
    // default case: T5ZZ binning
    int m1_max = 2125;
    int m1_min = 775;
    int m2_max = 2025;
    int m2_min = 75;
    int m1_div = 50;
    int m2_div = 50;

    if (model.find("TChiWZ") != std::string::npos) {
      m1_max = 712.5;
      m1_min = 87.5;
      m2_max = 305;
      m2_min = -5;
      m1_div = 25;
      m2_div = 10;
    }

    int nbinsx = (m1_max - m1_min)/m1_div;
    int nbinsy = (m2_max - m2_min)/m2_div;
    
    TH2F *hExp   = new TH2F("hExp",   "hExp"  , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    TH2F *hObs   = new TH2F("hObs",   "hObs"  , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    TH2F *hExp1m = new TH2F("hExp1m", "hExp1m", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    TH2F *hExp2m = new TH2F("hExp2m", "hExp2m", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    TH2F *hExp1p = new TH2F("hExp1p", "hExp1p", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    TH2F *hExp2p = new TH2F("hExp2p", "hExp2p", nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);
    //TH2F *hSig   = new TH2F("hSig",   "hSig"  , nbinsx, m1_min, m1_max, nbinsy, m2_min, m2_max);

    f_temp->Write();
    f_temp->Close();
  }
  delete f_temp;

  //This file is created earlier by running combine
  TFile *limit_file = new TFile(Form("limit_%s_%d_%d.root", model.c_str(), m1, m2), "READ");
  TTree *limit_tree = (TTree*)limit_file->Get("limit");

  //This file is created earlier by running combine
  //TFile *significance_file = new TFile(Form("significance_%s_%d_%d.root", model.c_str(), m1, m2), "READ");
  //TTree *significance_tree = (TTree*)significance_file->Get("limit");

  //This file will hold 2d histograms with limit r-values and discovery significance
  TFile *f = new TFile(Form("r-values_%s.root", model.c_str()), "UPDATE");

  double value = -1.0;
  limit_tree->SetBranchAddress("limit", &value);

  double limit   = -1.0; //observed limit
  double rm2s    = -1.0; //expected - 2 sigma
  double rm1s    = -1.0; //expected - 1 sigma
  double rmedian = -1.0; //expected limit
  double rp1s    = -1.0; //expected + 1 sigma
  double rp2s    = -1.0; //expected + 2 sigma

  for(int i=0; i< 6; i++){
    limit_tree->GetEntry(i);
    if(i==0)      rm2s    = value;
    else if(i==1) rm1s    = value;
    else if(i==2) rmedian = value;
    else if(i==3) rp1s    = value;
    else if(i==4) rp2s    = value;
    else if(i==5) limit   = value;
  }

  //double sig = -1.0;
  //significance_tree->SetBranchAddress("limit", &sig);
  //significance_tree->GetEntry(0);

  delete limit_tree;
  limit_file->Close();
  delete limit_file;

  //delete significance_tree;
  //significance_file->Close();
  //delete significance_file;

  // // not sure what this is protecting against..
  // if( value > 21263 ) exit(0);
  
  f->cd();

  TH2F *hExp   = (TH2F*)f->Get("hExp");
  TH2F *hObs   = (TH2F*)f->Get("hObs");
  TH2F *hExp1m = (TH2F*)f->Get("hExp1m");
  TH2F *hExp2m = (TH2F*)f->Get("hExp2m");
  TH2F *hExp1p = (TH2F*)f->Get("hExp1p");
  TH2F *hExp2p = (TH2F*)f->Get("hExp2p");
  //TH2F *hSig   = (TH2F*)f->Get("hSig");

  Fill2d(hExp  , rmedian, m1, m2);    
  Fill2d(hObs  , limit  , m1, m2);  
  Fill2d(hExp1m, rm1s   , m1, m2);  
  Fill2d(hExp2m, rm2s   , m1, m2);  
  Fill2d(hExp1p, rp1s   , m1, m2);  
  Fill2d(hExp2p, rp2s   , m1, m2);  
//.........这里部分代码省略.........

void M2gHist( UInt_t chan)
{
	UInt_t cbin;
	Double_t factor;
	TString name;

	TFile* file;

	cbin = chan+1;

	// Target FULL
	// Prompt
	TH2D *hP = (TH3D*)full.Get( "THR_M2gTaggP_v_TChanM2gP");
	hP->GetXaxis()->SetRange( cbin, cbin);
	TH1D *hPf_y = hP->ProjectionY( "M2gP");

	// Random
	TH3D *hR = (TH3D*)full.Get( "THR_M2gTaggR_v_TChanM2gR");
	hR->GetXaxis()->SetRange( cbin, cbin);
	TH1D *hRf_y = hR->ProjectionY( "M2gR");
	hRf_y->Scale( 0.5);

	// Subtracted
	hSf = (TH1D*) hPf_y->Clone( "FullSubt");
	hSf->Sumw2();
	hSf->Add( hRf_y, -1.0);

	// Target EMPTY
	// Prompt
	TH2D *hP = (TH3D*)empty.Get( "THR_M2gTaggP_v_TChanM2gP");
	hP->GetXaxis()->SetRange( cbin, cbin);
	TH1D *hPf_y = hP->ProjectionY( "M2gP");

	// Random
	TH3D *hR = (TH3D*)empty.Get( "THR_M2gTaggR_v_TChanM2gR");
	hR->GetXaxis()->SetRange( cbin, cbin);
	TH1D *hRe_y = hR->ProjectionY( "M2gR");
	hRe_y->Scale( 0.5);

	// Subtracted
	hSe = (TH1D*) hPe_y->Clone( "EmptySubt");
	hSe->Sumw2();
	hSe->Add( hRe_y, -1.0);

	// COMPLETELY SUBTRACTED

	// Dead-time correction
	file = full;
	f_dead_f = DeadTimeSF( file);
	file = empty;
	f_dead_e = DeadTimeSF( file);

	// Scale full to empty
	factor = HistSF( hesc, hfsc, cbin, cbin)*f_dead_f/f_dead_e;

	hSf->Scale( factor);
	hS = (TH1D*) hSf->Clone( "Subtracted");
	hS->Add( hSe, -1.0);

	fhist = (TH1D*) hSf->Clone( "Full");
	ehist = (TH1D*) hSe->Clone( "Empty");
	shist = (TH1D*) hS->Clone( "Subt");
}