C++ TChain::AddFriend方法代码示例

double* getPscls(std::string infile, int nFiles) {

    TChain *dataCH = NULL;
    dataCH = new TChain(Form("%sJetAnalyzer/t", algo.Data()));
    TChain *dataCH2 = new TChain("hltanalysis/HltTree");
    std::ifstream instr(infile.c_str(), std::ifstream::in);
    std::string filename;
    for(int ifile=0; ifile<nFiles; ifile++) {
        instr >> filename;
        dataCH->Add(filename.c_str());
        dataCH2->Add(filename.c_str());
    }
    dataCH->AddFriend(dataCH2, "hltanalysis/HltTree");
    double ov1, ov2, ov3, ov4;
    ov1 = dataCH->GetEntries("HLT_PAJet20_NoJetID_v1 && HLT_PAJet80_NoJetID_v1");
    ov2 = dataCH->GetEntries("HLT_PAJet40_NoJetID_v1 && HLT_PAJet80_NoJetID_v1");
    ov3 = dataCH->GetEntries("HLT_PAJet60_NoJetID_v1 && HLT_PAJet80_NoJetID_v1");
    ov4 = dataCH->GetEntries("HLT_PAJet80_NoJetID_v1");
    double *pscls = new double[4];
    pscls[0] = ov4/ov1;
    pscls[1] = ov4/ov2;
    pscls[2] = ov4/ov3;
    pscls[3] = 1.;
    return pscls;
}

void scanDiffRecoTracks(std::string eName, std::string eoName,
			std::string branchReg = "recoTracks_*"){
  gSystem->Load("libFWCoreFWLite");
  gROOT->ProcessLine("AutoLibraryLoader::enable();");
  TChain* e = new TChain("Events");
  e->SetScanField(0);

  e->Add(eName.c_str());

  TChain* eo = new TChain("Events");
  eo->SetScanField(0);

  eo->Add(eoName.c_str());
  e->AddFriend(eo, "eo");

  TRegexp regg(branchReg.c_str(), kTRUE);

  TChain* tc = e ;
  TObjArray* tl = tc->GetListOfBranches();
 
  Int_t nBr = tl->GetSize();
  for (int iB=0;iB<nBr;++iB){
    TBranch* br = (TBranch*)(tl->At(iB)); 
    TString ts(br->GetName()); 
    if(ts.Index(regg)>=0){
      std::cout<<ts.Data()<<std::endl;
      tc->Scan(Form("Sum$(%s.obj.pt()>0):Sum$(eo.%s.obj.pt()>0):Sum$(%s.obj.pt()):Sum$(%s.obj.pt())-Sum$(eo.%s.obj.pt())",
 		    ts.Data(), ts.Data(), ts.Data(), ts.Data(), ts.Data()),"", "");
    }
  } //> e.tecoTracks.recoT0.txt

}

TChain* load (std::string methodName)
{
    TString trainingFileName ("training");
    TString fname0 = pathToData + trainingFileName + TString (".root");
    TString fname1 = trainingFileName + TString ("_prediction__") + methodName + TString (".root");
    TChain* ch = new TChain ("data");
    ch->Add (fname0);
    TChain* fr = new TChain ("data");
    fr->Add (fname1);
    ch->AddFriend (fr,"s");
    ch->Print ("");
    return ch;
}

int main() 
{
   MySelector *A = new MySelector();

   TChain* ch = new TChain("AK4PFCHS/t");
	TChain* c1 = new TChain("event/t");

 	FILE *input;
	char filename[300];

	int i = 0;

	TString filenames [] = {"../lists/file_list_25ns_76X.txt",
									"../lists/file_list_25ns_76X_test.txt",
									"../lists/ptbinned/file_list_5to10_76X.txt",
									"../lists/ptbinned/file_list_10to15_76X.txt",
									"../lists/ptbinned/file_list_15to30_76X.txt",
									"../lists/ptbinned/file_list_30to50_76X.txt",
									"../lists/ptbinned/file_list_50to80_76X.txt",
									"../lists/ptbinned/file_list_80to120_76X.txt",
									"../lists/ptbinned/file_list_120to170_76X.txt",
									"../lists/ptbinned/file_list_170to300_76X.txt",
									"../lists/ptbinned/file_list_300to470_76X.txt",
									"../lists/ptbinned/file_list_470to600_76X.txt",
									"../lists/ptbinned/file_list_600to800_76X.txt",
									"../lists/ptbinned/file_list_800to1000_76X.txt",
									"../lists/ptbinned/file_list_1000to1400_76X.txt",
									"../lists/ptbinned/file_list_1400to1800_76X.txt",
									"../lists/ptbinned/file_list_1800to2400_76X.txt",
									"../lists/ptbinned/file_list_2400to3200_76X.txt",
									"../lists/ptbinned/file_list_3200toInf_76X.txt"};

	input = fopen( filenames[i], "r" );

	if (input != NULL)
	{ // lets read each line and get the filename from it
		while (fscanf(input,"%s\n",filename) != EOF) 
		{
			printf("%s\n",filename); 
			ch->Add(filename);
			ch->AddFriend(c1,filename);
		}
	}

//	std::cout << "number of events is " << ch->GetEntries() << std::endl;

	ch -> Process(A);

//	std::cout <<"			Analyzed events #" <<  A -> TotalEvents << std::endl;

}

