C++ HiForest::GetEntries方法代码示例

void skimTree(char *infname = "../JetSample/hiForest_Jet80or95_GR_R_53_LV6_12Mar2014_0000CET_Track8_Jet21_0.root")
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infname);
    c->hasHltTree=0;
   c->hasPFTree=0;
   c->hasPhotonTree=0;
   c->hasTowerTree=0;
   c->hasHbheTree=0;
   c->hasEbTree=0;
   c->hasGenpTree=0;
   c->hasGenParticleTree=0;   
   c->hasAk5CaloJetTree=0;
   c->hasAkPu2CaloJetTree=0;
   c->hasAkPu3CaloJetTree=0;
   c->hasAkPu4CaloJetTree=0;
   c->hasAkPu5CaloJetTree=0;
   c->hasAkPu2JetTree=0;
   c->hasAkPu3JetTree=0;
   c->hasAkPu4JetTree=0;
   c->hasAkPu5JetTree=0;
   c->hasAkVs2PFJetTree=0;
   c->hasAkVs3PFJetTree=0;
   c->hasAkVs4PFJetTree=0;
   c->hasAkVs5PFJetTree=0;
   c->SetOutputFile("skim_jet.root");

   int filtered=0;
   // Main loop
   for (int i=0;i<c->GetEntries();i++)
   {
      c->GetEntry(i);
      if (i%1000==0) cout <<filtered<<" "<<i<<" / "<<c->GetEntries()<<endl;
      //if (c->evt.hiBin>=20) continue;
      int flag=0;
      int flag2=0;
      for (int j=0;j<c->akVs3Calo.nref;j++) {
        if (fabs(c->akVs3Calo.jteta[j])>2) continue;
	if (c->akVs3Calo.jtpt[j]>120) flag=1;
	if (c->akVs3Calo.jtpt[j]>50) flag2++;
	if (flag>=1&&flag2>=2) break;
      }
        if (flag>=1&&flag2>=2) {
	   c->FillOutput(); // Write output forest
	   filtered++;
	}  
   }

   delete c;
}

void correlateEvent(char *infname = "0.root", char *mbname = "mb.root", char *outfname = "matched_Jet.root",int startEntry=0,int endEntry=0)
{
   // Define the input file and HiForest

   HiForest *d = new HiForest(mbname);
   turnOffBranches(d);
   d->SetOutputFile(outfname);
   cout <<"good good"<<endl;

   HiForest *c = new HiForest(infname);
   turnOffBranches(c);
   int filtered=0;
   int idx=0;
   // Main loop
   if (endEntry==0) endEntry=c->GetEntries();

   for (int i=startEntry;i<=endEntry;i++)
   {
      c->GetEntry(i);
      if (i%100==0) cout <<filtered<<" "<<i<<" / "<<c->GetEntries()<<endl;
      int nMatched=0;

      while (nMatched<20)
      {
         idx++;
         cout <<i<<" / " << c->GetEntries() << " / "<<idx<<" / "<<nMatched<<"\r";
         if (idx>=d->GetEntries()) idx=0;
         if (idx==i) continue;
         d->evtTree->GetEntry(idx);
         if ((fabs(c->evt.hiBin-d->evt.hiBin)/(double)c->evt.hiBin)>0.05) continue;
         if ((fabs(c->evt.vz-d->evt.vz)>1)) continue;

         d->skimTree->GetEntry(idx);
         if (!(d->skim.pcollisionEventSelection && d->skim.pHBHENoiseFilter )) continue;        
         
         nMatched++;
         d->GetEntry(idx);
         d->FillOutput(); // Write output forest 
      }   
         cout <<std::endl;
   }

   delete c;
   delete d;
}

void stdsort(int startline = 0, string flist = "")
{
  //! Define the input and output file and HiForest
  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile(flist.data());

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }
  cout<<" here"<<endl;
  HiForest *c = new HiForest(listoffiles[startline].data(),0,0,0,0,true);
  TFile * outf = new TFile(Form("sortedHiForest_%d.root",startline),"recreate");
  c->outf = outf;
  c->SetOutputFile("null",true);

  c->LoadNoTrees();
  c->hasEvtTree = true;
  //! loop through all the events once to construct the cent,vz pair array we'll be sorting over
  cout << "Constructing the cent:vz pair array..." << endl;
  for (int i=0;i<c->GetEntries();i++)
  {
    c->GetEntry(i);
    pair<int,double> centvz;
    centvz.first = c->evt->hiBin;
    centvz.second = c->evt->vz;
    evtCentVz.push_back(centvz);
    if (i%1000==0) cout <<i<<" / "<<c->GetEntries()<<" "<<c->setupOutput<<endl;
  }
  c->ResetBooleans();
  //! Make the index array which will get sorted on first centrality
  int evtindecies[c->GetEntries()];
  for (int i=0;i<c->GetEntries();i++)
  {
    evtindecies[i] = i;
  }
  
  cout << "Sorting the cent:vz pair array..." << " "<<c->setupOutput<<endl;
  //! Sort the index array first on the centrality bin, then on the vz of the entry at that index
  qsort (evtindecies, c->GetEntries(), sizeof(int), comparecentvz);
  
  //! Now fill the tree in the new order
  cout << "Filling the tree in the sorted order..." << " "<<c->setupOutput<<endl;
  for (int i=0;i<c->GetEntries();i++)
  {
    c->GetEntry(evtindecies[i]);
    c->FillOutput();
    if (i%1000==0) cout <<i<<" / "<<c->GetEntries()<<" "<<c->setupOutput<<endl;
  }

  delete c;
}

void noiseSkim(char *infname = "/d100/yjlee/hiForest/PromptReco2011/HIHighPt/skim_Photon35/merged_HIData2011_HIHighPt_highPtExercise_photonSkim35GeVEB.root")
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infname);
   c->SetOutputFile("skim_jet.root");

   // Main loop
   for (int i=0;i<c->GetEntries();i++)
   {
      c->GetEntry(i);
      if (i%1000==0) cout <<i<<" / "<<c->GetEntries()<<endl;
      
      if (c->akPu3PF.jtpt[0]>=98||c->akPu2PF.jtpt[0]>=98&&c->akPu4PF.jtpt[0]>98){
         c->FillOutput(); // Write output forest
      }
   }

   delete c;
}

