C++ TClonesArray::Clear方法代码示例

void decayAndFill(int const kf, TLorentzVector* b, double const weight, TClonesArray& daughters)
{
   pydecay->Decay(kf, b);
   pydecay->ImportParticles(&daughters);

   TLorentzVector p1Mom;
   TLorentzVector p2Mom;
   TVector3 v00;

   int nTrk = daughters.GetEntriesFast();
   for (int iTrk = 0; iTrk < nTrk; ++iTrk)
   {
     TParticle* ptl0 = (TParticle*)daughters.At(iTrk);

     switch (ptl0->GetPdgCode())
     {
       //Will only have pions
       case 211:
         ptl0->Momentum(p1Mom);
         v00.SetXYZ(ptl0->Vx() * 1000., ptl0->Vy() * 1000., ptl0->Vz() * 1000.); // converted to μm
         break;
       case -211:
         ptl0->Momentum(p2Mom);
         break;
       default:
         break;
     }
   }
   daughters.Clear();

   fill(kf, b, weight, p1Mom, p2Mom, v00);
}

int main(int argc, char* argv[])
{ 
  //Upload the file with the data
  TFile* file = TFile::Open("/Users/Fer/Documents/traajo/samples/NeroNtuples_9.root"); // TFile::Open() instead of a constructor since it works over xrootd etc.
  //Upload the tree with the event data
  TTree *tree=(TTree*)file->Get("nero/events");

  //Create the vector to store all the particle identifiers
  std::vector<Int_t> * lepPdgId;

  //Create a variable to store all the lepton event data
  TClonesArray *leptondata = new TClonesArray("leptondata");

  //Specify where all the lepton event data will be stores
  tree->SetBranchAddress("lepP4", &leptondata);
  //Specify where all the lepton identifiers will be stored
  tree->SetBranchAddress("lepPdgId", &lepPdgId);

  //Get how many events we have to loop through
  int nentries = tree->GetEntries();

  //Loop through all the events
  for(int ientry = 0; ientry < nentries; ientry++) 
  {
    //Reset the lepton data 
    leptondata->Clear();
    //This line stores the proper data both in "leptondata" and in "lepPdgId"
    tree->GetEntry(ientry);
    
    //Only if "leptondata" is not empty continue, this is to avoid segmentation errors
    if(leptondata->GetSize() == 0) continue;

    //Loop through all the entries in the current event 
    for(int j=0; j<leptondata->GetEntriesFast()-1; j++) 
    {
        //Only if the identifier of the particle is + or - 11 (electron or antielectron) store the data in electrondata
        if(abs(lepPdgId->at(j))==11) continue;
        //Store all the data of the electron in this variable
        TLorentzVector *electrondata = (TLorentzVector *)leptondata->At(j);
        //Get some specific property such as momentum, position or energy
        cout << electrondata->E() << endl;
    }
  }
  return 0; 
}

void tclread()
{
// read file generated by tclwrite
// loop on all entries.
// histogram center of lines
   TFile *f = new TFile("tcl.root");
   TTree *T = (TTree*)f->Get("T");
   TH2F *h2 = new TH2F("h2","center of lines",40,0,1,40,0,1);

   TClonesArray *arr = new TClonesArray("TLine");
   T->GetBranch("tcl")->SetAutoDelete(kFALSE);
   T->SetBranchAddress("tcl",&arr);
   Long64_t nentries = T->GetEntries();
   for (Long64_t ev=0;ev<nentries;ev++) {
      arr->Clear();
      T->GetEntry(ev);
      Int_t nlines = arr->GetEntriesFast();
      for (Int_t i=0;i<nlines;i++) {
         TLine *line = (TLine*)arr->At(i);
         h2->Fill(0.5*(line->GetX1()+line->GetX2()), 0.5*(line->GetY1()+line->GetY2()));
      }
   }
   h2->Draw("lego");
}