//_____________________________________________________________________________
TChain *CreateAODFriendChain(TString sDataset)
{ 
  if (!sDataset.EndsWith(".txt")) return 0;

  TChain *chain = new TChain("aodTree");
  TChain *cFrid = new TChain("aodTree");

  TString dataFile;
  ifstream dataList(sDataset.Data(), ios::in); 
  while (!dataList.eof()) {
    dataFile.ReadLine(dataList,kFALSE);
    if (!dataFile.EndsWith("AliAOD.root")) continue;
    if (!gSystem->AccessPathName(dataFile.Data())) chain->Add(dataFile.Data());

    dataFile.ReplaceAll("AliAOD.root","AliAOD.VertexingHF.root");
    if (!gSystem->AccessPathName(dataFile.Data())) cFrid->Add(dataFile.Data());
  } dataList.close();

  chain->AddFriend(cFrid);
  return chain;
}

TChain * compJOO(
    TString infile="hltana_Exp_151058.root",
    TString header="HLT_HIJet90U"
    )
{
  TChain * dj = new TChain("hltanalysis/HltTree");

  dj->Add(infile);
  dj->AddFriend("onl = hltanalysisOnl/HltTree",infile);
  cout << "Trig Total: " << dj->GetEntries() << endl;

  TCut evtSel("HLT_HIJet90U && hiBin*2.5<100");

  TH1D * hJEtIc5pu = plotJEt(dj,evtSel,"recoJetCalPt[0]","hJEtIc5pu",0);
  TH1D * hJEtIc5puOnl = plotJEt(dj,evtSel,"onl.recoJetCalPt[0]","hJEtIc5puOnl",0);

  TCut matEvtSel(evtSel);
  bool doMat = true;
  if (doMat) matEvtSel = evtSel && "abs(onl.recoJetCalEta[0]-recoJetCalEta[0])<0.2";
  TH2D * hJEtOnlIc5pu = plotJEtCorr(dj,matEvtSel,"onl.recoJetCalPt[0]/recoJetCalPt[0]:recoJetCalPt[0]","hJEtOnlIc5pu",";leading E_{T}^{Ic5pu} [GeV];E_{T}^{Ic5pu Online}/E_{T}^{Ic5pu}");

  TCanvas * cJEt = new TCanvas("cJEt","cJEt",500,500);
  CPlot cpJEt("JEt","JEt","leading E_{T}^{Jet}","#");
  cpJEt.SetLogy();
  cpJEt.AddHist1D(hJEtIc5puOnl,"Online Reco","hist",kGreen+2,kOpenSquare);
  cpJEt.AddHist1D(hJEtIc5pu,"Offline Reco","E",kBlack,kFullCircle);
  cpJEt.SetLegend(0.62,0.78,0.94,0.9);
  cpJEt.SetLegendHeader(header);
  cpJEt.Draw(cJEt,false);

  TCanvas * cDEta = new TCanvas("cDEta","cDEta",500,500);
  cDEta->SetLogy();
  dj->Draw("onl.recoJetCalEta[0]-recoJetCalEta[0]>>hDEta",evtSel,"E");
  TCanvas * cDPhi = new TCanvas("cDPhi","cDPhi",500,500);
  cDPhi->SetLogy();
  dj->Draw("onl.recoJetCalPhi[0]-recoJetCalPhi[0]>>hDPhi",evtSel,"E");

  return dj;
}

void BkgFluctuation(TString infname = "../ntv6_datambv3rand_djcalo_c0to30_offset0.root",
                    TCut myCut = "cent<30", char *title = "",bool drawLegend = false,
                    bool drawSys = true
		    )
{
  gStyle->SetHistLineWidth(2);
  // ===========================================================
  // Get Input
  // ===========================================================
  TChain* t = new TChain("ntjt");
  t->Add(infname);
  t->AddFriend("tcone",infname);
  cout << infname << ": " << t->GetEntries() << endl;
  t->SetAlias("cptsub0Merge0","cpt[0][1]+cpt[0][2]-(cptbg[0][1]+cptbg[0][2])");
  t->SetAlias("cptsub0Merge1","cpt[0][3]-(cptbg[0][3])");
  t->SetAlias("cptsub0Merge2","cpt[0][4]+cpt[0][5]-(cptbg[0][4]+cptbg[0][5])");
  t->SetAlias("cptsub0Merge3","(Sum$(cpt[0])-cpt[0][0])-(Sum$(cptbg[0])-cptbg[0][0])");

  // ===========================================================
  // Analysis Setup
  // ===========================================================
  TH1::SetDefaultSumw2();
   const int nBin = 3;
   double bins[nBin+1] = {1.,4,8,160};  
  TCut evtCut = "nljet>100&&abs(nljetacorr)<1.6&&aljet>50&&abs(aljetacorr)<1.6&&jdphi>2./3*TMath::Pi()&&!maskEvt";
  TCut exclusion = "abs(nljetacorr)>0.8"; //&&abs(aljetacorr)>0.8";
  evtCut = evtCut&&exclusion;
  //Int_t colors[7] = {kBlue, kCyan-3, kYellow+1, kOrange+1, kGreen+2, kMagenta, kRed };
  Int_t colors[7] = {kYellow+2, kGreen+2,kRed };
  
  // ===========================================================
  // Book Histograms
  // ===========================================================
  vector<TH1D*> vhBgSub;
  for (Int_t i=0; i<nBin+1; ++i) {
    vhBgSub.push_back(new TH1D(Form("hBgSub_%d"),";(Sig Cone #Sigma p_{T}) - (Bkg Cone #Sigma p_{T}) (GeV/c);fraction;",50,-300,300));
  }

  Float_t numEvt = t->GetEntries(evtCut&&myCut);
  cout << "Total Sel evt: " << numEvt << endl;
  for (Int_t i=0; i<nBin+1; ++i) {
    TString var(Form("cptsub0Merge%i",i));
    cout << "Sel evt: " << t->GetEntries(evtCut&&myCut&&var) << endl;
    cout << "var: " << t->GetAlias(var) << endl;
    t->Project(vhBgSub[i]->GetName(),var,(evtCut&&myCut&&var)*"weight");
    Float_t mean = vhBgSub[i]->GetMean();
    Float_t rms = vhBgSub[i]->GetRMS();
    cout << "bgsub mean: " << mean << " rms: " << rms << endl;
    vhBgSub[i]->Scale(1./numEvt);
    cout << "Bin " << i << " rms/sqrt(250): " << rms/sqrt(250)
      << " r=0: " << rms/sqrt(250.)*sqrt(fracArea(0))
      << " r=0.8: " << rms/sqrt(250.)*sqrt(fracArea(0.8)) << endl;
  }

  TCanvas * c2 = new TCanvas("c2","c2",600,600);
  c2->SetLogy();
  vhBgSub[nBin]->SetAxisRange(1e-5,5,"Y");
  vhBgSub[nBin]->Draw("E");
  for (Int_t i=0; i<nBin; ++i) {
    vhBgSub[i]->SetLineColor(colors[i]);
    vhBgSub[i]->Draw("same hist");
  }
  vhBgSub[nBin]->Draw("Esame");

  TLegend *leg = new TLegend(0.19,0.70,0.53,0.92);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->SetTextSize(0.035);
  leg->AddEntry(vhBgSub[nBin],"CMS 0-30%","");
  leg->AddEntry(vhBgSub[nBin],"Random Cone (MB Sample)","");
  leg->AddEntry(vhBgSub[nBin],Form("> %.1f GeV/c",bins[0]),"pl");
  for (int i=0;i<nBin;++i) {
    if (i!=nBin-1){
      leg->AddEntry(vhBgSub[i],Form("%.1f - %.1f GeV/c",bins[i],bins[i+1]),"l");
    } else {
      leg->AddEntry(vhBgSub[i],Form("> %.1f GeV/c",bins[i]),"l");
    }
  }
  leg->Draw();
}