int writentuple(char *ksp="ppJet40")
{

  timer.Start();

  LoadLib();

  TString inname="";
  if(strcmp(ksp,"ppJet40")==0)inname = "root://eoscms//eos/cms/store/group/phys_heavyions/yjlee/pp2013/promptReco/PP2013_HiForest_PromptReco_JSon_Jet40Jet60_ppTrack_forestv84.root";
  else if(strcmp(ksp,"ppJet80")==0)inname = "root://eoscms//eos/cms/store/caf/user/yjlee/pp2013/promptReco/PP2013_HiForest_PromptReco_JsonPP_Jet80_PPReco_forestv82.root";

  //! Load Lib
  //gSystem->Load("/afs/cern.ch/user/p/pawan/scratch0/CMSSW_6_2_0/src/work/pPb/HiForest/V3/hiForest_h.so");

  //! Define the input file and HiForest
  //! CMSSW_5_3_3
  HiForest *c = new HiForest(inname,Form("Forest%s",ksp),cPP);
  cout<<"Loaded the hiforest tree : "<<c->GetName()<<endl;
  ShutoffBranches(c);

  //! Output file
  //! HIHighPt
  TFile *fout = new TFile(Form("ntuple_2013_%s.root",ksp),"RECREATE");  

  std::cout<<"\t"<<std::endl;
  std::cout<<"\t"<<std::endl;
  std::cout<<"**************************************************** "<<std::endl;
  std::cout<<Form("Running for %s ",ksp)<<std::endl;
  std::cout<<Form("pT  cut for %0.3f ",kptrecocut)<<std::endl;
  std::cout<<Form("eta cut for %0.3f ",ketacut)<<std::endl;
  std::cout<<"My hiForest Tree : " <<c->GetName()<<"\t Entries "<<c->GetEntries()<<std::endl;
  std::cout<<"Output file  "<<fout->GetName()<<std::endl;
  std::cout<<"**************************************************** "<<std::endl;
  std::cout<<"\t"<<std::endl;
  std::cout<<"\t"<<std::endl;

  //! shut off jet trees
  c->hasAk2CaloJetTree=0;
  c->hasAk4CaloJetTree=0;
  c->hasAk3CaloJetTree=0;
  c->hasAk5CaloJetTree=0;

  c->hasAkPu2CaloJetTree=0;
  c->hasAkPu4CaloJetTree=0;
  c->hasAkPu3CaloJetTree=0;
  c->hasAkPu5CaloJetTree=0;

  c->hasAk2PFJetTree=0;
  c->hasAk4PFJetTree=0;
  c->hasAk5PFJetTree=0;

  c->hasAkPu2PFJetTree=0;
  c->hasAkPu4PFJetTree=0;
  c->hasAkPu5PFJetTree=0;

  c->hasTrackTree=0;

  //! For jets
  Jets *mJets=0;
  Long64_t nentries = c->GetEntries();
  std::cout<<Form("# of entries in TTree for %s : ",ksp)<<nentries<<std::endl;

  string jetVars = "";
  jetVars += "evt:run:vz:trig:jet40:jet60:jet80:jet100:ntrk:pt1:raw1:eta1:phi1:chMax1:chSum1:phSum1:neSum1:pt2:raw2:eta2:phi2:chMax2:chSum2:phSum2:neSum2:pt3:raw3:eta3:phi3:chMax3:chSum3:phSum3:neSum3";
  TNtuple *ntjet=0;
  ntjet = new TNtuple("ntjet","",jetVars.data());

  for (Long64_t ievt=0; ievt<nentries;ievt++) {//! event loop
  //for (Long64_t ievt=0; ievt<100;ievt++) {//! event loop
    //! load the hiForest event
    c->GetEntry(ievt);
    
    //! events with Single vertex
    bool evSel = false;
    float trig=-9;
    if(strcmp(ksp,"ppJet40")==0){
      evSel = fabs(c->evt.vz)<15. && c->skim.pHBHENoiseFilter  && c->skim.pPAcollisionEventSelectionPA && (c->hlt.HLT_PAJet40_NoJetID_v1 || c->hlt.HLT_PAJet60_NoJetID_v1); 
      trig=1;
    }
    else if(strcmp(ksp,"ppJet80")==0){
      evSel = fabs(c->evt.vz)<15. && c->skim.pHBHENoiseFilter  && c->skim.pPAcollisionEventSelectionPA && (c->hlt.HLT_PAJet80_NoJetID_v1 || c->hlt.HLT_PAJet100_NoJetID_v1); 
      trig=2;
    }
    if(!evSel)continue;

    float pt1 = -9, pt2    = -9, pt3    = -9,
      raw1    = -9, raw2   = -9, raw3   = -9,
      eta1    = -9, eta2   = -9, eta3   = -9,
      phi1    = -9, phi2   = -9, phi3   = -9,
      chMax1  = -9, chMax2 = -9, chMax3 = -9,
      chSum1  = -9, chSum2 = -9, chSum3 = -9,
      phSum1  = -9, phSum2 = -9, phSum3 = -9,
      neSum1  = -9, neSum2 = -9, neSum3 = -9;

    
    float run   = c->evt.run;
    float evt   = c->evt.evt;
    float vz    = c->evt.vz;
    
    float jet40  = c->hlt.HLT_PAJet40_NoJetID_v1;
//.........这里部分代码省略.........

void test(char * tag= "0", char *infName = "/d102/yjlee/hiForest2MC/Pythia80_HydjetDrum_mix01_HiForest2_v22_simTrack05.root")
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infName,"",cPPb);
   c->hasPFTree=0;
   c->hasPhotonTree=0;
   c->hasTowerTree=0;
   c->hasHbheTree=0;
   c->hasEbTree=0;
   c->hasGenpTree=0;
   c->hasGenParticleTree=0;   
   c->hasAkPu2CaloJetTree=0;
   c->hasAkPu3CaloJetTree=0;
   c->hasAkPu4CaloJetTree=0;
//   c->doTrackCorrections=1;
//   c->InitTree();
   
   // Output file
   TFile *output = new TFile(Form("output-%s.root",tag),"recreate");
   
   // Output
   TTree * t = new TTree("t","gammajet");
   
   JetData data(t,1);

   HistoData histos_MergedGeneralCalo("MergedGeneral");
   HistoData histos2_MergedGeneral("MergedGeneral2");   // phi dependent corr
   
   TH1D *hWeight = new TH1D("hWeight","",1000,0,100);
   TH1D *hWeight2 = new TH1D("hWeight2","",1000,0,100);
   TH1D *hPt = new TH1D("hPt","",100,0,100);
   TH1D *hNoWPt = new TH1D("hNoWPt","",100,0,100);


   TNtuple *nt = new TNtuple("nt","","m:eta:phi:pt:pt1:pt2:ch1:ch2:phi1:phi2:dxy1:dxy2:hiBin:N");
   TNtuple *ntEvt = new TNtuple("ntEvt","","N");
   nt->SetAutoFlush(30000);
   cout <<nt->GetAutoFlush()<<endl;
   TCanvas *cc = new TCanvas("cc","",600,600);   
//	 nt->SetCircular(1000);
   // Main loop

   TLorentzVector *v2 = new TLorentzVector;
   TLorentzVector *v = new TLorentzVector;
   TLorentzVector phi;

   for (int i=0;i<c->GetEntries()/1.;i++) {
      c->GetEntry(i);
      if (!c->selectEvent()) continue;
      if (i%1000==0){
         cout <<i<<" / "<<c->GetEntries()<<endl;
       } 
      int N=0;   
      for (int j=0;j<c->track.nTrk;j++) {
         if (!c->selectTrack(j)) continue;
         if (fabs(c->track.trkEta[j])>2.4) continue;
         if (fabs(c->track.trkPt[j])<0.4) continue;
         N++;
      }   
      ntEvt->Fill(N);
      for (int j=0;j<c->track.nTrk;j++) {
         if (!c->selectTrack(j)) continue;
         if (fabs(c->track.trkPt[j])<1) continue; 
//      if (fabs(c->track.trkDxy1[j]/c->track.trkDxyError1[j])<1) continue;
         for (int k=j+1;k<c->track.nTrk;k++) {
            if (j==k) continue;
            if (!c->selectTrack(k)) continue;
//            if (c->track.trkCharge[k]==c->track.trkCharge[j]) continue; 
     	    if (fabs(c->track.trkPt[k])<1) continue;
            v->SetPtEtaPhiM(c->track.trkPt[j],c->track.trkEta[j],c->track.trkPhi[j],0.493677);
            v2->SetPtEtaPhiM(c->track.trkPt[k],c->track.trkEta[k],c->track.trkPhi[k],0.493677);
//          v2->SetPtEtaPhiM(c->track.trkPt[k],c->track.trkEta[k],c->track.trkPhi[k],0.13957);
//          v->SetPtEtaPhiM(c->track.trkPt[j],c->track.trkEta[j],c->track.trkPhi[j],0.938272);
//          v->SetPtEtaPhiM(c->track.trkPt[j],c->track.trkEta[j],c->track.trkPhi[j],0.13957);
//          v2->SetPtEtaPhiM(c->track.trkPt[k],c->track.trkEta[k],c->track.trkPhi[k],0.13957);
  	    phi = (*v) + (*v2);
//         if ((phi.M())>5) {
            if ((phi.M())>1.2||phi.M()<0.0) {
//	       phi.Delete();
	       continue;
	    }    
	    nt->Fill(phi.M(),phi.Eta(),phi.Phi(),phi.Pt(),v->Pt(),v2->Pt(),c->track.trkCharge[j],c->track.trkCharge[k],v->Phi(),v2->Phi(),c->track.trkDxy1[j],c->track.trkDxy1[k],c->evt.hiBin,N);
//	       phi.Delete();
   }
      }

      //cout <<data.mpt<<endl;
      t->Fill();
   }
  // t->Write();
   histos_MergedGeneralCalo.calcEff();
   histos2_MergedGeneral.calcEff();
   output->Write();
   output->Close();
}

