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

void RunUEAnalysis(int nentries = 1E3,
		   const char* jetfile  = "/gpfs01/star/i_bnl/gdwebb/Run9/500GeVJets/Data/10103041/st_physics_10103041_raw_*.jets.root",
		   const char* skimfile = "/gpfs01/star/i_bnl/gdwebb/Run9/500GeVJets/Data/10103041/st_physics_10103041_raw_*.skim.root",
		   const char* uefile = "/gpfs01/star/i_bnl/gdwebb/Run9/500GeVJets/Data/10103041/st_physics_10103041_raw_*.ue.root",
		   const char* outfile = "test_jetQAData500Gev_10103041.root",
		   const Int_t  mEmbed = 0)
{
  cout << "jetfile  = " << jetfile  << endl;
  cout << "skimfile = " << skimfile << endl;
  cout << "uefile   = " << uefile << endl;
  cout << "outfile  = " << outfile << endl;
  const float area  =  4*TMath::Pi()/3;
  const float area2 =  2*TMath::Pi()/3;

  cout << "area: " << area << endl;

  // Open jet & skim files
  TChain* jetChain  = new TChain("jet");
  TChain* skimChain = new TChain("jetSkimTree");
  TChain* ueChain   = new TChain("ue");

  jetChain->Add(jetfile);
  skimChain->Add(skimfile);
  ueChain->Add(uefile);
  nentries = jetChain->GetEntriesFast();
  cout << "nentries = " << nentries << endl;
  
  // Set jet buffer
  StJetEvent* jetEvent = 0;
  jetChain->SetBranchAddress("AntiKtR060NHits12",&jetEvent);

  // Set skim buffer
  StJetSkimEvent* skimEvent = 0;
  skimChain->SetBranchAddress("skimEventBranch",&skimEvent); 

  StUeEvent* ueEvent_transP = 0;
  ueChain->SetBranchAddress("transP_AntiKtR060NHits12",&ueEvent_transP);

  StUeEvent* ueEvent_transM = 0;
  ueChain->SetBranchAddress("transM_AntiKtR060NHits12",&ueEvent_transM);

  StUeEvent* ueEvent_away = 0;
  ueChain->SetBranchAddress("away_AntiKtR060NHits12",&ueEvent_away);

  StUeEvent* ueEvent_toward = 0;
  ueChain->SetBranchAddress("toward_AntiKtR060NHits12",&ueEvent_toward);

  // Open output file for writing
  TFile* ofile = TFile::Open(outfile,"recreate");
  assert(ofile);

  //Book histograms
  TH1F* hleadingjetpTJP2 = new TH1F("hleadingjetpTJP2",";Leading jet p_{T} [GeV]",40,0,80);

  //Event Histograms
  TH2F* htransP_numTracksVsMaxJetpTJP2 = new TH2F("htransP_numTracksVsMaxJetpTJP2","transP; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransP_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransP_sumTrackpTVsMaxJetpTJP2","transP;  jet p_{T} [GeV]; <#Sigmap_{T,ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransP_numTowersVsMaxJetpTJP2 = new TH2F("htransP_numTowersVsMaxJetpTJP2","transP; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransP_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransP_sumTowerEtVsMaxJetpTJP2","transP;  jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransP_sumpTVsMaxJetpTJP2 = new TH2F("htransP_sumpTVsMaxJetpTJP2","transP; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransM_numTracksVsMaxJetpTJP2 = new TH2F("htransM_numTracksVsMaxJetpTJP2","transM; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransM_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransM_sumTrackpTVsMaxJetpTJP2","transM;  jet p_{T} [GeV]; #Sigmatrack p_{T}",40,0,80,320,0.,20.);

  TH2F* htransM_numTowersVsMaxJetpTJP2 = new TH2F("htransM_numTowersVsMaxJetpTJP2","transM; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransM_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransM_sumTowerEtVsMaxJetpTJP2","transM;  jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransM_sumpTVsMaxJetpTJP2 = new TH2F("htransM_sumpTVsMaxJetpTJP2","transM; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransverse_numTracksVsMaxJetpTJP2 = new TH2F("htransverse_numTracksVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransverse_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransverse_sumTrackpTVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; #Sigmatrack p_{T}",40,0,80,320,0.,20.);

  TH2F* htransverse_numTowersVsMaxJetpTJP2 = new TH2F("htransverse_numTowersVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransverse_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransverse_sumTowerEtVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransverse_sumpTVsMaxJetpTJP2 = new TH2F("htransverse_sumpTVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransverse_dPtVsMaxJetpTJP2 = new TH2F("htransverse_dPtVsMaxJetpTJP2","transverse; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)] * (jet area)",40,0,80,320,0.,20.);

  TH2F* htransDiff_numTracksVsMaxJetpTJP2 = new TH2F("htransDiff_numTracksVsMaxJetpTJP2","transDiff; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransDiff_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransDiff_sumTrackpTVsMaxJetpTJP2","transDiff; jet p_{T} [GeV]; #Sigmatrack p_{T}",40,0,80,320,0.,20.);
  TH2F* htransDiff_numTowersVsMaxJetpTJP2 = new TH2F("htransDiff_numTowersVsMaxJetpTJP2","transDiff; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransDiff_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransDiff_sumTowerEtVsMaxJetpTJP2","transDiff; jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransDiff_sumpTVsMaxJetpTJP2 = new TH2F("htransDiff_sumpTVsMaxJetpTJP2","transDiff; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransMax_numTracksVsMaxJetpTJP2 = new TH2F("htransMax_numTracksVsMaxJetpTJP2","transMax; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransMax_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransMax_sumTrackpTVsMaxJetpTJP2","transMax;  jet p_{T} [GeV]; <#Sigmap_{T,ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransMax_numTowersVsMaxJetpTJP2 = new TH2F("htransMax_numTowersVsMaxJetpTJP2","transMax; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransMax_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransMax_sumTowerEtVsMaxJetpTJP2","transMax;  jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransMin_numTracksVsMaxJetpTJP2 = new TH2F("htransMin_numTracksVsMaxJetpTJP2","transMin; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransMin_sumTrackpTVsMaxJetpTJP2 = new TH2F("htransMin_sumTrackpTVsMaxJetpTJP2","transMin;  jet p_{T} [GeV]; <#Sigmap_{T,ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransMin_numTowersVsMaxJetpTJP2 = new TH2F("htransMin_numTowersVsMaxJetpTJP2","transMin; jet p_{T} [GeV]; <N_{0}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransMin_sumTowerEtVsMaxJetpTJP2 = new TH2F("htransMin_sumTowerEtVsMaxJetpTJP2","transMin;  jet p_{T} [GeV]; <#SigmaE_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* htransMax_sumpTVsMaxJetpTJP2 = new TH2F("htransMax_sumpTVsMaxJetpTJP2", "transMax; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* htransMin_sumpTVsMaxJetpTJP2 = new TH2F("htransMin_sumpTVsMaxJetpTJP2", "transMin; jet p_{T} [GeV]; <#Sigmap_{T}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);

  TH2F* haway_numTracksVsMaxJetpTJP2 = new TH2F("haway_numTracksVsMaxJetpTJP2","away; jet p_{T} [GeV]; <N_{ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
  TH2F* haway_sumTrackpTVsMaxJetpTJP2 = new TH2F("haway_sumTrackpTVsMaxJetpTJP2","away;  jet p_{T} [GeV]; <#Sigmap_{T,ch}>/[#Delta#eta#Delta(#Delta#phi)]",40,0,80,320,0.,20.);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  p_pfJet->SetJetContName("akt4PF");
  p_pfJet->SetGenJetContName("");//akt4Gen");
  p_pfJet->SetEventObjects(fEventObjects);
  p_pfJet->SetRadius(0.4);
  
  lwJetFromForestProducer *p_caloJet = new lwJetFromForestProducer("lwJetForestProd");
  p_caloJet->SetInput(caloJetTree);
  p_caloJet->SetJetContName("akt4Calo");
  p_caloJet->SetGenJetContName("");
  p_caloJet->SetEventObjects(fEventObjects);
  p_caloJet->SetRadius(0.4);
  
  //---------------------------------------------------------------
  //analysis modules
  //

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

  anaJetMatching *matchingPFCaloJet = new anaJetMatching("matchingPFCaloJet","matchingPFCaloJet");
  matchingPFCaloJet->ConnectEventObject(fEventObjects);
  matchingPFCaloJet->SetHiEvtName("hiEventContainer");
  matchingPFCaloJet->DoPFJet80(true);
  matchingPFCaloJet->DoExcludePhoton30(true);
  matchingPFCaloJet->SetJetsNameBase("akt4PF");
  matchingPFCaloJet->SetJetsNameTag("akt4Calo");
  matchingPFCaloJet->SetNCentBins(1);
  handler->Add(matchingPFCaloJet);  
  
  anaPFvsCaloJet *anaPFCaloJet = new anaPFvsCaloJet("anaPFvsCaloJet","anaPFvsCaloJet");
  anaPFCaloJet->ConnectEventObject(fEventObjects);
  anaPFCaloJet->SetHiEvtName("hiEventContainer");
  anaPFCaloJet->DoPFJet80(true);
  anaPFCaloJet->DoExcludePhoton30(true);
  anaPFCaloJet->SetPFJetsName("akt4PF");
  anaPFCaloJet->SetCaloJetsName("akt4Calo");
  handler->Add(anaPFCaloJet);

  anaPFvsCaloJet *anaCaloPFJet = new anaPFvsCaloJet("anaCalovsPFJet","anaCalovsPFJet");
  anaCaloPFJet->ConnectEventObject(fEventObjects);
  anaCaloPFJet->SetHiEvtName("hiEventContainer");
  anaCaloPFJet->DoPFJet80(true);
  anaCaloPFJet->DoExcludePhoton30(true);
  anaCaloPFJet->SetPFJetsName("akt4Calo");
  anaCaloPFJet->SetCaloJetsName("akt4PF");
  handler->Add(anaCaloPFJet);

  /*
  //PF-GEN matching
  anaJetMatching *matchingGenPFJet = new anaJetMatching("matchingGenPFJet","matchingGenPFJet");
  matchingGenPFJet->ConnectEventObject(fEventObjects);
  matchingGenPFJet->SetHiEvtName("");
  matchingGenPFJet->SetJetsNameBase("akt4Gen");
  matchingGenPFJet->SetJetsNameTag("akt4PF");
  matchingGenPFJet->SetNCentBins(1);
  handler->Add(matchingGenPFJet);  
  
  anaPFvsCaloJet *anaGenPFJet = new anaPFvsCaloJet("anaGenVsPFJet","anaGenVsPFJet");
  anaGenPFJet->ConnectEventObject(fEventObjects);
  anaGenPFJet->SetHiEvtName("");
  anaGenPFJet->SetPFJetsName("akt4Gen");
  anaGenPFJet->SetCaloJetsName("akt4PF");
  handler->Add(anaGenPFJet);
  */
  
  //---------------------------------------------------------------
  //Event loop
  //---------------------------------------------------------------	  
  Long64_t entries_tot =  chain->GetEntriesFast(); //93064
  if(nentries<0) nentries = chain->GetEntries();
  // Long64_t nentries = 20;//chain->GetEntriesFast();
  Printf("nentries: %lld  tot: %lld",nentries,entries_tot);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    if (jentry%10000==0) Printf("Processing event %d  %d",(int)(jentry), (int)(nentries));
    //Run producers
    // Printf("produce hiEvent");
    p_evt->Run(jentry);   //hi event properties

    //Printf("produce pf particles");
    //  p_pf->Run(jentry);    //pf particles
    p_pfJet->Run(jentry); //jets
    p_caloJet->Run(jentry); //jets
    
    //Execute all analysis tasks
    handler->ExecuteTask();
  }
    
  fEventObjects->Print();

  TFile *out = new TFile(outname,"RECREATE");
  TList *tasks = handler->GetListOfTasks();
  TIter next(tasks);
  anaBaseTask *obj;
  while ((obj = dynamic_cast<anaBaseTask*>(next()) ))
    if(obj->GetOutput()) obj->GetOutput()->Write(obj->GetName(),TObject::kSingleKey);
  
  out->Write();
  out->Close();
  
}

//.........这里部分代码省略.........
        chain->Add(fname);
    } else if (process == "ZnnH145") {
        fname = dirMC + dijet + "ZH_ZToNuNu_HToBB_M-145_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZnnH150") {
        fname = dirMC + dijet + "ZH_ZToNuNu_HToBB_M-150_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);

    // WlnH
    } else if (process == "WlnH110") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-110_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH115") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-115_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH120") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-120_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH125") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-125_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH130") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-130_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH135") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-135_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH140") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-140_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH145") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-145_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "WlnH150") {
        fname = dirMC + dijet + "WH_WToLNu_HToBB_M-150_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);

    // ZllH
    } else if (process == "ZllH110") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-110_8TeV-powheg-herwigpp3" + suffix;  // NOTE: '3' in name
        chain->Add(fname);
    } else if (process == "ZllH115") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-115_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH120") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-120_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH125") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-125_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH130") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-130_8TeV-powheg-herwigpp3" + suffix;  // NOTE: '3' in name
        chain->Add(fname);
    } else if (process == "ZllH135") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-135_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH140") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-140_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH145") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-145_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);
    } else if (process == "ZllH150") {
        fname = dirMC + dijet + "ZH_ZToLL_HToBB_M-150_8TeV-powheg-herwigpp" + suffix;
        chain->Add(fname);

    } else {
        std::cout << "Error: Unrecognized process. I quit!" << std::endl;
        return;
    }

    TCut selection = regression.c_str();
    
    // Make output directory if it doesn't exist
    if (gSystem->AccessPathName(outdir))
        gSystem->mkdir(outdir);
    TString outname = outdir + prefix + Form("%s.root", process.Data());

    TFile* f1 = TFile::Open(outname, "RECREATE");
    TTree* t1 = (TTree*) chain->CopyTree(selection);

    if (t1->GetEntriesFast() > skimentries){
        TTree* t2 = t1->CopyTree("", "", skimentries);
        t2->SetName(TString(t1->GetName())+"_train");
        TTree* t3 = t1->CopyTree("", "", 1000000000, skimentries);
        t3->SetName(TString(t1->GetName())+"_test");
        t2->Write();
        t3->Write();
        assert(t1->GetEntriesFast() == t2->GetEntriesFast() + t3->GetEntriesFast());
        delete t1;
        std::clog << process << ": skimmed from " << chain->GetEntriesFast() << " to " << t2->GetEntriesFast() << " (training) and " << t3->GetEntriesFast() << " (test) entries." << std::endl;
    } else {
        t1->Write();
        std::clog << process << ": skimmed from " << chain->GetEntriesFast() << " to " << t1->GetEntriesFast() << " entries." << std::endl;
    }
    
    f1->Close();

    return;
}

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

  //anti-kt jet finder on reconstructed PUPPI particles ptmin=1
  LWJetProducer *lwjaktpt100 = new LWJetProducer("LWJetProducerAKTR030PuppiPtMin100","LWJetProducerAKTR030PuppiPtMin100");
  lwjaktpt100->ConnectEventObject(fEventObjects);
  lwjaktpt100->SetJetType(LWJetProducer::kAKT);
  lwjaktpt100->SetRadius(0.3);
  lwjaktpt100->SetGhostArea(0.005);
  lwjaktpt100->SetPtMinConst(1.);
  lwjaktpt100->SetParticlesName("puppiParticles");
  lwjaktpt100->SetJetContName("JetsAKTR030PuppiPtMin100");
  if(doPuppi) handler->Add(lwjaktpt100);

  //anti-kt jet finder on reconstructed PUPPI particles ptmin=2
  LWJetProducer *lwjaktpt200 = new LWJetProducer("LWJetProducerAKTR030PuppiPtMin200","LWJetProducerAKTR030PuppiPtMin200");
  lwjaktpt200->ConnectEventObject(fEventObjects);
  lwjaktpt200->SetJetType(LWJetProducer::kAKT);
  lwjaktpt200->SetRadius(0.3);
  lwjaktpt200->SetGhostArea(0.005);
  lwjaktpt200->SetPtMinConst(2.);
  lwjaktpt200->SetParticlesName("puppiParticles");
  lwjaktpt200->SetJetContName("JetsAKTR030PuppiPtMin200");
  if(doPuppi) handler->Add(lwjaktpt200);

  //anti-kt jet finder on generated particles
  LWJetProducer *lwjaktGen = new LWJetProducer("LWGenJetProducerAKTR030","LWGenJetProducerAKTR030");
  lwjaktGen->ConnectEventObject(fEventObjects);
  lwjaktGen->SetJetType(LWJetProducer::kAKT);
  lwjaktGen->SetRadius(0.3);
  lwjaktGen->SetGhostArea(0.005);
  lwjaktGen->SetPtMinConst(0.);
  lwjaktGen->SetParticlesName("genParticles");
  lwjaktGen->SetJetContName("GenJetsAKTR030");
  handler->Add(lwjaktGen);
 
  //Initialization of all analysis modules
  anaPuppiParticles *pupAna = new anaPuppiParticles("pupAna","pupAna");
  pupAna->ConnectEventObject(fEventObjects);
  pupAna->SetHiEvtName("hiEventContainer");
  pupAna->SetParticlesName("pfParticles");
  pupAna->SetJetsName("aktPUR030");
  if(doPuppi) handler->Add(pupAna);

  //particle-detector-level jet matching
  anaJetMatching *match = new anaJetMatching("jetMatching","jetMatching");
  match->ConnectEventObject(fEventObjects);
  match->SetHiEvtName("hiEventContainer");
  match->SetJetsNameBase("JetsAKTR030Puppi");
  match->SetJetsNameTag("GenJetsAKTR030");
  if(doPuppi) handler->Add(match);

  anaJetMatching *matchPt100 = new anaJetMatching("jetMatchingPtMin100","jetMatchingPtMin100");
  matchPt100->ConnectEventObject(fEventObjects);
  matchPt100->SetHiEvtName("hiEventContainer");
  matchPt100->SetJetsNameBase("JetsAKTR030PuppiPtMin100");
  matchPt100->SetJetsNameTag("GenJetsAKTR030");
  if(doPuppi) handler->Add(matchPt100);

  anaJetMatching *matchPt200 = new anaJetMatching("jetMatchingPtMin200","jetMatchingPtMin200");
  matchPt200->ConnectEventObject(fEventObjects);
  matchPt200->SetHiEvtName("hiEventContainer");
  matchPt200->SetJetsNameBase("JetsAKTR030PuppiPtMin200");
  matchPt200->SetJetsNameTag("GenJetsAKTR030");
  if(doPuppi) handler->Add(matchPt200);
  
  //---------------------------------------------------------------
  //Event loop
  //---------------------------------------------------------------	  
  Long64_t entries_tot =  chain->GetEntriesFast(); //93064
  if(nentries<0) nentries = chain->GetEntries();
  // Long64_t nentries = 20;//chain->GetEntriesFast();
  Printf("nentries: %lld  tot: %lld",nentries,entries_tot);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {

    //Run producers
    Printf("produce hiEvent");
    p_evt->Run(jentry);   //hi event properties
    Printf("produce pf particles");
    p_pf->Run(jentry);    //pf particles
    Printf("produce gen particles");
    p_gen->Run(jentry);   //generated particles
    Printf("produce PU jets");
    p_PUJet->Run(jentry); //forest jets
	    
    //Execute all analysis tasks
    // handler->ExecuteTask();
  }
    
  fEventObjects->Print();

  TFile *out = new TFile(outname,"RECREATE");
  TList *tasks = handler->GetListOfTasks();
  TIter next(tasks);
  anaBaseTask *obj;
  while ((obj = dynamic_cast<anaBaseTask*>(next()) ))
    if(obj->GetOutput()) obj->GetOutput()->Write(obj->GetName(),TObject::kSingleKey);
  
  out->Write();
  out->Close();
  
}