void RateEstimateVpt(TString infile="alldatacleaned.root") //Contains both cleaned and uncleaned data
{
  //TChain * tree = new TChain("icPu5JetAnalyzer_clean/t");
  TChain * tree = new TChain("icPu5JetAnalyzer/t");
  tree->Add(infile);
  //tree->AddFriend("hltanalysis_clean/HltTree",infile);
  tree->AddFriend("hltanalysis/HltTree",infile);
  //  cout << "Total: " << tree->GetEntries() << endl; 

  TCut evtSel("Run==152957&&HLT_HIMinBiasHfOrBSC&&hiBin<40&&abs(jteta[0])<2");


  //Canvas definition and choosing uncorrected and corrected jet pt variables 
  TCanvas *e1 = new TCanvas("e1","Rate Estimate vs. Threshold",200,10,700,500);
  e1->SetLogy();
  TH1D * RJetPt = new TH1D("RJetPt",";Leading Jet p_{T} (GeV/c);count;",300,0,300);
  TH1D * JetPt = new TH1D("JetPt",";Leading Jet p_{T} (GeV/c);count;",300,0,300);
  tree->Draw("rawpt[0]>>RJetPt",evtSel,"goff");
  tree->Draw("jtpt[0]>>JetPt",evtSel,"goff");
  cout << "Integral: " << JetPt->Integral(0,101) << endl;

  //Defines total # of entries for normalization and max rate
  TH1F * h4 = new TH1F("h4","",50,0,300);
  tree->Draw("LumiBlock>>h4",evtSel,"goff");
  Float_t R = 1 ; //Defines Maximum Rate
  Float_t TA = h4->GetEntries();

  //Computes Accept Fraction as a function of pt 
  const Int_t n = 61;
  Float_t x[n], y[n], z[n], w[n], t[n];
  for (Int_t i=0;i<n;i++) {
    x[i] = 300-i*5; 
    y[i] = RJetPt->Integral(x[i],301);
    z[i] = R*y[i]/TA;
    w[i] = JetPt->Integral(x[i],301);
    t[i] = R*w[i]/TA;
    printf(" i %i %f %f %f \n",i,x[i],y[i],z[i]);
  }

  //Defines graphs and axis. 
  gr = new TGraph(n,x,z);
  gs = new TGraph(n,x,t);
  gr->GetYaxis()->SetRangeUser(0.001,5000);

  // Attractive features for graphs
  gr->SetLineColor(kRed);
  gs->SetLineColor(kBlue);
  gr->SetMarkerColor(kRed);
  gr->SetMarkerStyle(kOpenCircle);
  gs->SetMarkerColor(kBlue);
  gs->SetMarkerStyle(kOpenCircle);
  gr->SetTitle("");
  gr->GetXaxis()->SetTitle("Jet Pt Threshold (GeV/c)");
  gr->GetYaxis()->SetTitle("Rate Estimate (Hz)");

  //Draw graphs
  gr->Draw("LAP");
  gs->Draw("LP");

  //Defines Legened
  TLegend *L = new TLegend(0.52292,0.675103,0.932471,0.919831);
  L->SetHeader("AllPhysics, |#eta| < 2, All Centrality" );
  L->SetBorderSize(0);
  L->SetFillStyle(0);
  L->AddEntry(gr,"Uncorrected p_{T}","pl");
  L->AddEntry(gs,"Corrected p_{T}","pl");
  L->Draw();

}

