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

void Init(){

  gSystem->Load("libSTAT");
  TStatToolkit toolkit;
  Double_t chi2;
  TVectorD fitParam, fitParam1,fitParam2;
  TVectorD errParam, errParam1,errParam2;
  TMatrixD covMatrix,covMatrix1,covMatrix2;
  
  Int_t npoints;
  //
  TFile f("driftitsTPC.root");
  TTree * tree = (TTree*)f.Get("Test");
  tree->SetAlias("side","(-1+(pTPC.fP[1]>0)*2)");     //side
  tree->SetAlias("z","(pTPC.fP[1]+0.0)");     //z position
  tree->SetAlias("dr","(1-abs(pTPC.fP[1])/250.)");    //norm drift length
  tree->SetAlias("tl","(pTPC.fP[3]+pITS.fP[3])*0.5"); //tan lampbda
  tree->SetAlias("sa","sin(pTPC.fAlpha+0.)");         //sin alpha
  tree->SetAlias("ca","cos(pTPC.fAlpha+0.)");         //cos alpha
  
  tree->SetAlias("dz","(pTPC.fP[1]-pITS.fP[1])");      //z delta
  tree->SetAlias("dy","(pTPC.fP[0]-pITS.fP[0])");      //z delta
  tree->SetAlias("dtl","(pTPC.fP[3]-pITS.fP[3])");     //delta tan lampbda
  tree->SetAlias("etl","sqrt(pITS.fC[9]+0.)");         //error tan lampbda
  tree->SetAlias("ez","sqrt(pITS.fC[2]+0.)");          //error z
}

void leadFrag(char *infname = "dj_HCPR-J50U-hiGoodMergedTracks_OfficialSelv2_Final0.root",
    TCut myCut = "cent<30",
    TString module="djcalo_tower")
{
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get(module+"/djTree");

   const int nBin = 5;
   double bins[nBin+1] = {0.5,1.5,4,8,20,1000};  
   double colors[nBin] = {kGray,kBlue-3,38,kOrange-8,kRed-6};

   // Selection cut
   TCut evtCut = "nljet>120&&abs(nljeta)<2&&aljet>50&&abs(aljeta)<2&&jdphi>2./3*TMath::Pi()";
   t->SetAlias("Aj","(nljet-aljet)/(nljet+aljet)");

   // histograms
   TH2D * hFrag_LJEt = new TH2D("hFrag_LJEt",";p_{T}^{J1};p_{T}^{Tower}/p_{T}^{J1}",50,120,300,50,0,1);
   hFrag_LJEt->GetXaxis()->CenterTitle();
   hFrag_LJEt->GetYaxis()->CenterTitle();
   TH2D * hLeadFrag_LJEt = new TH2D("hLeadFrag_LJEt",";p_{T}^{J1};p_{T}^{Highest Tower}/p_{T}^{J1}",50,120,300,50,0,1);
   hLeadFrag_LJEt->GetXaxis()->CenterTitle();
   hLeadFrag_LJEt->GetYaxis()->CenterTitle();

   // project tree
   t->Project("hLeadFrag_LJEt","lppt[0]/nljet:nljet",evtCut);
   t->Project("hFrag_LJEt","ppt/nljet:nljet",evtCut);


   // plot
   TCanvas *c2 = new TCanvas("c2","",600,600);
   c2->SetRightMargin(0.15);
   c2->SetLogz();
   hFrag_LJEt->Draw("colz");

   TCanvas *c3 = new TCanvas("c3","",600,600);
   c3->SetRightMargin(0.15);
   c3->SetLogz();
   hLeadFrag_LJEt->Draw("colz");

   TCanvas *c4 = new TCanvas("c4","",600,600);
   hLeadFrag_LJEt->ProfileX()->Draw("E");
   TH1D * hLeadFrag_LJEt_pfx = (TH1D*)gDirectory->Get("hLeadFrag_LJEt_pfx");
   hLeadFrag_LJEt_pfx->SetAxisRange(0,0.6,"Y");
   hFrag_LJEt->ProfileX()->Draw("same hist");
   TH1D * hFrag_LJEt_pfx = (TH1D*)gDirectory->Get("hFrag_LJEt_pfx");
}