//.........这里部分代码省略.........
   corr.load("trkCorr_HIN12017");
   corr.setOption1(true);
   corr.setOption2(true);
  }
*/
   TrackCorrector2D corr("/afs/cern.ch/work/y/ymao/analysis/AsymmetryPA/Corrections/trackCorrections_HIN12017v5_XSecWeighted.root");
   if(doTrackCorrections) corr.load("trkCorr_HIN12017");

    Evts * offSel = &(c->evt); 

    Skims * my_skim = &(c->skim); 

    Hlts * trigSel = &(c->hlt); 
    //jet tree
//    if(coll=="HI") 
        Jets * my_ct = &(c->akPu4PF); 
   //     Jets * my_ct = &(c->akPu4Calo); 
 //   else 
 //   Jets * my_ct = &(c->ak3PF);
    //  Jets * jetthres =  &(c->icPu5); 
   //track tree
    Tracks * my_tr = &(c->track);
    
    //GenParticle tree
    GenParticles * my_GenPart = &(c->genparticle);

    int curr_bin = nbin-1 ;
    int Nevents[nbin] = {0} ;
      Int_t Nevt_40_60[nbin] = {0} ;
      Int_t Nevt_60_75[nbin] = {0} ;
      Int_t Nevt_75_95[nbin] = {0} ;
      Int_t Nevt_95_120[nbin] = {0} ;
      Int_t Nevt_120[nbin] = {0} ;
    cout <<"Number of events ="<<c->GetEntries()<<endl ;
    for(int evi = 0; evi < c->GetEntries(); evi++) {
        c->GetEntry(evi);
        int noise_evt = my_skim->pHBHENoiseFilter ;
        //        int ecal_noise = my_skim->phiEcalRecHitSpikeFilter ;
        //        if(ecal_noise==0) continue ;
        
        double vz = offSel->vz ;
        int hiBin = offSel->hiBin ;
        weight = 1. ;      
 
         int pileup_Gplus ;
        double HFbin ;
        HFbin = (offSel->hiHFplusEta4)+(offSel->hiHFminusEta4);
    /*    if(coll=="PPb")
          HFbin = (offSel->hiHFplus);
        else 
         HFbin = (offSel->hiHFminus);
*/
       if(my_hists->IsMC!=kTRUE){
            int evt_sel ;
            pileup_Gplus = my_skim->pVertexFilterCutGplus ;
           if(coll=="PbPb"|| coll=="HI"|| coll=="PP2011")  evt_sel = my_skim->pcollisionEventSelection ;
            else  evt_sel = my_skim->pPAcollisionEventSelectionPA;
         //   if(evt_sel==0) continue ;
            if(!evt_sel) continue ;
           //for 0-90% selection using HF sum energy at |eta|>4
          //  if(HFbin<2.87) continue ; 
           //for 0-90% selection using HF sum energy at eta>4
         //   if(HFbin<2.66) continue ; 
        }
        if(my_hists->IsMC!=kTRUE){
          //  if(noise_evt==0) continue ;

//.........这里部分代码省略.........
            hDeltaRAll[nj][icen] = new TH1F(Form("hDeltaRAll%d_%d",nj,icen),Form("#DeltaR (all) for algorithm %s cent %d",cjets[nj],icen),100,0,1);
            hDeltaRSel[nj][icen] = new TH1F(Form("hDeltaRSel%d_%d",nj,icen),Form("#DeltaR (sel) for algorithm %s cent %d",cjets[nj],icen),100,0,1);
            
            for(int ir=0;ir<25;ir++){
                //! Response vs DeltaR
                hRspVsDeltaR[nj][icen][ir] = new TH1F(Form("hRspVsDeltaR%d_%d_%d",nj,icen,ir),Form(" <recopt/refpt> vs. #DeltaR (%d) algorithm %s cent %d",ir,cjets[nj],icen),rbins,rbinl,rbinh);
            }
        }//! icen
        //for eta clsoure in different jet pt bins
        for(int ipt=0;ipt<bins;ipt++){
            hratiocorrrefpt_etaptbin[nj][ipt]= new TH2F(Form("hratiocorrrefpt_etaptbin%d_%d",nj,ipt),Form("Gen matched jet Reco jet / Gen jet p_{T} (corr.) distribution jet vs eta in pt bin %d %s",ipt,cjets[nj]),
                                                  100,-5.0,5.0,rbins,rbinl,rbinh);
            hratiorawrefpt_etaptbin[nj][ipt]= new TH2F(Form("hratiorawrefpt_etaptbin%d_%d",nj,ipt),Form("Gen matched jet Reco jet RawPt/ Gen jet p_{T} (corr.) distribution jet vs eta in pt bin %d %s",ipt,cjets[nj]),
                                                  100,-5.0,5.0,rbins,rbinl,rbinh); 
        }
        //for different HFplusEta4 bins
        for(int ihf=0;ihf<nhfbin;ihf++){
            hratiocorrrefpt_genhfb[nj][ihf]= new TH2F(Form("hratiocorrrefpt_genhfbin%d_%d",nj,ihf),Form("Gen matched jet Reco jet / Gen jet p_{T} (corr.) distribution jet hiHF bin %d %s",ihf,cjets[nj]),
                                                      500,0,1000,rbins,rbinl,rbinh);
        }
        
    }//! nj
  std::cout<<"Initialized the histograms " <<std::endl;
  ///////////////////////////////////////////////////////////////////////////////////////// 


  //! Centrality reweighting function

  //! vertex z reweighting




  Long64_t nentries = c->GetEntries();
  std::cout<<Form("# of entries in TTree for %s : ",ksp)<<nentries<<std::endl;
  std::cout<<std::endl;
  hEvt->Fill(2,nentries);

  //! weight  for the merging of the samples for different pT hat bins
  Float_t wxs = xsection/(nentries/100000.);
  
  Int_t iEvent=0; 
  for (Long64_t ievt=0; ievt<nentries;ievt++) {//! event loop
    //for (Long64_t ievt=0; ievt<5000;ievt++) {//! event loop
    //! load the hiForest event
    c->GetEntry(ievt);
/*
    int hiBin       = c->evt.hiBin;
    float vz        = c->evt.vz;
    float hiHF      = c->evt.hiHF;
    int ntracks     = c->evt.hiNtracks;
    float HFplusEta = c->evt.hiHFplusEta4;
    float HFminusEta = c->evt.hiHFminusEta4 ; 
    float HFplusEta4 = HFplusEta+HFminusEta;
    //! testing
    //if(hiBin>4 && strcmp(ksp,"pbpb")==0)continue;

    //! apply vertex cut
    if(fabs(vz)>kVzcut)continue;


    //! Centrality bin
    if(hiBin<0 || hiBin>100)continue;
    int multb=GetMultBin(ntracks);
    int hiHFb=GetHFplusEta4Bin(HFplusEta4);  
*/      

void partonFlavorTree(double tag=0, char *infName = "/mnt/hadoop/cms/store/user/dgulhan/HIMC/Jet80/Track8_Jet21_STARTHI53_LV1/merged3/HiForest_Pythia_Hydjet_Jet80_Track8_Jet21_STARTHI53_LV1_merged_forest_0.root",collisionType cType = cPbPb)
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infName,"",cType);
   c->hasHltTree=0;
   c->hasPFTree=0;
   c->hasPhotonTree=0;
   c->hasTowerTree=0;
   c->hasHbheTree=0;
   c->hasEbTree=0;
   c->hasGenpTree=0;
   c->hasGenParticleTree=0;   
   c->hasAk5CaloJetTree=0;
   c->hasAkPu2CaloJetTree=0;
   c->hasAkPu3CaloJetTree=0;
   c->hasAkPu4CaloJetTree=0;
   c->hasAkPu5CaloJetTree=0;
   c->hasAkPu2JetTree=0;
   c->hasAkPu3JetTree=0;
   c->hasAkPu4JetTree=0;
   c->hasAkPu5JetTree=0;
   c->hasAkVs2PFJetTree=0;
   c->hasAkVs3PFJetTree=0;
   c->hasAkVs4PFJetTree=0;
   c->hasAkVs5PFJetTree=0;