void analyzePuppi(std::vector<std::string> urls, const char *outname = "eventObjects.root", Long64_t nentries = 20, Int_t firstF = -1, Int_t lastF = -1) {

  // std::vector<std::string> urls = CollectFiles(list);

  // Printf("anaFile: %d",anaFile);
  
  std::cout << "nfiles: " << urls.size() << std::endl;
  for (auto i = urls.begin(); i != urls.end(); ++i)
    std::cout << *i << std::endl;

  size_t firstFile = 0;
  size_t lastFile = urls.size();

  if(firstF>-1) {
    firstFile = (size_t)firstF;
    lastFile = (size_t)lastF;
  }
  std::cout << "firstFile: " << firstFile << "  lastFile: " << lastFile << std::endl;
  
  //add files to chain
  TChain *chain = NULL;
  chain = new TChain("hiEvtAnalyzer/HiTree");
  for(size_t i=firstFile; i<lastFile; i++) chain->Add(urls[i].c_str());
  Printf("hiTree done");
  
  TChain *pfTree = new TChain("pfcandAnalyzer/pfTree");
  for(size_t i=firstFile; i<lastFile; i++) pfTree->Add(urls[i].c_str());
  chain->AddFriend(pfTree);
  Printf("pfTree done");
  
  // TChain *muTree = new TChain("hltMuTree/HLTMuTree");
  // for(size_t i=firstFile; i<lastFile; i++) muTree->Add(urls[i].c_str());
  // chain->AddFriend(muTree);
  // Printf("muTree done");

  //TChain *jetTree = new TChain("akPu3PFJetAnalyzer/t");
  TChain *jetTree = new TChain("akPu2PFJetAnalyzer10/t");
  for(size_t i=firstFile; i<lastFile; i++) jetTree->Add(urls[i].c_str());
  chain->AddFriend(jetTree);
  Printf("jetTree done");

  TChain *genTree = new TChain("HiGenParticleAna/hi");
  for(size_t i=firstFile; i<lastFile; i++) genTree->Add(urls[i].c_str());
  chain->AddFriend(genTree);
  Printf("genTree done");
  
  TList *fEventObjects = new TList();


  //---------------------------------------------------------------
  // producers
  //
  hiEventProducer *p_evt = new hiEventProducer("hiEvtProd");
  p_evt->SetInput(chain);
  p_evt->SetHIEventContName("hiEventContainer");
  p_evt->SetEventObjects(fEventObjects);

  pfParticleProducer *p_pf = new pfParticleProducer("pfPartProd");
  p_pf->SetInput(chain);
  p_pf->SetpfParticlesName("pfParticles");
  p_pf->SetEventObjects(fEventObjects);

  genParticleProducer *p_gen = new genParticleProducer("genParticleProd");
  p_gen->SetInput(chain);
  p_gen->SetGenParticlesName("genParticles");
  p_gen->SetEventObjects(fEventObjects);

  lwJetFromForestProducer *p_PUJet = new lwJetFromForestProducer("lwJetForestProd");
  p_PUJet->SetInput(chain);
  p_PUJet->SetJetContName("aktPUR030");
  p_PUJet->SetEventObjects(fEventObjects);
  p_PUJet->SetRadius(0.3);
  
  //---------------------------------------------------------------
  //analysis modules
  //

  //handler to which all modules will be added
  anaBaseTask *handler = new anaBaseTask("handler","handler");

  //analysis modules which also produce
  anaPuppiProducer *pupProd = new anaPuppiProducer("pupProd","pupProd");
  pupProd->ConnectEventObject(fEventObjects);
  pupProd->SetHiEvtName("hiEventContainer");
  pupProd->SetPFPartName("pfParticles");
  pupProd->SetPuppiPartName("puppiParticles");
  pupProd->SetJetsName("aktPUR030");
  pupProd->SetAddMetricType(anaPuppiProducer::kMass);
  if(doPuppi) handler->Add(pupProd);

  //anti-kt jet finder on reconstructed PUPPI particles ptmin=0
  LWJetProducer *lwjakt = new LWJetProducer("LWJetProducerAKTR030Puppi","LWJetProducerAKTR030Puppi");
  lwjakt->ConnectEventObject(fEventObjects);
  lwjakt->SetJetType(LWJetProducer::kAKT);
  lwjakt->SetRadius(0.3);
  lwjakt->SetGhostArea(0.005);
  lwjakt->SetPtMinConst(0.);
  lwjakt->SetParticlesName("puppiParticles");
  lwjakt->SetJetContName("JetsAKTR030Puppi");
  if(doPuppi) handler->Add(lwjakt);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    Int_t count = 0;

    for (Int_t iDir=0; iDir<nDirs; ++iDir)
    {
      TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
      if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
        continue;

      if (offset > 0)
      {
        --offset;
        continue;
      }

      if (count++ == aRuns)
        break;

      TString presentDirName(aDataDir);
      presentDirName += "/";
      presentDirName += presentDir->GetName();

      chain->Add(presentDirName + "/AliESDs.root/esdTree");
    }
  }
  else
  {
    // Open the input stream
    ifstream in;
    in.open(aDataDir);

    ofstream outfile;
    if (check)
      outfile.open(check);

    Int_t count = 0;

    // Read the input list of files and add them to the chain
    TString line;
    while (in.good())
    {
      in >> line;

      if (line.Length() == 0)
        continue;

      if (offset > 0)
      {
        offset--;
        continue;
      }

      if (count++ == aRuns)
        break;

      TString esdFile(line);

      if (addFileName)
        esdFile += "/AliESDs.root";
        
      TString esdFileFriend(esdFile);
      esdFileFriend.ReplaceAll("AliESDs.root", "AliESDfriends.root");
        
      if (check)
      {
        TFile* file = TFile::Open(esdFile);
        if (!file)
          continue;
        file->Close();
        
        if (chainFriend)
        {
          TFile* file = TFile::Open(esdFileFriend);
          if (!file)
            continue;
          file->Close();
        }
        
        outfile << line.Data() << endl;
        printf("%s\n", line.Data());
      }        
        
        // add esd file
      chain->Add(esdFile);

        // add file
      if (chainFriend)
        chainFriend->Add(esdFileFriend);
    }

    in.close();
    
    if (check)
      outfile.close();
  }
  
  if (chainFriend)
    chain->AddFriend(chainFriend);

  return chain;
}

