C++ TTree::Draw方法代码示例

void csv(TString input="tmva.csvoutput.txt", TString par1="par2", TString par2="par3", TString par3="", TString value="eventEffScaled_5") {
  std::cout << "Usage:" << std::endl
            << ".x scripts/csv.C    with default arguments" << std::endl
            << ".x scripts/csv.C(filename, par1, par2, value)" << std::endl
            << std::endl
            << "  Optional arguments:" << std::endl
            << "    filename        path to CSV file" << std::endl
            << "    par1            name of X-parameter branch" << std::endl
            << "    par2            name of Y-parameter branch (if empty, efficiency is drawn as a function of par1)" << std::endl
            << "    value           name of result (efficiency) branch" << std::endl
            << std::endl;

  TTree *tree = new TTree("data", "data");
  tree->ReadFile(input);

  gStyle->SetPalette(1);
  gStyle->SetPadRightMargin(0.14);

  TCanvas *canvas = new TCanvas("csvoutput", "CSV Output", 1200, 900);

  tree->SetMarkerStyle(kFullDotMedium);
  tree->SetMarkerColor(kRed);
  if(par2.Length() > 0) {
    //tree->Draw(Form("%s:%s", par2.Data(), par1.Data()));
    if(par3.Length() > 0)
      tree->Draw(Form("%s:%s:%s:%s", par1.Data(), par2.Data(), par3.Data(), value.Data()), "", "COLZ"); //, "", "Z");
    else
      tree->Draw(Form("%s:%s:%s", par2.Data(), par1.Data(), value.Data()), "", "COLZ"); //, "", "Z");

    TH1 *histo = tree->GetHistogram();
    if(!histo)
      return;

    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(Form("Classifier parameter %s", par2.Data()));
    if(par3.Length() > 0)
      histo->GetZaxis()->SetTitle(Form("Classifier parameter %s", par3.Data()));
    else
      histo->GetZaxis()->SetTitle("");

    if(par3.Length() == 0) {
      float x = 0;
      float y = 0;
      float val = 0;
      double maxVal = tree->GetMaximum(value);
      double minVal = tree->GetMinimum(value);

      tree->SetBranchAddress(par1, &x);
      tree->SetBranchAddress(par2, &y);
      tree->SetBranchAddress(value, &val);
      TLatex l;
      l.SetTextSize(0.03);
    
      Long64_t nentries = tree->GetEntries();
      for(Long64_t entry=0; entry < nentries; ++entry) {
        tree->GetEntry(entry);
    
        l.SetTextColor(textColor(val, maxVal, minVal));
        l.DrawLatex(x, y, Form("%.3f", val*100));
      }
    }
  }
  else {
    tree->Draw(Form("%s:%s", value.Data(), par1.Data()));
    TH1 *histo = tree->GetHistogram();
    if(!histo) 
      return;
    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(value);
  }
}

int main (int argc, char **argv)
{
  const char* chInFile = "ws.root";
  const char* chOutFile = "ws_data.root";
  const char* chRootFile = "BDT20.root";
  const char* chCut = "";

  char option_char;
  while ( (option_char = getopt(argc,argv, "i:o:r:c:")) != EOF )
    switch (option_char)
      {
         case 'i': chInFile = optarg; break;
         case 'o': chOutFile = optarg; break;
         case 'r': chRootFile = optarg; break;
         case 'c': chCut = optarg; break;
         case '?': fprintf (stderr,
                            "usage: %s [i<input file> o<output file>]\n", argv[0]);
      }

  cout << "In  Ws = " << chInFile << endl;
  cout << "Out Ws = " << chOutFile << endl;
  cout << "Data From = " << chRootFile << endl;
  cout << "Extra Cut = " << chCut << endl;


   TFile* in_file = new TFile(chInFile);
   RooWorkspace *rws = (RooWorkspace*) in_file->Get("rws");

   TFile* tree_file = new TFile(chRootFile);
   TTree* tree = (TTree*) tree_file->Get("tree");

   TFile* out_file = new TFile(chOutFile, "RECREATE");

   RooArgSet allVars(*rws->var("m"),*rws->var("t"),*rws->var("et"),*rws->var("cpsi"),*rws->var("ctheta"),*rws->var("phi"),*rws->var("d"),*rws->cat("dilution"));
   RooDataSet* data = new RooDataSet("data","data",allVars);
   RooDataSet* dataBkg = new RooDataSet("dataBkg","dataBkg",allVars);

   //TCut* cut = new TCut("5.17<bs_mass && bs_mass<5.57 && bs_pdl>-0.044 && bs_pdl<0.3 ");
   TCut* cut = new TCut("5.17<bs_mass && bs_mass<5.57 && bs_epdl<0.025 && bs_pdl<0.4 && bs_pdl>-0.44");
   *cut += chCut;
   tree->Draw(">>entry_list", *cut, "entrylist");
   TEntryList* event_list = (TEntryList*) out_file->Get("entry_list");

   Double_t dM, dT, dEt, dCpsi, dCtheta, dPhi, dd;
   Int_t ddDefined;
   tree->SetBranchAddress("bs_mass", &dM);
   tree->SetBranchAddress("bs_pdl", &dT);
   tree->SetBranchAddress("bs_epdl", &dEt);
   tree->SetBranchAddress("bs_angle_cpsi", &dCpsi);
   tree->SetBranchAddress("bs_angle_ctheta", &dCtheta);
   tree->SetBranchAddress("bs_angle_phi", &dPhi);
   tree->SetBranchAddress("newtag_ost", &dd);
   tree->SetBranchAddress("newtag_ost_defined", &ddDefined);


   for (Long_t i=0; i<event_list->GetN(); i++){
     tree->GetEntry(event_list->GetEntry(i));

       *rws->var("m")=dM;
       *rws->var("t")=dT/0.0299792458;
       *rws->var("et")=dEt/0.0299792458;
       *rws->var("cpsi")=dCpsi;
       *rws->var("ctheta")=dCtheta;
       *rws->var("phi")=dPhi;

       *rws->var("d")=0;
       rws->cat("dilution")->setIndex(0);
       if ( ddDefined==1 ){
    	   rws->cat("dilution")->setIndex(1);
    	   *rws->var("d")=dd;
       }

       data->add(allVars);
       if (dM<5.29 || dM>5.44)
           dataBkg->add(allVars);
   }

   rws->import(*data);
   rws->import(*dataBkg);
   rws->Write("rws");
   out_file->Close();
   in_file->Close();
   tree_file->Close();

   cout << endl << "Done." << endl;
}