//________________________________________________________________________
void draw_tree()
{
   TFile* file = TFile::Open("example.root");
   // create a TTree   
   TTree *tree = (TTree*)file->Get("tree");
   
   tree->SetAlias("mult", "fOffset*fX.fElements");
   
   TCanvas* c1 = new TCanvas("c1");
   c1->Divide(2,3);
   c1->cd(1);
   tree->Draw("fOffset*fX.fElements");
   
   c1->cd(2);
   tree->Draw("mult");
   
   PlotAndCheck(tree,c1,3,"fX.fElements > 1.0");
   PlotAndCheck(tree,c1,5,"fY.fElements > 1.0");
   
}

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
		    )
{
   // ===========================================================
   // Get Input
   // ===========================================================
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get("ntjt");


   t->SetAlias("metxMergedAll","(metOutOfConex0+metOutOfConex1+metOutOfConex2+metOutOfConex3+metOutOfConex4+metOutOfConex5)");
   t->SetAlias("metxMerged0","(metOutOfConex0)");
   t->SetAlias("metxMerged1","(metOutOfConex1)");
   t->SetAlias("metxMerged2","(metOutOfConex2)");
   t->SetAlias("metxMerged3","(metOutOfConex3)");
   t->SetAlias("metxMerged4","(metOutOfConex4+metOutOfConex5)");

   // ===========================================================
   // Analysis Setup
   // ===========================================================
   const int nBin = 5;
   double bins[nBin+1] = {0.5,1.0,1.5,4,8,1000};
   int colors[5] = {kBlue-10,kYellow-7, kOrange-2,kGreen-5,kRed-3};
   
   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";
//   TCut evtCut = "nljet>120&&abs(nljetacorr)<2&&aljet>50&&abs(aljetacorr)<2&&jdphi>3&&!maskEvt";

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

   // ===========================================================
   // Find Average Weights
   // ===========================================================
   TH1D *hw[nBinAj];
   float meanWt[nBinAj];
   for ( int iaj = 0 ; iaj< nBinAj ; iaj++) {
      hw[iaj] = new TH1D(Form("hw_aj%d",iaj),"",1000,0,100);
      TCut ajCut = Form("Aj>%f && Aj<%f", ajBins[iaj],ajBins[iaj+1]);
      t->Draw(Form("weight>>hw_aj%d",iaj), evtCut&&myCut&&ajCut);
      meanWt[iaj] = hw[iaj]->GetMean();
      cout << " <Weight>: " << meanWt[iaj] << endl;
   }

   // ===========================================================
   // Draw Weighted Averages
   // ===========================================================
   TH1D *ppos[nBin];
   TH1D *pneg[nBin];
   TH1D *pe[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);
      TH1D *he[nBinAj];
      
      // =================================
      // Get Weighted Mean for each Aj bin
      // =================================
      h1->Sumw2();
      h2->Sumw2();
      t->Draw(Form("Aj>>h1%d",i), "weight"*(evtCut&&myCut));
      t->Draw(Form("Aj>>h2%d",i), Form("((-weight*metxMerged%d))",i)*(evtCut&&myCut));
      
      pe[i]=(TH1D*)h2->Clone();
      pe[i]->SetName(Form("p%d",i));     
      pe[i]->Divide(h1);
      ppos[i] = new TH1D(Form("ppos%d",i),"",nBinAj,ajBins);
      ppos[i]->SetLineColor(1);     
      ppos[i]->SetMarkerColor(colors[i]);
      ppos[i]->SetFillColor(colors[i]);
      ppos[i]->SetFillStyle(1001);
      pneg[i] = new TH1D(Form("pneg%d",i),"",nBinAj,ajBins);
      pneg[i]->SetLineColor(1);     
      pneg[i]->SetMarkerColor(colors[i]);
      pneg[i]->SetFillColor(colors[i]);
      pneg[i]->SetFillStyle(1001);

      // =================================
      // Caculated Stat Error of the Mean
      // =================================
      cout << "Stat Error for pt bin " << i << ": ";
      for ( int iaj = 0 ; iaj< nBinAj ; iaj++) {
         he[iaj] = new TH1D(Form("he%d_aj%d",i,iaj),"",100,-200,200);
         TCut ajCut = Form("Aj>%f && Aj<%f", ajBins[iaj],ajBins[iaj+1]);
         t->Draw(Form("((metxMerged%d))>>he%d_aj%d",i,i,iaj), "weight" * evtCut&&myCut&&ajCut);
         float theError = he[iaj]->GetRMS()/ (sqrt(he[iaj]->GetEntries()));
	 cout << theError << " ";
	 pe[i]->SetBinError(iaj+1, theError);
      }
      cout << endl;
   }

   // Stack
   for (int i=nBin-1;i>=0;i--)
   {
      for(int iaj = 0 ; iaj< nBinAj ; iaj++) {
//.........这里部分代码省略.........

void plotRBDphi( 
		double dphiCut=3.026,
		// double dphiCut=2.094,
	      double dphiCut2 = 3.026,
		 TString infname = "data.root",
		 TString pythia = "pythia.root",
		 TString mix = "mix.root",
		 bool useWeight = true,
		 bool drawXLabel = false,
		 bool drawLeg = true)
{		
  useWeight = 0;

  int threshold1 = 80;
  int threshold2 = 120;
  gStyle->SetErrorX(0); 
  TString cut1=Form("abs(eta1) < 2 && abs(eta2) < 2 && et1>%d&&et2>50",threshold1);
  TString cut2=Form("abs(eta1) < 2 && abs(eta2) < 2 && et1>%d&&et2>50",threshold2);
  cout <<cut1.Data()<<endl;
  cout <<cut2.Data()<<endl;

  TString trigcut = "";
  TString cstring = "";

  // 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");

  // open output
  TFile *outfile = new TFile("output.root","recreate");
  TNtuple *ntOut = new TNtuple("ntOut","","npart");


  nt->SetAlias("et1","pt1");
  nt->SetAlias("et2","pt2");

  ntMix->SetAlias("et1","pt1");
  ntMix->SetAlias("et2","pt2");
  ntPythia->SetAlias("et1","pt1");
  ntPythia->SetAlias("et2","pt2");

  //  ntMix->SetAlias("pt1","et1");
  //  ntMix->SetAlias("pt2","et2");

  //  ntPythia->SetAlias("pt1","et1");
  //  ntPythia->SetAlias("pt2","et2");

  ntPythia->SetAlias("weight","1");
  ntMix->SetAlias("weight","1");

  const int nBin = 6;
  double m[nBin+1] = {-1.5,-0.5,3.5,7.5,11.5,20.5,40.5};
  double npart[nBin] = {2,358.623,232.909,97.9521};
  double npart2[nBin] = {2,358.623,232.909,97.9521};
  
  double npartValue[40];
  npartValue[0] = 393.633;
  npartValue[1] = 368.819;
  npartValue[2] = 343.073;
  npartValue[3] = 317.625;
  npartValue[4] = 292.932;
  npartValue[5] = 271.917;
  npartValue[6] = 249.851;
  npartValue[7] = 230.72;
  npartValue[8] = 212.465;
  npartValue[9] = 194.752;
  npartValue[10] = 178.571;
  npartValue[11] = 163.23;
  npartValue[12] = 149.187;
  npartValue[13] = 136.011;
  npartValue[14] = 123.414;
  npartValue[15] = 111.7;
  npartValue[16] = 100.831;
  npartValue[17] = 90.7831;
  npartValue[18] = 80.9823;
  npartValue[19] = 72.6236;
  npartValue[20] = 64.1508;
  npartValue[21] = 56.6284;
  npartValue[22] = 49.9984;
  npartValue[23] = 43.3034;
  npartValue[24] = 37.8437;
  npartValue[25] = 32.6659;
  npartValue[26] = 27.83;
  npartValue[27] = 23.7892;
  npartValue[28] = 20.1745;
  npartValue[29] = 16.8453;
  npartValue[30] = 14.0322;
  npartValue[31] = 11.602;
  npartValue[32] = 9.52528;
  npartValue[33] = 7.6984;
  npartValue[34] = 6.446;
  npartValue[35] = 4.96683;
//.........这里部分代码省略.........

void doPostProcessing(TString infname, TString outfile, Int_t events, 
		      Float_t xsec, Float_t kfactor,
		      Float_t filt_eff, bool SortBasketsByEntry = false) {
  
  cout << "Processing File " << infname << endl;
  
  TFile *f = TFile::Open(infname.Data(), "READ");
  if (! f || f->IsZombie()) {
    cout << "File does not exist!" << endl;
    return;
  }
  
  TTree* t = (TTree*)f->Get("Events");
  if (! t || t->IsZombie()) {
    cout << "Tree does not exist!" << endl;
    return;
  }
        
  //-------------------------------------------------------------
  // Removes all non *_CMS2.* branches
  //-------------------------------------------------------------`
  t->SetBranchStatus("*", 0);
  t->SetBranchStatus("*_CMS2.*", 1);

  // Removes the branches (if they exist) that we want to replace
  //evt_xsec_excl
  TString bName = t->GetAlias("evt_xsec_excl");
  //cout << "evt_xsec_excl " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0); 
  }

  //evt_xsec_incl
  bName = t->GetAlias("evt_xsec_incl");
  //cout << "evt_xsec_incl " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0);   
  }
  
  //evt_kfactor
  bName = t->GetAlias("evt_kfactor");
  //cout << "evt_kfactor " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0); 
  }

  //evt_nEvts
  bName = t->GetAlias("evt_nEvts");
  //cout << "evt_nEvts " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0); 
  }

  //evt_filt_eff
  bName = t->GetAlias("evt_filt_eff");
  //cout << "evt_filt_eff " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0); 
  }

  //evt_scale1fb
  bName = t->GetAlias("evt_scale1fb");
  //cout << "evt_scale1fb " << bName << endl;
  if(bName != "") {
    bName.ReplaceAll(".obj", "*");
    t->SetBranchStatus(bName.Data(), 0); 
  }

  TFile *out = TFile::Open(outfile.Data(), "RECREATE");
  TTree *clone;
  if(SortBasketsByEntry)
    clone = t->CloneTree(-1, "fastSortBasketsByEntry");
  else 
    clone = t->CloneTree(-1, "fast");
   

  //-------------------------------------------------------------

  //Calculate scaling factor and put variables into tree 
  Float_t scale1fb = xsec*kfactor*1000*filt_eff/(Float_t)events;
  cout << "scale1fb: " << scale1fb << endl; 

  TBranch* b1 = clone->Branch("evtscale1fb", &scale1fb, "evt_scale1fb/F");
  TBranch* b2 = clone->Branch("evtxsecexcl", &xsec, "evt_xsec_excl/F");
  TBranch* b3 = clone->Branch("evtxsecincl", &xsec, "evt_xsec_incl/F");
  TBranch* b4 = clone->Branch("evtkfactor", &kfactor, "evt_kfactor/F");
  TBranch* b5 = clone->Branch("evtnEvts", &events, "evt_nEvts/I");
  TBranch* b6 = clone->Branch("evtfilteff", &filt_eff, "evt_filt_eff/F");
   
  clone->SetAlias("evt_scale1fb",  "evtscale1fb");
  clone->SetAlias("evt_xsec_excl", "evtxsecexcl");
  clone->SetAlias("evt_xsec_incl",  "evtxsecincl");
  clone->SetAlias("evt_kfactor",   "evtkfactor");
  clone->SetAlias("evt_nEvts",     "evtnEvts");
  clone->SetAlias("evt_filt_eff",     "evtfilteff");