void analyzeJES(std::vector<std::string> urls, const char *outname = "eventObjects.root", Long64_t nentries = 20, Int_t firstF = -1, Int_t lastF = -1, Int_t firstEvent = 0) {

  /*
    ptminType: minimum raw pt for particles used by puppi
    0 : 0 GeV
    1 : 1 GeV
    2 : 2 GeV

    jetSignalType: jets used to select jetty region in events
    0 : detector-level jets (akPu4PF)
    1 : particle-level jets (gen jets)
   */


  TString jetName = "aktPuppiR040";
  TString jetTreeName = "akPuppi4PFJetAnalyzer";
  jetName = "aktCsR040";
  jetTreeName = "akCs4PFJetAnalyzer";
  //jetName = "aktVsR040";
  //jetTreeName = "akVs4PFJetAnalyzer";

  std::cout << "analyzing JES for: " << jetName << " tree: " << jetTreeName << std::endl;
   
  std::cout << "nfiles: " << urls.size() << std::endl;
  for (auto i = urls.begin(); i != urls.end(); ++i)
    std::cout << *i << std::endl;

  size_t firstFile = 0;
  size_t lastFile = urls.size();

  if(firstF>-1) {
    firstFile = (size_t)firstF;
    lastFile = std::min((size_t)lastF,lastFile);
  }
  std::cout << "firstFile: " << firstFile << "  lastFile: " << lastFile << std::endl;

  Int_t lastEvent = nentries;
  if(firstEvent>0) {
    lastEvent = firstEvent + nentries;
  }
  std::cout << "firstEvent: " << firstEvent << std::endl;
  
  //add files to chain
  TChain *chain = NULL;
  chain = new TChain("hiEvtAnalyzer/HiTree");
  for(size_t i=firstFile; i<lastFile; i++) chain->Add(urls[i].c_str());
  Printf("hiTree done");
  
  TChain *jetTree = new TChain(Form("%s/t",jetTreeName.Data()));
  for(size_t i=firstFile; i<lastFile; i++) jetTree->Add(urls[i].c_str());
  chain->AddFriend(jetTree);
  Printf("jetTree done");

  TList *fEventObjects = new TList();

  //---------------------------------------------------------------
  // producers
  //
  hiEventProducer *p_evt = new hiEventProducer("hiEvtProd");
  p_evt->SetInput(chain);
  p_evt->SetHIEventContName("hiEventContainer");
  p_evt->SetEventObjects(fEventObjects);

  lwJetFromForestProducer *p_PUJet = new lwJetFromForestProducer("lwJetForestProd");
  p_PUJet->SetInput(chain);
  p_PUJet->SetJetContName(jetName);
  p_PUJet->SetGenJetContName("akt4Gen");
  p_PUJet->SetEventObjects(fEventObjects);
  p_PUJet->SetRadius(0.4);
  
  //---------------------------------------------------------------
  //analysis modules
  //

  //handler to which all modules will be added
  anaBaseTask *handler = new anaBaseTask("handler","handler");

  anaJetEnergyScale *anajesForest = new anaJetEnergyScale("anaJESForest","anaJESForest");
  anajesForest->ConnectEventObject(fEventObjects);
  anajesForest->SetHiEvtName("hiEventContainer");
  anajesForest->SetGenJetsName("akt4Gen");
  anajesForest->SetRecJetsName(jetName);
  anajesForest->SetNCentBins(4);
  anajesForest->SetUseForestMatching(true);
  //anajesForest->SetMaxDistance(0.2);
  handler->Add(anajesForest);

  anaJetEnergyScale *anajesForestQuarks = new anaJetEnergyScale("anaJESForestQuarks","anaJESForestQuarks");
  anajesForestQuarks->ConnectEventObject(fEventObjects);
  anajesForestQuarks->SetHiEvtName("hiEventContainer");
  anajesForestQuarks->SetGenJetsName("akt4Gen");
  anajesForestQuarks->SetRecJetsName(jetName);
  anajesForestQuarks->SetNCentBins(4);
  anajesForestQuarks->SetUseForestMatching(true);
  //anajesForestQuarks->SetMaxDistance(0.2);
  anajesForestQuarks->SetRefPartonFlavor(0,2);
  handler->Add(anajesForestQuarks);

  anaJetEnergyScale *anajesForestGluons = new anaJetEnergyScale("anaJESForestGluons","anaJESForestGluons");
  anajesForestGluons->ConnectEventObject(fEventObjects);
//.........这里部分代码省略.........

void analyzeZJetMCResponse(std::vector<std::string> urls, const char *outname = "eventObjects.root", Long64_t nentries = 20, Int_t firstF = -1, Int_t lastF = -1, Int_t firstEvent = 0) {

  TString jetName = "aktPuppiR040";
  TString jetTreeName = "akPuppi4PFJetAnalyzer";
  jetName = "aktPuR030"; //"aktCsR040";
  jetTreeName = "akPu3PFJetAnalyzer"; //"akCs4PFJetAnalyzer"; //akCs4PFJetAnalyzer

  std::cout << "analyzing Z-jet response for: " << jetName << " tree: " << jetTreeName << std::endl;
  
  std::cout << "nfiles: " << urls.size() << std::endl;
  for (auto i = urls.begin(); i != urls.end(); ++i)
    std::cout << *i << std::endl;

  size_t firstFile = 0;
  size_t lastFile = urls.size();

  if(firstF>-1) {
    firstFile = (size_t)firstF;
    lastFile = (size_t)lastF;
  }
  std::cout << "firstFile: " << firstFile << "  lastFile: " << lastFile << std::endl;

  Int_t lastEvent = nentries;
  if(firstEvent>0) {
    lastEvent = firstEvent + nentries;
  }
  
  //add files to chain
  TChain *chain = NULL;
  chain = new TChain("hiEvtAnalyzer/HiTree");
  for(size_t i=firstFile; i<lastFile; i++) chain->Add(urls[i].c_str());
  Printf("hiTree done");

  // TChain *hltTree = new TChain("hltanalysis/HltTree");
  // for(size_t i=firstFile; i<lastFile; i++) hltTree->Add(urls[i].c_str());
  // chain->AddFriend(hltTree);
  // Printf("hltTree done");

  // TChain *skimTree = new TChain("skimanalysis/HltTree");
  // for(size_t i=firstFile; i<lastFile; i++) skimTree->Add(urls[i].c_str());
  // chain->AddFriend(skimTree);
  // Printf("skimTree done");

  // TChain *muTree = new TChain("ggHiNtuplizer/EventTree");
  // for(size_t i=firstFile; i<lastFile; i++) muTree->Add(urls[i].c_str());
  // chain->AddFriend(muTree);
  // Printf("muTree done");

  TChain *genTree = new TChain("HiGenParticleAna/hi");
  for(size_t i=firstFile; i<lastFile; i++) genTree->Add(urls[i].c_str());
  //chain->AddFriend(genTree);
  Printf("genTree done: %d",(int)genTree->GetEntries());

  TChain *jetTree = new TChain(Form("%s/t",jetTreeName.Data()));
  for(size_t i=firstFile; i<lastFile; i++) jetTree->Add(urls[i].c_str());
  chain->AddFriend(jetTree);
  Printf("jetTree done");

  TList *fEventObjects = new TList();

  //---------------------------------------------------------------
  // producers
  //
  hiEventProducer *p_evt = new hiEventProducer("hiEvtProd");
  p_evt->SetInput(chain);
  p_evt->SetHIEventContName("hiEventContainer");
  p_evt->SetEventObjects(fEventObjects);

  lwMuonProducer *p_mu = new lwMuonProducer("lwMuonProd");
  p_mu->SetInput(chain);
  p_mu->SetlwMuonsRecoName("lwMuonsReco");
  p_mu->SetlwMuonsGeneName("");
  p_mu->SetEventObjects(fEventObjects);

  genParticleProducer *p_gen = new genParticleProducer("genParticleProd");
  p_gen->SetInput(genTree);
  p_gen->SetGenParticlesName("genParticles");
  p_gen->SetEventObjects(fEventObjects);

  lwJetFromForestProducer *p_PUJet = new lwJetFromForestProducer("lwJetForestProd");
  p_PUJet->SetInput(chain);
  p_PUJet->SetJetContName(jetName);
  p_PUJet->SetGenJetContName("");//akt4Gen");
  p_PUJet->SetEventObjects(fEventObjects);
  p_PUJet->SetRadius(0.4);

  //---------------------------------------------------------------
  //analysis modules
  //

  //handler to which all modules will be added
  anaBaseTask *handler = new anaBaseTask("handler","handler");

  //Z to mumu
  anaZToMuMu *ZToMuMu = new anaZToMuMu("ZToMuMu","ZToMuMu");
  ZToMuMu->ConnectEventObject(fEventObjects);
  ZToMuMu->SetHiEvtName("hiEventContainer");
  ZToMuMu->SetMuonsName("genParticles");
  ZToMuMu->SetCheckPid(true);
  ZToMuMu->SetZsName("zMuMuBosons");
//.........这里部分代码省略.........

void plotPFRandomCone(TString str = "forest.root", Int_t nRCPerEvent = 1, int maxEvents = -1) {

  // gStyle->SetOptStat(0000);
  // gStyle->SetOptTitle(0);

  double minEta = -1.5;
  double maxEta = 1.5;
  double minPhi = -TMath::Pi();
  double maxPhi = TMath::Pi();
  double radiusRC = 0.4;
  
  TChain *fChain = new TChain("hiEvtAnalyzer/HiTree");
  TChain *pfTree = new TChain("pfcandAnalyzer/pfTree");
  TChain *skimTree = new TChain("skimanalysis/HltTree");
  TChain *hltTree = new TChain("hltanalysis/HltTree");
  //  TFile *f = TFile::Open(str.Data());
  fChain->Add(str.Data());
  pfTree->Add(str.Data());
  skimTree->Add(str.Data());
  hltTree->Add(str.Data());
  fChain->AddFriend(pfTree);
  fChain->AddFriend(skimTree);
  fChain->AddFriend(hltTree);
  
  if(!fChain) {
    Printf("Couldn't find pfTree. Aborting!");
    return;
  }

  Int_t MinBiasTriggerBit;
  Int_t phfCoincFilter;
  
  Int_t           hiBin;
  TBranch        *b_hiBin;   //!
  fChain->SetBranchAddress("hiBin", &hiBin, &b_hiBin);
  fChain->SetBranchAddress("HLT_HIL1MinimumBiasHF1ANDExpress_v1",&MinBiasTriggerBit);
  fChain->SetBranchAddress("phfCoincFilter3",&phfCoincFilter);
  
  ForestPFs                    fPFs;              //!PFs in tree
  if (fChain->GetBranch("nPFpart"))
    fChain->SetBranchAddress("nPFpart", &fPFs.nPFpart, &fPFs.b_nPFpart);
  if (fChain->GetBranch("pfId"))
    fChain->SetBranchAddress("pfId", fPFs.pfId, &fPFs.b_pfId);
  if (fChain->GetBranch("pfPt"))
    fChain->SetBranchAddress("pfPt", fPFs.pfPt, &fPFs.b_pfPt);
  if (fChain->GetBranch("pfVsPtInitial"))
    fChain->SetBranchAddress("pfVsPtInitial", fPFs.pfVsPt, &fPFs.b_pfVsPt);
  if (fChain->GetBranch("pfEta"))
    fChain->SetBranchAddress("pfEta", fPFs.pfEta, &fPFs.b_pfEta);
  if (fChain->GetBranch("pfPhi"))
    fChain->SetBranchAddress("pfPhi", fPFs.pfPhi, &fPFs.b_pfPhi);

  //Printf("nentries: %d",(Int_t)fChain->GetEntries());

  TList *fOutput =  new TList();
  TH1::SetDefaultSumw2();
  
  TH3D *h2CentPtRCEta = new TH3D("h2CentPtRCEta","h2CentPtRCEta;centrality;p_{T,RC};#eta",100,0,100,250,-50.,200.,60,-6,6);
  fOutput->Add(h2CentPtRCEta);
  TH3D *h2CentPtRCEtaVS = new TH3D("h2CentPtRCEtaVS","h2CentPtRCEtaVS;centrality;p_{T,RC}^{VS};#eta",100,0,100,250,-50.,200.,60,-6,6);
  fOutput->Add(h2CentPtRCEtaVS);

  TH3D *h2MultPtRCEta = new TH3D("h2MultPtRCEta","h2MultPtRCEta;centrality;p_{T,RC};#eta",3000,0,6000,150,-50.,100.,60,-6,6);
  fOutput->Add(h2MultPtRCEta);
  TH3D *h2MultPtRCEtaVS = new TH3D("h2MultPtRCEtaVS","h2MultPtRCEtaVS;centrality;p_{T,RC}^{VS};#eta",3000,0,6000,150,-50.,100.,60,-6,6);
  fOutput->Add(h2MultPtRCEtaVS);
  Printf("histos defined");

  Int_t startEntry = 0;
  Int_t lastEntry = fChain->GetEntries();//100;
  Printf("events in chain: %d",lastEntry);
  if(maxEvents<lastEntry)
    lastEntry = maxEvents;
  Printf("lastEntry: %d",lastEntry);

  TRandom3 *rnd = new TRandom3();
 
  for (int j=startEntry; j<lastEntry; j++) {
    fChain->GetEntry(j);
    if(j%100==0) std::cout << "entry: "<< j << std::endl;
    //if(!MinBiasTriggerBit) continue;
    if(!phfCoincFilter) continue;
    
    double etaRC = rnd->Rndm() * (maxEta - minEta) + minEta;
    double phiRC = rnd->Rndm() * (maxPhi - minPhi) + minPhi;

    double ptRC = 0.;
    double ptRCVS = 0.;
    Int_t pfCount = 0;
    for(Int_t i = 0; i<fPFs.nPFpart; i++) {
      double pt = fPFs.pfPt[i];
      double ptVS = fPFs.pfVsPt[i];
      double phi = fPFs.pfPhi[i];
      double eta = fPFs.pfEta[i];

      double dr = deltaR(phi,phiRC,eta,etaRC);
      if(dr<radiusRC) {
        ptRC+=pt;
        ptRCVS+=ptVS;
      }
//.........这里部分代码省略.........

void analyzeZgHistos(std::vector<std::string> urls, const char *outname = "eventObjects.root", Long64_t nentries = 20, Int_t firstF = -1, Int_t lastF = -1, Int_t firstEvent = 0, int activateJetShift = 0) {

  /*
    ptminType: minimum raw pt for particles used by puppi
    0 : 0 GeV
    1 : 1 GeV
    2 : 2 GeV

    jetSignalType: jets used to select jetty region in events
    0 : detector-level jets (akPu4PF)
    1 : particle-level jets (gen jets)
   */

  TString jetSDName = "aktCs4PFSoftDrop";
  TString jetTreeSDName = "akCsSoftDrop4PFJetAnalyzer";//"akCs4PFSoftDropJetAnalyzer";//
  TString jetName = "aktPu4Calo";
  TString jetTreeName = "akPu4CaloJetAnalyzer";

  bool doDRVar = false;

  bool doJetShift = false;
  if(activateJetShift>0) doJetShift = true;
  float jetShift = 0.04;
  if(activateJetShift==1) jetShift = 0.04;
  if(activateJetShift==2) jetShift = -0.04;

  double minptjet = 30.;
  //bool doDijet = false;
  //if(!doDijet) minptjet = 80.;

  std::cout << "analyzing subjets for: " << jetName << " tree: " << jetTreeName << std::endl;
  std::cout << "analyzing subjets for: " << jetSDName << " tree: " << jetTreeSDName << std::endl;
   
  std::cout << "nfiles: " << urls.size() << std::endl;
  for (auto i = urls.begin(); i != urls.end(); ++i)
    std::cout << *i << std::endl;

  size_t firstFile = 0;
  size_t lastFile = urls.size();

  if(firstF>-1) {
    firstFile = (size_t)firstF;
    lastFile = std::min((size_t)lastF,lastFile);
  }
  std::cout << "firstFile: " << firstFile << "  lastFile: " << lastFile << std::endl;

  Int_t lastEvent = nentries;
  if(firstEvent>0) {
    lastEvent = firstEvent + nentries;
  }
  std::cout << "firstEvent: " << firstEvent << std::endl;
  
  //add files to chain
  TChain *chain = NULL;
  chain = new TChain("hiEvtAnalyzer/HiTree");
  for(size_t i=firstFile; i<lastFile; i++) chain->Add(urls[i].c_str());
  Printf("hiTree done");

  TChain *skimTree = new TChain("skimanalysis/HltTree");
  for(size_t i=firstFile; i<lastFile; i++) skimTree->Add(urls[i].c_str());
  chain->AddFriend(skimTree);
  Printf("skimTree done");

  TChain *hltTree = new TChain("hltanalysis/HltTree");
  for(size_t i=firstFile; i<lastFile; i++) hltTree->Add(urls[i].c_str());
  chain->AddFriend(hltTree);
  Printf("hltTree done");

  TChain *jetTreeSD = new TChain(Form("%s/t",jetTreeSDName.Data()));
  for(size_t i=firstFile; i<lastFile; i++) jetTreeSD->Add(urls[i].c_str());
  chain->AddFriend(jetTreeSD);
  Printf("jetTreeSD done");

  TChain *jetTree = new TChain(Form("%s/t",jetTreeName.Data()));
  for(size_t i=firstFile; i<lastFile; i++) jetTree->Add(urls[i].c_str());
  Printf("jetTree done");

  TChain *rhoTree = new TChain("hiFJRhoAnalyzer/t");
  for(size_t i=firstFile; i<lastFile; i++) rhoTree->Add(urls[i].c_str());
  chain->AddFriend(rhoTree);
  
  TList *fEventObjects = new TList();

  //---------------------------------------------------------------
  // producers
  //
  hiEventProducer *p_evt = new hiEventProducer("hiEvtProd");
  p_evt->SetInput(chain);
  p_evt->SetHIEventContName("hiEventContainer");
  p_evt->SetEventObjects(fEventObjects);

  triggerProducer *p_trg = new triggerProducer("trigProd");
  p_trg->SetInput(chain);
  p_trg->SetTriggerMapName("triggerMap");
  p_trg->AddTrigger("HLT_HIPuAK4CaloJet100_Eta5p1_v1");
  p_trg->AddTrigger("HLT_HIPuAK4CaloJet100_Eta5p1_Cent30_100_v1");
  p_trg->SetEventObjects(fEventObjects);

  lwJetFromForestProducer *p_SDJet = new lwJetFromForestProducer("lwJetForestProdSD");
  p_SDJet->SetInput(chain);
//.........这里部分代码省略.........

std::pair<TString,TString> TMVAClassification (
    TString infilename,
    AnalysisType analysisType = AnalysisType::DIRECT,
    TString additionalRootFileName = "")
{
    TMVA::Tools::Instance();

    std::string tmstr (now ());
    TString tmstmp (tmstr.c_str ());
   
  
    std::cout << "==> Start TMVAClassification" << std::endl;
    std::cout << "-------------------- open input file ---------------- " << std::endl;
    TString fname = infilename; //pathToData + infilename + TString (".root");
    if (analysisType != AnalysisType::TRANSFORMED)
        fname = pathToData + infilename + TString (".root");
    std::cout << "open file " << std::endl << fname.Data () << std::endl;


    std::cout << "-------------------- get tree ---------------- " << std::endl;
    TString treeName = "data";
    if (analysisType == AnalysisType::TRANSFORMED)
        treeName = "transformed";

    std::cout << "-------------------- create tchain with treeName ---------------- " << std::endl;
    std::cout << treeName << std::endl;
    TChain* tree = new TChain (treeName);
    std::cout << "add file" << std::endl;
    std::cout << fname << std::endl;
    tree->Add (fname);
    TChain* treeFriend (NULL);
    if (additionalRootFileName.Length () > 0)
    {
        std::cout << "-------------------- add additional input file ---------------- " << std::endl;
        std::cout << additionalRootFileName << std::endl;
        treeFriend = new TChain (treeName);
        treeFriend->Add (additionalRootFileName);
        tree->AddFriend (treeFriend,"p");
    }
//    tree->Draw ("mass:prediction");
//    return std::make_pair(TString("hallo"),TString ("nix"));
    TString outfileName;
    if (analysisType == AnalysisType::BACKGROUND)
    {
        outfileName = TString ("BACK_" + infilename) + tmstmp + TString (".root");
    }
    else
        outfileName += TString ( "TMVA__" ) + tmstmp + TString (".root");

    std::cout << "-------------------- open output file ---------------- " << std::endl;
    TFile* outputFile = TFile::Open( outfileName, "RECREATE" );

    std::cout << "-------------------- prepare factory ---------------- " << std::endl;
    TMVA::Factory *factory = new TMVA::Factory( "TMVAClassification", outputFile,
						"AnalysisType=Classification:Transformations=I:!V" );
    std::cout << "-------------------- add variables ---------------- " << std::endl;


    for (auto varname : variableNames)
    {
	factory->AddVariable (varname.c_str (), 'F');
    }

    for (auto varname : spectatorNames)
    {
	factory->AddSpectator (varname.c_str (), 'F');
    }
    
   
    std::cout << "-------------------- add trees ---------------- " << std::endl;
    TCut signalCut ("signal==1");
    TCut backgroundCut ("signal==0");
    if (analysisType == AnalysisType::TRANSFORMED)
    {
        signalCut = "(signal_original==1 && signal_in==0)";
        backgroundCut = "(signal_original==0 && signal_in==0)";
    }
    if (analysisType == AnalysisType::BACKGROUND)
    {
        signalCut     = TString("(signal==0) * (prediction > 0.7)");
        backgroundCut = TString("(signal==0) * (prediction < 0.4)");
    }
    //tree->Draw ("prediction",signalCut);
    //return std::make_pair(TString("hallo"),TString ("nix"));
    factory->AddTree(tree, "Signal", 1.0, baseCut + signalCut, "TrainingTesting");
    factory->AddTree(tree, "Background", 1.0, baseCut + backgroundCut, "TrainingTesting");


    
    TCut mycuts = ""; // for example: TCut mycuts = "abs(var1)<0.5 && abs(var2-0.5)<1";
    TCut mycutb = ""; // for example: TCut mycutb = "abs(var1)<0.5";

    /* // Set individual event weights (the variables must exist in the original TTree) */
    if (analysisType == AnalysisType::BACKGROUND)
    {
        factory->SetSignalWeightExpression ("prediction");
        factory->SetBackgroundWeightExpression ("1");
    }

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