C++ HiForest类代码示例

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 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 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;
}

void anaJetTrackRpA()
{
    std::cout << "start working\n";
     TFile *fcrel3 ; 
     TH1D *C_rel ;
   //for ak3PF jets, applied for pPb data only 
//    TFile* fcrel3 = TFile::Open("/afs/cern.ch/user/d/dgulhan/public/Corrections/Casym_pPb_double_hcalbins_algo_ak3PF_pt100_140_jet80_alphahigh_20_phicut250.root", "readonly");
    //for akPu3PF jets and pPb data
    if(coll=="PPb") fcrel3 = TFile::Open(Form("/afs/cern.ch/user/d/dgulhan/public/Corrections/Casym_pPb_double_hcalbins_algo_%s_pt100_140_jet80_alphahigh_20_phicut250.root", algo.Data()), "readonly");
    if(coll=="PbP") fcrel3 = TFile::Open(Form("/afs/cern.ch/user/d/dgulhan/public/Corrections/Casym_Pbp_double_hcalbins_algo_%s_pt100_140_jet80_alphahigh_20_phicut250.root", algo.Data()), "readonly");
    if(fcrel3)  C_rel=(TH1D*)fcrel3->Get("C_asym");
  //  fcrel3->Close();

   TrigName=getenv("TRIG"); 
   hist_class *my_hists = new hist_class("pfjet");
    cout <<my_hists->IsMC<<endl ;
   cout <<"analysis trig = " <<TrigName <<endl ;

  //this function is to correct for UE subtraction (we are using akPu3PF algorithm)
   TF1 * fUE ;
   TF1 * fgaus ;
    fgaus=new TF1("fgaus","gaus(0)",-20,20);
   if(my_hists->IsMC==kTRUE){
    fUE = new TF1("fUE","[0]/pow(x,[1])");
  //   //for ak3PF jets
 //    fUE->SetParameters(0.4183,0.4244);
   //  //for akPu3PF jets
    if(algo=="akPu3PF") fUE->SetParameters(1.052,0.5261);
     else  fUE->SetParameters(0.4183,0.4244);
    //for gauss smearing
    fgaus->SetParameters(1,0,1);
    }
   else {
       fUE = new TF1("fUE","1-[0]/pow(x,[1])",20,300);
        //   //for ak3PF jets
  //   fUE->SetParameters(0.8648,0.8167);
   //  //for akPu3PF jets
     if(algo=="akPu3PF") fUE->SetParameters(0.3015,0.8913);
     else  fUE->SetParameters(0.8648,0.8167); 
    }

   TF1 * fVz = new TF1("fVx","[0]+[1]*x+[2]*TMath::Power(x,2)+[3]*TMath::Power(x,3)+[4]*TMath::Power(x,4)", -15., 15.);
    fVz->SetParameters(1.60182e+00,1.08425e-03,-1.29156e-02,-7.24899e-06,2.80750e-05);

    if(my_hists->IsMC==kTRUE){
      pthat=atoi(getenv("PTHAT")) ;
      ptmax=atoi(getenv("PTMAX")) ;
       cout <<"pthat = " <<pthat <<"  pthatmax =" <<ptmax <<endl ;
     }
    if(my_hists->IsMC!=kTRUE){ 
      if(coll=="HI")
        dataPath="/net/hidsk0001/d00/scratch/yjlee/merge/pbpbDijet_v20" ;//mit PbPb data path
    else {
    /*  if(TrigName=="Jet20")   
          dataPath="root://eoscms//eos/cms/store/group/phys_heavyions/krajczar/inbound/mnt/hadoop/cms/store/user/krajczar/pPb_Jet20_Full_v1" ;
   else if(TrigName=="Jet40" || TrigName=="Jet60")
          dataPath="root://eoscms//eos/cms/store/group/phys_heavyions/krajczar/inbound/mnt/hadoop/cms/store/user/krajczar/pPb_Jet40Jet60_Full_v1" ;
    else  
          dataPath="root://eoscms//eos/cms/store/group/phys_heavyions/yjlee/pPb2013/promptReco";
      */
     if(TrigName=="Jet80" || TrigName=="Jet100")
        dataPath="root://eoscms//eos/cms/store/group/phys_heavyions/yjlee/pPb2013/promptReco";
     else  
        dataPath="root://eoscms//eos/cms/store/caf/user/ymao";
    }
  } //data Path
    else { //MC analysis
        if(coll=="HI") {
          if(pthat==50||pthat==80||pthat==100||pthat==170)
             dataPath= Form("/mnt/hadoop/cms/store/user/yenjie/HiForest_v27/"); //MIT MC normial
           else 
                dataPath= Form("/mnt/hadoop/cms/store/user/yenjie/HiForest_v28/"); //MIT MC normial
       }
       else if(coll=="PPb")
      //  dataPath=Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod16/Hijing_Pythia_pt%d/HiForest_v77_merged01", pthat);
        dataPath=Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod16/Signal_Pythia_pt%d/HiForest_v77_v2_merged01", pthat);
       else if(coll=="PbP")
        dataPath=Form("/mnt/hadoop/cms/store/user/dgulhan/Pbp/HP05/prod24/Hijing_Pythia_pt%d/HiForest_v84_merged02", pthat);
        else if(coll=="PP2011")
         dataPath= Form("/net/hisrv0001/home/zhukova/scratch/HIHighPt/forest/pthat%d", pthat); //lxplus path for pp
       else {       
       if(pthat==220)
         dataPath= Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod16/Signal_Pythia_pt%d/HiForest_v77_v2_merged02", pthat); //2013 pp tracking for 5.02TeV
      else 
       dataPath= Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod16/Signal_Pythia_pt%d/HiForest_v77_v2_merged01", pthat); //2013 pp tracking for 5.02TeV
        }
    }
    if(my_hists->IsMC!=kTRUE){  //real data analysis
        if(coll=="HI")             
            intputFile="promptskim-hihighpt-hltjet80-pt90-v20.root" ; //full dataset
	else if(coll=="PP2011")  
             intputFile="hiForest2_pp_ppreco_415_90percent.root";  //! 2011 pp data rereco
        else if(coll=="PbPb")
           intputFile="PbPHiForest2_PbPbPAHighPtJet80_cent50-100_pprereco.root"; 	
        else if(coll=="PPb"){
         if(TrigName=="Jet20")
       //     intputFile="mergedJet20_KK.root" ;
            intputFile="mergedJet20_pPb_Jet20_Full_UsingKKForest_v1.root" ;
         else    if(TrigName=="Jet40" || TrigName=="Jet60")
         //   intputFile="mergedJet40Jet60_KK.root" ;
//.........这里部分代码省略.........

int IndResponse5TeVOnFlyForest(double kPt=30,const char *ksp="ppb",int ifile=0)
{

  timer.Start();

  //! Load Lib
  gSystem->Load("../HiForest_V3/hiForest_h.so");
  
  //! Load the hiforest input root file
  HiForest *c = 0;
  //! Latest
//  if(strcmp(ksp,"ppb")==0)c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod13/Hijing_Pythia_pt%.f/HiForest_v72_v01_merged01/pt%.f_prod13_v72_merged_forest_%d.root",kPt,kPt,ifile),"PythiaHijing",cPPb); 
//  if(strcmp(ksp,"ppb")==0)c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod16/Hijing_Pythia_pt%.f/HiForest_v77_merged01/pt%.f_HP04_prod16_v77_merged_forest_%d.root",kPt,kPt,ifile),"PythiaHijing",cPPb); 
  if(strcmp(ksp,"ppb")==0)c = new HiForest(Form("/mnt/hadoop/cms/store/user/kjung/pPb_Full_BTagForestMerged/pPb_Full_BTagForest%.f_Fix3output.root", kPt),"PythiaHijing",cPPb);
 // if(strcmp(ksp,"ppb")==0)c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP04/prod25/HiForest_v85_merged01/pt%.f_HP04_prod25_v85_merged_forest_%d.root", kPt,ifile),"PythiaHijing",cPPb);
// else if(strcmp(ksp,"pbp")==0)c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/Pbp/HP05/prod24/Hijing_Pythia_pt%.f/HiForest_v84_merged02/pt%.f_HP05_prod24_v84_merged_forest_%d.root",kPt,kPt,ifile),"HijingPythia",cPPb); 
// else if(strcmp(ksp,"pbp")==0)   c = new HiForest(Form("/mnt/hadoop/cms/store/user/kjung/PbpMC/HP06/mergedPthat%.f/PbpJEC_pthat%.f_outputMerged.root",kPt,kPt),"HijingPythia",cPPb);
// else if(strcmp(ksp,"pbp")==0)    c = new HiForest(Form("/mnt/hadoop/cms/store/user/kjung/PbpMC/HP08/MergedFiles/pthat%.foutputMerged.root",kPt),"HijingPythia",cPPb);
 else if(strcmp(ksp,"pbp")==0)    c = new HiForest(Form("/net/hidsk0001/d00/scratch/maoyx/pPb/CMSSW_5_3_8_HI_patch2/src/MNguyen/combinePtHatBins/CorrPbpJet%.f.root",kPt),"HijingPythia",cPPb);
  else {
    c = new HiForest(Form("/mnt/hadoop/cms/store/user/kjung/pPbForest/Signal_Pythia_pt%.f/merged/pt%.f_HP04_prod16_v77_merged_forest_Sum.root",kPt,kPt),"Pythia",cPPb);
//    if(kPt==280)c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP03/prod10/Signal_Pythia_pt%0.0f/HiForest_v63_v02_merged03/pt%0.0f_HP03_prod09_merged_forest_%d.root",kPt,kPt,ifile),"Pythia",cPPb);
//    else c = new HiForest(Form("/mnt/hadoop/cms/store/user/dgulhan/pPb/HP03/prod10/Signal_Pythia_pt%0.0f/HiForest_v63_merged03/pt%0.0f_HP03_prod09_merged_forest_%d.root",kPt,kPt,ifile),"Pythia",cPPb);
  }
  double xsection=0;
  double xup=0;
  double xsub=0;
  double maxpthat=9999;

  if(kPt==15){
    maxpthat=30;
    xup =5.335e-01;
    xsub=3.378e-02;
  }
  else if(kPt==30){
    maxpthat=50;
    xup =3.378e-02;
    xsub=3.778e-03;
  }
  else if(kPt==50){
    maxpthat=80;
  //  maxpthat=9999;
    xup =3.778e-03;
  //  xsub=0;
    xsub=4.412e-04;
  }
  else if(kPt==80){
    maxpthat=120;
    xup =4.412e-04;
    xsub=6.147e-05;
  }
  else if(kPt==120){
    maxpthat=170;
    xup=6.147e-05;
    xsub=1.018e-05;
  }else if(kPt==170){
    maxpthat=220;
    xup=1.018e-05;
    xsub=2.477e-06;
  }else if(kPt==220){
    maxpthat=280;
    xup =2.477e-06;
    xsub=6.160e-07;
  }else if(kPt==280){
    maxpthat=9999;
    xup =6.160e-07;
    xsub=0;
  }
/*    maxpthat=370;
    xup =6.160e-07;
    xsub=1.088e-07;
    }else if(kPt==370){
    maxpthat=9999;
    xup=1.088e-07; 
    xsub=0;
    }
*/
  xsection = xup-xsub;

  std::cout<<std::endl;
  std::cout<<std::endl;

  
  /*
  const int knj = 25;
  const char *cjets[knj] = {"icPu5",
			    "ak1PF","ak2PF","ak3PF","ak4PF","ak5PF","ak6PF",
			    "akPu1PF","akPu2PF","akPu3PF","akPu4PF","akPu5PF","akPu6PF",			    
			    "ak1Calo","ak2Calo","ak3Calo","ak4Calo","ak5Calo","ak6Calo",
			    "akPu1Calo","akPu2Calo","akPu3Calo","akPu4Calo","akPu5Calo","akPu6Calo"};

  const char *calgo[knj] = {"icPu5",
			    "ak1PF","ak2PF","ak3PF","ak4PF","ak5PF","ak6PF",
			    "akPu1PF","akPu2PF","akPu3PF","akPu4PF","akPu5PF","akPu6PF",			    
			    "ak1Calo","ak2Calo","ak3Calo","ak4Calo","ak5Calo","ak6Calo",
			    "akPu1Calo","akPu2Calo","akPu3Calo","akPu4Calo","akPu5Calo","akPu6Calo"};
  */
  const int knj = 12;
  const char *cjets[knj] = {"ak3PF","ak4PF","ak5PF",
			    "akPu3PF","akPu4PF","akPu5PF",
//.........这里部分代码省略.........

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 sortforestCentVz(double etCut=40, char *infname = "/mnt/hadoop/cms/store/user/velicanu/forest/HiForestTrack_v3.root")
// void sortforestCentVz(double etCut=40, char *infname = "/mnt/hadoop/cms/store/user/velicanu/HIHighPt/HIRun2011_hiHighPtTrack_PromptSkim_forest_v0/9902eec616cc8b0649a7a8bb69754615/HiForest_998_1_rJ1.root")
// void sortforestCentVz(double etCut=40, char *infname = "/net/hisrv0001/home/dav2105/hdir/HIHighPt/HIRun2011_hiHighPtTrack_PromptSkim_forest_v0/9902eec616cc8b0649a7a8bb69754615/HiForest_1000_1_8wo.root")
// void sortforestCentVz(double etCut=40, char *infname = "/net/hisrv0001/home/dav2105/hdir/HIHighPt/HIRun2011_hiHighPtTrack_PromptSkim_forest_v0/9902eec616cc8b0649a7a8bb69754615/HiForest_1001_1_TRm.root")
void frankforestsort(int startline = 0)
{
  //! Define the input and output file and HiForest
  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  // ifstream infile("/net/hidsk0001/d00/scratch/dav2105/forest/jobsort/jet80sortlist.txt");
  ifstream infile("/net/hidsk0001/d00/scratch/dav2105/forest/jobsort/frankjet80sortlist.txt");
  // ifstream infile("/net/hidsk0001/d00/scratch/dav2105/forest/jobsort/sortlist.txt");

  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);
  // c->SetOutputFile(Form("sortedHiForest_%d.root",startline));
  TFile * outf = new TFile(Form("sortedHiForest_%d.root",startline),"recreate");
  c->outf = outf;
  c->SetOutputFile("null",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;
    // if (i > 1000) break;
  }
  
  //! 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;
    // if (i > 1000) break;
	// cout << "In original loop | " << evtindecies[i] << " : (" << evtCentVz[evtindecies[i]].first <<" , "<<evtCentVz[evtindecies[i]].second << endl;
  }
  
  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);
  // qsort (evtindecies, 1000, 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();
	// cout << "In fill loop | " << evtindecies[i] << " : (" << evtCentVz[evtindecies[i]].first <<" , "<<evtCentVz[evtindecies[i]].second << endl;
    if (i%1000==0) cout <<i<<" / "<<c->GetEntries()<<" "<<c->setupOutput<<endl;
    // if (i > 1000) break;
  }

  delete c;
}

void correlate(
    const char* infname = "/data_CMS/cms/yilmaz/HiForest_HYDJET_Track8_Jet21_STARTHI53_LV1_merged_forest_0.root",
    const char* outname = "histograms_01.root",
    bool MC = 1,
    bool PbPb = 1,
    double jetEtaMax = 1.6,
    int centIndex = 0, int etaBin = 0, int leadJetPtBin = 0, int trackPtBin = 0
) {

    cout<<"Begin"<<endl;

    bool mini = 1;

    int Nevents = 1000;

    bool usePF = 1;
    int R = 3;

    double etaCOM = -0.465;

    bool doFlow = 0;
    bool doInclusiveJets = 1;
    bool doTracks = 1;
    bool doGenParticles = 1;

    bool fillTracks = 0;

    if(mini) {
        fillTracks = 0;
        doTracks = 0;
    }

    double pi = TMath::Pi();

    double trkMin = 0.5;

    double leadPtMin = 120;
    double subleadPtMin = 50;
    double dphiMin = 2.*pi/3.;

    int frame = 1; // Dijet frame for z

    int analysisId = 0;

    double ptSubLeadMin = 30;
    double ptLeadMin = 120;

    double etLeadMin[10] = {100, 100,120,150,180,200  };
    double etLeadMax[10] = {1000,120,150,180,200,1000 };

    double tkMin[10] = {4.,  4.,5.,7.,10., 20.};
    double tkMax[10] = {100.,5.,7.,10.,20.,1000.};

    TF1* gaus = new TF1("gaus","gaus",-5,5);
    gaus->SetParameter(0,1);
    gaus->SetParameter(1,0);
    gaus->SetParameter(2,1);

    TRandom* engin = new TRandom();

    //  double fitMin[10] = {100,80,60,50,100};
    //  double fitMax[10] = {250,200,150,100,250};

    double fitMin[10] = {100,90,80,70,100};
    double fitMax[10] = {300,300,300,300,300};

    TH1::SetDefaultSumw2();
    TH2::SetDefaultSumw2();
    cout<<"x"<<endl;

    string name[10] = {"c0to10","c10to20","c20to30","c30to50","c50to100","c0to30","c30to100","c0to100"};
    cout<<"x"<<endl;

    TFile* outf = new TFile(outname,"recreate");

    TNtuple *nt;
    nt = new TNtuple("nt","nt","x");

    TH3D* hAxisLead = new TH3D("hAxisLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hAxisSubLead = new TH3D("hAxisSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH3D* hCorrLead = new TH3D("hCorrLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hCorrSubLead = new TH3D("hCorrSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH3D* hGenParticleLead = new TH3D("hGenParticleLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hGenParticleSubLead = new TH3D("hGenParticleSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH1D* hPtLead = new TH1D("hPtLead","",500,0,500);
    TH1D* hPtSubLead = new TH1D("hPtSubLead","",500,0,500);

    bool pp = 0;

    double ajMin[10] = {0,    0.11, 0.22, 0.33, 0. , 0, 0};
    double ajMax[10] = {0.11, 0.22, 0.33, 1.,   1. , 0 ,0};

    int nEta = 4;
    double etaMin[10] = {0,  0,    0.5,  1.};
    double etaMax[10] = {5., 0.5,  1.,   2.};

    outf->cd();
//.........这里部分代码省略.........

void drawTrkCorrPtvCent(
                           TString outdir="fig/trkcorrv14"
)
{
   TH1::SetDefaultSumw2();
   gSystem->mkdir(outdir,kTRUE);
   float xmin=1,xmax=179.9;
   TString title="Iterative Tracking";

   /////////////////////////////////////////////////////////////////////////////////////
   // Load Histograms
   /////////////////////////////////////////////////////////////////////////////////////
   HiForest * c = new HiForest("/net/hidsk0001/d00/scratch/yjlee/merge/v27/pthat170/Dijet170_HydjetDrum_v27_mergedV1.root");
   c->doTrackCorrections = true;
   c->InitTree();
   
   TrackingCorrections * trkCorr = c->trackCorrections[0];
      
   cout << endl << "========= plot =========" << endl;
   Int_t etaPM=2; // 7 +2,-3 for |eta|<1.2, 7 =5,-6 for full eta
   Float_t jetPtMin=0;
   Int_t jetBegBin = trkCorr->jetBin_->FindBin(jetPtMin);
   Int_t jetEndBin = trkCorr->numJEtBins_;
   cout << Form("jet pt %.0f bin: ",jetPtMin) << jetBegBin << " to " << jetEndBin << endl;
   cout << "========================" << endl;
   bool doTestCorr = true;
   
   string infpath=trkCorr->sample_[0]->GetName();
   TString src=infpath.substr(infpath.find_last_of('/')+1);
   src.ReplaceAll(".root","");
   TString tag = src+"_"+trkCorr->trkCorrModule_+Form("_vs_Pt_jet%.0f_ieta%d_wts%d",jetPtMin,etaPM,trkCorr->weightSamples_);

   /////////////////////////////////////////////////////////////////////////////////////
   // Inspect Projection
   /////////////////////////////////////////////////////////////////////////////////////
   // Get Eff/fake histograms
   int numCentBin=trkCorr->numCentBins_;
   TH1D * vhCorrPtRef[2][5], *vhCorrPt[2][5];
   Int_t colors[10] = {kBlack,kRed,kYellow+2,kGreen+2,kBlue};
   Int_t styles[2] = {kFullCircle,kOpenCircle};
   for (Int_t lv=0; lv<2; ++lv) {
     for (Int_t c=0; c<numCentBin; ++c) {
       vhCorrPt[lv][c] = (TH1D*) trkCorr->InspectCorr(lv,c,c,jetBegBin,jetEndBin,2,7-etaPM-1,7+etaPM);
       vhCorrPt[lv][c]->SetMarkerStyle(styles[lv]);
       handsomeTH1(vhCorrPt[lv][c],colors[c]);
       vhCorrPt[lv][c]->SetAxisRange(xmin,xmax,"X");
     }
   }
   
   // Draw Histograms
   TCanvas * cEff = new TCanvas("cEff","cEff",500,500);
   cEff->SetLogx();
   vhCorrPt[0][0]->SetAxisRange(0,1,"Y");
   vhCorrPt[0][0]->SetTitle(";Track p_{T} (GeV/c);A #times #epsilon");
   vhCorrPt[0][0]->SetTitleOffset(1.2);
   vhCorrPt[0][0]->SetTitleSize(0.055);
   vhCorrPt[0][0]->Draw("E");
   vhCorrPt[1][0]->Draw("sameE");
   for (Int_t lv=0; lv<2; ++lv) {
     for (Int_t c=numCentBin-1; c>=0; --c) {
       vhCorrPt[lv][c]->Draw("sameE");
     }
   }
   TLegend *leg0 = new TLegend(0.16,0.786,0.46,0.92);
   leg0->SetFillStyle(0);
   leg0->SetBorderSize(0);
   leg0->SetTextSize(0.04);
   leg0->AddEntry(vhCorrPt[0][0],"PYTHIA+HYDJET","");
   if (jetPtMin >= 40) leg0->AddEntry(vhCorrPt[0][0],Form("Jet p_{T} #geq %.0f GeV/c",jetPtMin),"");
   leg0->AddEntry(vhCorrPt[0][0],Form("Track %.1f < #eta < %.1f",trkCorr->etaBin_->GetBinLowEdge(7-etaPM-1), trkCorr->etaBin_->GetBinLowEdge(7+etaPM+1)),"");
   leg0->Draw();
   TLine * l = new TLine(xmin,1,xmax,1);
   l->SetLineStyle(2);
   l->Draw();
	
   TLegend *leg = new TLegend(0.34,0.25,0.56,0.55);
   leg->SetFillStyle(0);
   leg->SetBorderSize(0);
   leg->SetTextSize(0.035);
   leg->AddEntry(vhCorrPt[0][0],title,"");
// 	leg->AddEntry(vhCorrPt[0][0],"0-5%","p");
// 	leg->AddEntry(vhCorrPt[0][1],"5-10%","p");
// 	leg->AddEntry(vhCorrPt[0][2],"10-20%","p");
// 	leg->AddEntry(vhCorrPt[0][3],"30-50%","p");
// 	leg->AddEntry(vhCorrPt[0][4],"50-90%","p");
	leg->AddEntry(vhCorrPt[0][0],"0-30%","p");
	leg->AddEntry(vhCorrPt[0][1],"30-100%","p");
   leg->Draw();
   
	drawText("CMS Simulation",0.64,0.89);
	drawText("Fake Rate",0.69,0.26);
   
   cEff->Print(outdir+"/"+tag+".gif");
   cEff->Print(outdir+"/"+tag+".pdf");

   /////////////////////////////////////////////////////////////////////////////////////
   // Inspect Events
   /////////////////////////////////////////////////////////////////////////////////////
	TCanvas * cPtHat = new TCanvas("cPtHat","cPtHat",1000,500);
	cPtHat->Divide(2,1);
//.........这里部分代码省略.........

int IndResponse(double kPt=80,const char *ksp="pbpb")
{

  timer.Start();

  //! Load Lib
  gSystem->Load("hiForest_h.so");
  
  //! Load the hiforest input root file
  HiForest *c = 0;
  if(strcmp("pp",ksp)==0)c = new HiForest(Form("/net/hisrv0001/home/icali/hadoop/Pythia/Z2/ppDijet_merged/pp276Dijet%0.0f_merged.root",kPt),"pp2012",1,1);
  else  {
    if(kPt==30)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia30_HydjetDrum_mix01_HiForest2_v19.root","pbpb2012",0,1);
    else if(kPt==50)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia50_HydjetDrum_mix01_HiForest2_v19.root","pbpb2012",0,1);
    else if(kPt==80)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia80_HydjetDrum_mix01_HiForest2_v20.root","pbpb2012",0,1);
    else if(kPt==120)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia120_HydjetDrum_mix01_HiForest2_v21_ivan.root","pbpb2012",0,1);
    else if(kPt==170)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia170_HydjetDrum_mix01_HiForest2_v19.root","pbpb2012",0,1);
    else if(kPt==200)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia200_HydjetDrum_mix01_HiForest2_v21_ivan.root","pbpb2012",0,1);
    else if(kPt==250)c = new HiForest("/d102/yjlee/hiForest2MC/Pythia250_HydjetDrum_mix01_HiForest2_v21_ivan.root","pbpb2012",0,1);
    //c = new HiForest(Form("/net/hisrv0001/home/icali/hadoop/Hydjet1.8/Z2/Dijet_merged/Dijet%0.0f_merged.root",kPt),"pbpb2012",0,1);
  }


  double xsection=0;
  double xup=0;
  double xsub=0;
  double maxpthat=9999;
  if(kPt==15){
    maxpthat=30;
    xup=1.566e-01;
    xsub=1.079e-02;
  }
  else if(kPt==30){
    maxpthat=50;
    xup=1.079e-02;
    xsub=1.021e-03;
  }
  else if(kPt==50){
    maxpthat=80;
    xup=1.021e-03;
    xsub=9.913e-05;
  }
  else if(kPt==80){
    maxpthat=120;
    xup=9.913e-05;
    xsub=1.128e-05;
  }
  else if(kPt==120){
    maxpthat=170;
    xup=1.128e-05;
    xsub=1.470e-06;
  }
  else if(kPt==170){
    maxpthat=200;
    xup=1.470e-06;
    xsub=5.310e-07;
  }
  else if(kPt==200){
    maxpthat=250;
    xup=5.310e-07;
    xsub=1.192e-07;
  }
  else if(kPt==250){
    maxpthat=300;
    xup=1.192e-07;
    xsub=3.176e-08;
  }
  else if(kPt==300){
    maxpthat=9999;
    xup=3.176e-08;
    xsub=0;
  }
  xsection = xup-xsub;


  std::cout<<std::endl;
  std::cout<<std::endl;


  //! Don't want to loop over trees which is not used in the analysis
  //! event trees
  //c->hasEvtTree=0;

  //! jet trees
  //! Switch on only the jet trees which you  require
  const char *cjets[7] = {"icPu5","ak2PF","ak3PF","ak4PF","akPu2PF","akPu3PF","akPu4PF"};
  c->SelectJetAlgo(cjets,7);
  

  //! photon tree
  //c->hasPhotonTree=0;


  //! Select only the jet branches which you are going to use
  //! This increases the speed for running over trees with lot of branches.
  //! This is currently for only jet algos and evtTree  uniformly applied to all the 

  //! Event Tree 
  const char *evlist[]={"hiBin","vz","hiHF"};
  const int kevbr = sizeof(evlist)/sizeof(const char *);
//.........这里部分代码省略.........

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

void analyzeDiJetMPT(
   TString jetAlgo = "akPu3PF",
   TString trkCol = "anaTrack",
   TString inname="/mnt/hadoop/cms/store/user/yenjie/HiForest_v27/Dijet${pthat}_HydjetDrum_v27_mergedV1.root",
   TString outname="output.root",
   bool isPP=false,
   int dataSrcType = 1, // 0 mc, 1 hi, 2 pp 2.76 TeV, 3 pp 7TeV
   double samplePtHat=0,
   double ptHatMax=9999,
   double sampleWeight = 1, // data: 1, mc: s = 0.62, b = 0.38
   double vzMax = 15,
   int maxEntries = -1,
   double leadingJetPtMin=-1,
   double subleadingJetPtMin=-1,
   double sigDPhi=-1,
   bool genJetMode=false,
   int smearCentBin = 0,
   bool doCentReWeight=false,
   TString mcfname="",
   TString datafname="output-data-Photon-v7_v30.root"
                      )
{
   ///////////////////////////////////////////////////
   // Setup Analysis
   ///////////////////////////////////////////////////
   int saveTracks = 1; // 0=none, 1=all, 10=cone
   double cutjetPt = 10;
   double cutjetEta = 2;
   double cutPtTrk=1.;
   double cutEtaTrk = 2.4;
   double trkJetAssoR=0.3;

   TString tag=Form("%s_%.0f_%.0f_%.0f_saveTrk%d_jmin%.0f_tmin%.0f_genj%d_sm%d",jetAlgo.Data(),leadingJetPtMin,subleadingJetPtMin,sigDPhi*1000,saveTracks,cutjetPt,cutPtTrk,genJetMode,smearCentBin);
   outname.ReplaceAll(".root",Form("_%s.root",tag.Data()));
   cout << "Input: " << inname << " isPP: " << isPP << endl;
   cout << "Sample pthat = " << samplePtHat << " ptHatMax = " << ptHatMax << endl;
   cout << "Track pt min = " << cutPtTrk << endl;
   cout << "skim: leading Jet > " << leadingJetPtMin << " subleading > " << subleadingJetPtMin << " dphi > " << sigDPhi  << endl;
   cout << "Genjet mode: " << genJetMode << endl;
   cout << "SmearCentBin  = " << smearCentBin << endl;
   cout << "Output: " << outname << endl;

   // Centrality reweiting
   CentralityReWeight cw(datafname,mcfname,"offlSel&&pt1>100&&pt2>0&&acos(cos(phi2-phi1))>2./3*3.14159");

   // Define the input file and HiForest
   HiForest * c = new HiForest(inname,jetAlgo.Data(),isPP);
   c->doTrackCorrections = true;
   c->doTrackingSeparateLeadingSubleading = false;
   c->InitTree();

   // intialize jet variables
   Jets* anajet = &(c->akPu3PF) ;
   
   // Output file
   cout << "Output: " << outname << endl;
   TFile *output = new TFile(outname,"recreate");
   if (doCentReWeight&&mcfname!="") {
      cw.Init(); //cw.hCentData->Draw(); cw.hCentMc->Draw("same");
   }

   // basics
   output->cd();
   TH1D * hCent = new TH1D("hCent","",40,0,40);
   TH1D * hVz = new TH1D("hVz","",60,-30,30);
   TH1D * hPtHatBeforeSel = new TH1D("hPtHatBeforeSel","",200,0,1000);
   TH1D * hPtHat = new TH1D("hPtHat","",200,0,1000);
   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
   ///////////////////////////////////////////////////
//.........这里部分代码省略.........

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 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();
}