Double_t CalibTree::Loop(int loop, TFile *fout, bool useweight, int nMin,
			 bool inverse, double rmin, double rmax, int ietaMax,
			 int applyL1Cut, double l1Cut, bool last, 
			 double fraction, bool writeHisto, bool debug) {

  if (fChain == 0) return 0;
  Long64_t nbytes(0), nb(0), kprint(0);
  Long64_t nentryTot = fChain->GetEntriesFast();
  Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? 
    (Long64_t)(fraction*nentryTot) : nentryTot;
  if (detIds.size() == 0) {
    for (Long64_t jentry=0; jentry<nentries; jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // Find DetIds contributing to the track
      bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) :
		     ((t_Run < runlo_) || (t_Run > runhi_)));
      if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_)) {
	bool isItRBX = (cSelect_ && exclude_ && cSelect_->isItRBX(t_DetIds));
	++kprint;
	if (!isItRBX) {
	  for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { 
	    if (selectPhi((*t_DetIds)[idet])) {
	      unsigned int detid = truncateId((*t_DetIds)[idet],
					      truncateFlag_,debug);
	      if (debug && (kprint<=10)) {
		std::cout << "DetId[" << idet << "] Original " << std::hex 
			  << (*t_DetIds)[idet] << " truncated " << detid 
			  << std::dec;
	      }
	      if (std::find(detIds.begin(),detIds.end(),detid) == detIds.end()){
		detIds.push_back(detid);
		if (debug && (kprint<=10)) std::cout << " new";
	      }
	      if (debug && (kprint<=10)) std::cout << std::endl;
	    }
	  }
	  // Also look at the neighbouring cells if available
	  if (t_DetIds3 != 0) {
	    for (unsigned int idet=0; idet<(*t_DetIds3).size(); idet++) { 
	      if (selectPhi((*t_DetIds3)[idet])) {
		unsigned int detid = truncateId((*t_DetIds3)[idet],
						truncateFlag_,debug);
		if (std::find(detIds.begin(),detIds.end(),detid)==detIds.end()){
		  detIds.push_back(detid);
		}
	      }
	    }
	  }
	}
      }
    }
  }
  if (debug) {
    std::cout << "Total of " << detIds.size() << " detIds and " 
	      << histos.size() << " histos found" << std::endl;
    // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
    for (unsigned int k=0; k<detIds.size(); ++k) {
      int subdet, depth, zside, ieta, iphi;
      unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
      std::cout << "DetId[" << k << "] " << subdet << ":" << zside*ieta << ":"
		<< depth << ":" << iphi << "  " << std::hex << detIds[k] 
		<< std::dec << std::endl;
    }
  }
  unsigned int k(0);
  for (std::map<unsigned int, TH1D*>::const_iterator itr = histos.begin();
       itr != histos.end(); ++itr,++k) {
    if (debug) {
      std::cout << "histos[" << k << "] " << std::hex << itr->first 
		<< std::dec << " " << itr->second;
      if (itr->second != 0) std::cout << " " << itr->second->GetTitle();
      std::cout << std::endl;
    }
    if (itr->second != 0) itr->second->Delete();
  }

  for (unsigned int k=0; k<detIds.size(); ++k) {
    char name[20], title[100];
    sprintf (name, "Hist%d_%d", detIds[k], loop);
    int subdet, depth, zside, ieta, iphi;
    unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
    sprintf (title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside*ieta, depth, loop);
    TH1D* hist = new TH1D(name,title,100, 0.0, 5.0);
    hist->Sumw2();
    if (debug) {
      std::cout << "Book Histo " << k << " " << title << std::endl;
    }
    histos[detIds[k]] = hist;
  }
  std::cout << "Total of " << detIds.size() << " detIds and " << histos.size() 
	    << " found in " << nentries << std::endl;

  nbytes = nb = 0;
  std::map<unsigned int, myEntry > SumW;
  std::map<unsigned int, double  > nTrks;

  int ntkgood(0);
  for (Long64_t jentry=0; jentry<nentries; jentry++) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   nVtx_pfMetT = 0;
   nVtx_pfMetX = 0;
   nVtx_pfMetY = 0;
   sumPt_catagory = 0;
   sumPt_sumPt = 0;
   sumPt_MetX = 0;
   sumPt_MetY = 0;
   sumPt_pfMetT = 0;
   sumPt_pfMetX = 0;
   sumPt_pfMetY = 0;

   fChain->SetBranchAddress("Count_catagory", &Count_catagory, &b_Count_catagory);
   fChain->SetBranchAddress("Count_counts", &Count_counts, &b_Count_counts);
   fChain->SetBranchAddress("Count_MetX", &Count_MetX, &b_Count_MetX);
   fChain->SetBranchAddress("Count_MetY", &Count_MetY, &b_Count_MetY);
   fChain->SetBranchAddress("Count_pfMetT", &Count_pfMetT, &b_Count_pfMetT);
   fChain->SetBranchAddress("Count_pfMetX", &Count_pfMetX, &b_Count_pfMetX);
   fChain->SetBranchAddress("Count_pfMetY", &Count_pfMetY, &b_Count_pfMetY);
   fChain->SetBranchAddress("nVtx_catagory", &nVtx_catagory, &b_nVtx_catagory);
   fChain->SetBranchAddress("nVtx_nVtx", &nVtx_nVtx, &b_nVtx_nVtx);
   fChain->SetBranchAddress("nVtx_MetX", &nVtx_MetX, &b_nVtx_MetX);
   fChain->SetBranchAddress("nVtx_MetY", &nVtx_MetY, &b_nVtx_MetY);
   fChain->SetBranchAddress("nVtx_pfMetT", &nVtx_pfMetT, &b_nVtx_pfMetT);
   fChain->SetBranchAddress("nVtx_pfMetX", &nVtx_pfMetX, &b_nVtx_pfMetX);
   fChain->SetBranchAddress("nVtx_pfMetY", &nVtx_pfMetY, &b_nVtx_pfMetY);
   fChain->SetBranchAddress("sumPt_catagory", &sumPt_catagory, &b_sumPt_catagory);
   fChain->SetBranchAddress("sumPt_sumPt", &sumPt_sumPt, &b_sumPt_sumPt);
   fChain->SetBranchAddress("sumPt_MetX", &sumPt_MetX, &b_sumPt_MetX);
   fChain->SetBranchAddress("sumPt_MetY", &sumPt_MetY, &b_sumPt_MetY);
   fChain->SetBranchAddress("sumPt_pfMetT", &sumPt_pfMetT, &b_sumPt_pfMetT);
   fChain->SetBranchAddress("sumPt_pfMetX", &sumPt_pfMetX, &b_sumPt_pfMetX);
   fChain->SetBranchAddress("sumPt_pfMetY", &sumPt_pfMetY, &b_sumPt_pfMetY);

   Long64_t nentries = fChain->GetEntriesFast();
   //nentries = 30;

   Long64_t nbytes = 0, nb = 0;

   TCanvas *myCan = new TCanvas("myCan","myCan");
   TProfile *pr_MEx_vtx         = new TProfile("pr_MEx_vtx","pr_MEx_vtx",20,0,40);
   TProfile *pr_MEy_vtx         = new TProfile("pr_MEy_vtx","pr_MEy_vtx",20,0,40);
   TProfile *pr_MEx_pfMetT_vtx[vtxBin];
   TProfile *pr_MEx_vtx_pfMetT[metBin];
   TProfile *pr_MEy_vtx_pfMetT[metBin];
   TProfile *pr_MEx_pfMetT = new TProfile("pr_MEx_pfMetT","pr_MEx_pfMetT",30,0,90);
   TString name;
   for(int i(0);i<vtxBin;i++){
     name = Form("pr_MEx_pfMetT_vtx%d", (i+1)*metBin);
     pr_MEx_pfMetT_vtx[i]  = new TProfile(name,name,30,0,90);
   }
   for(int i(0); i<metBin;i++){
     name = Form("pr_MEx_vtx_pfMetT%d", (i+1)*metBin);
     pr_MEx_vtx_pfMetT[i]  = new TProfile(name,name,20,0,40);
     pr_MEx_vtx_pfMetT[i]->SetMaximum(5);
     pr_MEx_vtx_pfMetT[i]->SetMinimum(-3);

     name = Form("pr_MEy_vtx_pfMetT%d", (i+1)*metBin);
     pr_MEy_vtx_pfMetT[i]  = new TProfile(name,name,20,0,40);
     pr_MEy_vtx_pfMetT[i]->SetMaximum(5);
     pr_MEy_vtx_pfMetT[i]->SetMinimum(-3);
   }
   TH1D *h_pfMetT                = new TH1D("h_pfMetT","h_pfMetT",50,0,100);

   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = fChain->LoadTree(jentry);
      if (ientry < 0) break;

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

    TCut selection = baseline.c_str();
    // Different baseline for dimuons
    // MET filters
    selection += metfilter.c_str();

    // JSON & trigger
    if (fname.Contains("Data")) {
        TCut trigger = mettrigger.c_str();
        selection += trigger;
        //selection.Print();
      }
      /*
    } else if (process == "Data_METBTag_R" || process == "Data_METBTag_P") {
        TCut trigger = metbtagtrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleMu_R" || process == "Data_SingleMu_P") {
        TCut trigger = singlemutrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleEl_R" || process == "Data_SingleEl_P") {
        TCut trigger = singleeltrigger.c_str();
        selection += trigger;
        //selection.Print();
    }
      */


    chain->Add(fname);

    // Sum Count, CountWithPU, CountWithPUP, CountWithPUM
    TObjArray * files = chain->GetListOfFiles();
    TIter next(files);
    TChainElement * chainElem = 0;
    TFile * f2 = 0;
    TH1D * h1 = new TH1D("Count", ";Counts", 16, 0, 16);
    TH1F * htemp = 0;

    while ((chainElem = (TChainElement*) next())) {
      //#ifndef XROOTD
      //      f2 = TFile::Open("dcache:" + TString(chainElem->GetTitle()));
      //#else
      std::cout << "chainElem->GetTitle() " << chainElem->GetTitle() << std::endl; 
        f2 = TFile::Open( TString(chainElem->GetTitle()));
	//#endif
        htemp = (TH1F*) f2->Get("Count");
        h1->SetBinContent(1, h1->GetBinContent(1)+htemp->GetBinContent(1));
	/*
        htemp = (TH1F*) f2->Get("CountWithPU");
        h1->SetBinContent(2, h1->GetBinContent(2)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUP");
        h1->SetBinContent(3, h1->GetBinContent(3)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUM");
        h1->SetBinContent(4, h1->GetBinContent(4)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithMCProd");
        h1->SetBinContent(5, h1->GetBinContent(5)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUMCProd");
        h1->SetBinContent(6, h1->GetBinContent(6)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("countWithSignalQCDcorrections");
        h1->SetBinContent(7, h1->GetBinContent(7)+htemp->GetBinContent(1));
        */
        std::clog << fname << ": skimmed from " << chainElem->GetTitle() << std::endl;
    }
    
    // LHE Count
    
    TH1D * h2 = new TH1D("LHECount", ";LHE Counts", 16, 0, 16);
    TString process_lhe = fname;
    if (process_lhe.BeginsWith("WJets"))
        process_lhe = "WJets";
    else if (process_lhe.BeginsWith("ZJets"))
        process_lhe = "ZJets";
    else 
        process_lhe = "";
    const std::vector<std::string>& lhecuts = GetLHECuts(process_lhe.Data());
    for (unsigned int i=0; i < lhecuts.size(); i++) {
        TCut cut2 = lhecuts.at(i).c_str();
        h2->SetBinContent(i+1, chain->GetEntries(cut2));
    }
    
    

    // Make output directory if it doesn't exist
    if (gSystem->AccessPathName(outdir))
        gSystem->mkdir(outdir);
    TString outname = outdir + prefix + Form("%s.root", outputname.Data());
    std::cout << "outname is " << outname << std::endl;

    TFile* f1 = TFile::Open(outname, "RECREATE");
    TTree* t1 = (TTree*) chain->CopyTree(selection);
    std::clog << fname << ": skimmed from " << chain->GetEntriesFast() << " to " << t1->GetEntriesFast() << " entries." << std::endl;
    
    t1->Write();
    h1->Write();
     h2->Write();
    f1->Close();

    return;
}