//   c->doTrackCorrections=1;
//   c->InitTree();
   
   // Output file
   TFile *output = new TFile(Form("output-%.0f.root",tag),"recreate");
   
   // Output
   TTree * t = new TTree("t","gammajet");
   
   JetData data(t,1);

   TH1D *hWeight = new TH1D("hWeight","",1000,0,100);
   TH1D *hWeight2 = new TH1D("hWeight2","",1000,0,100);
   TH1D *hPt = new TH1D("hPt","",150,0,150);
   TH1D *hGenPt = new TH1D("hGenPt","",150,0,150);
   TH1D *hNoWPt = new TH1D("hNoWPt","",150,0,150);
   TH1D *hRmin = new TH1D("hRmin","",100,0,10);
   // Main loop
   for (int i=0;i<c->GetEntries();i++) {
      c->GetEntry(i);
      data.hiBin = c->evt.hiBin;
      if (i % 1000 == 0) cout <<i<<" / "<<c->GetEntries()<<endl;
      
      // Event selection
//      if (fabs(c->evt.vz)>15) continue;
//      if (!c->selectEvent()) continue;
      // Select leading and subleading jet
      
      data.nJet=0;
      for (int j=0;j<c->akVs3Calo.nref;j++) {
         if (fabs(c->akVs3Calo.jteta[j])>2) continue;
	 if (c->akVs3Calo.jtpt[j]<100) continue;
	 data.jetPt[data.nJet]=c->akVs3Calo.jtpt[j];
	 data.jetEta[data.nJet]=c->akVs3Calo.jteta[j];
	 data.jetPhi[data.nJet]=c->akVs3Calo.jtphi[j];
	 data.jetFlavor[data.nJet]=c->akVs3Calo.refparton_flavor[j];
	 data.refPt[data.nJet]=c->akVs3Calo.refpt[j];
//	 cout <<data.nJet<<endl;
	 
	 data.jetTrkMult1[data.nJet]=0;
	 data.jetTrkMult2[data.nJet]=0;
	 data.jetTrkMult3[data.nJet]=0;
	 data.jetTrkMult4[data.nJet]=0;
	 data.jetTrkMult5[data.nJet]=0;
	 data.jetTrkMult6[data.nJet]=0;
	 data.jetTrkMult7[data.nJet]=0;
	 data.jetTrkMult8[data.nJet]=0;
	 data.jetTrkSum1[data.nJet]=0;
	 data.jetTrkSum2[data.nJet]=0;
	 data.jetTrkSum3[data.nJet]=0;
	 data.jetTrkSum4[data.nJet]=0;
	 data.jetTrkSum5[data.nJet]=0;
	 data.jetTrkSum6[data.nJet]=0;
	 data.jetTrkSum7[data.nJet]=0;
	 data.jetTrkSum8[data.nJet]=0;
	 data.jetTrkW1[data.nJet]=0;
	 data.jetTrkW2[data.nJet]=0;
	 data.jetTrkW3[data.nJet]=0;
	 data.jetTrkW4[data.nJet]=0;
	 data.jetTrkW5[data.nJet]=0;
	 data.jetTrkW6[data.nJet]=0;
	 data.jetTrkW7[data.nJet]=0;
	 data.jetTrkW8[data.nJet]=0;
	 
	 
         for (int k=0;k<c->track.nTrk;k++) {
            if (fabs(c->track.trkEta[k])>2.4) continue;
   	    if ((c->track.trkPt[k])<4) continue;
            if (!(c->track.highPurity[k] &&
              fabs(c->track.trkDxy1[k]/c->track.trkDxyError1[k])<3.0 &&
              fabs(c->track.trkDz1[k]/c->track.trkDzError1[k])<3.0 && 
              (c->track.trkPtError[k]/c->track.trkPt[k])<0.1)) continue;

            
            double dR = c->getDR( c->track.trkEta[k], c->track.trkPhi[k], data.jetEta[data.nJet], data.jetPhi[data.nJet]);
//.........这里部分代码省略.........

void dijetTrack(double tag=0, char *infName = "/d102/yjlee/hiForest2MC/Pythia80_HydjetDrum_mix01_HiForest2_v22_simTrack05.root")
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infName);

   // Turn off the non-necessary trees
   c->hasPFTree=0;
   c->hasPhotonTree=0;
   c->hasTowerTree=0;
   c->hasHbheTree=0;
   c->hasEbTree=0;
   c->hasGenpTree=0;
   c->hasGenParticleTree=0;   
   c->hasAkPu2CaloJetTree=0;
   c->hasAkPu3CaloJetTree=0;
   c->hasAkPu4CaloJetTree=0;
   c->doTrackCorrections=1;
   c->InitTree();
   
   // Output file
   TFile *output = new TFile(Form("output-%.0f.root",tag),"recreate");
   
   // Book a output tree
   TTree * t = new TTree("t","gammajet");

   // My tree data format   
   JetData data(t,1);

   // Book histograms
   TH1F *hJetTrackDphi    = new TH1F("hJetTrackDphi"   , ";#Delta #phi;Weighted Entries", 100,  0, TMath::Pi());
   TH1F *hJetTrackDeta    = new TH1F("hJetTrackDeta"   , ";#Delta #eta;Weighted Entries", 100, -2, 2);
   TH2F *hJetTrack2D      = new TH2F("hJetTrack2D"     , ";#Delta #eta; #Delta #phi"    , 100, -2, 2,100,0, TMath::Pi());
   TH1F *hGenJetTrackDphi = new TH1F("hGenJetTrackDphi", ";#Delta #phi;Weighted Entries", 100,  0, TMath::Pi());
   TH1F *hGenJetTrackDeta = new TH1F("hGenJetTrackDeta", ";#Delta #eta;Weighted Entries", 100, -2, 2);
   TH2F *hGenJetTrack2D   = new TH2F("hGenJetTrack2D"  , ";#Delta #eta; #Delta #phi"    , 100, -2, 2,100,0, TMath::Pi());
   
   // Main loop
   for (int i=0;i<c->GetEntries();i++) {
      c->GetEntry(i);
      data.hiBin = c->evt.hiBin;
      if (i % 1000 == 0) cout <<i<<" / "<<c->GetEntries()<<endl;
      data.leadingJetPt = -1;
      data.subleadingJetPt = -1;
      data.leadingJetIt = -1;
      data.subleadingJetIt = -1;
      data.genleadingJetPt = -1;
      data.gensubleadingJetPt = -1;
      
      // =================================
      // Select leading and subleading jet
      // =================================
      for (int j=0;j<c->icPu5.nref;j++) {
         if (fabs(c->icPu5.jteta[j])>2) continue;
         if (c->icPu5.jtpt[j]>data.leadingJetPt) {
	    data.leadingJetPt = c->icPu5.jtpt[j];
	    data.leadingJetEta = c->icPu5.jteta[j];
	    data.leadingJetPhi = c->icPu5.jtphi[j];
	    data.leadingJetIt = j;
	 }   
	 if (c->icPu5.jtpt[j]>data.subleadingJetPt && c->icPu5.jtpt[j] < data.leadingJetPt) {
	    data.subleadingJetPt = c->icPu5.jtpt[j];
	    data.subleadingJetEta = c->icPu5.jteta[j];
	    data.subleadingJetPhi = c->icPu5.jtphi[j];
	    data.subleadingJetIt = j;
         }
	 if (c->icPu5.jtpt[j]<data.subleadingJetPt) break;	 
      } 

      // =================================================
      // Select generator level leading and subleading jet
      // =================================================
      for (int j=0;j<c->icPu5.ngen;j++) {
         if (fabs(c->icPu5.geneta[j])>2) continue;
         if (c->icPu5.genpt[j]>data.genleadingJetPt) {
	    data.genleadingJetPt = c->icPu5.genpt[j];
	    data.genleadingJetEta = c->icPu5.geneta[j];
	    data.genleadingJetPhi = c->icPu5.genphi[j];
	 }   
	 if (c->icPu5.genpt[j]>data.gensubleadingJetPt && c->icPu5.genpt[j] < data.genleadingJetPt) {
	    data.gensubleadingJetPt  = c->icPu5.genpt[j];
	    data.gensubleadingJetEta = c->icPu5.geneta[j];
	    data.gensubleadingJetPhi = c->icPu5.genphi[j];
         }
	 if (c->icPu5.genpt[j]<data.gensubleadingJetPt) break;	 
      } 
      
      // ====================================================================
      // Now we have leading and subleading jet from generator and reco level
      // Loop over tracks in reco level
      // ====================================================================
      for (int j=0;j<c->track.nTrk;j++) {
         if (fabs(c->track.trkEta[j])>2.4) continue;
	 if (fabs(c->track.trkPt[j]) <0.5) continue;
	 double dphi1 = acos(cos(c->track.trkPhi[j]-data.leadingJetPhi));
         double deta1 = (c->track.trkEta[j]-data.leadingJetEta);

         // Fill the histogram, with track pT as the weight of the entry.
	 hJetTrackDphi->Fill(dphi1       , c->track.trkPt[j]);
	 hJetTrackDeta->Fill(deta1       , c->track.trkPt[j]);
	 hJetTrack2D  ->Fill(deta1, dphi1, c->track.trkPt[j]);
//.........这里部分代码省略.........

void randomcone(int condor_iter, string flist, string tag, int centmin, int centmax)
{
  using namespace std;
  
  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile(flist.data());
  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }
  
  HiForest * c = new HiForest(listoffiles[condor_iter].data(),"forest",cPPb,0);
  
  // TH2D * etaphi = new TH2D();
  int nevents = c->GetEntries();
  const int netabins = 48;
  
  TRandom3 r;
  float rConeEta = 0;
  float rConePhi = 0;
  float dR = 0.3;
  float thisdR = 0;
  float pfPtSum;
  float pfVsPtSum;
  float pfVsPtInitialSum;
  
  TFile * outf = new TFile(Form("randcone_%s_centmin%d_centmax%d_%d.root",tag.data(),centmin,centmax,condor_iter),"recreate");
  TH1D * hpfPtSum = new TH1D("hpfPtSum","pfPtSum;#Sigma p_{T}",240,0,120);
  TH1D * hpfVsPtSum = new TH1D("hpfVsPtSum","pfVsPtSum;#Sigma p_{T}",200,-50,50);
  TH1D * hpfVsPtInitialSum = new TH1D("hpfVsPtInitialSum","hpfVsPtInitialSum;#Sigma p_{T}",200,-50,50);
  
  TH1D * hpfPtSumEta[netabins];
  TH1D * hpfVsPtSumEta[netabins];
  TH1D * hpfVsPtInitialSumEta[netabins];
  for(int i = 0 ; i < netabins ; ++i)
  {
    hpfPtSumEta[i] = new TH1D(Form("hpfPtSumEta_%d",i),Form("pfPtSum_%d;#Sigma p_{T}",i),200,-100,100);
    hpfVsPtSumEta[i] = new TH1D(Form("hpfVsPtSumEta_%d",i),Form("pfVsPtSum_%d;#Sigma p_{T}",i),200,-100,100);
    hpfVsPtInitialSumEta[i] = new TH1D(Form("hpfVsPtInitialSumEta_%d",i),Form("hpfVsPtInitialSum_%d;#Sigma p_{T}",i),200,-100,100);
  }
  
  TH1D * hConeEta = new TH1D("hConeEta","hConeEta;#eta",netabins,-2.4,2.4);
  
  for(int i = 0 ; i < nevents ; ++i)
  {
    if(i%1000==0) cout<<i<<"/"<<nevents<<endl;
    c->GetEntry(i);
    if(c->evt.hiBin < centmin || c->evt.hiBin >= centmax) continue;
    rConeEta = 4.8*(r.Rndm()-0.5); // uniform dist. +- 2.4 eta
    rConePhi = 2*TMath::Pi()*(r.Rndm()-0.5); // uniform dist. +- pi
    hConeEta->Fill(rConeEta);
    int whichbin = hConeEta->FindBin(rConeEta);
    pfPtSum = 0.0;
    pfVsPtSum = 0.0;
    pfVsPtInitialSum = 0.0;
    for(int j = 0 ; j < c->pf.nPFpart ; ++j)
    {
      thisdR = getdR(rConeEta,c->pf.pfEta[j],rConePhi,c->pf.pfPhi[j]);
      if(thisdR < dR)
      {
        pfPtSum+=c->pf.pfPt[j];
        pfVsPtSum+=c->pf.pfVsPt[j];
        pfVsPtInitialSum+=c->pf.pfVsPtInitial[j];
      }
    }
    hpfPtSum->Fill(pfPtSum);
    hpfVsPtSum->Fill(pfVsPtSum);
    hpfVsPtInitialSum->Fill(pfVsPtInitialSum);
    hpfPtSumEta[whichbin-1]->Fill(pfPtSum);
    hpfVsPtSumEta[whichbin-1]->Fill(pfVsPtSum);
    hpfVsPtInitialSumEta[whichbin-1]->Fill(pfVsPtInitialSum);
  }
  TH1D * hmeanpfPtSumEta = new TH1D("hmeanpfPtSumEta","hmeanpfPtSumEta;#eta;mean #Sigma p_{T}",netabins,-2.4,2.4);
  TH1D * hmeanpfVsPtSumEta = new TH1D("hmeanpfVsPtSumEta","hmeanpfVsPtSumEta;#eta;mean #Sigma p_{T}",netabins,-2.4,2.4);
  TH1D * hmeanpfVsPtInitialSumEta = new TH1D("hmeanpfVsPtInitialSumEta","hmeanpfVsPtInitialSumEta;#eta;mean #Sigma p_{T}",netabins,-2.4,2.4);
  for(int i = 1 ; i <= netabins ; ++i)
  {
    hmeanpfPtSumEta->SetBinContent(i,hpfPtSumEta[i-1]->GetMean());
    hmeanpfPtSumEta->SetBinError(i,hpfPtSumEta[i-1]->GetMeanError());
    hmeanpfVsPtSumEta->SetBinContent(i,hpfVsPtSumEta[i-1]->GetMean());
    hmeanpfVsPtSumEta->SetBinError(i,hpfVsPtSumEta[i-1]->GetMeanError());
    hmeanpfVsPtInitialSumEta->SetBinContent(i,hpfVsPtInitialSumEta[i-1]->GetMean());
    hmeanpfVsPtInitialSumEta->SetBinError(i,hpfVsPtInitialSumEta[i-1]->GetMeanError());
  }
  TH1D * one = new TH1D("one","",1,0,1);
  one->Fill(0.5);
  outf->Write();
  outf->Close();
}

//.........这里部分代码省略.........
      hEtaAll[nj][icen] = new TH1F(Form("hEtaAll%d_%d",nj,icen),Form("Denominator eta  for algorithm %s cent %d",cjets[nj],icen),20,-2.0,2.0);
      hPhiAll[nj][icen] = new TH1F(Form("hPhiAll%d_%d",nj,icen),Form("Denominator  phi  for algorithm %s cent %d",cjets[nj],icen),20,-pi,pi);
      
      hPtSel [nj][icen] = new TH1F(Form("hPtSel%d_%d",nj,icen),Form("Numerator pT for algorithm %s cent %d",cjets[nj],icen),40,30,430);
      hEtaSel[nj][icen] = new TH1F(Form("hEtaSel%d_%d",nj,icen),Form("Numerator eta  for algorithm %s cent %d",cjets[nj],icen),20,-2.0,2.0);
      hPhiSel[nj][icen] = new TH1F(Form("hPhiSel%d_%d",nj,icen),Form("Numerator  phi  for algorithm %s cent %d",cjets[nj],icen),20,-pi,pi);

      hDeltaR[nj][icen]    = new TH1F(Form("hDeltaR%d_%d",nj,icen),Form("#DeltaR for algorithm %s cent %d",cjets[nj],icen),100,0,1);
      hDeltaRAll[nj][icen] = new TH1F(Form("hDeltaRAll%d_%d",nj,icen),Form("#DeltaR (all) for algorithm %s cent %d",cjets[nj],icen),100,0,1);
      hDeltaRSel[nj][icen] = new TH1F(Form("hDeltaRSel%d_%d",nj,icen),Form("#DeltaR (sel) for algorithm %s cent %d",cjets[nj],icen),100,0,1);

      for(int ir=0;ir<25;ir++){
	//! Response vs DeltaR
        hRspVsDeltaR[nj][icen][ir] = new TH1F(Form("hRspVsDeltaR%d_%d_%d",nj,icen,ir),Form(" <recopt/refpt> vs. #DeltaR (%d) algorithm %s cent %d",ir,cjets[nj],icen),rbins,rbinl,rbinh);
      }
    }//! icen
  }//! nj
  std::cout<<"Initialized the histograms " <<std::endl;
  ///////////////////////////////////////////////////////////////////////////////////////// 


  //! Centrality reweighting function
  //TF1* fcen = new TF1("fcen","exp(-1.0*pow(x+1.11957e+01,2)/pow(1.34120e+01,2)/2)",0,40);
  TF1 *fcen = new TF1("fcen","[0]*exp([1]+[2]*x+[3]*x*x+[4]*x*x*x)",0,40);
  fcen->SetParameters(2.10653e-02,5.61607,-1.41493e-01,1.00586e-03,-1.32625e-04);

  //! vertex z reweighting
  TF1 *fVz = new TF1("fVz","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x");
  if(strcmp(ksp,"pp")==0) fVz->SetParameters(8.41684e-01,-2.58609e-02,4.86550e-03,-3.10581e-04,2.07918e-05);
  else  fVz->SetParameters(7.62788e-01,-1.13228e-02,5.85199e-03,-3.04550e-04,4.43440e-05);


  Long64_t nb = 0;
  Long64_t nentries = c->GetEntries();
  std::cout<<Form("# of entries in TTree for %s : ",ksp)<<nentries<<std::endl;
  std::cout<<std::endl;
  hEvt->Fill(2,nentries);

  //! weight  for the merging of the samples for different pT hat bins
  Float_t wxs = xsection/(nentries/100000.);
  
  Int_t iEvent=0; 
  for (Long64_t ievt=0; ievt<nentries;ievt++) {//! event loop
  //for (Long64_t ievt=0; ievt<1;ievt++) {//! event loop
    //! load the hiForest event
    nb += c->GetEntry(ievt);

    int hiBin       = c->evt.hiBin;
    float vz        = c->evt.vz;
    float hiHF      = c->evt.hiHF;

    //! testing
    //if(hiBin>4 && strcmp(ksp,"pbpb")==0)continue;

    //! apply vertex cut
    if(fabs(vz)>kVzcut)continue;


    //! Centrality bin
    if(hiBin<0 || hiBin>39)continue;
    int centb=-1;
    double wcen=1;
    double wvz=fVz->Eval(vz);
    if(strcmp(ksp,"pbpb")==0){
      centb=GetCentBin(hiBin);
      wcen = fcen->Eval(hiBin);

void photon_JEC()
{
    TH1::SetDefaultSumw2();

    Double_t alphabins[] = {0.075, 0.125, 0.175, .25, .35};

    TH1D *RdRmcalpha[nptbins];
    TProfile *Ratio[2][nptbins];

    for(int i = 0; i < nptbins; i++)
    {
        TString name;
        name = "RdRmcalpha_";
        name += ptbins[i];
        RdRmcalpha[i] = new TH1D(name,name+";#alpha;R_{data}/R_{MC}",4,alphabins);
        name = "Rd_";
        name += ptbins[i];
        Ratio[0][i] = new TProfile(name,"",4,alphabins);
        name = "Rmc_";
        name += ptbins[i];
        Ratio[1][i] = new TProfile(name,"",4,alphabins);
    }

    TH1D *alphas = new TH1D("alphas",";#alpha = p_{T}^{Jet3}/p_{T}^{#gamma}",100,0,0.8);

    TH1D *hPhotonPt[2];
    hPhotonPt[0] = new TH1D("hPhotonPt_data",";photon p_{T}",200,0,400);
    hPhotonPt[1] = (TH1D*)hPhotonPt[0]->Clone("hPhotonPt_mc");
    TH1D *hPhotonEta[2];
    hPhotonEta[0] = new TH1D("hPhotonEta_data",";photon #eta",26,-3,3);
    hPhotonEta[1] = (TH1D*)hPhotonEta[0]->Clone("hPhotonEta_mc");

    TH1D *hJet2Pt[2];
    hJet2Pt[0] = new TH1D("hJet2Pt_data",";jet 2 p_{T}",200,0,400);
    hJet2Pt[1] = (TH1D*)hJet2Pt[0]->Clone("hJet2Pt_mc");
    TH1D *hJet2Eta[2];
    hJet2Eta[0] = new TH1D("hJet2Eta_data",";jet 2 #eta",26,-3,3);
    hJet2Eta[1] = (TH1D*)hJet2Eta[0]->Clone("hJet2Eta_mc");

    TH1D *hJet3Pt[2];
    hJet3Pt[0] = new TH1D("hJet3Pt_data",";jet 3 p_{T}",200,0,400);
    hJet3Pt[1] = (TH1D*)hJet3Pt[0]->Clone("hJet3Pt_mc");
    TH1D *hJet3Eta[2];
    hJet3Eta[0] = new TH1D("hJet3Eta_data",";jet 3 #eta",26,-3,3);
    hJet3Eta[1] = (TH1D*)hJet3Eta[0]->Clone("hJet3Eta_mc");


    TString files[2];

    files[0] = "/mnt/hadoop/cms/store/user/luck/pA2013_forests/PA2013_HiForest_PromptReco_HLT_Photon40.root";
    files[1] = "/mnt/hadoop/cms/store/user/luck/pA2013_MC/HiForest2_QCDPhoton30_5020GeV_100k.root";

    bool montecarlo[2] = {false, true};

    //loop over files
    //do the smaller MC file first, for some reason I have segfaults when I process the MC second.
    for(int ii = 1; ii > -1; ii--)
    {
        HiForest *c = new HiForest(files[ii], "Forest", cPPb, montecarlo[ii]);
        c->InitTree();

        //loop over events in each file
        int nentries = c->GetEntries();
        for(int jentry = 0; jentry<nentries; jentry++)
        {
            if (jentry% 1000 == 0)  {
                printf("%d / %d\n",jentry,nentries);
            }

            c->GetEntry(jentry);

            //collision selection
            if( !((montecarlo[ii] || c->skim.pHBHENoiseFilter) && c->skim.phfPosFilter1 && c->skim.phfNegFilter1 && c->skim.phltPixelClusterShapeFilter && c->skim.pprimaryvertexFilter) )
                continue;

            if(c->photon.nPhotons == 0)
                continue;

            //loop over photons in the event
            Float_t leadingPt = 40; //minPt is 40GeV
            Int_t leadingIndex = -1;
            for(int i = 0; i<c->photon.nPhotons; i++)
            {
                if(c->photon.pt[i] > leadingPt)
                {
                    leadingPt = c->photon.pt[i];
                    leadingIndex = i;
                }
            }

            if(leadingIndex == -1) // no photons above minPt
                continue;
            if ( TMath::Abs(c->photon.eta[leadingIndex]) > 1.3 ) // barrel photons only
                continue;
            if(c->photon.ecalRecHitSumEtConeDR04[leadingIndex] > TMath::Min(3., 0.05*c->photon.energy[leadingIndex]) )
                continue;
            if(c->photon.hcalTowerSumEtConeDR04[leadingIndex] > TMath::Min(2.4, 0.05*c->photon.energy[leadingIndex]) )
                continue;
            //if(c->photon.trkSumPtHollowConeDR04[leadingIndex] > 0.10*leadingPt )
            //continue;
//.........这里部分代码省略.........

void ffStudy(double tag=0, char *infName = "/d102/yjlee/hiForest2MC/Pythia80_HydjetDrum_mix01_HiForest2_v22_simTrack05.root")
{
   // Define the input file and HiForest
   HiForest *c = new HiForest(infName);
   c->hasPFTree=0;
   c->hasPhotonTree=0;
   c->hasTowerTree=0;
   c->hasHbheTree=0;
   c->hasEbTree=0;
   c->hasGenpTree=0;
   c->hasGenParticleTree=0;   
   c->hasAkPu2CaloJetTree=0;
   c->hasAkPu3CaloJetTree=0;
   c->hasAkPu4CaloJetTree=0;
   c->doTrackCorrections=1;
   c->InitTree();
   
   // Output file
   TFile *output = new TFile(Form("output-%.0f.root",tag),"recreate");
   
   // Output
   TTree * t = new TTree("t","gammajet");
   
   JetData data(t,1);

   HistoData histos_MergedGeneralCalo("MergedGeneral");
   HistoData histos2_MergedGeneral("MergedGeneral2");   // phi dependent corr
   
   TH1D *hWeight = new TH1D("hWeight","",1000,0,100);
   TH1D *hWeight2 = new TH1D("hWeight2","",1000,0,100);
   TH1D *hPt = new TH1D("hPt","",100,0,100);
   TH1D *hNoWPt = new TH1D("hNoWPt","",100,0,100);
   
   // Main loop
   for (int i=0;i<c->GetEntries();i++) {
      c->GetEntry(i);
      data.hiBin = c->evt.hiBin;
      if (i % 1000 == 0) cout <<i<<" / "<<c->GetEntries()<<endl;
      data.leadingJetPt = -1;
      data.subleadingJetPt = -1;
      data.leadingJetIt = -1;
      data.subleadingJetIt = -1;
      data.genleadingJetPt = -1;
      data.gensubleadingJetPt = -1;
      
      // Event selection
      
      
      // Select leading and subleading jet
      for (int j=0;j<c->icPu5.nref;j++) {
         if (fabs(c->icPu5.jteta[j])>2) continue;
         if (c->icPu5.jtpt[j]>data.leadingJetPt) {
	    data.leadingJetPt = c->icPu5.jtpt[j];
	    data.leadingJetEta = c->icPu5.jteta[j];
	    data.leadingJetPhi = c->icPu5.jtphi[j];
	    data.leadingJetIt = j;
	 }   
	 if (c->icPu5.jtpt[j]>data.subleadingJetPt && c->icPu5.jtpt[j] < data.leadingJetPt) {
	    data.subleadingJetPt = c->icPu5.jtpt[j];
	    data.subleadingJetEta = c->icPu5.jteta[j];
	    data.subleadingJetPhi = c->icPu5.jtphi[j];
	    data.subleadingJetIt = j;
         }
	 if (c->icPu5.jtpt[j]<data.subleadingJetPt) break;	 
      } 

      // Select generator level leading and subleading jet
      for (int j=0;j<c->icPu5.ngen;j++) {
         if (fabs(c->icPu5.geneta[j])>2) continue;
         if (c->icPu5.genpt[j]>data.genleadingJetPt) {
	    data.genleadingJetPt = c->icPu5.genpt[j];
	    data.genleadingJetEta = c->icPu5.geneta[j];
	    data.genleadingJetPhi = c->icPu5.genphi[j];
	 }   
	 if (c->icPu5.genpt[j]>data.gensubleadingJetPt && c->icPu5.genpt[j] < data.genleadingJetPt) {
	    data.gensubleadingJetPt = c->icPu5.genpt[j];
	    data.gensubleadingJetEta = c->icPu5.geneta[j];
	    data.gensubleadingJetPhi = c->icPu5.genphi[j];
         }
	 if (c->icPu5.genpt[j]<data.gensubleadingJetPt) break;	 
      } 
      
      //if (data.subleadingJetPt<50||data.subleadingJetPt>80) continue;
      
      // MPT calculation
      data.mpt = 0;
      data.cormpt = 0;
      data.cormpt2 = 0;
      data.genMpt = 0;
      data.genPMpt = 0;
      
      for (int j=0;j<c->track.nTrk;j++) {
         if (fabs(c->track.trkEta[j])>2.4) continue;
	 if (fabs(c->track.trkPt[j])<0.5) continue;
	 double mptTrk = -c->track.trkPt[j]*cos(c->track.trkPhi[j]-data.leadingJetPhi);
	 data.mpt+=mptTrk;
      }

      for (int j=0;j<c->track.nTrk;j++) {
         if (fabs(c->track.trkEta[j])>2.4) continue;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   TH2D * hJetPt2D = new TH2D("hJetPt2D","",100,0,500,100,0,500);
   TH1D * hJDPhi = new TH1D("hJDPhi","",40,0,Pi());
   TH1D * hAj = new TH1D("hAj","",32,0,0.8);   
   TH1D * hTrkPt = new TH1D("hTrkPt",";p_{T} (GeV/c)",200,0,200);
   TH1D * hTrkCorrPt = new TH1D("hTrkCorrPt",";p_{T} (GeV/c)",200,0,200);
   TH1D * hGenpPt = new TH1D("hGenpPt",";p_{T} (GeV/c)",200,0,200);
   TH1D* smearingHist=0;
   if ( smearCentBin > 0 ) {
      smearingHist = new TH1D("smearingH","",100,-2,2);
   }
   EvtSel evt;
   DiJet gj;
   TTree * tgj = new TTree("tgj","dijet jet tree");
   BookGJBranches(tgj,evt,gj);

   // pp triggers
   int HLT_Jet40_v1=0;
   if (dataSrcType==2) {
      c->hltTree->SetBranchAddress("HLT_Jet40_v1",&HLT_Jet40_v1);
   } else if (dataSrcType==3) {
      c->hltTree->SetBranchAddress("HLT_Jet40_v1",&HLT_Jet40_v1);
   }

   ////////////////////////
   // Smearing Setup
   ////////////////////////
   if (smearCentBin>0) {
     LoadParameters();
   }

   ///////////////////////////////////////////////////
   // Main loop
   ///////////////////////////////////////////////////
   if (maxEntries<0||maxEntries > c->GetEntries()) maxEntries = c->GetEntries();
   for (int i=0;i<maxEntries;i++) {
      c->GetEntry(i);
      // Event Info
      evt.run = c->hlt.Run;
      evt.evt = c->hlt.Event;
      evt.cBin = c->evt.hiBin;
      if (isPP) evt.cBin = 39;
      evt.evtPlane = c->evt.hiEvtPlanes[21];
      evt.nJ = anajet->nref;
      evt.nT = c->track.nTrk;
      evt.trig = (c->hlt.HLT_HIJet80_v1 > 0);
      evt.offlSel = (c->skim.pcollisionEventSelection > 0);
      if (!c->hasSkimTree) evt.offlSel = (c->evt.hiNtracks>0 && c->evt.hiHFplus>=4 && c->evt.hiHFminus>=4);
      evt.noiseFilt = (c->skim.pHBHENoiseFilter > 0);
      if (dataSrcType>1) {
         if (dataSrcType==2) {
            evt.trig = (HLT_Jet40_v1>0);
         } else if (dataSrcType==3) {
            evt.trig = (HLT_Jet40_v1>0);
         }
//          evt.offlSel = (c->skim.phfCoincFilter && c->skim.ppurityFractionFilter);
      }
      if (!c->hasHltTree) evt.trig = true;
      evt.anaEvtSel = evt.offlSel;
      if (dataSrcType>0) evt.anaEvtSel = evt.anaEvtSel && evt.trig && evt.noiseFilt;
      evt.vz = c->track.vz[1];
      // Get Centrality Weight
      if (doCentReWeight) evt.weight = cw.GetWeight(evt.cBin);
      else evt.weight = 1;
      evt.npart = getNpart(evt.cBin);
      evt.ncoll = getNcoll(evt.cBin);
      evt.sampleWeight = sampleWeight/maxEntries; // for different mc sample, 1 for data