//.........这里部分代码省略.........

void analyzePuRejection(string var = "etaDijet"){


  TH1::SetDefaultSumw2();


  TString infname = "/d101/yetkin/analysis/d0204/ntuple_data_pPb_akPu3PF_forest71_20130204_01.root";
  TString refname = "/d101/yetkin/analysis/d0204/ntuple_data_PbPb_akPu3PF_forest71_20130204_01.root";
  TString mixname = "/d101/yetkin/analysis/d0204/hijing.root";

  TFile *inf = new TFile(infname.Data());
  TTree *nt =(TTree*)inf->FindObjectAny("ntdijet");
  TTree *ntevt =(TTree*)inf->FindObjectAny("ntevt");
  nt->AddFriend(ntevt);

  TFile *infMix = new TFile(mixname.Data());
  TTree *ntMix =(TTree*)infMix->FindObjectAny("ntdijet");
  TTree *ntevtMix =(TTree*)infMix->FindObjectAny("ntevt");
  TTree *ntw =(TTree*)infMix->FindObjectAny("ntw");
  ntMix->AddFriend(ntw);
  ntMix->AddFriend(ntevtMix);

  nt->SetAlias("pt1","jtpt1");
  nt->SetAlias("pt2","jtpt2");
  nt->SetAlias("eta1","jteta1");
  nt->SetAlias("eta2","jteta2");
  nt->SetAlias("phi1","jtphi1");
  nt->SetAlias("phi2","jtphi2");

  ntMix->SetAlias("pt1","jtpt1");
  ntMix->SetAlias("pt2","jtpt2");
  ntMix->SetAlias("eta1","jteta1");
  ntMix->SetAlias("eta2","jteta2");
  ntMix->SetAlias("phi1","jtphi1");
  ntMix->SetAlias("phi2","jtphi2");

  nt->SetAlias("cent","hfp");
  ntMix->SetAlias("cent","hfp");

  nt->SetAlias("etaDijet","(eta1+eta2)/2");
  ntMix->SetAlias("etaDijet","(eta1+eta2)/2");

  nt->SetAlias("x","pt2/pt1");
  ntMix->SetAlias("x","pt2/pt1");

  nt->SetAlias("dphi","acos(cos(phi1-phi2))");
  ntMix->SetAlias("dphi","acos(cos(phi1-phi2))");


  nt->SetAlias("var",var.data());
  ntMix->SetAlias("var",var.data());

  nt->SetAlias("w","pt1> -999");
  ntMix->SetAlias("w","weight*(pthat > 50)");

  TCut noises[2] = {"noise == -1","noise == 1"};

  TTree* t[2] = {nt,ntMix};

  TH2D* h[2];
  TH1D* hMean[3];
  TH1D* hSigma[3];

  double yMax = 250;
  double yMin = 0;

  string varName = "N_{trk}^{offline}";

  if(var == "etaDijet"){
    yMax = 2;
    yMin = -2;
    varName = "#eta_{dijet}=(#eta_{1}+#eta_{2})/2";
  }

  if(var == "x"){
    yMax = 1;
    yMin = 0;
    varName = "p_{T,2}/p_{T,1}";
  }

  if(var == "dphi"){
    yMax = pi;
    yMin = pi-pi/3;
    varName = "#Delta#phi_{1,2}";
  }


  TCanvas* c1 = new TCanvas("c1","",1200,600);
  c1->Divide(2,1);

  for(int i = 0; i < 2; ++i){
    c1->cd(i+1)->SetLogz();
    h[i] = new TH2D(Form("h%d",i),Form(";E_{T}^{HF[#eta > 4]};%s",varName.data()),14,0,70,250,yMin,yMax);
    h[i]->GetXaxis()->CenterTitle();
    h[i]->GetYaxis()->CenterTitle();

    t[0]->Draw(Form("var:cent>>%s",h[i]->GetName()),dijet&&deltaPhi&&jetID&&vtx&&noises[i],"colz");
 
    fitSlices(h[i]);

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

void plotFeedDown(int ntest=1, int centL=0,int centH=100)
{
   // B cross-section
   TFile *inf = new TFile("output_pp_Bmeson_5TeV_y1.root");
//   TFile *inf = new TFile("outputBplus_D_pp_rap24.root");
//    TFile *inf = new TFile("outputBplus_pp.root");
   TH1D *hBPtMax = (TH1D*)inf->Get("hmaxall");
   TH1D *hBPtMin = (TH1D*)inf->Get("hminall");
   TH1D *hBPt = (TH1D*)inf->Get("hpt");
   hBPt->SetName("hBPt");
   hBPtMax->SetName("hBPtMax");
   hBPtMin->SetName("hBPtMin");

   TH1D *hBMaxRatio = (TH1D*)hBPt->Clone("hBMaxRatio");
   hBMaxRatio->Divide(hBPtMax);

   TH1D *hBMinRatio = (TH1D*)hBPt->Clone("hBMinRatio");
   hBMinRatio->Divide(hBPtMin);

   
   hBPt->Rebin(1); 
   
   // D cross-section
//   TFile *infD = new TFile("outputD0_D_pp.root");
   TFile *infD = new TFile("output_pp_d0meson_5TeV_y1.root");
   TH1D *hDPtMax = (TH1D*)infD->Get("hmaxall");
   TH1D *hDPtMin = (TH1D*)infD->Get("hminall");
   TH1D *hDPt = (TH1D*)infD->Get("hpt");
   hDPt->SetName("hDPt");
   hDPtMax->SetName("hDPtMax");
   hDPtMin->SetName("hDPtMin");
   hDPt->Rebin(1); 

   // ratio of B->D0: not correct85% from PYTHIA
   //hBPt->Scale(0.85);
   hBPt->Scale(0.598);
   
   // c->D (55.7%)
   hDPt->Scale(0.557);

   
   TFile *inf2 = new TFile("/data/HeavyFlavourRun2/BtoDPythia/treefile_merged.root");
//   TFile *inf2 = new TFile("test.root");
   TTree *hi = (TTree*) inf2->Get("ana/hi");

   hi->SetAlias("yD","log((sqrt(1.86484*1.86484+pt*pt*cosh(eta)*cosh(eta))+pt*sinh(eta))/sqrt(1.86484*1.86484+pt*pt))");			    
   hi->SetAlias("yB","log((sqrt(5.3*5.3+pt*pt*cosh(eta)*cosh(eta))+pt*sinh(eta))/sqrt(5.3*5.3+pt*pt))");			    
   hi->SetAlias("yJ","log((sqrt(3.09692*3.09692+pt*pt*cosh(eta)*cosh(eta))+pt*sinh(eta))/sqrt(3.09692*3.09692+pt*pt))");			    


   // 6.5, 8, 10, 13, 30

/*
   TH1D *hBNoCut = (TH1D*)hBPt->Clone("hBNoCut");
   TH1D *hBHasD  = (TH1D*)hBPt->Clone("hBHasD");
   
   hi->Draw("pt>>hBHasD","(abs(pdg)>500&&abs(pdg)<600&&abs(yB)<2.4)&&Sum$(abs(pdg)==421&&abs(yD)<2)>0");
   hi->Draw("pt>>hBNoCut","(abs(pdg)>500&&abs(pdg)<600&&abs(yB)<2.4)");
;

   hBNoCut->Divide(hBHasD);
   hBPt->Divide(hBNoCut);
  */  

// 0-100%
   int npoint = 7;
   	
   double ptBins_npjpsi[8] = {1,3,6.5,8,10,13,30,300};
  
   double raa_npjpsi[7];//      = {1,0.6, 0.52,0.43,0.43,0.34,0.5};  
   double raaStat_npjpsi[7];//  = {1,0.4,0.12,0.08,0.09,0.07,0.5};
   double raaSyst_npjpsi[7];//  = {0,0,0.06,0.05,0.05,0.04,0};

/*
0-10, 10-20, 20-30, 30-40, 40-50, 50-100
double nonPromptJpsiRAA_2012[]           = {0.,0.38,0.43,0.48,0.52,0.65,0.69};
double nonPromptJpsiRAAError_2012[]      = {0.,0.02,0.03,0.03,0.04,0.06,0.07};
double nonPromptJpsiRAAErrorSyst_2012[]  = {0.,0.04,0.05,0.05,0.06,0.07,0.07};
*/

   if (centL==0&&centH==100) {
      raa_npjpsi[0]=1.0;	raaStat_npjpsi[0]=0.0;		raaSyst_npjpsi[0]=1.0;	// no measurement
      raa_npjpsi[1]=0.6;	raaStat_npjpsi[1]=0.0;		raaSyst_npjpsi[1]=0.4;	// prelim
      raa_npjpsi[2]=0.52;	raaStat_npjpsi[2]=0.12;		raaSyst_npjpsi[2]=0.06;	// np jpsi pas
      raa_npjpsi[3]=0.43;	raaStat_npjpsi[3]=0.08;		raaSyst_npjpsi[3]=0.05;	// np jpsi pas
      raa_npjpsi[4]=0.43;	raaStat_npjpsi[4]=0.09;		raaSyst_npjpsi[4]=0.05;	// np jpsi pas
      raa_npjpsi[5]=0.34;	raaStat_npjpsi[5]=0.07;		raaSyst_npjpsi[5]=0.04;	// np jpsi pas
      raa_npjpsi[6]=0.5;	raaStat_npjpsi[6]=0.0;		raaSyst_npjpsi[6]=0.25;	// b-jet
   }
   
   if (centL==0&&centH==10) {
      raa_npjpsi[0]=1.0;	raaStat_npjpsi[0]=0.0;		raaSyst_npjpsi[0]=1.0;	// no measurement
      raa_npjpsi[1]=1.0;	raaStat_npjpsi[1]=0.0;		raaSyst_npjpsi[1]=1.0;	// no measurement
      raa_npjpsi[2]=0.38;	raaStat_npjpsi[2]=0.02;		raaSyst_npjpsi[2]=0.04;	// np jpsi pas
      raa_npjpsi[3]=0.38;	raaStat_npjpsi[3]=0.02;		raaSyst_npjpsi[3]=0.04;	// np jpsi pas
      raa_npjpsi[4]=0.38;	raaStat_npjpsi[4]=0.02;		raaSyst_npjpsi[4]=0.04;	// np jpsi pas
      raa_npjpsi[5]=0.38;	raaStat_npjpsi[5]=0.02;		raaSyst_npjpsi[5]=0.04;	// np jpsi pas
      raa_npjpsi[6]=0.39;	raaStat_npjpsi[6]=0.0;		raaSyst_npjpsi[6]=0.20;	// b-jet
   
   }
//.........这里部分代码省略.........

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 DoAnalysisWithTree(TFile* file, int run, int x, int y)
{

  gout << 0 << " " << run << " " << x << " " << y << " ";

  TF1* flandau = new TF1("flandau","landau",0,400);

  TCanvas *c1 = new TCanvas("c1","",1280,960);
  c1->Divide(4,2);

  TCanvas *c2 = new TCanvas("c2","",1280,960);
  c2->Divide(4,2);

  TCanvas *c3 = new TCanvas("c3","",1280,960);
  c3->Divide(4,2);

  TCanvas *c4 = new TCanvas("c4","",1280,960);
  c4->Divide(4,2);

  TTree* tree = (TTree*)file->Get("T");

  tree->SetAlias("Average_HODO_HORIZONTAL","Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
  tree->SetAlias("Average_HODO_VERTICAL","Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
  tree->SetAlias("Valid_HODO_HORIZONTAL","Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) > 0");
  tree->SetAlias("Valid_HODO_VERTICAL","Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) > 0");

  ofstream fout((const char*)Form("DataRunByRun/mapper_run%d.dat",run));

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

      TString drawstring1 = "TOWER_CALIB_TILE_MAPPER[";
      drawstring1 += i*2;
      drawstring1 += "].energy>>hs1_";
      drawstring1 += i+1;
      drawstring1 += "(100,0,400)";

      c1->cd(i+1);
      tree->Fit("flandau",drawstring1,"","","");

      TString drawstring2 = "TOWER_CALIB_TILE_MAPPER[";
      drawstring2 += i*2;
      drawstring2 += "].energy>>hs2_";
      drawstring2 += i+1;
      drawstring2 += "(100,0,400)";

      c2->cd(i+1);
      tree->Fit("flandau",drawstring2,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL","","");
      fout << run << " " << 9999 << " " << 9999 << " " << flandau->GetParameter(1) << " " << flandau->GetParError(1) << endl;
      gout << flandau->GetParameter(1) << " " << flandau->GetParError(1) << " ";

      TString drawstring3 = "TOWER_CALIB_TILE_MAPPER[";
      drawstring3 += i*2;
      drawstring3 += "].energy>>hs3_";
      drawstring3 += i+1;
      drawstring3 += "(100,0,400)";

      c3->cd(i+1);
      tree->Fit("flandau",drawstring3,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL && abs(TOWER_CALIB_C2[3].energy)<200","","");

      TString drawstring4 = "TOWER_CALIB_TILE_MAPPER[";
      drawstring4 += i*2;
      drawstring4 += "].energy>>hs4_";
      drawstring4 += i+1;
      drawstring4 += "(100,0,400)";

      c4->cd(i+1);
      tree->Fit("flandau",drawstring4,"Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL && abs(TOWER_CALIB_C2[3].energy)>200","","");

    }

  c1->Print(Form("FigsRunByRun/mapper_run%d_step1.png",run));
  c2->Print(Form("FigsRunByRun/mapper_run%d_step2.png",run));
  c3->Print(Form("FigsRunByRun/mapper_run%d_step3.png",run));
  c4->Print(Form("FigsRunByRun/mapper_run%d_step4.png",run));
  c1->Print(Form("FigsRunByRun/mapper_run%d_step1.pdf",run));
  c2->Print(Form("FigsRunByRun/mapper_run%d_step2.pdf",run));
  c3->Print(Form("FigsRunByRun/mapper_run%d_step3.pdf",run));
  c4->Print(Form("FigsRunByRun/mapper_run%d_step4.pdf",run));

  delete c1;
  delete c2;
  delete c3;
  delete c4;

  gout << endl;

  return;

  // --- let's try to have a look at the hodoscope positions

  TCanvas* c5 = new TCanvas("c5","",800,800);
  TH2D* th2d_hodo_fine = new TH2D("th2d_hodo_fine","",40,-0.5,7.5,40,-0.5,7.5);
  TH2D* th2d_hodo_coarse = new TH2D("th2d_hodo_coarse","",8,-0.5,7.5,8,-0.5,7.5);
  th2d_hodo_fine->GetXaxis()->SetTitle("Horizontal Hodoscope");
  th2d_hodo_fine->GetYaxis()->SetTitle("Vertical Hodoscope");
  th2d_hodo_coarse->GetXaxis()->SetTitle("Horizontal Hodoscope");
  th2d_hodo_coarse->GetYaxis()->SetTitle("Vertical Hodoscope");
  tree->Draw("Average_HODO_HORIZONTAL:Average_HODO_VERTICAL>>th2d_hodo_fine","Valid_HODO_VERTICAL && Valid_HODO_HORIZONTAL","colz");
  c5->Print(Form("FigsRunByRun/hodoscope_fine_run%d.png",run));
//.........这里部分代码省略.........

void plotRJ(  double ajCut=0.24,
	      int cbin = 0,
		 TString infname = "data.root",
		 TString pythia = "pythia.root",
		 TString mix = "mix.root",
		 bool useWeight = true,
		 bool drawXLabel = false,
		 bool drawLeg = true)
{		
  gStyle->SetErrorX(0); 
  TString cut="et1>120 && et2>50 && abs(dphi)>2.5 && abs(eta1) < 2 && abs(eta2) < 2";
  TString cutpp="et1>120 && et2>50 && abs(dphi)>2.5 && abs(eta1) < 2 && abs(eta2) < 2";
  TString trigcut = "";
  TString cstring = "";

  // 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");

  nt->SetAlias("et1","pt1");
  nt->SetAlias("et2","pt2");

  ntPythia->SetAlias("pt1","et1");
  ntPythia->SetAlias("pt2","et2");

  ntMix->SetAlias("pt1","et1");
  ntMix->SetAlias("pt2","et2");

  // open output
  TFile *outfile = new TFile("output.root","recreate");
  TNtuple *ntOut = new TNtuple("ntOut","","a:b");

  const int nBin = 4;
  double m[nBin+1] = {-1.5,-0.5,3.5,11.5, 40.5};
  double npart0[nBin] = {2,355.788,163.34,46.6685}; // Flat Averages
  double npart[nBin] = {2,358.623,232.909,97.9521};
  
  
 /*
  double npart[40];
  npart[0] = 393.633;
  npart[1] = 368.819;
  npart[2] = 343.073;
  npart[3] = 317.625;
  npart[4] = 292.932;
  npart[5] = 271.917;
  npart[6] = 249.851;
  npart[7] = 230.72;
  npart[8] = 212.465;
  npart[9] = 194.752;
  npart[10] = 178.571;
  npart[11] = 163.23;
  npart[12] = 149.187;
  npart[13] = 136.011;
  npart[14] = 123.414;
  npart[15] = 111.7;
  npart[16] = 100.831;
  npart[17] = 90.7831;
  npart[18] = 80.9823;
  npart[19] = 72.6236;
  npart[20] = 64.1508;
  npart[21] = 56.6284;
  npart[22] = 49.9984;
  npart[23] = 43.3034;
  npart[24] = 37.8437;
  npart[25] = 32.6659;
  npart[26] = 27.83;
  npart[27] = 23.7892;
  npart[28] = 20.1745;
  npart[29] = 16.8453;
  npart[30] = 14.0322;
  npart[31] = 11.602;
  npart[32] = 9.52528;
  npart[33] = 7.6984;
  npart[34] = 6.446;
  npart[35] = 4.96683;
  npart[36] = 4.23649;
  npart[37] = 3.50147;
  npart[38] = 3.16107;
  npart[39] = 2.7877;
  */

  TH1D *hTmp = new TH1D("hTmp","",100,0,400);
  TH1D *h = new TH1D("h","",nBin,m);
  TH1D *h2 = new TH1D("h2","",nBin,m);
  
  nt->Draw("bin>>h",Form("(et1-et2)/(et1+et2)>%f&&%s",ajCut,cut.Data()));
  nt->Draw("bin>>h2",Form("%s",cut.Data()));
  TGraphAsymmErrors *g = calcEff(h2,h,npart);

  ntPythia->Draw("bin>>h",Form("(et1-et2)/(et1+et2)>%f&&%s",ajCut,cut.Data()));
  ntPythia->Draw("bin>>h2",Form("%s",cut.Data()));
//.........这里部分代码省略.........

void MakeATGCRatioHistograms() {
  char* dir = "/uscms_data/d2/andersj/Wjj/2012/data/Moriond2013/ReducedTrees/";
  TFile* fSM = new TFile((dir + string("Lambda00_Kappa00_G00.root")).c_str());
  TFile* fSM2 = new TFile((string("Lambda00_Kappa00_G00.root")).c_str());
  TFile* f_L01K00G00 = new TFile((string("Lambda01_Kappa00_G00.root")).c_str());
  TFile* f_L03K00G00 = new TFile((string("Lambda03_Kappa00_G00.root")).c_str());
  TFile* f_L05K00G00 = new TFile((string("Lambda05_Kappa00_G00.root")).c_str());
  TFile* f_L07K00G00 = new TFile((dir + string("Lambda07_Kappa00_G00.root")).c_str());
  TFile* f_L09K00G00 = new TFile((dir + string("Lambda09_Kappa00_G00.root")).c_str());
  TFile* f_L11K00G00 = new TFile((dir + string("Lambda11_Kappa00_G00.root")).c_str());
  TFile* f_L00K05G00 = new TFile((string("Lambda00_Kappa05_G00.root")).c_str());
  TFile* f_L00K11G00 = new TFile((dir + string("Lambda00_Kappa11_G00.root")).c_str());
  TFile* f_L00K16G00 = new TFile((dir + string("Lambda00_Kappa16_G00.root")).c_str());
  TFile* f_L00K20G00 = new TFile((dir + string("Lambda00_Kappa20_G00.root")).c_str());
  TFile* f_L00K00G11 = new TFile((dir + string("Lambda00_Kappa00_G11.root")).c_str());
  TFile* f_L00K00G40 = new TFile((dir + string("Lambda00_Kappa00_G40.root")).c_str());
  TFile* f_L00K00G60 = new TFile((dir + string("Lambda00_Kappa00_G60.root")).c_str());


  TTree* trSM = (TTree*) fSM->Get("WJet");
  TTree* trSM2 = (TTree*) fSM2->Get("WJet");
  TTree* tr_L01K00G00 = (TTree*) f_L01K00G00->Get("WJet");
  TTree* tr_L03K00G00 = (TTree*) f_L03K00G00->Get("WJet");
  TTree* tr_L05K00G00 = (TTree*) f_L05K00G00->Get("WJet");
  TTree* tr_L07K00G00 = (TTree*) f_L07K00G00->Get("WJet");
  TTree* tr_L09K00G00 = (TTree*) f_L09K00G00->Get("WJet");
  TTree* tr_L11K00G00 = (TTree*) f_L11K00G00->Get("WJet");
  TTree* tr_L00K05G00 = (TTree*) f_L00K05G00->Get("WJet");
  TTree* tr_L00K11G00 = (TTree*) f_L00K11G00->Get("WJet");
  TTree* tr_L00K16G00 = (TTree*) f_L00K16G00->Get("WJet");
  TTree* tr_L00K20G00 = (TTree*) f_L00K20G00->Get("WJet");
  TTree* tr_L00K00G11 = (TTree*) f_L00K00G11->Get("WJet");
  TTree* tr_L00K00G40 = (TTree*) f_L00K00G40->Get("WJet");
  TTree* tr_L00K00G60 = (TTree*) f_L00K00G60->Get("WJet");


  char* dijetPt = "sqrt(JetPFCor_Pt[0]**2+JetPFCor_Pt[1]**2+2*JetPFCor_Pt[0]*JetPFCor_Pt[1]*cos(JetPFCor_Phi[0]-JetPFCor_Phi[1]))";
  trSM->SetAlias("dijetPt", dijetPt);
  trSM2->SetAlias("dijetPt", dijetPt);
  tr_L01K00G00->SetAlias("dijetPt", dijetPt);
  tr_L03K00G00->SetAlias("dijetPt", dijetPt);
  tr_L05K00G00->SetAlias("dijetPt", dijetPt);
  tr_L07K00G00->SetAlias("dijetPt", dijetPt);
  tr_L09K00G00->SetAlias("dijetPt", dijetPt);
  tr_L11K00G00->SetAlias("dijetPt", dijetPt);
  tr_L00K05G00->SetAlias("dijetPt", dijetPt);
  tr_L00K11G00->SetAlias("dijetPt", dijetPt);
  tr_L00K16G00->SetAlias("dijetPt", dijetPt);
  tr_L00K20G00->SetAlias("dijetPt", dijetPt);
  tr_L00K00G11->SetAlias("dijetPt", dijetPt);
  tr_L00K00G40->SetAlias("dijetPt", dijetPt);
  tr_L00K00G60->SetAlias("dijetPt", dijetPt);


//   char* cutsNoBoost = "(effwt*puwt)*((W_pt<200.) && (dijetPt>70.) && (event_met_pfmet >25) &&(abs(W_muon_eta)<2.1) && (W_muon_pt>25.) && (abs(JetPFCor_Eta[0])<2.4) && (abs(JetPFCor_Eta[1])<2.4) && (abs(JetPFCor_Eta[0]-JetPFCor_Eta[1])<1.5) &&(abs(JetPFCor_dphiMET[0])>0.4) &&(W_mt>30.) &&(JetPFCor_Pt[0]>40.) &&(JetPFCor_Pt[1]>35.) &&(JetPFCor_Pt[2]<30.) &&(JetPFCor_bDiscriminatorCSV[0]<0.244) &&(JetPFCor_bDiscriminatorCSV[1]<0.244) && (Mass2j_PFCor>70. && Mass2j_PFCor<100.))";

//  char* cutsBoosted = "(effwt*puwt)*((W_pt>200.) && (event_met_pfmet >50) &&(abs(W_muon_eta)<2.1) && (W_muon_pt>30.) && (GroomedJet_CA8_pt[0]>200.) &&(abs(GroomedJet_CA8_eta[0])<2.4) &&(JetPFCor_bDiscriminatorCSV[0]<0.244) &&(GroomedJet_numberbjets <1)&&(ggdboostedWevt==1) && (GroomedJet_CA8_tau2tau1[0]<0.55) && (GroomedJet_CA8_mass[0]>70. && GroomedJet_CA8_mass[0]<100.))";



  char* cutsNoBoost = "(effwt*puwt)*((W_pt<200.) && (dijetPt>70.) && (event_met_pfmet >25) && (abs(JetPFCor_Eta[0])<2.4) && (abs(JetPFCor_Eta[1])<2.4) && (abs(JetPFCor_Eta[0]-JetPFCor_Eta[1])<1.5) &&(abs(JetPFCor_dphiMET[0])>0.4) &&(W_mt>30.) &&(JetPFCor_Pt[0]>40.) &&(JetPFCor_Pt[1]>35.) &&(JetPFCor_Pt[2]<30.) )";

 char* cutsBoosted = "(effwt*puwt)*((W_pt>200.) && (event_met_pfmet >50) && (GroomedJet_CA8_pt[0]>200.) &&(abs(GroomedJet_CA8_eta[0])<2.4) )";


  char* cutsNoBoostWt = "(effwt*puwt*weight)*((W_pt<200.) && (dijetPt>70.) && (event_met_pfmet >25) && (abs(JetPFCor_Eta[0])<2.4) && (abs(JetPFCor_Eta[1])<2.4) && (abs(JetPFCor_Eta[0]-JetPFCor_Eta[1])<1.5) &&(abs(JetPFCor_dphiMET[0])>0.4) &&(W_mt>30.) &&(JetPFCor_Pt[0]>40.) &&(JetPFCor_Pt[1]>35.) &&(JetPFCor_Pt[2]<30.) )";

 char* cutsBoostedWt = "(effwt*puwt*weight)*((W_pt>200.) && (event_met_pfmet >50) && (GroomedJet_CA8_pt[0]>200.) &&(abs(GroomedJet_CA8_eta[0])<2.4) )";


  ///// ----- Create all the histograms ----------
 double xbins[] = {100., 110., 120., 130., 140., 150., 160., 180., 200., 220., 240., 260., 300., 340., 
		   380., 420., 500., 600.};
  TH1D* hSMNoBoost = new TH1D("hSMNoBoost", "", sizeof(xbins)/ sizeof(xbins[0])-1, xbins);
  hSMNoBoost->Sumw2();
  hSMNoBoost->SetLineWidth(2);
  hSMNoBoost->GetXaxis()->SetTitle("dijet p_{T} (GeV)");
  hSMNoBoost->GetYaxis()->SetTitle("Ratio over SM");

  TH1D* hSMNoBoost2 = hSMNoBoost->Clone("hSMNoBoost2");
  TH1D* h_L01K00G00_NoBoost = hSMNoBoost->Clone("h_L01K00G00_NoBoost");
  TH1D* h_L03K00G00_NoBoost = hSMNoBoost->Clone("h_L03K00G00_NoBoost");
  TH1D* h_L05K00G00_NoBoost = hSMNoBoost->Clone("h_L05K00G00_NoBoost");
  TH1D* h_L07K00G00_NoBoost = hSMNoBoost->Clone("h_L07K00G00_NoBoost");
  TH1D* h_L09K00G00_NoBoost = hSMNoBoost->Clone("h_L09K00G00_NoBoost");
  TH1D* h_L11K00G00_NoBoost = hSMNoBoost->Clone("h_L11K00G00_NoBoost");
  TH1D* h_L00K05G00_NoBoost = hSMNoBoost->Clone("h_L00K05G00_NoBoost");
  TH1D* h_L00K11G00_NoBoost = hSMNoBoost->Clone("h_L00K11G00_NoBoost");
  TH1D* h_L00K16G00_NoBoost = hSMNoBoost->Clone("h_L00K16G00_NoBoost");
  TH1D* h_L00K20G00_NoBoost = hSMNoBoost->Clone("h_L00K20G00_NoBoost");
  TH1D* h_L00K00G11_NoBoost = hSMNoBoost->Clone("h_L00K00G11_NoBoost");
  TH1D* h_L00K00G40_NoBoost = hSMNoBoost->Clone("h_L00K00G40_NoBoost");
  TH1D* h_L00K00G60_NoBoost = hSMNoBoost->Clone("h_L00K00G60_NoBoost");

  TH1D* hSMBoosted = hSMNoBoost->Clone("hSMBoosted");
  TH1D* hSMBoosted2 = hSMNoBoost->Clone("hSMBoosted2");
  TH1D* h_L01K00G00_Boosted = hSMNoBoost->Clone("h_L01K00G00_Boosted");
  TH1D* h_L03K00G00_Boosted = hSMNoBoost->Clone("h_L03K00G00_Boosted");
  TH1D* h_L05K00G00_Boosted = hSMNoBoost->Clone("h_L05K00G00_Boosted");
  TH1D* h_L07K00G00_Boosted = hSMNoBoost->Clone("h_L07K00G00_Boosted");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    max = 150;
  }

  if(iPlot == 101){
    Nbin = 20;
    max = 2;
    min = -2;
  }

  // projection histogram
  TH1D *h = new TH1D(Form("h",cbin),"",Nbin,min,max);
  TH1D *hReference = new TH1D(Form("hReference",cbin),"",Nbin,min,max);
  TH1D *hDataMix = new TH1D(Form("hDataMix",cbin),"",Nbin,min,max);

  TH1D *hB = new TH1D(Form("hB",cbin),"",Nbin,min,max);
  TH1D *hReferenceB = new TH1D(Form("hReferenceB",cbin),"",Nbin,min,max);
  TH1D *hDataMixB = new TH1D(Form("hDataMixB",cbin),"",Nbin,min,max);

  TH1D *hFull = new TH1D("hFull","",Nbin,min,max);
  TH1D *hReferenceFull = new TH1D("hReferenceFull","",Nbin,min,max);
  TH1D *hDataMixFull = new TH1D("hDataMixFull","",Nbin,min,max);

  TH1D* hNorm = new TH1D("hNorm","",1000,0,1000);
  TH1D* hNormReference = new TH1D("hNormReference","",1000,0,1000);
  TH1D* hNormDataMix = new TH1D("hNormDataMix","",1000,0,1000);

  hB->SetLineStyle(2);
  hReferenceB->SetLineStyle(2);
  hDataMixB->SetLineStyle(2);

  //  ntReference->SetAlias("pt1","et1");
  //  ntReference->SetAlias("pt2","et2");

  nt->SetAlias("pt1","jtpt1");
  nt->SetAlias("pt2","jtpt2");
  nt->SetAlias("eta1","jteta1");
  nt->SetAlias("eta2","jteta2");
  nt->SetAlias("phi1","jtphi1");
  nt->SetAlias("phi2","jtphi2");
  ntMix->SetAlias("pt1","jtpt1");
  ntMix->SetAlias("pt2","jtpt2");
  ntMix->SetAlias("eta1","jteta1");
  ntMix->SetAlias("eta2","jteta2");
  ntMix->SetAlias("phi1","jtphi1");
  ntMix->SetAlias("phi2","jtphi2");
  ntReference->SetAlias("pt1","jtpt1");
  ntReference->SetAlias("pt2","jtpt2");
  ntReference->SetAlias("eta1","jteta1");
  ntReference->SetAlias("eta2","jteta2");
  ntReference->SetAlias("phi1","jtphi1");
  ntReference->SetAlias("phi2","jtphi2");

  nt->SetAlias("cent","hfp");
  ntMix->SetAlias("cent","hfp");
  ntReference->SetAlias("cent","hfp");

  nt->SetAlias("dphi","acos(cos(phi1-phi2))");
  ntMix->SetAlias("dphi","acos(cos(phi1-phi2))");
  ntReference->SetAlias("dphi","acos(cos(phi1-phi2))");

  if(iPlot == 0){
    nt->SetAlias("var","pt2/pt1");
    ntReference->SetAlias("var","pt2/pt1");
    ntMix->SetAlias("var","pt2/pt1");
  }

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
		    )
{
   // ===========================================================
   // Get Input
   // ===========================================================
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get("ntjt");
   t->SetAlias("metxMergedAll","(metx0+metx1+metx2+metx3+metx4+metx5)");
   t->SetAlias("metxMerged0","(metx0+metx1+metx2+metx3)");
   t->SetAlias("metxMerged1","(metx1+metx2+metx3)");
   t->SetAlias("metxMerged2","(metx2+metx3)");
   t->SetAlias("metxMerged3","(metx3)");
   t->SetAlias("metxMerged4","(metx4+metx5)");
   //TString name2(infname);
   //name2.ReplaceAll("ntv2deadforpix_","");
   //name2.ReplaceAll("ntv2_","");
   //t->AddFriend("dj=djTree",name2);

   // ===========================================================
   // Analysis Setup
   // ===========================================================
   Int_t plotLayer=10;
   const int nBin = 5;
   double bins[nBin+1] = {0.5,1.0,2,4,8,1000};  
   
   const int nBinAj = 4;
   double ajBins[nBinAj+1] = {0.0001,0.13,0.24,0.35,0.69999};
   // Selection cut
   TCut evtCut = "nljet>100&&abs(nljetacorr)<1.6&&aljet>40&&abs(aljetacorr)<1.6&&jdphi>2./3*TMath::Pi()&&!maskEvt";

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

   // ===========================================================
   // Find Average Weights
   // ===========================================================
   TH1D *hw[nBinAj];
   float meanWt[nBinAj];
   for ( int iaj = 0 ; iaj< nBinAj ; iaj++) {
      hw[iaj] = new TH1D(Form("hw_aj%d",iaj),"",1000,0,100);
      TCut ajCut = Form("Aj>%f && Aj<%f", ajBins[iaj],ajBins[iaj+1]);
      t->Draw(Form("weight>>hw_aj%d",iaj), evtCut&&myCut&&ajCut);
      meanWt[iaj] = hw[iaj]->GetMean();
      cout << " <Weight>: " << meanWt[iaj] << endl;
   }

   // ===========================================================
   // Draw Weighted Averages
   // ===========================================================
   TH1D *p[nBin];
   TH1D *pe[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);
      TH1D *he[nBinAj];
      
      // =================================
      // Get Weighted Mean for each Aj bin
      // =================================
      h1->Sumw2();
      h2->Sumw2();
      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);

      // =================================
      // Caculated Stat Error of the Mean
      // =================================
      cout << "Stat Error for pt bin " << i << ": ";
      for ( int iaj = 0 ; iaj< nBinAj ; iaj++) {
         he[iaj] = new TH1D(Form("he%d_aj%d",i,iaj),"",100,-200,200);
         TCut ajCut = Form("Aj>%f && Aj<%f", ajBins[iaj],ajBins[iaj+1]);
         t->Draw(Form("((metxMerged%d))>>he%d_aj%d",i,i,iaj), "weight" * evtCut&&myCut&&ajCut);
         float theError = he[iaj]->GetRMS()/ (sqrt(he[iaj]->GetEntries()));
	 cout << theError << " ";
	 p[i]->SetBinError(iaj+1, theError);
      }
      cout << endl;
   }
   // correct the error bars..  cummulative -> it's own.
   subtractError(p[0],p[1]);
   subtractError(p[1],p[2]);
   subtractError(p[2],p[3]);
   
   TH1D *pall;
   TH1D *pallE;
   TH1D *h1 = new TH1D(Form("hAll1"),"",nBinAj,ajBins);
   TH1D *h2 = new TH1D(Form("hAll2"),"",nBinAj,ajBins);
//.........这里部分代码省略.........