void CalibTree::makeplots(double rmin, double rmax, int ietaMax,
			  bool useweight, double fraction, bool debug) {

  if (fChain == 0) return;
  Long64_t nentryTot = fChain->GetEntriesFast();
  Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? 
    (Long64_t)(fraction*nentryTot) : nentryTot;

  // Book the histograms
  std::map<int,std::pair<TH1D*,TH1D*> > histos;
  for (int ieta=-ietaMax; ieta<=ietaMax; ++ieta) {
    char name[20], title[100];
    sprintf(name,"begin%d",ieta);
    if (ieta==0) sprintf(title,"Ratio at start");
    else         sprintf(title,"Ratio at start for i#eta=%d",ieta);
    TH1D* h1 = new TH1D(name,title,50,rmin,rmax);
    h1->Sumw2();
    sprintf(name,"end%d",ieta);
    if (ieta==0) sprintf(title,"Ratio at the end");
    else         sprintf(title,"Ratio at the end for i#eta=%d",ieta);
    TH1D* h2 = new TH1D(name,title,50,rmin,rmax);
    h2->Sumw2();
    histos[ieta] = std::pair<TH1D*,TH1D*>(h1,h2);
  }
  //Fill the histograms
  Long64_t nbytes(0), nb(0);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    if (std::find(entries.begin(), entries.end(), jentry) != entries.end()) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if (goodTrack()) {
      double Etot(0);
      for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { 
	double hitEn(0);
	unsigned int id    = (*t_DetIds)[idet];
        unsigned int detid = truncateId(id,truncateFlag_,false);
	if (Cprev.find(detid) != Cprev.end()) 
	  hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
	else 
	  hitEn = (*t_HitEnergies)[idet];
	if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id);
	Etot += hitEn;
      }
      double evWt   = (useweight) ? t_EventWeight : 1.0; 
      double pmom   = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
      double ratioi = t_eHcal/(pmom-t_eMipDR);
      double ratiof = Etot/(pmom-t_eMipDR);
      if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
	if (ratioi>=rmin && ratioi<=rmax) {
	  histos[0].first->Fill(ratioi,evWt);
	  histos[t_ieta].first->Fill(ratioi,evWt);
	}
	if (ratiof>=rmin && ratiof<=rmax) {
	  histos[0].second->Fill(ratiof,evWt);
	  histos[t_ieta].second->Fill(ratiof,evWt);
	}
      }
    }
  }

  //Fit the histograms
  TH1D *hbef1 = new TH1D("Eta1Bf","Mean vs i#eta",2*ietaMax,-ietaMax,ietaMax);
  TH1D *hbef2 = new TH1D("Eta2Bf","Median vs i#eta",2*ietaMax,-ietaMax,ietaMax);
  TH1D *haft1 = new TH1D("Eta1Af","Mean vs i#eta",2*ietaMax,-ietaMax,ietaMax);
  TH1D *haft2 = new TH1D("Eta2Af","Median vs i#eta",2*ietaMax,-ietaMax,ietaMax);
  for (int ieta=-ietaMax; ieta<=ietaMax; ++ieta) {
    int    bin   = (ieta < 0) ? (ieta+ietaMax+1) : (ieta+ietaMax);
    TH1D*  h1    = histos[ieta].first;
    double mean1 = h1->GetMean();
    double err1  = h1->GetMeanError();
    std::pair<double,double> fit1 = fitMean(h1,1);
    if (debug) {
      std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " 
		<< err1 << " and " << fit1.first <<" +- " << fit1.second 
		<< std::endl;
    }
    if (ieta != 0) {
      hbef1->SetBinContent(bin,mean1);      hbef1->SetBinError(bin,err1);
      hbef2->SetBinContent(bin,fit1.first); hbef2->SetBinError(bin,fit1.second);
    }
    h1->Write();
    TH1D* h2 = histos[ieta].second;
    double mean2 = h2->GetMean();
    double err2  = h2->GetMeanError();
    std::pair<double,double> fit2 = fitMean(h2,1);
    if (debug) { 
      std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " 
		<< err2 << " and " << fit2.first <<" +- " << fit2.second 
		<< std::endl;
    }
    if (ieta != 0) {
      haft1->SetBinContent(bin,mean2);      haft1->SetBinError(bin,err2);
      haft2->SetBinContent(bin,fit2.first); haft2->SetBinError(bin,fit2.second);
    }
    h2->Write();
  }
  fitPol0(hbef1,debug); fitPol0(hbef2,debug);
  fitPol0(haft1,debug); fitPol0(haft2,debug);
}