void TMVAClassificationElecTau(std::string ordering_ = "Pt", std::string bkg_ = "qqH115vsWZttQCD") {

    TMVA::Tools::Instance();

    TString outfileName( "TMVAElecTau"+ordering_+"Ord_"+bkg_+".root" );
    TFile* outputFile = TFile::Open( outfileName, "RECREATE" );

    TMVA::Factory *factory = new TMVA::Factory( "TMVAClassificationElecTau"+ordering_+"Ord_"+bkg_, outputFile,
            "!V:!Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D" );
    factory->AddVariable( "pt1", "pT-tag1", "GeV/c"         , 'F'  );
    factory->AddVariable( "pt2", "pT-tag2", "GeV/c"         , 'F'  );
    factory->AddVariable( "Deta","|y-tag1 - y-tag2|",""     , 'F'  );
    //factory->AddVariable( "opposite:=abs(eta1*eta2)/eta1/eta2","sign1*sign2",""             , 'F'  );
    //factory->AddVariable( "Dphi", "#Delta#phi" ,""             , 'F'  );
    factory->AddVariable( "Mjj", "M(tag1,tag2)", "GeV/c^{2}"  , 'F'  );

    factory->AddSpectator( "eta1",  "#eta_{tag1}" , 'F' );
    factory->AddSpectator( "eta2",  "#eta_{tag2}" , 'F' );

    factory->SetWeightExpression( "sampleWeight" );

    TString fSignalName              = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/ElecTauStream2011/nTupleVBFH115-powheg-PUS1_Open_ElecTauStream.root";
    TString fBackgroundNameDYJets    = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/ElecTauStream2011/nTupleZjets-alpgen-PUS1_Open_ElecTauStream.root";
    TString fBackgroundNameWJets     = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/ElecTauStream2011/nTupleWJets-madgraph-PUS1_Open_ElecTauStream.root";
    TString fBackgroundNameQCD       = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/ElecTauStream2011/nTupleQCD_Open_ElecTauStream.root";
    TString fBackgroundNameTTbar     = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/ElecTauStream2011/nTupleTTJets-madgraph-PUS1_Open_ElecTauStream.root";


    TFile *fSignal(0);
    TFile *fBackgroundDYJets(0);
    TFile *fBackgroundWJets(0);
    TFile *fBackgroundQCD(0);
    TFile *fBackgroundTTbar(0);

    fSignal           = TFile::Open( fSignalName );
    fBackgroundDYJets = TFile::Open( fBackgroundNameDYJets );
    fBackgroundWJets  = TFile::Open( fBackgroundNameWJets );
    fBackgroundQCD    = TFile::Open( fBackgroundNameQCD );
    fBackgroundTTbar  = TFile::Open( fBackgroundNameTTbar );

    if(!fSignal || !fBackgroundDYJets || !fBackgroundWJets || !fBackgroundQCD || !fBackgroundTTbar) {
        std::cout << "ERROR: could not open files" << std::endl;
        exit(1);
    }

    TString tree = "outTree"+ordering_+"Ord";

    TCut mycuts = "";
    TCut mycutb = "";

    TCut cutA  = "pt1>0 && tightestHPSWP>0";
    TCut cutB  = "pt1>0 && combRelIsoLeg1<0.1";
    TCut cutBl = "pt1>0 && combRelIsoLeg1<0.3";
    TCut cutC  = "pt1>0 && diTauCharge==0";
    TCut cutD  = "pt1>0 && MtLeg1<40";

    // select events for training
    TFile* dummy = new TFile("dummy.root","RECREATE");
    TH1F* allEvents = new TH1F("allEvents","",1,-10,10);
    float totalEvents, cutEvents;

    // signal: all
    TTree *signal           = ((TTree*)(fSignal->Get(tree)))->CopyTree(cutA&&cutB&&cutC&&cutD);
    cout << "Copied signal tree with full selection: " << ((TTree*)(fSignal->Get(tree)))->GetEntries() << " --> "  << signal->GetEntries()  << endl;
    allEvents->Reset();
    signal->Draw("eta1>>allEvents","sampleWeight");
    cutEvents  = allEvents->Integral();
    Double_t signalWeight =   1.0;
    cout << "Signal: expected yield " << cutEvents << " -- weight " << signalWeight << endl;

    // Z+jets: all
    TTree *backgroundDYJets = ((TTree*)(fBackgroundDYJets->Get(tree)))->CopyTree(cutA&&cutB&&cutC&&cutD);
    cout << "Copied DYJets tree with full selection: " << ((TTree*)(fBackgroundDYJets->Get(tree)))->GetEntries() << " --> "  << backgroundDYJets->GetEntries()  << endl;
    allEvents->Reset();
    backgroundDYJets->Draw("eta1>>allEvents","sampleWeight");
    cutEvents  = allEvents->Integral();
    Double_t backgroundDYJetsWeight = 1.0;
    cout << "ZJets: expected yield " << cutEvents << " -- weight " << backgroundDYJetsWeight << endl;

    // W+jets: iso+Mt
    TTree *backgroundWJets  = ((TTree*)(fBackgroundWJets->Get(tree)))->CopyTree(cutB&&cutD);
    cout << "Copied WJets tree with iso+Mt selection: " << ((TTree*)(fBackgroundWJets->Get(tree)))->GetEntries() << " --> "  << backgroundWJets->GetEntries()  << endl;
    allEvents->Reset();
    backgroundWJets->Draw("eta1>>allEvents","sampleWeight");
    totalEvents  = allEvents->Integral();
    allEvents->Reset();
    backgroundWJets->Draw("eta1>>allEvents","sampleWeight*(tightestHPSWP>0 && diTauCharge==0)");
    cutEvents  = allEvents->Integral();
    Double_t backgroundWJetsWeight  =  cutEvents / totalEvents;
    cout << "WJets: expected yield " << cutEvents  << " -- weight " << backgroundWJetsWeight << endl;

    // QCD: Mt+loose iso
    TTree *backgroundQCD    = ((TTree*)(fBackgroundQCD->Get(tree)))->CopyTree(cutD&&cutBl);
    cout << "Copied QCD tree with Mt selection: " << ((TTree*)(fBackgroundQCD->Get(tree)))->GetEntries() << " --> "  << backgroundQCD->GetEntries()  << endl;
    allEvents->Reset();
    backgroundQCD->Draw("eta1>>allEvents","sampleWeight");
    totalEvents  = allEvents->Integral();
    allEvents->Reset();
    backgroundQCD->Draw("eta1>>allEvents","sampleWeight*(tightestHPSWP>0 && diTauCharge==0 && combRelIsoLeg1<0.1)");
    cutEvents  = allEvents->Integral();
//.........这里部分代码省略.........

  void draw_sigvsback(std::string const & cname,
		      std::string const & dr,
		      std::string const & binning,
		      std::string const & we_sig,
		      std::string const & sig_label,
		      std::string const & we_back,		     
		      std::string const & back_label,		     
		      std::string const & op = "",
		      std::string const & title = "",
		      std::string const & xtitle = "",
		      std::string const & ytitle = "",
		      bool const significance = false,
		      bool const efficiency = false,
		      bool const forwards = false) {
    
    TString hname_sig = "h_sig";
    TString draw_str_sig = "";
    draw_str_sig += dr;
    draw_str_sig += ">>";
    draw_str_sig += hname_sig;
    draw_str_sig += binning;
    
    TString hname_back = "h_back";
    TString draw_str_back = "";
    draw_str_back += dr;
    draw_str_back += ">>";
    draw_str_back += hname_back;
    draw_str_back += binning;
    
    TCanvas * canvas = new TCanvas("temp");
    if(tree2)
      tree2->Draw(draw_str_sig.Data(),
		  we_sig.c_str(),
		  op.c_str());
    else
      tree->Draw(draw_str_sig.Data(),
		 we_sig.c_str(),
		 op.c_str());
    delete canvas;
    
    canvas = new TCanvas("temp");
    tree->Draw(draw_str_back.Data(),
	       we_back.c_str(),
	       op.c_str());
    delete canvas;
    
    TH1F * hist_sig = (TH1F*)gDirectory->Get(hname_sig.Data());
    TH1F * hist_back = (TH1F*)gDirectory->Get(hname_back.Data()); 
    if(hist_sig->GetEntries() == 0) 
      std::cout << "draw_sigvsback: No signal for selection: \"" << we_sig << "\"\n";
    if(hist_back->GetEntries() == 0) 
      std::cout << "draw_sigvsback: No background for selection: \"" << we_back << "\"\n";    
    
    hist_back->SetStats(0);
    hist_back->SetTitle(title.c_str());
    hist_back->GetXaxis()->SetTitle(xtitle.c_str());
    hist_back->GetXaxis()->CenterTitle(); 
    hist_back->GetYaxis()->SetTitle(ytitle.c_str());
    hist_back->GetYaxis()->CenterTitle();
    hist_sig->SetLineColor(kRed);
    hist_back->SetLineColor(kBlue);

    TLegend * legend = new TLegend(0.6, 0.9, 0.9, 0.6);
    /*
    legend->SetHeader(("Total: "+to_string_with_precision(hist_sig->Integral()+hist_back->Integral())+" events").c_str());
    ((TLegendEntry*)legend->GetListOfPrimitives()->First())->SetTextAlign(22);
    
    legend->AddEntry(hist_sig, (sig_label+": "+to_string_with_precision(hist_sig->Integral())+" events").c_str());
    legend->AddEntry(hist_back, (back_label+": "+to_string_with_precision(hist_back->Integral())+" events").c_str());
    */

    legend->AddEntry(hist_sig, (sig_label).c_str());
    legend->AddEntry(hist_back, (back_label).c_str());

    canvas = new TCanvas(cname.c_str());   
    TPad * pad_base = new TPad();
    pad_base->Draw();
    pad_base->SetFillColor(0);
    double ymin_base = hist_sig->GetYaxis()->GetXmin();
    double ymax_base = hist_sig->GetYaxis()->GetXmax();
    double dy_base = (ymax_base-ymin_base)/0.8;
    double xmin_base = hist_sig->GetXaxis()->GetXmin();
    double xmax_base = hist_sig->GetXaxis()->GetXmax();
    double dx_base = (xmax_base-xmin_base)/0.8;
    pad_base->Range(xmin_base-0.1*dx_base,ymin_base-0.1*dy_base,xmax_base+0.1*dx_base,ymax_base+0.1*dy_base);
    pad_base->cd();
    hist_back->Draw();
    hist_sig->Draw("same");
    if(!significance && !efficiency) legend->Draw();

    //////    
    hist_sig->Scale(run_pot / signal_pot);
    hist_back->Scale(run_pot / background_pot);
    //////

    TGraph * graph_sig = nullptr;
    TPad * pad_sig = nullptr;
    if(significance) {
      graph_sig = getgraphsig(hist_sig, hist_back, forwards);
      canvas->cd();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  
  cout<<"starting canvas loop..."<<endl;
  for(int ci =0; ci<ncanvases; ci++){ //loop over canvases
    //create indicators and strings:
    int year;
    if(ci<(ncanvases/2)){//2011 for 1st half canvases, 2015 for 2nd half
      year = 2011;
    }else year = 2015;
    TString yearstring = Form("%d",year);
    TString filetype;
    if(ci%2==0){//true every other canvas
      filetype = "data";
    }else filetype = "MC";
    int file_num=-1000;//assign file_num based on combinations of data and MC
    if(year==2011){
      if(filetype=="data") file_num=0;
      if(filetype=="MC") file_num=1;
    }
    if(year==2015){
      if(filetype=="data") file_num=2;
      if(filetype=="MC") file_num=3;
    }
    int B0;//these will be assigned within the histogram loop
    TString Bdecay;
    TString branch;
    TString  tracktype = "nLongTracks";
    TString cistring = Form("%d",ci);
    cout<<"loop "<<ci<<" indicators and strings created:"<<endl<<"year: "<<year<<endl;
    cout<<"file: "<<file_num<<endl;

    //create the plots
    placeholder = "c"+cistring;
    placeholder2=yearstring+filetype+tracktype;
    c[ci] = new TCanvas(placeholder,placeholder2,1200,800); //create the canvases
    c[ci]->cd();
    gStyle->SetOptStat("");
    leg[ci] = new TLegend(0.7, 0.7, .97, .93);//create legend
    placeholder = "hs"+cistring;
    hs[ci] = new THStack(placeholder,placeholder2); //create the stack to hold the histograms
    cout<<"starting histogram loop..."<<endl;
    for(int hi=0;hi<nhpc;hi++){ //loop over histograms in a canvas; B0=0 and B0=1
      if((int)floor(hi)%2 ==0){//every other histogram
        B0=0;
        Bdecay = "B^{-}->(D^{0}->K^{-} #pi^{+})#mu^{-}";
        if((filetype=="data")&&(year==2011)){
          branch="tupleb2D0Mu/tupleb2D0Mu";
        }else branch="Tuple_b2D0MuX/DecayTree";
      } else{
        B0=1;
        Bdecay = "B^{0}->(D^{-}->K^{+} #pi^{-} #pi^{-})#mu^{+}";
        if((filetype=="data")&&(year==2011)){
          branch="tupleb2DpMu/tupleb2DpMu";
        }else branch="Tuple_b2DpMuX/DecayTree";
      }
      //create convenient strings
      TString histring = Form("%d",hi);
      TString hname = "h"+cistring+histring;
      
      //create histograms
      if(tracktype=="nTracks") h[ci][hi] = new TH1F(hname,tracktype,104,0,1144);
      if(tracktype=="nLongTracks") h[ci][hi] = new TH1F(hname,tracktype,131,0,262);
      cout<<"histogram loop "<<hi<<" strings and histograms created"<<endl;

      //navigate files
      cout<<"navigating file..."<<endl;
      TTree *MyTree;
      f[file_num]->GetObject(branch,MyTree);
      //draw histograms
      cout<<"drawing histogram "<<hi<<"..."<<endl;
      h[ci][hi]->SetLineColor(hi+1);
      placeholder = tracktype+">>"+hname;
      MyTree->Draw(placeholder);
      cout<<"stacking histogram "<<hi<<"..."<<endl;
      hs[ci]->Add(h[ci][hi]);//stack histograms
      leg[ci]->AddEntry(h[ci][hi],Bdecay,"l");//fill legend
      if(ci==0) legf->AddEntry(h[ci][hi],Bdecay,"l");//fill combined legend; all the same, so using the first one is fine
    }
    //draw stacked histograms
    cout<<"drawing stack "<<ci<<"..."<<endl;
    placeholder = yearstring+": "+filetype+": "+tracktype;
    hs[ci]->SetTitle(placeholder);
    hs[ci]->Draw("nostack");
    leg[ci]->Draw();
    cf->cd(ci+1);
    hs[ci]->Draw("nostack");//not using c[ci]->DrawClonePad(); because too many legends
    //save stuff:
    cout<<"saving files..."<<endl;
    placeholder = outputlocation+plotfilename+"(";//the closing page is added after the loop
    c[ci]->Print(placeholder);
    placeholder = outputlocation+"c"+cistring+".C";
    c[ci]->SaveAs(placeholder);
    cout<<endl;
  }
  cf->cd();
  legf->Draw();
  placeholder = outputlocation+plotfilename+")";
  cf->Print(placeholder);
  gROOT->SetBatch(kFALSE);
  cout<<"done"<<endl;
}

void balanceMetVsAj(TString infname = "dj_HCPR-J50U-hiGoodMergedTracks_OfficialSelv2_Final0_120_50.root",
                    TCut myCut = "cent<30", char *title = "",bool drawLegend = false,
                    bool drawSys = true
                   )
{
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get("ntjt");
  
   t->SetAlias("metxMerged0","metx0+metx1+metx2");
   t->SetAlias("metxMerged1","metx1+metx2");
   t->SetAlias("metxMerged2","metx2");
   t->SetAlias("metxMerged3","metx3+metx4");
   const int nBin = 4;
   double bins[nBin+1] = {0.5,1.5,4,8,1000};  
   double colors[nBin] = {38, kOrange-8,kBlue-3,kRed};

   const int nBinAj = 4;
   double ajBins[nBinAj+1] = {0.0001,0.11,0.22,0.33,0.49999};
   // Selection cut
   TCut evtCut = "nljet>120&&abs(nljetacorr)<2&&aljet>50&&abs(aljetacorr)<2&&jdphi>2./3*TMath::Pi()&&!maskEvt";

   cout << "Sel evt: " << t->GetEntries(evtCut&&myCut) << endl;

   TH1D *p[nBin];

   for (int i=0;i<nBin;i++)
   {
      TH1D *h1 = new TH1D(Form("h1%d",i),"",nBinAj,ajBins);
      TH1D *h2 = new TH1D(Form("h2%d",i),"",nBinAj,ajBins);
      h1->Sumw2();
      h2->Sumw2();
  //    t->Project(Form("h%d",i),"Aj", "1"*(evtCut&&myCut));
//      t->Project(Form("h2%d",i),"Aj", Form("((-1*metxMerged%d))",i)*(evtCut&&myCut));
      t->Draw(Form("Aj>>h1%d",i), "weight"*(evtCut&&myCut));
      t->Draw(Form("Aj>>h2%d",i), Form("((-weight*metxMerged%d))",i)*(evtCut&&myCut));
      p[i]=(TH1D*)h2->Clone();
      p[i]->SetName(Form("p%d",i));     
      p[i]->Divide(h1);
      p[i]->SetLineColor(1);     
      p[i]->SetMarkerColor(colors[i]);
      p[i]->SetFillColor(colors[i]);
//      p[i]->SetFillStyle(3004+fabs(i-1));
      p[i]->SetFillStyle(1001);
     
   }

   TH1D *pall;
   TH1D *h1 = new TH1D(Form("hAll1"),"",nBinAj,ajBins);
   TH1D *h2 = new TH1D(Form("hAll2"),"",nBinAj,ajBins);
   h1->Sumw2();
   h2->Sumw2();
   t->Draw(Form("Aj>>hAll1"), "weight"*(evtCut&&myCut));
   t->Draw(Form("Aj>>hAll2"), Form("((-weight*metx))")*(evtCut&&myCut));

   pall=(TH1D*)h2->Clone();
   pall->SetName("pall");
   pall->Divide(h1);
   pall->SetXTitle("A_{J}");
   pall->SetYTitle("<#slash{p}_{T}^{#parallel}> (GeV/c)");
   pall->GetXaxis()->CenterTitle();
   pall->GetYaxis()->CenterTitle();
   pall->GetXaxis()->SetLabelSize(22);
   pall->GetXaxis()->SetLabelFont(43);
   pall->GetXaxis()->SetTitleSize(24);
   pall->GetXaxis()->SetTitleFont(43);
   pall->GetYaxis()->SetLabelSize(22);
   pall->GetYaxis()->SetLabelFont(43);
   pall->GetYaxis()->SetTitleSize(24);
   pall->GetYaxis()->SetTitleFont(43);
   pall->GetXaxis()->SetTitleOffset(1.8);
   pall->GetYaxis()->SetTitleOffset(2.4);



   pall->SetNdivisions(505);
   pall->SetAxisRange(-59.9,59.9,"Y");
   pall->SetMarkerSize(1);
   pall->Draw("E");
   for (int i=0;i<nBin;++i) {
      p[i]->SetLineWidth(1);
      p[i]->Draw("hist same");
   }
   pall->Draw("E same");

   if (drawSys == 1) {
      for(int i = 0; i < nBinAj; ++i){
        double x = pall->GetBinCenter(i+1);
        double y = pall->GetBinContent(i+1);
        // Quote the difference between GEN and RECO in >8 Bin (20%) before adjusting eff as systematics
        double err = -p[nBin-1]->GetBinContent(i+1)*0.2;
 
        DrawTick(y,err,err,x,1,0.02,1);
      }
   }

   // Legend
   TLegend *leg = new TLegend(0.10,0.68,0.70,0.96);
   leg->SetFillStyle(0);
   leg->SetBorderSize(0);
   leg->SetTextFont(63);
//.........这里部分代码省略.........

void hfCentrality(char* fname = "r151878.root",  char* trg="anaPixelHitJet50U")
{
   
   TFile *f1=new TFile(Form("/d101/kimy/macro/pixelTrees/rootFiles/%s",fname),"r");
   TTree *trPixMB = (TTree*)f1->Get("anaPixelHitMB/PixelTree");
   TTree *trPixJet = (TTree*)f1->Get("anaPixelHitJet50U/PixelTree");
   
   //   const int nscEtBin = 10;
   //   double scetBin[nscEtBin+1] = { 15.1,18,21,24,28,32,36,40,48,55,70};
   
   TH1D* hhfMB  = new TH1D(Form("hhfMB",trg),";Sum HF Energy (TeV);Fraction of minimum bias events",70,0,160);
   TH1D* hhfMB1  = (TH1D*)hhfMB->Clone("hhfMB1");
   TH1D* hhfMB2  = (TH1D*)hhfMB->Clone("hhfMB2");
   TH1D* hhfMB3  = (TH1D*)hhfMB->Clone("hhfMB3");
   TH1D* hhfJET  = (TH1D*)hhfMB->Clone("hhfJET");

   TCanvas* c1 = new TCanvas(Form("c1_%s",trg),"",400,400);

   trPixMB->Draw("hf/1000.>>hhfMB");
   trPixMB->Draw("hf/1000.>>hhfMB1","cBin<=4");
   trPixMB->Draw("hf/1000.>>hhfMB2","cBin>=4 && cBin<=12");
   trPixMB->Draw("hf/1000.>>hhfMB3","cBin>=12 && cBin<=35");
   
   trPixJet->Draw("hf/1000.>>hhfJET");

   
   handsomeTH1(hhfMB1,1);
   handsomeTH1(hhfMB2,2);
   handsomeTH1(hhfMB3,4);
   hhfMB1->SetFillColor(1);
   hhfMB2->SetFillColor(2);
   hhfMB3->SetFillColor(4);

   
   hhfMB->Scale(1./hhfMB->GetEntries());
   hhfMB->SetAxisRange(1e-5,1,"Y");
   hhfJET->Scale(1./hhfJET->GetEntries());
   hhfJET->SetAxisRange(1e-5,1,"Y");

   hhfMB->SetLineWidth(2);
   handsomeTH1(hhfMB);
   hhfMB->DrawCopy();
   
   //   TLine* t1 = new TLine(.511,1e-5,.511,hhfMB->GetBinContent(hhfMB->FindBin(0.511)));
   TLine* t1 = new TLine(0,1e-5,.511,hhfMB->GetBinContent(hhfMB->FindBin(0.511)));
   TLine* t2 = new TLine(35.397,1e-5,35.397,hhfMB->GetBinContent(hhfMB->FindBin(35.397)));
   TLine* t3 = new TLine(79.370,1e-5,79.370,hhfMB->GetBinContent(hhfMB->FindBin(79.370)));
   t1->SetLineWidth(1);
   t2->SetLineWidth(1);
   t3->SetLineWidth(1);
   t1->SetLineStyle(7);
   t2->SetLineStyle(7);
   t3->SetLineStyle(7);
   t1->Draw();
   t2->Draw();
   t3->Draw();
   
   //  hhfMB3->DrawCopy();
   //  hhfMB2->DrawCopy("same");
   //  hhfMB1->DrawCopy("same");
   gPad->SetLogy();

   TLegend* leg0 = new TLegend(0.2813131,0.7115054,0.9,0.8439785,NULL,"brNDC");
   TLegend* leg0b =new TLegend(0.3813131,0.7115054,1,0.8439785,NULL,"brNDC");
   TLegend* leg1 = new TLegend(0.1767677,0.3826344,0.4467677,0.5875591,NULL,"brNDC");
   TLegend* leg2 = new TLegend(0.3611111,0.3876344,0.6111111,0.5865591,NULL,"brNDC");
   TLegend* leg3 = new TLegend(0.5606061,0.3876344,0.8106061,0.5865591,NULL,"brNDC");

   easyLeg(leg0b,"Centrality, HLT_HIJet50U");

   easyLeg(leg1,"30%-100%");
   easyLeg(leg2,"10%-30% ");
   easyLeg(leg3," 0%-10% ");

   //   TCanvas* c2 = new TCanvas(Form("c2_%s",trg),"",400,400);
   
   float triggerRatio = 0.29 / 49.;
   
   hhfJET->SetLineWidth(2);
   hhfJET->SetLineColor(2);
   hhfJET->SetFillColor(2);
   hhfJET->SetFillStyle(3544);
   hhfJET->Scale(triggerRatio);
   hhfJET->DrawCopy("same");
   
   easyLeg(leg0,"Centrality");
   leg0->AddEntry(hhfMB,"HLT_MinBiasHForBSC_Core","l");
   leg0->AddEntry(hhfJET,"HLT_HiJet50U","l");
   leg0->Draw();
   //   leg1->Draw();
   //   leg2->Draw();
   //   leg3->Draw();
     
   TLatex *bint = new TLatex(0.1867677,0.4876344,"30%-100%");
   bint->SetTextFont(63);
   bint->SetTextSize(13);
   bint->SetNDC();
   bint->Draw();

   bint = new TLatex(0.4011111,0.4876344,"10%-30% ");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
				totalCutBkgPbp[in][ir][ipt] = trigCut && recoCut && massRangeBkg && ctauRange && accRecoPlus && accRecoMinus && rapRangePbp[in][ir][ipt] && ptRange[in][ir][ipt];
				totalCutBkgpPb[in][ir][ipt] = trigCut && recoCut && massRangeBkg && ctauRange && accRecoPlus && accRecoMinus && rapRangepPb[in][ir][ipt] && ptRange[in][ir][ipt];
			}
		}
	}


	/////////////////////////////////////////////////////////////	
	/////////////////////////////////////////////////////////////	
	// define 1D hist
	TH1D *hMeanPt[nEt][nRap][nPt];
	TH1D *hMeanPtBkg[nEt][nRap][nPt];
	double meanVal[nEt][nRap][nPt];
	double maxVal[nEt][nRap][nPt];
	double meanValBkg[nEt][nRap][nPt];
	double maxValBkg[nEt][nRap][nPt];

	//TCanvas* c1 = new TCanvas("c1","c1",600,600);
	TCanvas* cmulti[nEt][nRap];
	TLegend *legUR = new TLegend(0.58,0.68,0.90,0.90,NULL,"brNDC");
	TLegend *legUM = new TLegend(0.30,0.68,0.65,0.90,NULL,"brNDC");
	TLegend *legUL = new TLegend(0.17,0.68,0.51,0.90,NULL,"brNDC");
	TLegend *legBM = new TLegend(0.30,0.20,0.65,0.42,NULL,"brNDC");
	SetLegendStyle(legUR);
	SetLegendStyle(legUM);
	SetLegendStyle(legUL);
	SetLegendStyle(legBM);
	TLatex* latex = new TLatex();
	latex->SetNDC();
	latex->SetTextAlign(12);
	latex->SetTextSize(0.04);

	for (int in=0; in<nEt; in++){
		for (int ir=0; ir<nRap; ir++ ) {
	 		cmulti[in][ir]= new TCanvas(Form("cmulti_%d_%d",in,ir),Form("%s",rapstrArr[ir].c_str()),1500,600);
			cmulti[in][ir]->Divide(5,2);
			for (int ipt=0; ipt<nPt; ipt++) {
				hMeanPt[in][ir][ipt] = new TH1D(Form("hMeanPt_%d_%d_%d",in, ir, ipt),";p_{T} (GeV/c);events", nbin, ptArr[ipt], ptArr[ipt+1]);
				hMeanPt[in][ir][ipt] ->Sumw2();
				SetHistStyle(hMeanPt[in][ir][ipt],4,0);
				hMeanPtBkg[in][ir][ipt] = new TH1D(Form("hMeanPtBkg_%d_%d_%d",in, ir, ipt),";p_{T} (GeV/c);events", nbin, ptArr[ipt], ptArr[ipt+1]);
				hMeanPtBkg[in][ir][ipt] ->Sumw2();
				SetHistStyle(hMeanPtBkg[in][ir][ipt],3,10);
				cmulti[in][ir]->cd(ipt+1);
				gPad->SetLogy(1);
/*				
				treeDataPbp1->Draw(Form("Reco_QQ_4mom.Pt()>>%s",hMeanPt[in][ir][ipt]->GetName()),(totalCutPbp[in][ir][ipt]&&(!runCut)),"",nEntry);
				treeDataPbp2->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPt[in][ir][ipt]->GetName()),(totalCutPbp[in][ir][ipt]&&runCut),"",nEntry);
				treeDatapPb->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPt[in][ir][ipt]->GetName()),totalCutpPb[in][ir][ipt],"",nEntry);
				treeDataPbp1->Draw(Form("Reco_QQ_4mom.Pt()>>%s",hMeanPtBkg[in][ir][ipt]->GetName()),(totalCutBkgPbp[in][ir][ipt]&&(!runCut)),"pe same",nEntry);
				treeDataPbp2->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPtBkg[in][ir][ipt]->GetName()),(totalCutBkgPbp[in][ir][ipt]&&runCut),"pe same",nEntry);
				treeDatapPb->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPtBkg[in][ir][ipt]->GetName()),totalCutBkgpPb[in][ir][ipt],"pe same",nEntry);
*/

				treeDataPbp1->Draw(Form("Reco_QQ_4mom.Pt()>>%s",hMeanPt[in][ir][ipt]->GetName()),(totalCutPbp[in][ir][ipt]&&(!runCut)),"");
				treeDataPbp2->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPt[in][ir][ipt]->GetName()),(totalCutPbp[in][ir][ipt]&&runCut),"");
				treeDatapPb->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPt[in][ir][ipt]->GetName()),totalCutpPb[in][ir][ipt],"");
				treeDataPbp1->Draw(Form("Reco_QQ_4mom.Pt()>>%s",hMeanPtBkg[in][ir][ipt]->GetName()),(totalCutBkgPbp[in][ir][ipt]&&(!runCut)),"pe same");
				treeDataPbp2->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPtBkg[in][ir][ipt]->GetName()),(totalCutBkgPbp[in][ir][ipt]&&runCut),"pe same");
				treeDatapPb->Draw(Form("Reco_QQ_4mom.Pt()>>+%s",hMeanPtBkg[in][ir][ipt]->GetName()),totalCutBkgpPb[in][ir][ipt],"pe same");

				meanVal[in][ir][ipt] =  hMeanPt[in][ir][ipt]->GetMean(1);
				maxVal[in][ir][ipt] = 50*hMeanPt[in][ir][ipt]->GetMaximum();
				meanValBkg[in][ir][ipt] =  hMeanPtBkg[in][ir][ipt]->GetMean(1);
				maxValBkg[in][ir][ipt] = 50*hMeanPtBkg[in][ir][ipt]->GetMaximum();
				hMeanPt[in][ir][ipt]->GetYaxis()->SetRangeUser(0.1,maxVal[in][ir][ipt]);
				hMeanPtBkg[in][ir][ipt]->GetYaxis()->SetRangeUser(0.1,maxValBkg[in][ir][ipt]);
				std::cout << hMeanPt[in][ir][ipt]->GetName() <<" : entries = " << hMeanPt[in][ir][ipt]->GetEntries() << std::endl;
				cout << " Cand <p_{T}> = " << meanVal[in][ir][ipt] << endl;			
				cout << " Bkg <p_{T}> = " << meanValBkg[in][ir][ipt] << endl;			
				latex->DrawLatex(0.54, 0.87, rapstrArr[ir].c_str());
				latex->DrawLatex(0.54, 0.80, Form("Cand <p_{T}> = %.2f (GeV/c)",meanVal[in][ir][ipt]));
				latex->DrawLatex(0.54, 0.72, Form("Bkg <p_{T}> = %.2f (GeV/c)",meanValBkg[in][ir][ipt]));
				//c1->Update();
				//c1->SaveAs(Form("%s.png",hMeanPt[in][ir][ipt]->GetName()));
				//c1->Clear();
			}
			cmulti[in][ir]->Update();
			cmulti[in][ir]->SaveAs(Form("dir_meanPt/meanPt_%s_%d_%d_wBkg.pdf",strBin,in,ir));
			//cmulti[in][ir]->Clear();
		}
	}
	
	/// Save as a root file
	TFile *outFile = new TFile(Form("meanPt_%s_wBkg.root",strBin),"RECREATE");
	std::cout << "strBin : " << strBin << std::endl;
	outFile->cd();
	for (int in=0; in<nEt; in++){
		for (int ir=0; ir<nRap; ir++ ) {
			for (int ipt=0; ipt<nPt; ipt++) {
				hMeanPt[in][ir][ipt]->Write();
				hMeanPtBkg[in][ir][ipt]->Write();
			}
		}
	}
	outFile->Close();

	return 0;	

}