//  -------------------------------------------------------------------------
//
//   ----- General Macro for R3B CALIFA S438 Data Display
//         Author: Hector Alvarez <[email protected]>
//         Last Update: 30/09/14
//         Comments:
//         MAIN DISPLAY OF S438 CALIFA CRYSTALHIT
//	
//  -------------------------------------------------------------------------
//
//   Usage: 
//      > root -l 
//      ROOT> .L plotAll.C
//      ROOT> plot("inputFile")
//     
//     where inputFile is the input file :) 
//  -------------------------------------------------------------------------
void plot(TString inputFile) {

  // CONFIGURATION
  ifstream input1;
  input1.open("./angles/petals_angels_angles.txt"); 
  Double_t polar[128]; Double_t azimuthal[128]; Int_t number=0;
  for(Int_t i=0;i<128;i++) input1 >> number >> polar[i] >> azimuthal[i]; 
  //

  gROOT->SetStyle("Default");
  //gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);

  TFile *file1 = TFile::Open(inputFile);

  //SETTINGS: reduce memory by disabling those canvas that you do not want!
  Bool_t plotMultiplicities=kTRUE; Bool_t plotEnergyvsAngle=kTRUE; Bool_t plotProtonCorr=kTRUE; Bool_t plotAddBack=kTRUE;

  //HISTOGRAMS DEFINITION FOR CRYSTALHITS
  TH1F* hMult; TH1F* hMult2; TH2F* hEnergyvsId; TH2F* hToT_EnergyvsId;
  if(plotMultiplicities) {
    hMult = new TH1F("hMult","Total CrystalHits vs crystalId",128,0,127);
    hMult2 = new TH1F("hMult2","Number of CrystalHit with max. E of the event vs crystalId",128,0,127);
    hEnergyvsId = new TH2F("hEnergyvsId","Energy vs crystalId",128,0,127,3000,-100,29900);
    hToT_EnergyvsId = new TH2F("hToT_EnergyvsId","ToT_Energy vs crystalId",128,0,127,3000,-1000,299000);
  }  
  TH2F* hP1EnergyvsPolar; TH2F* hP2EnergyvsPolar; TH2F* hP1ToT_EnergyvsPolar; TH2F* hP2ToT_EnergyvsPolar;
  if(plotEnergyvsAngle) {
    hP1EnergyvsPolar = new TH2F("hP1EnergyvsPolar","PETAL1: Energy vs polar",50,22,68,3000,-100,29900);
    hP2EnergyvsPolar = new TH2F("hP2EnergyvsPolar","PETAL2: Energy vs polar",50,22,68,3000,-100,29900);
    hP1ToT_EnergyvsPolar = new TH2F("hP1ToT_EnergyvsPolar","PETAL1: ToT_Energy vs polar",50,22,68,3000,-1000,299000);
    hP2ToT_EnergyvsPolar = new TH2F("hP2ToT_EnergyvsPolar","PETAL2: ToT_Energy vs polar",50,22,68,3000,-1000,299000);
  }
  TH2F* hToT_EnergyCorr; TH2F* hAngleCorr; TH2F* hToT_EnergyCorrLarge; TH2F* hAngleCorrLarge;
  if(plotProtonCorr) {
    hToT_EnergyCorr = new TH2F("hToT_EnergyCorr","ToT_Energy Petal1 vs ToT_Energy Petal2 (FULL COMBINATORIAL)",3000,-100,29900,3000,-100,29900);
    hEnergyCorr = new TH2F("hEnergyCorr","Energy Petal1 vs Energy Petal2 (FULL COMBINATORIAL)",3000,-100,29900,3000,-100,29900);
    hAngleCorr = new TH2F("hAngleCorr","Polar Petal1 vs polar Petal2 (FULL COMBINATORIAL)",50,22,68,50,22,68);
    hToT_EnergyCorrLarge = new TH2F("hToT_EnergyCorrLarge","ToT_Energy Petal1 vs ToT_Energy Petal2 (ONLY LARGEST ENERGY HIT)",3000,-100,29900,3000,-100,29900);
    hEnergyCorrLarge = new TH2F("hEnergyCorrLarge","Energy Petal1 vs Energy Petal2 (ONLY LARGEST ENERGY HIT)",3000,-100,29900,3000,-100,29900);
    hAngleCorrLarge = new TH2F("hAngleCorrLarge","Polar Petal1 vs polar Petal2 (ONLY LARGEST ENERGY HIT)",50,22,68,50,22,68);
  }  
  TH2F* hP1ABEnergy; TH2F* hP2ABEnergy; TH2F* hP1ABToT_Energy; TH2F* hP2ABToT_Energy;
  if(plotAddBack) {
    hP1ABEnergy= new TH2F("hP1ABEnergy","AddBack Energy in Petal1: second largest Energy vs largest Energy",3000,-100,29900,3000,-100,29900);
    hP2ABEnergy = new TH2F("hP2ABEnergy","AddBack Energy in Petal2: second largest Energy vs largest Energy",3000,-100,29900,3000,-100,29900);
    hP1ABToT_Energy = new TH2F("hP1ABToT_Energy","AddBack ToT_Energy in Petal1: second largest ToT_Energy vs largest ToT_Energy",3000,-1000,299000,3000,-1000,299000);
    hP2ABToT_Energy = new TH2F("hP2ABToT_Energy","AddBack ToT_Energy in Petal2: second largest ToT_Energy vs largest ToT_Energy",3000,-1000,299000,3000,-1000,299000);
  }

  TTree* caloTree = (TTree*)file1->Get("cbmsim");
 
  //Crystal Hits
  TClonesArray* crystalHitCA;  
  R3BCaloCrystalHit** crystalHit;
  crystalHitCA = new TClonesArray("R3BCaloCrystalHit",5);
  TBranch *branchCrystalHit = caloTree->GetBranch("CaloCrystalHit");
  if(branchCrystalHit) branchCrystalHit->SetAddress(&crystalHitCA);
  
  Int_t* crystalId; 
  Double_t* energies;
  Double_t* ToT_energies;
  Double_t maxEnergyP1=0; Double_t maxEnergyP2=0; Double_t maxToT_EnergyP1=0; Double_t maxToT_EnergyP2=0; 
  Int_t winnerEnergyP1=0; Int_t winnerEnergyP2=0; Int_t winnerToT_EnergyP1=0; Int_t winnerToT_EnergyP2=0;
  Int_t secondEnergyP1=0; Int_t secondEnergyP2=0; Int_t secondToT_EnergyP1=0; Int_t secondToT_EnergyP2=0;

  Long64_t nevents = caloTree->GetEntries();
  Int_t crystalHitsPerEvent =0;

  for(Int_t i=0;i<nevents;i++){
    if(i%100000 == 0) printf("Event:%i\n",i);
    crystalHitCA->Clear();
    caloTree->GetEvent(i);
    crystalHitsPerEvent = crystalHitCA->GetEntries(); 
    maxEnergyP1=-1; maxEnergyP2=-1; maxToT_EnergyP1=-1; maxToT_EnergyP2=-1; 
    winnerEnergyP1=-1; winnerEnergyP2=-1; winnerToT_EnergyP1=-1; winnerToT_EnergyP2=-1;
    secondEnergyP1=-1; secondEnergyP2=-1; secondToT_EnergyP1=-1; secondToT_EnergyP2=-1;

    if(crystalHitsPerEvent>0) {
      crystalHit = new R3BCaloCrystalHit*[crystalHitsPerEvent];
      crystalId = new Int_t[crystalHitsPerEvent]; 
      energies = new Double_t[crystalHitsPerEvent];
//.........这里部分代码省略.........

void PHPythia8::fillPythiaNode(Pythia8::Pythia *pythia)
{

  Int_t numpart = pythia->event.size() - 1;
  TClonesArray *fParticles = new TClonesArray("TParticle",1000);
  fParticles->Clear();
  TClonesArray &a = *((TClonesArray*)fParticles);
  
  if(verbosity > 2) cout << "PHPythia8::process_event - numparticles: " << numpart << endl;
  for (Int_t i = 1; i <= numpart; i++) 
    {
      /* This would require changing all of the container
	 new(a[i]) TParticle(pythia->event[i].id(),
	 pythia->event[i].isFinal(),
	 pythia->event[i].mother1()  - 1,
	 pythia->event[i].mother2()  - 1,
	 pythia->event[i].daughter1() - 1,
	 pythia->event[i].daughter2() - 1,
	 pythia->event[i].px(),       // [GeV/c]
	 pythia->event[i].py(),       // [GeV/c] 
	 pythia->event[i].pz(),       // [GeV/c]
	 pythia->event[i].e(),        // [GeV]
	 pythia->event[i].xProd(),    // [mm] 
	 pythia->event[i].yProd(),    // [mm]
	 pythia->event[i].zProd(),    // [mm] 
	 pythia->event[i].tProd());   // [mm/c]
      */
      new(a[i-1]) TMCParticle(pythia->event[i].status(),
			      pythia->event[i].id(),
			      pythia->event[i].mother1(), //TMC is old and has no space for mother2...
			      pythia->event[i].daughter1(),
			      pythia->event[i].daughter2(),
			      pythia->event[i].px(),       // [GeV/c]
			      pythia->event[i].py(),       // [GeV/c] 
			      pythia->event[i].pz(),       // [GeV/c]
			      pythia->event[i].e(),        // [GeV]
			      pythia->event[i].m(),
			      pythia->event[i].xProd(),    // [mm]
			      pythia->event[i].yProd(),    // [mm]
			      pythia->event[i].zProd(),    // [mm] 
			      pythia->event[i].tProd(),
			      pythia->event[i].tau());   // [mm/c]
      
      if (verbosity > 8) cout << "particle: " << i << "  " << pythia->event[i].id() << endl;
      
    }
  
  if (verbosity > 3) cout << "PHPythia8::process_event - finished filling TClonesArray" << endl;
  
  Int_t numParticles = a.GetLast() + 1;	
  int nstable = 0;
  for (Int_t ipart=0; ipart<numParticles; ipart++)
    { 
      // get the particle information
      TMCParticle *particle = (TMCParticle*)a.At(ipart);	
      phpythia->addParticle(*particle);
      if ( phpythia->isStable(ipart) ) ++nstable;
    }
  
  if (verbosity > 3) cout << "PHPythia8::process_event - nstable: " << nstable << endl;
  
  phpythiaheader->SetEvt(eventcount);		// Event number
  phpythiaheader->SetNpart(numpart);	// Number of particle entries in entire history
  phpythiaheader->SetProcessid(pythia->info.code());	// Process ID
  phpythiaheader->SetParId_1(pythia->info.id1());	// KF codes for partons participating in hard scattering
  phpythiaheader->SetParId_2(pythia->info.id2());	// KF codes for partons participating in hard scattering
  phpythiaheader->SetX1(pythia->info.x1());		// Bjorken x1,x2
  phpythiaheader->SetX2(pythia->info.x2());
  phpythiaheader->SetSvar(pythia->info.sHat());	// partonic s,t,u
  phpythiaheader->SetTvar(pythia->info.tHat());
  phpythiaheader->SetUvar(pythia->info.uHat());
  phpythiaheader->SetQsqr(pythia->info.Q2Ren());	// Q squared
  phpythiaheader->SetPt(pythia->info.pTHat());		// transverse momentum
  
  // primary vertex information. Assume position is at zero by default
  phpythiaheader->SetPrimaryVertexX(0);
  phpythiaheader->SetPrimaryVertexY(0);
  phpythiaheader->SetPrimaryVertexZ(0);
  
  // rapidity of parton parton com frame
  phpythiaheader->SetYcom(pythia->info.y());	
   
  
}

int main(int argc, char* argv[])
{ 
  //Essentials

  //Upload the file with the data, make sure the adress of the file matches the one in your computer
  TFile* file = TFile::Open("/Users/Fer/Documents/traajo/samples/NeroNtuples_9.root"); // TFile::Open() instead of a constructor since it works over xrootd etc. =D
  
  //Upload the tree with the event data
  TTree *tree=(TTree*)file->Get("nero/events");

/////////////////////////////////////////////////////

  //Lepton criteria

  //Create a variable to store all the lepton event data
  TClonesArray *leptondata = new TClonesArray("leptondata");

  //Create the vector to store all the particle identifiers
  std::vector<Int_t> * lepPdgId= 0;

  //Specify where all the lepton event data will be stored
  tree->SetBranchAddress("lepP4", &leptondata);

  //Specify where all the lepton identifiers will be stored
  tree->SetBranchAddress("lepPdgId", &lepPdgId);

  //Histogram to plot the distribution of lepton mass 
  TH1F *lepmass = new TH1F("lepmass", "Lepton mass", 50, 0, 150);

/////////////////////////////////////////////////////

  //MET criteria

  //Create a variable to store all the "met" data
  TClonesArray *metdata = new TClonesArray("metdata");

  //Specify where all the "met" data will be stored
  tree->SetBranchAddress("metP4", &metdata);

  //Histogram to plot the distribution of the transverse mass 
  TH1F *metmass = new TH1F("metmass", "Missing transverse mass", 50, 0, 150);

/////////////////////////////////////////////////////

  //Eta criteria

  //Create the variable for Eta
  Double_t eta;

/////////////////////////////////////////////////////

  //Jet criteria

  //Create a variable to store all the jet event data
  TClonesArray *jetdata = new TClonesArray("jetdata");

  //Specify where all the jet event data will be stored
  tree->SetBranchAddress("jetP4", &jetdata);

  //Histogram to plot the distribution of jet mass 
  TH1F *jetmass = new TH1F("jetmass", "Jet mass", 50, 0, 150);

  //Variable to store the amount of jets
  Double_t size;

/////////////////////////////////////////////////////

  //Histogram to plot the distribution of the whole mass 
  TH1F *wholemass = new TH1F("wholemass", "Whole mass", 50, 0, 150);  

/////////////////////////////////////////////////////

  //Histogram variables

  //Create the canvas were the histograms will be ploted
  TCanvas* c1 = new TCanvas("c1", "Masses", 600, 600);

  //Divide that canvas to plot all histograms together
  c1->Divide(2,2);

/////////////////////////////////////////////////////

  //Variables for the for loop

  //Get how many events we have to loop through
  int nentries = tree->GetEntries();

  //Create a variable to store the mass values
  Double_t mass;

  //Loop through all the events
  for(int ientry = 0; ientry < nentries; ientry++) 
  {

    //Variable of the whole data
    TLorentzVector addable_lorentz_wholedata;

    //Reset the lepton data 
    leptondata->Clear();

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

//.........这里部分代码省略.........
      if (!(isnan(info->RhoKt6PFJets) || 
            isinf(info->RhoKt6PFJets))) {
        rhoEleIso = info->RhoKt6PFJets;
      }
      EleEAEra = kDataEra_2012_Data;
    }

    //use only odd numbered events to evaluate efficiency for data. even numbered events were used for training
    //if (info->evtNum % 2 == 0 && !matchGen) continue;


    // trigger requirement               
    Bool_t passTrigger = kFALSE;
    if(dataType == 0) {
      if ((info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) == kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) == kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) == kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) passTrigger = kTRUE;
    } else if(dataType == 1) {
      if(DataEraInput == 2) {
        if(info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50)  continue;
        if(info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30)  continue;
        if(info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17)                         continue;
      }
        
      if ((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) == kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) passTrigger = kTRUE;
    }
    if(dataType != -1 && !passTrigger) continue;     
      
    // good vertex requirement
    if(!(info->hasGoodPV)) continue;

    if(matchGen) genBr->GetEntry(ientry);
      
    electronArr->Clear();
    muonArr->Clear(); 
    pfcandidateArr->Clear(); 
    electronBr->GetEntry(ientry);
    muonBr->GetEntry(ientry);
    pfcandidateBr->GetEntry(ientry);

    //********************************************************
    //Low Met Requirement
    //********************************************************
    TVector3 met;        
    if(info->pfMEx!=0 || info->pfMEy!=0) {       
      met.SetXYZ(info->pfMEx, info->pfMEy, 0);
    }
    if (met.Pt() > 25) continue;
      

    //********************************************************
    //Loop over TAG electrons
    //********************************************************
    for(Int_t i=0; i<electronArr->GetEntriesFast(); i++) {

      const higgsana::TElectron *tag = (higgsana::TElectron*)((*electronArr)[i]);
	
      if(matchGen) {
        Bool_t match1 = (higgsana::deltaR(tag->eta, tag->phi, gen->eta_1, gen->phi_1) < 0.5);
        Bool_t match2 = (higgsana::deltaR(tag->eta, tag->phi, gen->eta_2, gen->phi_2) < 0.5);
        if(!match1 && !match2)
          continue;
      }

      if(tag->pt          < 20)  continue;
      if(fabs(tag->scEta) > 2.5) continue;

//.........这里部分代码省略.........
  vector<Double_t> nSelCorrv, nSelCorrVarv;
  vector<Double_t> accv, accCorrv;
  vector<Double_t> accErrv, accErrCorrv;

  const baconhep::TTrigger triggerMenu("../../BaconAna/DataFormats/data/HLT_50nsGRun");

  //
  // loop through files
  //
  for(UInt_t ifile=0; ifile<fnamev.size(); ifile++) {  

    // Read input file and get the TTrees
    cout << "Processing " << fnamev[ifile] << " ..." << endl;
    infile = TFile::Open(fnamev[ifile]); 
    assert(infile);

    eventTree = (TTree*)infile->Get("Events"); assert(eventTree);  
    eventTree->SetBranchAddress("Info",              &info); TBranch *infoBr     = eventTree->GetBranch("Info");
    eventTree->SetBranchAddress("GenEvtInfo",         &gen); TBranch *genBr      = eventTree->GetBranch("GenEvtInfo");
    eventTree->SetBranchAddress("GenParticle", &genPartArr); TBranch *genPartBr  = eventTree->GetBranch("GenParticle");
    eventTree->SetBranchAddress("Electron",   &electronArr); TBranch *electronBr = eventTree->GetBranch("Electron");
    eventTree->SetBranchAddress("PV",   &vertexArr); TBranch *vertexBr = eventTree->GetBranch("PV");

    nEvtsv.push_back(0);
    nSelv.push_back(0);
    nSelCorrv.push_back(0);
    nSelCorrVarv.push_back(0);

    //
    // loop over events
    //
    for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
      genBr->GetEntry(ientry);
      genPartArr->Clear(); genPartBr->GetEntry(ientry);
      infoBr->GetEntry(ientry);

      Int_t glepq1=-99;
      Int_t glepq2=-99;

      if (fabs(toolbox::flavor(genPartArr, BOSON_ID))!=LEPTON_ID) continue;
      TLorentzVector *vec=new TLorentzVector(0,0,0,0);
      TLorentzVector *lep1=new TLorentzVector(0,0,0,0);
      TLorentzVector *lep2=new TLorentzVector(0,0,0,0);
      toolbox::fillGen(genPartArr, BOSON_ID, vec, lep1, lep2,&glepq1,&glepq2,1);
      if(vec->M()<MASS_LOW || vec->M()>MASS_HIGH) continue;      
      delete vec; delete lep1; delete lep2;

      vertexArr->Clear();
      vertexBr->GetEntry(ientry);
      double npv  = vertexArr->GetEntries();
      Double_t weight=gen->weight;
      if(doPU>0) weight*=h_rw->GetBinContent(h_rw->FindBin(info->nPUmean));

      nEvtsv[ifile]+=weight;        

      // trigger requirement               
      if (!isEleTrigger(triggerMenu, info->triggerBits, kFALSE)) continue;
      
      // good vertex requirement
      if(!(info->hasGoodPV)) continue;
      
      electronArr->Clear();
      electronBr->GetEntry(ientry);

      for(Int_t i1=0; i1<electronArr->GetEntriesFast(); i1++) {
  	const baconhep::TElectron *ele1 = (baconhep::TElectron*)((*electronArr)[i1]);

void compare(TString filenameA, TString filenameB)
{
    TFile *fileA = TFile::Open(filenameA);
    TTree *tree = (TTree *)fileA->Get("cbmsim");
    tree->AddFriend("cbmsimB = cbmsim", filenameB);

    TClonesArray *digisA = new TClonesArray("R3BLandDigi");
    tree->GetBranch("LandDigi")->SetAutoDelete(kFALSE);
    tree->SetBranchAddress("LandDigi", &digisA);

    TClonesArray *digisB = new TClonesArray("R3BLandDigi");
    tree->GetBranch("cbmsimB.LandDigi")->SetAutoDelete(kFALSE);
    tree->SetBranchAddress("cbmsimB.LandDigi", &digisB);



    TH1D *hNumLandDigis = new TH1D("hNumLandDigis", "N(LandDigis, A) - N(LandDigis, B)", 51, -25., 25.);
    TH1D *hEtot  = new TH1D("hEtot", "E(tot, A) - E(tot, B)", 201, -100., 100.);
    TH1D *hEtotA = new TH1D("hEtotA", "E(tot, A)", 2000, 1., 2001.);
    TH1D *hEtotB = new TH1D("hEtotB", "E(tot, B)", 2000, 1., 2001.);
    TH1D *hTmin  = new TH1D("hTmin", "T(min, A) - T(min, B)", 201, -10., 10.);

    UInt_t nentries = tree->GetEntries();
    for (UInt_t ev = 0; ev < nentries; ev++) {
        digisA->Clear();
        digisB->Clear();

        tree->GetEntry(ev);


        int nA = digisA->GetEntries();
        int nB = digisB->GetEntries();
        hNumLandDigis->Fill(nA - nB);


        R3BLandDigi *digi;
        double EtotA = 0.;
        for (int i = 0; i < nA; i++) {
            digi = (R3BLandDigi *)digisA->At(i);
            EtotA += digi->GetQdc();
        }
        hEtotA->Fill(EtotA);
        double EtotB = 0.;
        for (int i = 0; i < nB; i++) {
            digi = (R3BLandDigi *)digisB->At(i);
            EtotB += digi->GetQdc();
        }
        hEtotB->Fill(EtotB);
        hEtot->Fill(EtotA - EtotB);


        double TminA = 1000;
        for (int i = 0; i < nA; i++) {
            digi = (R3BLandDigi *)digisA->At(i);
            if (digi->GetTdc() < TminA) {
                TminA = digi->GetTdc();
            }
        }
        double TminB = 1000;
        for (int i = 0; i < nB; i++) {
            digi = (R3BLandDigi *)digisB->At(i);
            if (digi->GetTdc() < TminB) {
                TminB = digi->GetTdc();
            }
        }
        hTmin->Fill(TminA - TminB);
    }


    TCanvas *c1 = new TCanvas("c1", "", 900, 700);
    c1->Divide(2, 2);
    c1->cd(1);
    hNumLandDigis->Draw();
    c1->cd(2);
    hEtot->Draw();
    c1->cd(3);
    hTmin->Draw();
    c1->cd(4);
    hEtotA->Draw();
    hEtotB->SetLineColor(kRed);
    hEtotB->Draw("same");

}

//.........这里部分代码省略.........
      }
      EleEAEra = kDataEra_2012_Data;
    }

    //use only odd numbered events to evaluate efficiency for data. even numbered events were used for training
    //Don't need this for zee gamma
    //if (info->evtNum % 2 == 0 && !matchGen) continue;

    // trigger requirement               
    Bool_t passTrigger = kFALSE;
    if(dataType == 0 ) {
      if ((info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) == kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) == kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) == kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) passTrigger = kTRUE;
    } else if(dataType == 1 ) {
      if(info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50)  continue;
      if(info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30)  continue;
      if(info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17)                         continue;
 
      if ((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) == kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) passTrigger = kTRUE;
    } else if (dataType < 0) {
      if ((info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) == kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) == kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) == kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) passTrigger = kTRUE;
      if ((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) == kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) passTrigger = kTRUE;
    }
    if(!passTrigger) continue;     



    // good vertex requirement
    if(!(info->hasGoodPV)) continue;

    electronArr->Clear();
    muonArr->Clear(); 
    pfcandidateArr->Clear(); 
    genparticleArr->Clear(); 
    electronBr->GetEntry(ientry);
    photonBr->GetEntry(ientry);
    muonBr->GetEntry(ientry);
    pfcandidateBr->GetEntry(ientry);
    if(matchGen) {
      genparticleBr->GetEntry(ientry);
    }  


    //********************************************************
    //Loop over TAG electrons
    //********************************************************
    vector<Int_t> probeAlreadyUsed;
    for(Int_t i=0; i<electronArr->GetEntriesFast(); ++i) {
      probeAlreadyUsed.push_back(kFALSE);
    }

    for(Int_t i=0; i<electronArr->GetEntriesFast(); ++i) {

      const higgsana::TElectron *tag = (higgsana::TElectron*)((*electronArr)[i]);
	
      Bool_t TagIsEle = MatchedToStatus1Ele(tag, genparticleArr);

      if(tag->pt          < 20)  continue;
      if(fabs(tag->scEta) > 2.5) continue;
//       if(dataType == 1) {
//         if(tag->pt          < 30)  continue;
//       }

void HHToBBGGSelectionCCOneFakePhoton(const string inputfile,          // input file
                       const string outputfile,         // output directory
                       Int_t SampleType = 1
  ) {
  
  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //============================================================================================================== 
  bool printdebug = false;


  //*****************************************************************************************
  //Setup Jet Energy Corrections
  //*****************************************************************************************
  std::vector<cmsana::JetCorrectorParameters> correctionParameters;
  correctionParameters.push_back(cmsana::JetCorrectorParameters( ( getenv("CMSSW_BASE") + string("/src/CMSAna/JetEnergyCorrections/data/GR_R_52_V9_L1FastJet_AK5PF.txt")).c_str()));
  correctionParameters.push_back(cmsana::JetCorrectorParameters( ( getenv("CMSSW_BASE") + string("/src/CMSAna/JetEnergyCorrections/data/GR_R_52_V9_L2Relative_AK5PF.txt")).c_str()));
  correctionParameters.push_back(cmsana::JetCorrectorParameters( ( getenv("CMSSW_BASE") + string("/src/CMSAna/JetEnergyCorrections/data/GR_R_52_V9_L3Absolute_AK5PF.txt")).c_str()));
  cmsana::FactorizedJetCorrector *JetCorrector = new cmsana::FactorizedJetCorrector(correctionParameters);


  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================  
  Double_t nEvents = 0;

  
  //*****************************************************************************************
  // Set up output ntuple
  //*****************************************************************************************
  TFile *outFile = new TFile(outputfile.c_str(),"RECREATE"); 
  TH1F *NEvents =  new TH1F("NEvents",";;",1,-0.5,0.5);
  TH1F *N2bjets = new TH1F("N2bjets","Number of Events w/ 2 bjets",1,-0.5,0.5); // Count the number of events that have two real cjets
  TH1F *N4genjets = new TH1F("N4genjets", "Number of Events w/ at 2 bjets and >=2 genjets", 1, -0.5, 0.5); // Count the number of events that have two real bjets + at least two other jets
  TH1F *NPUMean = new TH1F("NPUMean",";NPUMean;Number of Events", 100, -0.5, 99.5);
  cmsana::HHToBBGGEventTree *outputEventTree = new cmsana::HHToBBGGEventTree;
  outputEventTree->CreateTree();


  //*****************************************************************************************
  // Set up input
  //*****************************************************************************************
  TFile *infile=0;
  TTree *eventTree=0;
  
  // Data structures to store info from TTrees
  cmsana::TEventInfo *info  = new cmsana::TEventInfo();
  TClonesArray *genparticleArr = new TClonesArray("cmsana::TGenParticle");
  TClonesArray *genjetArr = new TClonesArray("cmsana::TGenJet");
  TClonesArray *photonArr = new TClonesArray("cmsana::TPhoton");
  TClonesArray *muonArr = new TClonesArray("cmsana::TMuon");
  TClonesArray *electronArr = new TClonesArray("cmsana::TElectron");
  TClonesArray *jetArr = new TClonesArray("cmsana::TJet");
  TClonesArray *pfcandidateArr = new TClonesArray("cmsana::TPFCandidate");

  // Read input file and get the TTrees
  cout << "Processing " << inputfile << "..." << endl;
  infile = TFile::Open(inputfile.c_str(),"read");
  assert(infile);

    
  eventTree = (TTree*)infile->Get("Events"); assert(eventTree);  
  eventTree->SetBranchAddress("Info",     &info);        TBranch *infoBr     = eventTree->GetBranch("Info");
  eventTree->SetBranchAddress("GenParticle", &genparticleArr); TBranch *genparticleBr = eventTree->GetBranch("GenParticle");
  eventTree->SetBranchAddress("GenJet", &genjetArr); TBranch *genjetBr = eventTree->GetBranch("GenJet");
  cout << "NEvents = " << eventTree->GetEntries() << endl;

  // Read efficiency file for Fake Photon Rate
  TFile *glu = TFile::Open("/afs/cern.ch/work/v/vlambert/public/releases/CMSSW_5_3_9_patch3/src/PhotonEfficiencies/Efficiencies/PhotonMistagRate_gluon.root");
  TH2F *gluon_efficiencies = (TH2F*)glu->Get("MistagRate_CSVMedium_PtEta");  // access 2D histogram with efficiency weights for gluons
  TFile *qk = TFile::Open("/afs/cern.ch/work/v/vlambert/public/releases/CMSSW_5_3_9_patch3/src/PhotonEfficiencies/Efficiencies/PhotonMistagRate_quark.root");
  TH2F *quark_efficiencies = (TH2F*)qk->Get("MistagRate_CSVMedium_PtEta");  // access 2D histogram with efficiency weights for quarks
  
  // Efficiencies for charm jets faking bjets
  TFile *file4Weight = TFile::Open("/afs/cern.ch/work/v/vlambert/public/releases/CMSSW_5_3_9_patch3/src/PhotonEfficiencies/Efficiencies/BJetMistagRate_type4_nocuts.root");
  TH2F *type4Weight = (TH2F*)file4Weight->Get("MistagRate_CSVMedium_Pt_Eta");

  // Efficiencies for Photons
  TFile *FakePhoton = TFile::Open("/afs/cern.ch/work/v/vlambert/public/releases/CMSSW_5_3_9_patch3/src/PhotonEfficiencies/Efficiencies/PhotonEfficiency_PromptPhoton.root");
  TH2F *Photon_efficiency = (TH2F*)FakePhoton->Get("Efficiency_PtEta");
  
  // null vector for default four vectors
  cmsana::FourVector null(0.0,0.0,0.0,0.0);

  // loop over events
  for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
    if (ientry % 1000 == 0) cout << "Processed Event " << ientry << endl;
    infoBr->GetEntry(ientry);

    NEvents->Fill(0);
    NPUMean->Fill(info->nPUMean);

    //***********************************************************
    // Definition of Pileup Energy density
    //***********************************************************
    Double_t rho = info->RhoKt6PFJets;

    genparticleArr->Clear(); 
    genjetArr->Clear(); 

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

//.........这里部分代码省略.........
    nEvtsv.push_back(0);
    nSelv.push_back(0);
    nSelBv.push_back(0);
    nSelEv.push_back(0);
    nSelCorrv.push_back(0);
    nSelBCorrv.push_back(0);
    nSelECorrv.push_back(0);
    nSelCorrVarv.push_back(0);
    nSelBCorrVarv.push_back(0);
    nSelECorrVarv.push_back(0);
    
    for(Int_t iy=0; iy<=hHLTErr->GetNbinsY(); iy++) {
      for(Int_t ix=0; ix<=hHLTErr->GetNbinsX(); ix++) {
        hHLTErr ->SetBinContent(ix,iy,0);
        hHLTErrB->SetBinContent(ix,iy,0);
        hHLTErrE->SetBinContent(ix,iy,0);
      }
    }
    for(Int_t iy=0; iy<=hGsfSelErr->GetNbinsY(); iy++) {
      for(Int_t ix=0; ix<=hGsfSelErr->GetNbinsX(); ix++) {
        hGsfSelErr ->SetBinContent(ix,iy,0);
        hGsfSelErrB->SetBinContent(ix,iy,0);
        hGsfSelErrE->SetBinContent(ix,iy,0);
      }
    }

    //
    // loop over events
    //
    for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
//      if(ientry==10000) break;
      infoBr->GetEntry(ientry);
      genBr->GetEntry(ientry);      
      genPartArr->Clear(); genPartBr->GetEntry(ientry);
  
      if (charge==-1 && toolbox::flavor(genPartArr, -BOSON_ID)!=LEPTON_ID) continue;
      if (charge==1 && toolbox::flavor(genPartArr, BOSON_ID)!=-LEPTON_ID) continue;
      if (charge==0 && fabs(toolbox::flavor(genPartArr, BOSON_ID))!=LEPTON_ID) continue;
    
      vertexArr->Clear();
      vertexBr->GetEntry(ientry);
      double npv  = vertexArr->GetEntries();
      Double_t weight=gen->weight;
      if(doPU>0) weight*=h_rw->GetBinContent(h_rw->FindBin(info->nPUmean));

      nEvtsv[ifile]+=weight;
      
      // trigger requirement                
      if (!isEleTrigger(triggerMenu, info->triggerBits, kFALSE)) continue;
      
      // good vertex requirement
      if(!(info->hasGoodPV)) continue;
      
      electronArr->Clear();
      electronBr->GetEntry(ientry);
      Int_t nLooseLep=0;
      const baconhep::TElectron *goodEle=0;
      TLorentzVector vEle(0,0,0,0);
      TLorentzVector vElefinal(0,0,0,0);
      Bool_t passSel=kFALSE;
     
      for(Int_t i=0; i<electronArr->GetEntriesFast(); i++) {
        const baconhep::TElectron *ele = (baconhep::TElectron*)((*electronArr)[i]);
	vEle.SetPtEtaPhiM(ele->pt, ele->eta, ele->phi, ELE_MASS);
	
	//double ele_pt = gRandom->Gaus(ele->scEt*getEleScaleCorr(ele->scEta,0), getEleResCorr(ele->scEta,0));

int getSimulatedDigits(TString digitdir)
{

    AliRunLoader *rl = AliRunLoader::Open(digitdir+TString("/galice.root"));

    AliPHOSLoader *prl = (AliPHOSLoader*)rl->GetDetectorLoader("PHOS");

    prl->LoadSDigits();
    prl->LoadDigits();

    Int_t nSimEvents = rl->GetNumberOfEvents();
    TFile *mydigitsFile = new TFile("mydigits.root", "RECREATE");
    TTree *mydigitTree = new TTree("mydigitTree", "mydigitTree");
    TClonesArray * mydigits = new TClonesArray(AliPHOSDigit::Class());
    mydigitTree->Branch("Digits", &mydigits);

    for (Int_t ev = 0; ev < nSimEvents; ev++)
    {
        rl->GetEvent(ev);

        Int_t nPhosDigits = prl->SDigits()->GetEntries();


        Int_t nDigsFound = 0;

        for (Int_t iDig = 0; iDig < nPhosDigits; iDig++)
        {
            const AliPHOSDigit *digit = prl->SDigit(iDig);

            Int_t id = digit->GetId();
            if (id > 3*geom->GetNPhi()*geom->GetNZ()) continue;

            for (Int_t n = 0; n < nDigsFound; n++)
            {
                AliPHOSDigit *tmpDig = (AliPHOSDigit*)mydigits->At(n);
                if (id == tmpDig->GetId())
                {
                    *tmpDig += *digit;
                    digit = 0;
                    break;
                }
            }
            if (digit)
            {
                AliPHOSDigit *newDig = new((*mydigits)[nDigsFound]) AliPHOSDigit(-1, id, (float)0.0, (float)0.0);
                *newDig += *digit;

                nDigsFound++;
            }
        }
        for(int i = 0; i <nDigsFound; i++)
        {

            AliPHOSDigit *tmpDig = (AliPHOSDigit*)mydigits->At(i);
            Int_t relId[4];
            geom->AbsToRelNumbering(tmpDig->GetId(), relId);

            // Some necessary ugly hacks needed at the moment, need to get these numbers from OCDB
            tmpDig->SetEnergy((float)((int)(tmpDig->GetEnergy()/phosCalibData->GetADCchannelEmc(relId[0], relId[3], relId[2]))));
            //std::cout << "Sample time step: " << phosCalibData->GetSampleTimeStep() << " " << phosCalibData->GetTimeShiftEmc(relId[0], relId[3], relId[2]) << geom << std::endl;
            tmpDig->SetTime(1.5e-8);
            tmpDig->SetTimeR(1.5e-8);
        }
        mydigitTree->Fill();
        mydigits->Clear("C");
    }
    mydigitsFile->Write();
    simTree = mydigitTree;
    //mydigit
    //mydigitsFile->Close();
    return 0;
}

// Main macro function
//--------------------------------------------------------------------------------------------------
void SkimTightPlusRecoPerFile(string inputFilename, string outputFilename) 
{
  gBenchmark->Start("SkimNtuples");
    
  TTree::SetMaxTreeSize(kMaxLong64);
  

  mithep::MuonIDMVA *muonIDMVA = new mithep::MuonIDMVA();
  muonIDMVA->Initialize("BDTG method",
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml", 
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml", 
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml", 
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml", 
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml", 
                        "/home/sixie/CMSSW_analysis/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml",
                        mithep::MuonIDMVA::kIDIsoCombinedDetIso);
  
  mithep::ElectronIDMVA *electronIDMVA = new mithep::ElectronIDMVA();
  electronIDMVA->Initialize("BDTG method",
                            "MitPhysics/data/ElectronMVAWeights/Subdet0LowPt_IDIsoCombined_BDTG.weights.xml", 
                            "MitPhysics/data/ElectronMVAWeights/Subdet1LowPt_IDIsoCombined_BDTG.weights.xml", 
                            "MitPhysics/data/ElectronMVAWeights/Subdet2LowPt_IDIsoCombined_BDTG.weights.xml", 
                            "MitPhysics/data/ElectronMVAWeights/Subdet0HighPt_IDIsoCombined_BDTG.weights.xml", 
                            "MitPhysics/data/ElectronMVAWeights/Subdet1HighPt_IDIsoCombined_BDTG.weights.xml", 
                            "MitPhysics/data/ElectronMVAWeights/Subdet2HighPt_IDIsoCombined_BDTG.weights.xml",
                            mithep::ElectronIDMVA::kIDIsoCombined);



  // Don't write TObject part of the objects
  mithep::TEventInfo::Class()->IgnoreTObjectStreamer();
  mithep::TMuon::Class()->IgnoreTObjectStreamer();
  mithep::TElectron::Class()->IgnoreTObjectStreamer();
  mithep::TJet::Class()->IgnoreTObjectStreamer();
  mithep::TPhoton::Class()->IgnoreTObjectStreamer();

  // Data structures to store info from TTrees
  mithep::TEventInfo *info  = new mithep::TEventInfo();
  TClonesArray *muonArr     = new TClonesArray("mithep::TMuon");
  TClonesArray *electronArr     = new TClonesArray("mithep::TElectron");
  TClonesArray *jetArr    = new TClonesArray("mithep::TJet");
  TClonesArray *photonArr     = new TClonesArray("mithep::TPhoton");
    
  UInt_t nInputEvts = 0;
  UInt_t nPassEvts  = 0;
  UInt_t nEventsTotal = 0;

  TFile* outfile = new TFile(outputFilename.c_str(), "RECREATE");
  
  //
  // Initialize data trees and structs
  // 
  TTree *outEventTree = new TTree("Events","Events"); 
  outEventTree->Branch("Info",     &info);
  outEventTree->Branch("Muon",     &muonArr);
  outEventTree->Branch("Electron", &electronArr);
  outEventTree->Branch("PFJet",    &jetArr);
  outEventTree->Branch("Photon",   &photonArr);

  cout << "Skimming " << inputFilename << "..." << endl;
  TTree *eventTree = getTreeFromFile(inputFilename.c_str(),"Events"); 
  nEventsTotal += getNEvents(inputFilename.c_str()); 
  assert(eventTree);
    
  // Set branch address to structures that will store the info  
  eventTree->SetBranchAddress("Info",     &info);        TBranch *infoBr     = eventTree->GetBranch("Info");
  eventTree->SetBranchAddress("Muon",     &muonArr);     TBranch *muonBr     = eventTree->GetBranch("Muon");
  eventTree->SetBranchAddress("Electron", &electronArr); TBranch *electronBr = eventTree->GetBranch("Electron");
  eventTree->SetBranchAddress("PFJet",    &jetArr);      TBranch *jetBr      = eventTree->GetBranch("PFJet");
  eventTree->SetBranchAddress("Photon",    &photonArr);  TBranch *photonBr   = eventTree->GetBranch("Photon");
     
  for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) { 
    infoBr->GetEntry(ientry);
    muonArr->Clear();     muonBr->GetEntry(ientry);
    electronArr->Clear(); electronBr->GetEntry(ientry);      
    jetArr->Clear(); jetBr->GetEntry(ientry);
    photonArr->Clear(); photonBr->GetEntry(ientry);
     
    if (ientry % 100000 == 0) cout << "Events: " << ientry << endl;

    nInputEvts++;
      
    Bool_t keep = kTRUE;
    
    Double_t rho = 0;
    if (!(TMath::IsNaN(info->PileupEnergyDensity) || isinf(info->PileupEnergyDensity))) rho = info->PileupEnergyDensity;

    Int_t NTightMuons = 0;
    Int_t NRecoMuons = 0;
    for(Int_t i=0; i<muonArr->GetEntries(); i++) {
      const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
      if ( fabs(mu->eta) < 2.4
           && 
           mu->pt > 10.0
        ) {
        NRecoMuons++;

//         if ( passMuonCuts(mu)) NTightMuons++;
//.........这里部分代码省略.........

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

  //********************************************************
  // Get Tree
  //********************************************************
  eventTree = getTreeFromFile(inputFilename.c_str(),"Events"); 
  TBranch *infoBr;
  TBranch *electronBr;
  TBranch *muonBr;
  TBranch *jetBr;


  //*****************************************************************************************
  //Loop over muon Data Tree
  //*****************************************************************************************
  // Set branch address to structures that will store the info  
  eventTree->SetBranchAddress("Info",       &info);      infoBr       = eventTree->GetBranch("Info");
  eventTree->SetBranchAddress("Electron", &electronArr); electronBr = eventTree->GetBranch("Electron");
  eventTree->SetBranchAddress("Muon", &muonArr);         muonBr = eventTree->GetBranch("Muon");
  eventTree->SetBranchAddress("PFJet", &jetArr);         jetBr = eventTree->GetBranch("PFJet");
  
  for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {       	
    infoBr->GetEntry(ientry);
    if (ientry % 100000 == 0) cout << "Event " << ientry << endl;
		
    mithep::RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
    if(hasJSON && !rlrm.HasRunLumi(rl)) continue;  // not certified run? Skip to next event...
     
    //for the skimmed input, I already required the HLT bits.
    //    if (!passHLT(info->triggerBits, info->runNum, kTRUE)) continue;

    //********************************************************
    // Load the branches
    //********************************************************
    electronArr->Clear(); 
    muonArr->Clear(); 
    jetArr->Clear(); 
    electronBr->GetEntry(ientry);
    muonBr->GetEntry(ientry);
    jetBr->GetEntry(ientry);

    //event weight
    Double_t eventweight = info->eventweight * lumi;

    //********************************************************
    // TcMet
    //********************************************************
    TVector3 met;        
    if(info->tcMEx!=0 || info->tcMEy!=0) {       
      met.SetXYZ(info->tcMEx, info->tcMEy, 0);
    }

    //********************************************************
    // TcMet
    //********************************************************

    Int_t NSoftMuons = 0;
    Int_t NLeptons = 0;
    vector<Int_t> leptonType;
    vector<Int_t> leptonIndex;
    vector<Double_t> leptonPt;
    vector<Double_t> leptonEta;
    vector<Double_t> leptonPhi;
    vector<Int_t> leptonCharge;

    Int_t NJets = 0;
    const mithep::TJet *leadingJet = 0;