//.........这里部分代码省略.........
  tree->SetBranchAddress("genVV_jet1_phi", &jetGen1_phi);
  tree->SetBranchAddress("genVV_jet2_phi", &jetGen2_phi);
  tree->SetBranchAddress("genVV_jet3_phi", &jetGen3_phi);
  tree->SetBranchAddress("genVV_jet4_phi", &jetGen4_phi);  
  tree->SetBranchAddress("genVV_jet5_phi", &jetGen5_phi);

// 

 

 // Set the channel
  //----------------------------------------------------------------------------
  Float_t SelectedChannel = -999;

  /*  if      (flavorChannel == "MuMu") SelectedChannel =  0;
  else if (flavorChannel == "EE"  ) SelectedChannel =  1;
  else if (flavorChannel == "EMu" ) SelectedChannel =  2;
  else if (flavorChannel == "MuE" ) SelectedChannel =  3;
  else if (flavorChannel == "All" ) SelectedChannel = -1;
  */

  int kk = 0;

 //----------------------------------------------------------------------------
  // Loop
  //----------------------------------------------------------------------------
  

  // Float_t Nentries = 179545; 
  //Float_t Nentries = 187525;
  Float_t Nentries = 10745; // with powheg 
  //Float_t Nentries = 9792; // with madgraph 

  for (int ievent=0; ievent<tree->GetEntriesFast(); ievent++) {
  //for (int ievent=0; ievent<Nentries; ievent++) {
  //for (int ievent=Nentries; ievent<tree->GetEntriesFast(); ievent++) {
   
    tree->GetEntry(ievent);

    //Double_t mybaseW =  5984.0/1933235;//5812.3/1933235; // madgraph (1933232)
    //Double_t mybaseW = 5984.0/999864;//  5812.3/999864; // powheg (999860)
    Double_t mybaseW = 182.852 /109986; // GGWW 
    //Double_t mybaseW = 5984.0/539594; //5812.3/539594; // mcnlo



    Float_t luminosity = 19.365;

    Double_t efficiencyW =  effW * triggW * puW;
    // Double_t totalW = effW * triggW * baseW * efficiencyW * luminosity;

    Double_t totalW = puW * mybaseW * luminosity;//efficiencyW;

    Double_t totalWGen = nllW *  mybaseW * luminosity ; //Removing PU from gen level
    
    Double_t totalWReco =  effW * triggW * puW * nllW *  mybaseW * luminosity;// * effW * triggW ;//puW *  mybaseW * luminosity;//puW * effW * triggW * mybaseW * luminosity;

  
    // The GEN selection begins here
    //--------------------------------------------------------------------------
    
    /// ---> 1) Need status 1 leptons to define the same fiducial region
    /// ---> 2) Count how many GEN leptons we have in each bin, applying the fidual region cuts
    /// ---> 3) Apply also, OF, jetbin and opposite-charged cuts.
    
    

void Hyperon(const char* dataset="collection.xml")
{
  /* $Id$ */
  TStopwatch timer;
  timer.Start();


    TStringToken libs("Core,Tree,Geom,VMC,Physics,Minuit,Gui,XMLParser,Minuit2,Proof,STEERBase,ESD,AOD,OADB,ANALYSIS,ANALYSISalice,CDB,RAWDatabase,STEER,CORRFW,PHOSUtils,PHOSbase,PHOSpi0Calib,PHOSrec,PHOSshuttle,PHOSsim", ",");
  while( libs.NextToken() )
    gSystem->Load( Form("lib%s", libs.Data()) );

  gSystem->Load("libTree");
  gSystem->Load("libGeom");
  gSystem->Load("libVMC");
  gSystem->Load("libPhysics");
  gSystem->Load("libPWGGAGammaConv");
  gSystem->Load("libPWGGAHyperon");
  
  //load analysis framework
  gSystem->Load("libANALYSIS");
  gSystem->Load("libANALYSISalice"); //AliAnalysisTaskSE
  
  // Connect to alien
  TString token = gSystem->Getenv("GRID_TOKEN") ;
  if (1) // token == "OK" ) 
    TGrid::Connect("alien://");
  else 
    AliInfo("You are not connected to the GRID") ; 

  cout << "Pi0Analysis: processing collection " << dataset << endl;

  // Create the chain
  TChain* chain = new TChain("aodTree");

  TGridCollection * collection = dynamic_cast<TGridCollection*>(TAlienCollection::Open(dataset));
  
  TAlienResult* result = collection->GetGridResult("",0 ,0);
  TList* rawFileList = result->GetFileInfoList();
  
  for (Int_t counter=0 ; counter < rawFileList->GetEntries() ; counter++) {
    TFileInfo * fi =  static_cast<TFileInfo*>(rawFileList->At(counter)) ; 
    const char * rawFile = fi->GetCurrentUrl()->GetUrl() ;  
    printf("Processing %s\n", rawFile) ;
    chain->Add(rawFile);
    printf("Chain: %d entries.\n",chain->GetEntriesFast()); 
  }

  // Make the analysis manager
  AliAnalysisManager *mgr = new AliAnalysisManager("Hyperon");
  mgr->SetCommonFileName("histos.root");
  
  // AOD input handler
  AliAODInputHandler* aodH = new AliAODInputHandler();
  mgr->SetInputEventHandler(aodH);

  
  // Debug level
  mgr->SetDebugLevel(2);

   gROOT->LoadMacro("$ALICE_PHYSICS/OADB/macros/AddTaskCentrality.C"); 
   AliCentralitySelectionTask* taskCentrality = AddTaskCentrality(); 
   if (analysisMC) 
    taskCentrality->SetMCInput(); 

  // // Update it for Hyperon (YK 22.01.2015)
  // gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskEventplane.C");
  // AliEPSelectionTask *taskEP = AddTaskEventplane() ; 

  // gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C");
  // AliVZEROEPSelectionTask *selTask = AddTaskVZEROEPSelection();  

  // Add my task
  AliAnalysisHyperon* task = new AliAnalysisHyperon();
  // gROOT->LoadMacro("$ALICE_ROOT/PWGGA/PHOSTasks/PHOS_PbPb/AddTaskPHOSPi0Flow.C"); // Update it for Hyperon (YK 22.01.2015)
  // task = AddTaskPHOSPi0Flow();

  
  // // Create containers for input/output
  AliAnalysisDataContainer *cinput   = mgr->GetCommonInputContainer(); 
  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("Hyperon",TList::Class(),AliAnalysisManager::kOutputContainer,"HyperonHist.root");
  
  // // Connect input/output
  mgr->ConnectInput(task , 0, cinput);
  mgr->ConnectOutput(task, 1, coutput1);
 
  if (mgr->InitAnalysis()) {
    mgr->PrintStatus();
    mgr->StartAnalysis("local", chain);
  }

  timer.Stop();
  timer.Print();
}

void readTree(TString filelist="/global/project/projectdirs/star/pwg/starhf/zamiller/minBiasTemplates101015/Trees.list",TString outFile="/global/project/projectdirs/star/pwg/starhf/zamiller/minBiasTemplates101015/readTreeOutSubmitTEST.root")
{

  //TString filelist=Form("%s/%s",localDir,fileIn);
  //TString outFile=Form("%s/%s",outDir,fileOut);

  TFile fout(outFile,"RECREATE");
  TH1F::SetDefaultSumw2(); 
  TH1D::SetDefaultSumw2(); 

  for(Int_t c=0; c< numPtBins; c++)
  {
    lowpt[c] = anaConst::lpt[c];
    highpt[c] = anaConst::hpt[c];
  }
  for(Int_t ptbin = 0; ptbin < numPtBins; ptbin++)
  {
    hdEtaRawce[ptbin] = new TH1D(Form("hdEtaRawce_%i",ptbin),"hdEtaRawce",100, -5,5);
    hdEtaRawbe[ptbin] = new TH1D(Form("hdEtaRawbe_%i",ptbin),"hdEtaRawbe",100, -5,5);
    hEventTallyce[ptbin] = new TH1D(Form("ceEventTally_%i",ptbin),"ceEvent Tally",1,0,1);
    hEventTallybe[ptbin] = new TH1D(Form("beEventTally_%i",ptbin),"beEvent Tally",1,0,1);
    hEventTallybce[ptbin] = new TH1D(Form("bceEventTally_%i",ptbin),"bceEvent Tally",1,0,1);
    hdPhiRawce[ptbin] = new TH1D(Form("hdPhiRawce_%i",ptbin),"hdPhiRawce",Phibin, -10,10);
    hdPhiRawbe[ptbin] = new TH1D(Form("hdPhiRawbe_%i",ptbin),"hdPhiRawbe",Phibin, -10,10);
    hdPhiRawbce[ptbin] = new TH1D(Form("hdPhiRawbce_%i",ptbin),"hdPhiRawbce",Phibin, -10,10);
    hept[ptbin] = new TH1D(Form("hept_%i",ptbin),"hept",200,0.,20.);
  }

  // Make Chain
  TChain* chain = new TChain("tree");
  int nfile = 0;
  /*  nfile += chain->Add("Oct30_set13/pythia_tree_Oct30_13.root");
      nfile += chain->Add("Oct30_set14/pythia_tree_Oct30_14.root");
      nfile += chain->Add("Oct30_set15/pythia_tree_Oct30_15.root");
      nfile += chain->Add("Oct30_set16/pythia_tree_Oct30_16.root");
      nfile += chain->Add("Oct30_set17/pythia_tree_Oct30_17.root");
  //  nfile += chain->Add("Oct30_set12/pythia_tree_Oct30_12.root");
  //  nfile += chain->Add("Oct26_set3/pythia_tree_Oct26_3.root");
  // nfile += chain->Add("Oct18_set4/pythia_tree_Oct18_4.root");
  //nfile += chain->Add("liweiTemplate_part2.root");*/
  char filename[1000];
  ifstream fstream(filelist.Data());
  int ifile = 0;
  while (fstream >> filename)
  {
    nfile+= chain->Add(filename);
    std::cout << "Got File: " << filename << std::endl;
  }

  std::cout <<"Added "<<nfile<<" files"<<std::endl;
  std::cout<<"# entries in chain: "<<chain->GetEntries()<<std::endl;
  if (chain == 0) return;

  int ceNtrigger=0;
  int beNtrigger=0;
  int bceNtrigger=0;
  int ptbin,maxptbin;

  //define variables
  Int_t   Event, numberofcElectrons, numberofbElectrons,numberofbcElectrons, numberofHadrons, noTracks;   //
  Float_t celectron_id,celectron_status,celectron_pt,celectron_pz,celectron_phi,celectron_eta,celectron_y;                //track info
  Float_t belectron_id,belectron_status,belectron_pt,belectron_pz,belectron_phi,belectron_eta,belectron_y;
  Float_t bcelectron_id,bcelectron_status,bcelectron_pt,bcelectron_pz,bcelectron_phi,bcelectron_eta,bcelectron_y;
  Float_t assoh_id,assoh_status,assoh_pt,assoh_pz,assoh_phi,assoh_eta,assoh_y;

  int goodEvent = 0;


  Long64_t nentries = chain->GetEntriesFast();
  int x = 0; int n = nentries; int w = 25;

  for(int i = 0; i < nentries; i++){
    Long64_t ientry = chain->LoadTree(i);
    if (ientry < 0) break;
    TBranch *b_destep = chain->GetBranch("hf2eDecay");

    TLeaf* leaf_Event_id            = b_destep->GetLeaf("Event_id");
    TLeaf* leaf_refMult             = b_destep->GetLeaf("refMult");
    TLeaf* leaf_numberofcElectrons  = b_destep->GetLeaf("numberofcElectrons");
    TLeaf* leaf_numberofbElectrons  = b_destep->GetLeaf("numberofbElectrons");
    TLeaf* leaf_numberofbcElectrons = b_destep->GetLeaf("numberofbcElectrons");
    TLeaf* leaf_numberofHadrons     = b_destep->GetLeaf("numberofHadrons");
    TLeaf* leaf_noTracks            = b_destep->GetLeaf("noTracks");

    TLeaf* leaf_ce_id               = b_destep->GetLeaf("ce_id");
    TLeaf* leaf_ce_status           = b_destep->GetLeaf("ce_status");
    TLeaf* leaf_ce_pt               = b_destep->GetLeaf("ce_pt");
    TLeaf* leaf_ce_pz               = b_destep->GetLeaf("ce_pz");
    TLeaf* leaf_ce_phi              = b_destep->GetLeaf("ce_phi");
    TLeaf* leaf_ce_eta              = b_destep->GetLeaf("ce_eta");
    TLeaf* leaf_ce_y                = b_destep->GetLeaf("ce_y");

    TLeaf* leaf_be_id               = b_destep->GetLeaf("be_id");
    TLeaf* leaf_be_status           = b_destep->GetLeaf("be_status");
    TLeaf* leaf_be_pt               = b_destep->GetLeaf("be_pt");
    TLeaf* leaf_be_pz               = b_destep->GetLeaf("be_pz");
    TLeaf* leaf_be_phi              = b_destep->GetLeaf("be_phi");
    TLeaf* leaf_be_eta              = b_destep->GetLeaf("be_eta");
    TLeaf* leaf_be_y                = b_destep->GetLeaf("be_y");

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

//.........这里部分代码省略.........
  
  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);

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

  trackProducer *p_trk = new trackProducer("trackProd");
  if(doTracks) {
    p_trk->SetInput(chain);
    p_trk->SetTracksName("tracks");
    p_trk->SetEventObjects(fEventObjects);
  }
  
  //---------------------------------------------------------------
  //analysis modules
  //

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

  anaPFCandidates *anaPFCandJet = new anaPFCandidates("pfCandWithPFJets","pfCandWithPFJets");
  anaPFCandJet->ConnectEventObject(fEventObjects);
//  anaPFCandJet->SetHiEvtName("hiEventContainer");
  anaPFCandJet->SetParticlesName("pfParticles");
  anaPFCandJet->SetJetsName("akt4PF");
  handler->Add(anaPFCandJet);

  anaPFCandidates *anaPF = new anaPFCandidates("pfCand","pfCand");
  anaPF->ConnectEventObject(fEventObjects);
  //anaPF->SetHiEvtName("hiEventContainer");
  anaPF->SetParticlesName("pfParticles");
  anaPF->SetJetsName("akt4Gen");
  handler->Add(anaPF);

  anaPFCandidates *anaTrkCaloJet = new anaPFCandidates("tracksWithCaloJets","tracksWithCaloJets");
  anaTrkCaloJet->ConnectEventObject(fEventObjects);
  //anaTrkCaloJet->SetHiEvtName("hiEventContainer");
  anaTrkCaloJet->SetParticlesName("tracks");
  anaTrkCaloJet->SetJetsName("aktPu4Calo");
  if(doTracks) handler->Add(anaTrkCaloJet);

  anaPFCandidates *anaTrk = new anaPFCandidates("tracksWithGenJet","tracksWithGenJet");
  anaTrk->ConnectEventObject(fEventObjects);
  //anaTrk->SetHiEvtName("hiEventContainer");
  anaTrk->SetParticlesName("tracks");
  anaTrk->SetJetsName("akt4Gen");
  if(doTracks) handler->Add(anaTrk);
  
  //---------------------------------------------------------------
  //Event loop
  //---------------------------------------------------------------	  
  Long64_t entries_tot =  chain->GetEntriesFast(); //93064
  if(nentries<0) nentries = chain->GetEntries();
  // Long64_t nentries = 20;//chain->GetEntriesFast();
  Printf("nentries: %lld  tot: %lld",nentries,entries_tot);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {

    //Run producers
    // Printf("produce hiEvent");
    // p_evt->Run(jentry);   //hi event properties
    //Printf("produce pf particles");
    p_pf->Run(jentry);    //pf particles
    if(doTracks) p_trk->Run(jentry);    //tracks
    p_PUJet->Run(jentry); //jets
    
    //Execute all analysis tasks
    handler->ExecuteTask();
  }
    
  fEventObjects->Print();

  TFile *out = new TFile(outname,"RECREATE");
  TList *tasks = handler->GetListOfTasks();
  TIter next(tasks);
  anaBaseTask *obj;
  while ((obj = dynamic_cast<anaBaseTask*>(next()) ))
    if(obj->GetOutput()) obj->GetOutput()->Write(obj->GetName(),TObject::kSingleKey);
  
  out->Write();
  out->Close();
  
}

//.........这里部分代码省略.........
  // 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");
  handler->Add(ZToMuMu);

  anaZJetMCResponse *ZResp = new anaZJetMCResponse("ZJetMCResponse","ZJetMCResponse");
  ZResp->ConnectEventObject(fEventObjects);
  ZResp->SetHiEvtName("hiEventContainer");
  ZResp->SetZsName("zMuMuBosons");
  ZResp->SetJetsName(jetName);
  ZResp->SetNCentBins(3);
  handler->Add(ZResp);


 //---------------------------------------------------------------
  //Event loop
  //---------------------------------------------------------------	  
  Long64_t entries_tot =  chain->GetEntriesFast(); //93064
  if(nentries<0) lastEvent = chain->GetEntries();
  // Long64_t nentries = 20;//chain->GetEntriesFast();
  Printf("nentries: %lld  tot: %lld",nentries,entries_tot);
  for (Long64_t jentry=firstEvent; jentry<lastEvent; ++jentry) {
     if (jentry%10000==0) Printf("Processing event %d  %d",(int)(jentry), (int)(lastEvent));
    //Run producers
    // Printf("produce hiEvent");
    p_evt->Run(jentry);   //hi event properties

    //Printf("produce pf particles");
    p_gen->Run(jentry);   //gen particles
    p_PUJet->Run(jentry); //jets
    
    //Execute all analysis tasks
    handler->ExecuteTask();
  }
    
  fEventObjects->Print();

  TFile *out = new TFile(outname,"RECREATE");
  TList *tasks = handler->GetListOfTasks();
  TIter next(tasks);
  anaBaseTask *obj;
  while ((obj = dynamic_cast<anaBaseTask*>(next()) ))
    if(obj->GetOutput()) obj->GetOutput()->Write(obj->GetName(),TObject::kSingleKey);
  
  out->Write();
  out->Close();
  
}