//.........这里部分代码省略.........
      if (category == "wGwoGSF") CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf<0.5";
      if (category == "wGwGSFwoPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput<-0.1";
      if (category == "wGwGSFwPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput>-0.1";
    }

    if(discriminator == "-AntiEMed"){
      if (category == "NoEleMatch") CategorySelection = "Tau_GsfEleMatch<0.5"; 
      if (category == "woG") CategorySelection = "Tau_NumGammaCands<0.5"; 
      if (category == "wGwoGSF") CategorySelection = "Tau_NumGammaCands>0.5 && (Tau_HasGsf<0.5 || (Tau_HasGsf>0.5 && Elec_PFMvaOutput>-0.1))";
      if (category == "wGwGSFwoPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput<-0.1";
    }

  TCut ElecSelection = CategorySelection && PUSelection && ElecPtSelection && ElecAbsEtaSelection && ElecMatchSelection ;
  TCut TauSelection = CategorySelection && PUSelection && TauPtSelection && TauAbsEtaSelection && TauMatchSelection ;
  TCut Selection;
  if (variable.find("Elec")!=std::string::npos)Selection = ElecSelection;
  if (variable.find("Tau")!=std::string::npos)Selection = TauSelection;
  

  TH1F* hVariable   = new TH1F( "hVariable" ,"" , nBins ,xMin, xMax);
  hVariable->SetXTitle(Form("%s",variable.data()));

  if (matching->find("EleMatch")!=std::string::npos){
//     hVariable->SetFillColor(kRed);
//     hVariable->SetFillStyle(3345);
    hVariable->SetLineColor(kRed);
    hVariable->SetLineWidth(2);
  }
  if (matching->find("HadMatch")!=std::string::npos){
//     hVariable->SetFillColor(kBlue);
//     hVariable->SetFillStyle(3354);
    hVariable->SetLineColor(kBlue);
    hVariable->SetLineWidth(2);
  }  
  inputTree->Draw(Form("%s>>hVariable",variable.data()));

  cout<<"Variable plotted : "<<variable<<endl;
  cout<<"Matching applied : "<<matching->data()<<endl;
  cout<<"  Total number of Candidates : "<<hVariable->GetEntries()<<endl;
  inputTree->Draw(Form("%s>>hVariable",variable.data()),Selection);
  cout<<"  Number of Cantidates after selection: "<<hVariable->GetEntries()<<endl;
  hVariable->Scale(1./hVariable->Integral());
  leg->AddEntry(hVariable,Form("%s",matching->data()));

  histograms.push_back(hVariable);
  c1->Clear();
  }
//   double yMin = +1.e+6;
//   double yMax = -1.e+6;
  TH1* refHistogram = histograms.front();
  refHistogram->SetStats(false);
  refHistogram->SetTitle("");
//   refHistogram->SetMinimum(yMin);
//   refHistogram->SetMaximum(yMax);


  if (xAxisTitle == "HoHplusE" ) {
    refHistogram->SetMaximum(1.0);
    refHistogram->SetMinimum(0.01);
    c1->SetLogy();
  }

  if(xAxisTitle == "E_{#gamma}/(P_{in}-P_{out})" ){
    refHistogram->SetMaximum(0.03);
    refHistogram->SetMinimum(0.0);
  }

  if(xAxisTitle == "HadrMva(#tau)" ){
    refHistogram->SetMaximum(0.25);
    refHistogram->SetMinimum(0.0);
  }

  TAxis* xAxis = refHistogram->GetXaxis();
  xAxis->SetTitle(xAxisTitle.Data());
  xAxis->SetTitleOffset(1.15);
  //if(variable.find("AbsEta")!=std::string::npos)xAxis->SetLimits(AbsEtaMin, AbsEtaMax);
  TAxis* yAxis = refHistogram->GetYaxis();
  yAxis->SetTitle(yAxisTitle.Data());
  yAxis->SetTitleOffset(1.30);

  int numHistograms = histograms.size();
  float YMax = 0;
  for ( int iHistogram = 0; iHistogram < numHistograms; ++iHistogram ) {
    TH1* histogram = histograms[iHistogram];
    if(histogram->GetMaximum()>YMax) YMax = histogram->GetMaximum();
  }
  for ( int iHistogram = 0; iHistogram < numHistograms; ++iHistogram ) {
    TH1* histogram = histograms[iHistogram];
    yAxis->SetRangeUser(0.,YMax+0.10*YMax);
    std::string drawOption = "hist";
    if ( iHistogram > 0 ) drawOption.append("same");
    histogram->Draw(drawOption.data());
    leg->Draw();

  }//loop matchings
  string outputName = Form("plots/plotVariablesAntiEMVA/%s/plotVariablesAntiEMVA_v4_%s_%s_%s",category.Data(),category.Data(),variable.data(),Region.Data());
  c1->Print(std::string(outputName).append(".png").data());
  c1->Print(std::string(outputName).append(".pdf").data());

}

void test_msv(std::string var="m_sv",int nbins=25, double xmin=0, double xmax=250,std::string xtitle="m_sv", std::string ytitle="Events")
{
  SetStyle(); gStyle->SetLineStyleString(11,"20 10");
  TH1::SetDefaultSumw2(1);

  std::string dir = "/afs/cern.ch/work/a/arapyan/public/svfitsamples/";
  double sigscale = 10;
  double sigscale1 = 10; 
  std::stringstream scale; scale << sigscale;
  std::stringstream scale1; scale1 << sigscale1;

  //Cut definitions
  double luminosity = 3000;
  std::stringstream lumi; lumi << luminosity;
  std::string objcut = "(ptTau1>30 && ptTau2>30 && abs(etaTau1) <4.0 && abs(etaTau1)<4.0 )";
  //std::string jetcut = objcut+"*(ptJet1>30 && ptJet2>30 && abs(etaJet1) <4.7 && abs(etaJet2) <4.7 )";
  std::string mthcut = objcut+"*( ( (tauCat1==3 && tauCat2==2) || (tauCat2==3 && tauCat1==2) ) )";

  //signal region
  //std::string vbfcut = ththcut+"*eventWeight*(eventType==0)*"+lumi.str();
  //std::string zttcut = ththcut+"*eventWeight*(eventType==1)*"+lumi.str();
  //std::string ttbarcut = ththcut+"*eventWeight*(eventType==3)*"+lumi.str();
  //std::string ewkcut = ththcut+"*eventWeight*(eventType==2)*"+lumi.str();
  //std::string othercut = ththcut+"*eventWeight*(eventType==4 || eventType==2)*"+lumi.str();
  
  //--------------------------------------------------------------------------
  
  //Get the trees
  TTree *tree = load(dir+"hh.root"); 
  TTree *tree2 = load(dir+"Bjj-vbf.root"); 
  TTree *tree3 = load(dir+"vbf_bgd.root"); 

  //-------------------------------------------------------------------------
  
  //Get histograms
  TCanvas *canv0 = MakeCanvas("canv", "histograms", 600, 600);
  canv0->cd();
  std::string vardraw;
  TH1F *vbf = new TH1F("VBFH","",nbins,xmin,xmax);
  vardraw = var+">>"+"VBFH";
  tree->Draw(vardraw.c_str(),mthcut.c_str());
  InitSignal(vbf);
  vbf->SetLineColor(kBlack);
  //TH1F *ttbar = new TH1F("TTbar","",nbins,xmin,xmax);
  //vardraw = var+">>"+"TTbar";
  //tree->Draw(vardraw.c_str(),ttbarcut.c_str());
  //InitHist(ttbar, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(155,152,204), 1001);
  TH1F *ztt = new TH1F("Ztt","",nbins,xmin,xmax);
  vardraw = var+">>"+"Ztt";
  tree2->Draw(vardraw.c_str(),mthcut.c_str());
  InitHist(ztt, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(248,206,104), 1001);
  //TH1F *ewk = new TH1F("Ewk","",nbins,xmin,xmax);
  //vardraw = var+">>"+"Ewk";
  //tree->Draw(vardraw.c_str(),ewkcut.c_str());
  //InitHist(ewk, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  TH1F *other = new TH1F("Other","",nbins,xmin,xmax);
  vardraw = var+">>"+"Other";
  tree3->Draw(vardraw.c_str(),mthcut.c_str());
  InitHist(other, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  
  delete canv0;

  //----------------------------------------------------------------------------
  //Print out the yields
  /*  Double_t error=0.0;
  ofstream outfile;
  outfile.open("yields_test.txt");
  outfile << "Yields for the signal region." << std::endl;
  outfile << "VBF   "  << vbf->IntegralAndError(0,vbf->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "TTbar   "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(0,ztt->GetNbinsX(),error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "other   "  << other->IntegralAndError(0,other->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "S/sqrt(B)    "  << vbf->Integral()/(other->Integral()+ztt->Integral()+ttbar->Integral()) << endl;
  //--------------------------------------------------------------------------
  //continue outputing
  //outfile << "Ewk total    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << endl << endl << endl;
  outfile << "In the signal region (100,150GeV)  " <<endl;
  outfile << "VBF   "  << vbf->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "TTbar    "  << ttbar->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(5,11,error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "other    "  << other->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile.close();*/
  //-----------------------------------------------------------------------
  //Draw the histograms
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);
  canv->cd();
  //ewk->Add(other); ttbar->Add(ewk); 
  //ztt->Add(ttbar); vbf->Add(ztt);
  other->Add(ztt);
  //Error band stat
  TH1F* errorBand = (TH1F*)vbf ->Clone("errorBand");
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(13);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  //  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
  //     if(errorBand->GetBinContent(idx)>0){
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    histXaxis[5] = (char*) "N_{jets}";
    histbins[5] = 5;
    histmin[5] = 4;
    histmax[5] = 9;

    histnames[6] = (char*) "BDTG_ge4je2t";
    histXaxis[6] = (char*) "BDT output";
    histbins[6] = 12;
    histmin[6] = -0.95;
    histmax[6] = 0.85;

  }
  else std::cout << "Error2, wrong category name" << std::endl;
  ////
  TString cuts = "(oppositeLepCharge == 1) && (dR_leplep > 0.2) && (mass_leplep > 12) && isCleanEvent && PassZmask==1 && ";
  cuts += jettagcut;

  //// sample info
  TString htitle[21] = {"_singlet_s","_singlet_tW","_singlet_t","_singletbar_s","_singletbar_tW","_singletbar_t","_ttbarW","_ttbarZ","_ttbar_cc","_ttbar_bb","_ttbar_b","_ttbar","_wjets","_zjets","_zjets_lowmass","_ww","_wz","_zz","_MuEG","_DoubleElectron","_DoubleMu",};
//     Float_t lumi = 19450;

  ////////
  //SKIPPING HISTOGRAM INITIALIZATION FOR OTHER 49 VARIABLES (USE ABOVE THREE EXAMPLES FOR YOUR USE)
  ////////

  TDirectory *currentDir = gDirectory; 

  for (int i = 0; i < numhists; i++){

    std::cout << "-->Figuring out binning for " << histnames[i] << std::endl;

    //Set the binning on the x-axis by looking at the data and combining bins (if necessary) to make sure there are no poorly populated bins
    //Define poorly populated arbitrarily as N < 100
    currentDir->cd(); //This is because the "Draw" command only finds
    TString histName = Form("tempData_%s",histnames[i]);
    TH1 * tempHist = new TH1D(histName,"",histbins[i],histmin[i],histmax[i]);

    //Only look at data
    for (int isam=18; isam<21; isam++){
      TString sample = htitle[isam];
      TString fileName = "/afs/crc.nd.edu/user/w/wluo1/LHCP_2013/CMSSW_5_3_8_patch1/src/BEAN/DrawPlots/bin/treeFiles/dilSummaryTrees" +sample + "_2012_53x_July5th_all.root";
      std::cout << "  -->Including " << fileName << std::endl;

      TFile * tmpfile = new TFile(fileName);
      TTree * tmpTree = (TTree*)tmpfile->Get("summaryTree");

      currentDir->cd(); //This is because the "Draw" command only finds

      
      TString selection =  "(" + cuts + "&&";
      if (isam == 18) selection+= "MuonEle && isMuEGTriggerPass";
      else if (isam == 19) selection += "TwoEle && isDoubleElectronTriggerPass";
      else if (isam == 20) selection += "TwoMuon && isDoubleMuTriggerPass";
      selection += ")";

      TString var = string(histnames[i])+">>+"+histName ;
      tmpTree->Draw(var, selection, "goff");


      tmpfile->Close();
      delete tmpfile;
    }

    std::cout << "tempHist name = " << tempHist->GetName() << std::endl;
    std::cout << "tempHist: Entries = " << tempHist->GetEntries() << ", Integral = " << tempHist->Integral() << std::endl;
    double threshold = (tempHist->Integral())/20;

// draw the same thing but after reco
void genPlots02(std::string fullPath, int nOverlay = 500, bool custBinning = false)
{
    const int fVerbose(1);
    setTDRStyle();
    gStyle->SetOptStat(112211);
    gStyle->SetPalette(1);
    // Canvas
    c = new TCanvas("c2","c2",1000,600);
    const unsigned int nPadX = 1;
    const unsigned int nPadY = 1;
    c->Divide(nPadX,nPadY);
    const unsigned int nPads=nPadX*nPadY;
    for(unsigned int i=1; i<=nPads; i++)
    {
	TPad* pad= (TPad*)c->cd(i);
	pad->SetTopMargin(0.10);
	pad->SetRightMargin(0.20);
	pad->SetLeftMargin(0.15);
    }
    // Open file
    TFile *f = TFile::Open(fullPath.c_str());
    if (f==0)
    {
	cout << "File " << fullPath << " not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0)
	cout << "Succesfully opened file " << fullPath << endl;
    // Get TTree
    TTree* t = (TTree*) f->Get("events");
    if(fVerbose>0) cout << "Got TTree with " << t->GetEntries() << " entries" << endl;
    // Do a cut, if needed
    //t->Draw(">>lst","chi2lb>.1&&mlb>5.61&&mlb<5.63");
    //t->Draw(">>lst","chi2lb>.1&&isSig==1");
    t->Draw(">>lst","(rid1m&4)==4&&(rid2m&4)==4&&mjp>2.895&&mjp<3.295&&prob1m>0.1&&prob2m>0.1&&ptjp>2&&probjp>0.005&&ml0>1.101&&ml0<1.129&&probpr>0.02&&probpi>0.02&&rptpr>rptpi&&ptl0>3&&rptpr>1&&rptpi>0.5&&probl0>0.02&&alphal0<0.3&&d3l0>1&&d3l0/d3El0>10&&problb>0.001&&alphalb<0.3");
    TEventList *lst;
    lst = (TEventList*)gDirectory->Get("lst");
    t->SetEventList(lst);
    if(fVerbose>0) cout << "Got TTree with " << t->GetEntries() << " entries" << endl;

    // Do plots
    c->cd(1);
    //doPlot2d(t,"hrzL0vtx", "vrl0:TMath::Abs(vzl0)",30,0,300,30,0,120,"Tit","|z|","r","cm","cm");
    if (custBinning)
    {
	double newbinsX[]={0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50};
	const int newbinsX_size = sizeof(newbinsX)/sizeof(double);
	std::vector<double> binvecX(newbinsX,newbinsX+newbinsX_size);
	//double newbinsY[]={0,1,2,3,4,5,6,7,8,9,10};
	double newbinsY[]={0,0.5,1,2,4,8,16,32};
	const int newbinsY_size = sizeof(newbinsY)/sizeof(double);
	std::vector<double> binvecY(newbinsY,newbinsY+newbinsY_size);
	doPlot2d(t,"hrzL0vtxreco", "vrl0:TMath::Abs(vzl0)",binvecX, binvecY,"#Lambda vertices","|z|","r","cm","cm");
    }
    else
    {
	doPlot2d(t,"hrzL0vtxreco", "vrl0:TMath::Abs(vzl0)",30,0,50,30,0,30,"#Lambda vertices","|z|","r","cm","cm");
    }

    // add tracker
    TPad* pad;
    pad = (TPad*)c->cd(1);
    pad->Modified();
    pad->Update();
    repositionPalette("hrzL0vtxreco");
    pad->Update();
    pad->SetLogz();
    drawTracker(pad);

    if (nOverlay<=0) return;
    int maxN = nOverlay;
    if (maxN > t->GetEntries()) maxN = t->GetEntries();
    double vrl0,vzl0,ppr,ppi,etapr,etapi;
    t->SetBranchAddress("vrl0",&vrl0);
    t->SetBranchAddress("vzl0",&vzl0);
    t->SetBranchAddress("ppr",&ppr);
    t->SetBranchAddress("etapr",&etapr);
    t->SetBranchAddress("ppi",&ppi);
    t->SetBranchAddress("etapi",&etapi);
    double scalepr = 4;
    double scalepi = 8;

    { // reference indicator
	const double x1pr = 0; const double y1pr = -3;
	const double x2pr = scalepr; const double y2pr = y1pr;
	const double versatz = 14;
	const double x1pi = x1pr+versatz; const double y1pi = -3;
	const double x2pi = x2pr+versatz+scalepi; const double y2pi = y1pi;
	TArrow *a;
	a = new TArrow(x1pr,y1pr,x2pr,y2pr,.01,">");
	a->SetLineColor(24);
	a->Draw();
	TLatex tl;
	tl.SetTextSize(20);
	tl.SetTextFont(4);
	tl.DrawLatex(x1pr,y2pr-1.2,"p(p) / 1 GeV");
	a = new TArrow(x1pi,y1pi,x2pi,y2pi,.01,">");
	a->SetLineColor(20);
	a->Draw();
//.........这里部分代码省略.........

void plotDeltaPhi(int cbin,
		  TString infname,
		  TString pythia,
		  TString mix,
		  bool useWeight,
		  bool drawXLabel,
		  bool drawLeg)
{

  TString cut="pt1>100 && pt2>40";
  TString cutpp="pt1>100 && et2>40";
  TString cstring = "";
  if(cbin==0) {
    cstring = "0-10%";
    cut+=" && bin>=0 && bin<4";
    cut+=" && ((pt1-pt2)/(pt1+pt2) < 0.3)";
  } else if (cbin==1) {
    cstring = "10-30%";
    cut+=" && bin>=4 && bin<12";
    cut+=" && ((et1-et2)/(et1+et2) < 0.3)";
  } else {
    cstring = "30-100%";
    cut+=" && bin>=12 && bin<40";
    cut+=" && ((et1-et2)/(et1+et2) < 0.3)";
  }

  // open the data file
  TFile *inf = new TFile(infname.Data());
  TTree *nt =(TTree*)inf->FindObjectAny("nt");

  // open the pythia (MC) file
  TFile *infPythia = new TFile(pythia.Data());
  TTree *ntPythia = (TTree*) infPythia->FindObjectAny("nt");

  // open the datamix file
  TFile *infMix = new TFile(mix.Data());
  TTree *ntMix =(TTree*)infMix->FindObjectAny("nt");

  // projection histogram
  TH1D *h = new TH1D("h","",20,-2.4,2.4);
  TH1D *hPythia = new TH1D("hPythia","",20,-2.4,2.4);
  TH1D *hDataMix = new TH1D("hDataMix","",20,-2.4,2.4);
  
  nt->Draw("eta1>>h",Form("(%s)",cut.Data())); 
   
  if (useWeight) {
    // use the weight value caluculated by Matt's analysis macro
    ntPythia->Draw("eta1>>hPythia",Form("(%s)",cutpp.Data())); 
    ntMix->Draw("eta1>>hDataMix",Form("(%s)*weight",cut.Data())); 
  } else {
    // ignore centrality reweighting
     ntPythia->Draw("eta1>>hPythia",Form("(%s)",cutpp.Data()));
     ntMix->Draw("eta1>>hDataMix",Form("(%s)",cut.Data()));  
  }
  
  // calculate the statistical error and normalize
  h->Sumw2();
  h->Scale(1./h->GetEntries());
  h->SetMarkerStyle(20);

  hPythia->Scale(1./hPythia->Integral(0,20));
  hPythia->SetLineColor(kBlue);
  hPythia->SetFillColor(kAzure-8);
  hPythia->SetFillStyle(3005);
   
  hPythia->SetStats(0);
  hPythia->Draw("hist");

  if(drawXLabel) hPythia->SetXTitle("Leading jet #eta");

  hPythia->GetXaxis()->SetLabelSize(20);
  hPythia->GetXaxis()->SetLabelFont(43);
  hPythia->GetXaxis()->SetTitleSize(22);
  hPythia->GetXaxis()->SetTitleFont(43);
  hPythia->GetXaxis()->SetTitleOffset(1.5);
  hPythia->GetXaxis()->CenterTitle();

  //hPythia->GetXaxis()->SetNdivisions(905,true);

  hPythia->SetYTitle("Event Fraction");

  hPythia->GetYaxis()->SetLabelSize(20);
  hPythia->GetYaxis()->SetLabelFont(43);
  hPythia->GetYaxis()->SetTitleSize(20);
  hPythia->GetYaxis()->SetTitleFont(43);
  hPythia->GetYaxis()->SetTitleOffset(2.5);
  hPythia->GetYaxis()->CenterTitle();

  hPythia->SetAxisRange(9E-4,5.9,"Y");
  hDataMix->SetAxisRange(9E-4,5.9,"Y");
  h->SetAxisRange(9E-4,1.9,"Y");

  h->Draw("same");

  hDataMix->Scale(1./hDataMix->Integral(0,20));
  hDataMix->SetLineColor(kRed);
  hDataMix->SetFillColor(kRed-9);
  hDataMix->SetFillStyle(3004);
  hDataMix->Draw("same");

//.........这里部分代码省略.........

void MakePID(){

  FILE* aoq_in = fopen("cutParameters/PID/Sn132/peak_aoq_in.txt","w");
  FILE* aoq_out = fopen("cutParameters/PID/Sn132/peak_aoq_out.txt","w");
  FILE* zet_in = fopen("cutParameters/PID/Sn132/peak_zet_in.txt","w");
  FILE* zet_out = fopen("cutParameters/PID/Sn132/peak_zet_out.txt","w");

  TH1F* h1_AOQ[64][2];
  TH1F* h1_ZET[64][2];
  TCanvas* can_AOQ[4][2];
  TCanvas* can_ZET[4][2];
  for(Int_t ii = 0; ii<4;ii++){
    can_AOQ[ii][0] = new TCanvas(Form("can_AOQ_0_%i",ii), "AOQ 0");
    can_AOQ[ii][0] ->Divide(4,4);
    can_AOQ[ii][1] = new TCanvas(Form("can_AOQ_2_%i",ii), "AOQ 2");
    can_AOQ[ii][1] ->Divide(4,4);
    can_ZET[ii][0] = new TCanvas(Form("can_ZET_0_%i",ii), "ZET 0");
    can_ZET[ii][0] ->Divide(4,4);
    can_ZET[ii][1] = new TCanvas(Form("can_ZET_2_%i",ii), "ZET 2");
    can_ZET[ii][1] ->Divide(4,4);
  }

  Int_t ican = 0, ipad = 1, nentries = (int)3e5;//run = 9, 
  Float_t aoq = 2.64, zet = 50.;
  TF1* f1_aoq[64][2],*f1_zet[64][2];
  for(Int_t run = 9; run<64;run++){

    for(short ii = 0; ii < 2; ii++){     
      f1_aoq[run][ii] = new TF1(Form("f1_aoq_%i_%02d",ii,run),"gaus(0)+pol0(3)",2.63,2.665);
      f1_zet[run][ii] = new TF1(Form("f1_zet_%i_%02d",ii,run),"gaus(0)+pol0(3)",48,52);
      f1_zet[run][ii]->SetParLimits(0,10,1e6);
      f1_zet[run][ii]->SetParameter(1,12e3);
      f1_zet[run][ii]->SetParameter(1,zet);
      f1_zet[run][ii]->SetParLimits(1,zet-1.5,zet+4.5);
      f1_zet[run][ii]->SetParameter(2,0.5);
      f1_zet[run][ii]->SetParLimits(2,0,1.);

      f1_aoq[run][ii]->SetParLimits(0,10,1e6);
      f1_aoq[run][ii]->SetParameter(1,12e3);
      f1_aoq[run][ii]->SetParameter(1,aoq);
      f1_aoq[run][ii]->SetParLimits(1,aoq-0.1,aoq+0.1);
      f1_aoq[run][ii]->SetParameter(2,2e-3);
      f1_aoq[run][ii]->SetParLimits(2,0,5e-3);

    }//for (ii)

    TFile* _file = TFile::Open(Form("data/rootfiles/new/run%04d.root",run));
    TTree* tree = (TTree*)_file->Get("tree");
    h1_AOQ[run][0] = new TH1F(Form("h1_AOQ_%i_0",run),Form("run %i, AOQ 0", run),400,2.5,2.7);
    h1_AOQ[run][1] = new TH1F(Form("h1_AOQ_%i_2",run),Form("run %i, AOQ 2", run),400,2.5,2.7);
    h1_ZET[run][0] = new TH1F(Form("h1_ZET_%i_0",run),Form("run %i, ZET 0", run),200,45,55);
    h1_ZET[run][1] = new TH1F(Form("h1_ZET_%i_2",run),Form("run %i, ZET 2", run),200,45,55);

    tree->Draw(Form("AOQ[0]>>h1_AOQ_%i_0",run),"TMath::Abs(ZET[0]-50)<0.5","",nentries,0);
    tree->Draw(Form("ZET[0]>>h1_ZET_%i_0",run),"TMath::Abs(AOQ[0]-2.64)<0.005","",nentries,0);
    tree->Draw(Form("AOQ[2]>>h1_AOQ_%i_2",run),"TMath::Abs(ZET[0]-50)<0.5 && TMath::Abs(AOQ[0]-2.64)<0.005 && TMath::Abs(delta[2]-delta[4])<0.5","",nentries,0);
    tree->Draw(Form("ZET[2]>>h1_ZET_%i_2",run),"TMath::Abs(ZET[0]-50)<0.5 && TMath::Abs(AOQ[0]-2.64)<0.005 && TMath::Abs(AOQ[2]-2.64)<0.005","",nentries,0);

    can_AOQ[ican][0]->cd(ipad);
    h1_AOQ[run][0]->Draw();
    h1_AOQ[run][0]->Fit(f1_aoq[run][0]);
    can_AOQ[ican][1]->cd(ipad);
    h1_AOQ[run][1]->Draw();
    h1_AOQ[run][1]->Fit(f1_aoq[run][1]);

    can_ZET[ican][0]->cd(ipad);
    h1_ZET[run][0]->Draw();
    h1_ZET[run][0]->Fit(f1_zet[run][0]);
    can_ZET[ican][1]->cd(ipad);
    h1_ZET[run][1]->Draw();
    h1_ZET[run][1]->Fit(f1_zet[run][1]);

    if(ipad==16){
      ipad = 1;
      ican++;
    }else   ipad++;
    fprintf(aoq_in,"%i %f %f \n", run, f1_aoq[run][0]->GetParameter(1),f1_aoq[run][0]->GetParameter(2));
    fprintf(aoq_out,"%i %f %f \n", run, f1_aoq[run][1]->GetParameter(1),f1_aoq[run][1]->GetParameter(2));
    fprintf(zet_in,"%i %f %f \n", run, f1_zet[run][0]->GetParameter(1),f1_zet[run][0]->GetParameter(2));
    fprintf(zet_out,"%i %f %f \n", run, f1_zet[run][1]->GetParameter(1),f1_zet[run][1]->GetParameter(2));



    //delete _file;
    //delete tree;
  }

  fclose(aoq_in);
  fclose(aoq_out);
  fclose(zet_in);
  fclose(zet_out);

}

void efficiency(TString inputname, TString outputname="hists_efficiencies.root",TString mvacut="0.8",TString friendfile=""){


  TString INFV="(abs(nuvtxx_truth)<360-50&&abs(nuvtxy_truth)<600-50&&nuvtxz_truth>50&&nuvtxz_truth<1394-150)";
  TString SIGNAL_DEF="(abs(beamPdg)==14 && abs(neu)==12 && cc==1 &&"+ INFV + ")";
  TString WEIGHTING="projected_weight";


  gROOT->SetStyle("T2K");


  //Figure out how we are going to select
  //Normal selection mode
  TString SELECTION="mvaresult>"+mvacut;
  //Selecting with a friend tree
  if (friendfile!=""){
    SELECTION="mvatree.newmvaresult>"+mvacut;
  }

  //Get the file
  TFile *file = new TFile(inputname,"READ");
  //Get the tree
  TTree *MVASelection = (TTree*) file->Get("mvaselect/MVASelection");
  //If asked for, get the friend tree
  if (friendfile!=""){
    MVASelection->AddFriend("mvatree",friendfile);
  }

  //ENu
  TH1F *h_enu_sig = new TH1F("h_enu_sig","h_enu_sig",25,0,10);
  h_enu_sig->Sumw2();
  h_enu_sig->GetXaxis()->SetTitle("E_{#nu} (GeV)");
  h_enu_sig->GetYaxis()->SetTitle("No. #nu");
  TH1F *h_enu_sel_sig = new TH1F("h_enu_sel_sig","h_enu_sel_sig",25,0,10);
  h_enu_sel_sig->Sumw2();
  h_enu_sel_sig->SetLineColor(4);
  h_enu_sel_sig->SetLineColor(4);
  h_enu_sel_sig->SetFillColor(4);
  h_enu_sel_sig->SetFillStyle(3003);
  h_enu_sel_sig->GetXaxis()->SetTitle("E_{#nu} (GeV)");
  h_enu_sel_sig->GetYaxis()->SetTitle("No. #nu");

  MVASelection->Draw("Ev >> h_enu_sig",WEIGHTING+"*("+SIGNAL_DEF+")");
  MVASelection->Draw("Ev >> h_enu_sel_sig",WEIGHTING+"*("+SIGNAL_DEF+"&&"+SELECTION+")");

  TEfficiency *h_eff_enu = new TEfficiency(*h_enu_sel_sig,*h_enu_sig);
  h_eff_enu->SetTitle(TString("eff;")+h_enu_sig->GetXaxis()->GetTitle()+TString(";Efficiency (No units)"));
  //h_eff_enu->Draw();


  //Q2
  TH1F *h_q2_sig = new TH1F("h_q2_sig","h_q2_sig",25,0,10);
  h_q2_sig->Sumw2();
  h_q2_sig->SetFillStyle(3003);
  h_q2_sig->SetFillColor(1);
  h_q2_sig->GetXaxis()->SetTitle("Q^{2} (GeV^{2})");
  h_q2_sig->GetYaxis()->SetTitle("No. #nu");
  TH1F *h_q2_sel_sig = new TH1F("h_q2_sel_sig","h_q2_sel_sig",25,0,10);
  h_q2_sel_sig->Sumw2();
  h_q2_sel_sig->SetLineColor(4);
  h_q2_sel_sig->SetFillColor(4);
  h_q2_sel_sig->SetFillStyle(3003);
  h_q2_sel_sig->GetXaxis()->SetTitle("Q^{2} (GeV^{2})");
  h_q2_sel_sig->GetYaxis()->SetTitle("No. #nu");

  MVASelection->Draw("Q2 >> h_q2_sig",WEIGHTING+"*("+SIGNAL_DEF+")");
  MVASelection->Draw("Q2 >> h_q2_sel_sig",WEIGHTING+"*("+SIGNAL_DEF+"&&"+SELECTION+")");

  TEfficiency *h_eff_q2 = new TEfficiency(*h_q2_sel_sig,*h_q2_sig);
  h_eff_q2->SetTitle(TString("eff;")+h_q2_sig->GetXaxis()->GetTitle()+TString(";Efficiency (No units)"));



  //Lepton angle
  TH1F *h_lepangle_sig = new TH1F("h_lepangle_sig","h_lepangle_sig",25,0,180);
  h_lepangle_sig->Sumw2();
  h_lepangle_sig->SetFillStyle(3003);
  h_lepangle_sig->SetFillColor(1);
  h_lepangle_sig->GetXaxis()->SetTitle("Lepton angle to neutrino (degrees)");
  h_lepangle_sig->GetYaxis()->SetTitle("No. #nu");
  TH1F *h_lepangle_sel_sig = new TH1F("h_lepangle_sel_sig","h_lepangle_sel_sig",25,0,180);
  h_lepangle_sel_sig->Sumw2();
  h_lepangle_sel_sig->SetLineColor(4);
  h_lepangle_sel_sig->SetFillColor(4);
  h_lepangle_sel_sig->SetFillStyle(3003);
  h_lepangle_sel_sig->GetXaxis()->SetTitle("Lepton angle to neutrino (degrees)");
  h_lepangle_sel_sig->GetYaxis()->SetTitle("No. #nu");

  MVASelection->Draw("LepNuAngle*180/3.14159265359 >> h_lepangle_sig",WEIGHTING+"*("+SIGNAL_DEF+")");
  MVASelection->Draw("LepNuAngle*180/3.14159265359 >> h_lepangle_sel_sig",WEIGHTING+"*("+SIGNAL_DEF+"&&"+SELECTION+")");

  TEfficiency *h_eff_lepangle = new TEfficiency(*h_lepangle_sel_sig,*h_lepangle_sig);
  h_eff_lepangle->SetTitle(TString("eff;")+h_lepangle_sig->GetXaxis()->GetTitle()+TString(";Efficiency (No units)"));

  //Lepton momentum
  TH1F *h_lepmom_sig = new TH1F("h_lepmom_sig","h_lepmom_sig",25,0,10);
  h_lepmom_sig->Sumw2();
  h_lepmom_sig->SetFillStyle(3003);
  h_lepmom_sig->SetFillColor(1);
  h_lepmom_sig->GetXaxis()->SetTitle("Lepton momentum (GeV)");
//.........这里部分代码省略.........