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

void AnalyseEvents(ExRootTreeReader *treeReader, MyPlots *plots)
{
  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");

  Long64_t allEntries = treeReader->GetEntries();

  cout << "** Chain contains " << allEntries << " events" << endl;

  Jet *jet[2];
  MissingET *met;
  Electron *electron;

  Long64_t entry;

  Int_t i;

  // Loop over all events
  for(entry = 0; entry < allEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);

    // Analyse two leading jets
    if(branchJet->GetEntriesFast() >= 2)
    {
      jet[0] = (Jet*) branchJet->At(0);
      jet[1] = (Jet*) branchJet->At(1);

      plots->fJetPT[0]->Fill(jet[0]->PT);
      plots->fJetPT[1]->Fill(jet[1]->PT);
    }

    // Analyse missing ET
    if(branchMissingET->GetEntriesFast() > 0)
    {
      met = (MissingET*) branchMissingET->At(0);
      plots->fMissingET->Fill(met->MET);
    }

    // Loop over all electrons in event
    for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
    {
      electron = (Electron*) branchElectron->At(i);
      plots->fElectronPT->Fill(electron->PT);
    }
  }
}

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

//______________________________________________________________________________
void alice_esd_read()
{
   // Read tracks and associated clusters from current event.

   AliESDRun    *esdrun = (AliESDRun*)    esd->fESDObjects->FindObject("AliESDRun");
   TClonesArray *tracks = (TClonesArray*) esd->fESDObjects->FindObject("Tracks");

   // This needs further investigation. Clusters not shown.
   // AliESDfriend *frnd   = (AliESDfriend*) esd->fESDObjects->FindObject("AliESDfriend");
   // printf("Friend %p, n_tracks:%d\n", frnd, frnd->fTracks.GetEntries());

   if (track_list == 0) {
      track_list = new TEveTrackList("ESD Tracks"); 
      track_list->SetMainColor(6);
      //track_list->SetLineWidth(2);
      track_list->SetMarkerColor(kYellow);
      track_list->SetMarkerStyle(4);
      track_list->SetMarkerSize(0.5);

      gEve->AddElement(track_list);
   }

   TEveTrackPropagator* trkProp = track_list->GetPropagator();
   trkProp->SetMagField( 0.1 * esdrun->fMagneticField ); // kGaus to Tesla

   gProgress->Reset();
   gProgress->SetMax(tracks->GetEntriesFast());
   for (Int_t n=0; n<tracks->GetEntriesFast(); ++n)
   {
      AliESDtrack* at = (AliESDtrack*) tracks->At(n);

      // If ITS refit failed, take track parameters at inner TPC radius.
      AliExternalTrackParam* tp = at;
      if (! trackIsOn(at, kITSrefit)) {
         tp = at->fIp;
      }

      TEveTrack* track = esd_make_track(trkProp, n, at, tp);
      track->SetAttLineAttMarker(track_list);
      track_list->AddElement(track);

      // This needs further investigation. Clusters not shown.
      // if (frnd)
      // {
      //     AliESDfriendTrack* ft = (AliESDfriendTrack*) frnd->fTracks->At(n);
      //     printf("%d friend = %p\n", ft);
      // }
      gProgress->Increment(1);
   }

   track_list->MakeTracks();
}

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 digitsPHOS(Int_t nevents, Int_t nfiles)
{

 TH1F * hadc = new TH1F ("hadc", "hadc", 100, -10, 200);
 TH1F * hadcLog = new TH1F ("hadclog", "hadclog", 100, -2, 4);
 IlcRunLoader* runLoader = IlcRunLoader::Open("gilc.root","Event","READ");
 IlcPHOSLoader * phosLoader = dynamic_cast<IlcPHOSLoader*>(runLoader->GetLoader("PHOSLoader"));
 
 for (Int_t ievent=0; ievent <nevents; ievent++) {
  // for (Int_t ievent = 0; ievent < runLoader->GetNumberOfEvents(); ievent++) {
   runLoader->GetEvent(ievent) ;
   phosLoader->CleanDigits() ; 
   phosLoader->LoadDigits("READ") ;
   TClonesArray * digits    = phosLoader->Digits() ;
   printf("Event %d contains %d digits\n",ievent,digits->GetEntriesFast());
   
   
   for (Int_t j = 0; j < digits->GetEntries(); j++) {
     
     IlcPHOSDigit* dig = dynamic_cast<IlcPHOSDigit*> (digits->At(j));
     //cout << dig->GetEnergy() << endl;
     hadc->Fill(dig->GetEnergy());
     if(dig->GetEnergy()>0)
       hadcLog->Fill(TMath::Log10(dig->GetEnergy()));
     
   }
   
 }
 
 TFile fc("digits.PHOS.root","RECREATE");
 fc.cd();
 hadc->Write();
 hadcLog->Write();
 fc.Close();


}

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

void pythia8_susy() {

  Int_t maxEvts = 100; // Maximo numero de eventos

  char* path = gSystem->ExpandPathName("$PYTHIA8DATA");
  if (gSystem->AccessPathName(path)) {
    Warning("pythia8.C", 
            "Environment variable PYTHIA8DATA must contain path to pythi8100/xmldoc directory !");
    return;
  }
  
  // Load libraries
  gSystem->Load("$PYTHIA8/lib/libpythia8");
  gSystem->Load("$PYTHIA8/lib/liblhapdfdummy");
  
  gSystem->Load("libEG");
  gSystem->Load("libEGPythia8");
  
  //Definir archivo de salida
  TFile * outfile = new TFile("eventos_pythia8_SUSY.root","RECREATE");
  
  // Array of particles
  TClonesArray* particles = new TClonesArray("TParticle", 5000);
  
  //Definir el TTree
  TTree*tree= new TTree("tree","Arbol con particulas segun Pythia8");
  tree->Branch("particles",&particles);
  
  // Create pythia8 object
  TPythia8* pythia8 = new TPythia8();
  
  //*Configurar: Aqui seleccione el proceso que quiere simular    
  pythia8->ReadString("SUSY:all = on"); //Todos los procesos susy posibles
  //pythia8->ReadString("SUSY:qqbar2chi+-chi0 = on"); //Un proceso en especial
  
  //Importante: pasar a Pythia8 el nombre del archivo SLHA
  pythia8->ReadString("SLHA:file = SUSY_LM2_sftsht.slha"); //insertar aqui el nombre del archivo SLHA

  // Initialize 
  
  pythia8->Initialize(2212 /* p */, 2212 /* p */, 7000. /* TeV */);
  
  int iev = 0;
  
  // Event loop
  
  while( iev < maxEvts ) {
    
    pythia8->GenerateEvent();
    if (iev < 1) pythia8->EventListing();
    pythia8->ImportParticles(particles,"All");
    
    Int_t np = particles->GetEntriesFast();
    
    // Particle loop
    
    for (Int_t ip = 0; ip < np; ip++) {
      
      TParticle* part = (TParticle*) particles->At(ip);
      Int_t ist = part->GetStatusCode();
      Int_t pdg = part->GetPdgCode();
      
    }

   tree->Fill();
   ++iev;
     	 
  }
  
  pythia8->PrintStatistics();
  
  outfile->Write();
  outfile->Close();
  
   
}

Bool_t monojet::Process(Long64_t entry)
{
    GetEntry(entry);

    if( entry % 100000 == 0 ) cout << "Processing event number: " << entry << endl;
    //cout << "Processing event number: " << entry << endl;

    // To make the processing fast, apply a very looose selection
    if (((TLorentzVector*)((*metP4)[0]))->Pt() < 40. or jetP4->GetEntries() < 1) return kTRUE;

    //this is the type tree
    tm->SetValue("run",runNum);
    tm->SetValue("event",eventNum);
    tm->SetValue("lumi",lumiNum);

    float dR = 0.;
    TClonesArray *tightLep;
    TClonesArray *cleanJet;
    TClonesArray *cleanTau;

    tightLep = new TClonesArray("TLorentzVector",20);
    cleanJet = new TClonesArray("TLorentzVector",20);
    cleanTau = new TClonesArray("TLorentzVector",20);

    std::vector<bool>  jetMonojetId_clean;
    jetMonojetId_clean.clear();
    std::vector<bool>  jetMonojetIdLoose_clean;
    jetMonojetIdLoose_clean.clear();
    //std::vector<float> jetPuId_clean;
    //jetPuId_clean.clear();

    std::vector<float>  tauId_clean;
    tauId_clean.clear();
    std::vector<float>  tauIso_clean;
    tauIso_clean.clear();

    int n_tightlep = 0;

    // ********* Leptons ********** //
    for(int lepton = 0; lepton < lepP4->GetEntries(); lepton++) {
        TLorentzVector* Lepton = (TLorentzVector*) lepP4->At(lepton);
        // check if this is a tight lep, and check the overlap

        //iso_1 = divide(input_tree.lepIso[0],input_tree.lepP4[0].Pt())
        //if (input_tree.lepTightId[0]==0 or iso_1 > 0.12): continue

        if (Lepton->Pt() > 20. && (*lepTightId)[lepton] > 1) {
            n_tightlep +=1;
            new ( (*tightLep)[tightLep->GetEntriesFast()]) TLorentzVector(Lepton->Px(), Lepton->Py(), Lepton->Pz(), Lepton->Energy());

            //check overlap with jets
            for(int j = 0; j < jetP4->GetEntries(); j++) {
                TLorentzVector* Jet = (TLorentzVector*) jetP4->At(j);
                dR = deltaR(Lepton,Jet);
                if (dR > dR_cut) {
                    new ( (*cleanJet)[cleanJet->GetEntriesFast()]) TLorentzVector(Jet->Px(), Jet->Py(), Jet->Pz(), Jet->Energy());
                    jetMonojetId_clean.push_back((*jetMonojetId)[j]);
                    jetMonojetIdLoose_clean.push_back((*jetMonojetIdLoose)[j]);
                    //jetPuId_clean.push_back((*jetPuId)[j]);
                }
            }
            //check overlap with taus
            for(int tau = 0; tau < tauP4->GetEntries(); tau++) {
                TLorentzVector* Tau = (TLorentzVector*) tauP4->At(tau);
                dR = deltaR(Lepton,Tau);
                if (dR > dR_cut) new ( (*cleanTau)[cleanTau->GetEntriesFast()]) TLorentzVector(Tau->Px(), Tau->Py(), Tau->Pz(), Tau->Energy());
                tauId_clean.push_back((*tauId)[tau]);
                tauIso_clean.push_back((*tauIso)[tau]);
            } // tau overlap
        } // tight lepton selection
    }//lepton loop

    tm->SetValue("n_tightlep",n_tightlep);

    TLorentzVector fakeMET;

    // Z Selection
    TLorentzVector Z;
    if(lepP4->GetEntries() == 2 && n_tightlep > 0) {
        if (((*lepPdgId)[0]+(*lepPdgId)[1])==0 ) {
            Z = *((TLorentzVector*)((*lepP4)[0])) + *((TLorentzVector*)((*lepP4)[1]));
            fakeMET = *((TLorentzVector*)((*metP4)[0])) + Z;
        }
    }

    float MT = 0.0;
    //// W Selection
    if(lepP4->GetEntries() == 1 && n_tightlep == 1) {
        fakeMET = *((TLorentzVector*)((*metP4)[0])) + *((TLorentzVector*)((*lepP4)[0])) ;
        MT = transverseMass( ((TLorentzVector*)((*lepP4)[0]))->Pt(), ((TLorentzVector*)((*lepP4)[0]))->Phi(), ((TLorentzVector*)((*metP4)[0]))->Pt(), ((TLorentzVector*)((*metP4)[0]))->Phi());
    }

    tm->SetValue("mt",MT);

    // ********* Jets ********** //
    for(int jet = 0; jet < jetP4->GetEntries(); jet++) {
        TLorentzVector* Jet = (TLorentzVector*) jetP4->At(jet);
        //cout << (*jetMonojetId)[0] <<endl;
        //cout << Jet->Pt()<<endl;
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      genBr        = eventTree->GetBranch("Gen");
      if( !genBr ){
	printf("MC info is not found in signal MC file\n");
	assert(0);
      }
    }
    eventTree->SetBranchAddress("Electron",   &electronArr);   TBranch *electronBr   = eventTree->GetBranch("Electron");
    eventTree->SetBranchAddress("Dielectron", &dielectronArr); TBranch *dielectronBr = eventTree->GetBranch("Dielectron");
    eventTree->SetBranchAddress("Muon",       &muonArr);       TBranch *muonBr       = eventTree->GetBranch("Muon");
    eventTree->SetBranchAddress("PFJet",      &pfJetArr);      TBranch *pfJetBr      = eventTree->GetBranch("PFJet");
    eventTree->SetBranchAddress("Photon",     &photonArr);     TBranch *photonBr     = eventTree->GetBranch("Photon");
    eventTree->SetBranchAddress("PV",         &pvArr);         TBranch *pvBr         = eventTree->GetBranch("PV");
    
     for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) { 
//      for(UInt_t ientry=0; ientry< 100000; ientry++) { // For testing
      infoBr->GetEntry(ientry);
      if( isGenPresent)
	genBr->GetEntry(ientry);

      electronArr->Clear();   
      dielectronArr->Clear(); 
      muonArr->Clear();       
      pfJetArr->Clear();      
      photonArr->Clear();     
      pvArr->Clear();         
      
      nInputEvts++;
            
      Bool_t keep = kFALSE;

      dielectronBr->GetEntry(ientry);

      // Require at least one dielectron
      if(dielectronArr->GetEntriesFast() > 0) {
	// Apply cuts on the dielectron
	int nDielectronsPass = 0;
	for(Int_t i=0; i<dielectronArr->GetEntriesFast(); i++) {
	  mithep::TDielectron *dielectron = (mithep::TDielectron*)((*dielectronArr)[i]);
	  // Require at least one dielectron above 19 GeV and the other above 9 GeV
	  bool etCut = false;
	  if( (dielectron->scEt_1 > 19 && dielectron->scEt_2 > 9) ||
	      (dielectron->scEt_1 > 9 && dielectron->scEt_2 > 19) )
	    etCut = true;
	  // Require at least one dielectron to pass full ID
	  bool idCut = false;
// 	  if( passSmurf(extractElectron(dielectron,1)) ||
// 	      passSmurf(extractElectron(dielectron,2)) )
	  if( DYTools::energy8TeV == 1 ){
	    if( passEGMID2012( DYTools::extractElectron(dielectron,1), WP_MEDIUM, info->rhoLowEta)
		|| passEGMID2012(DYTools::extractElectron(dielectron,2), WP_MEDIUM, info->rhoLowEta) )
	      idCut = true;
	  }else{
	    if( passEGMID2011( DYTools::extractElectron(dielectron,1), WP_MEDIUM, info->rhoLowEta)
		|| passEGMID2011(DYTools::extractElectron(dielectron,2), WP_MEDIUM, info->rhoLowEta) )
	      idCut = true;
	  }
	  if( etCut && idCut )
	    nDielectronsPass++;
	}
	if(nDielectronsPass > 0)
	  keep = kTRUE;
      }
      
      if(keep) {
	// Some of the objects are dropped for skimmed events
// 	electronBr->GetEntry(ientry);

void renderHLTClusters(TTree* clustersTree)
{
    
    /**************************
     *  Visualization of Clusters
     **************************/
    const Int_t kMaxCl=100*160;
    Int_t fNColorBins = 5;
    
    TEvePointSet* clusters = new TEvePointSet(kMaxCl);
    clusters->SetOwnIds(kTRUE);
    
    TEvePointSetArray * cc = new TEvePointSetArray("TPC Clusters Colorized");
    cc->SetMainColor(kRed);
    cc->SetMarkerStyle(8);
    cc->SetMarkerSize(1.0);
    cc->InitBins("Cluster Charge", fNColorBins, 0., fNColorBins*60.);
    
    cc->GetBin(0)->SetMainColor(kGray);
    cc->GetBin(0)->SetMarkerSize(1.0);
    cc->GetBin(1)->SetMainColor(kBlue);
    cc->GetBin(1)->SetMarkerSize(1.0);
    cc->GetBin(2)->SetMainColor(kCyan);
    cc->GetBin(2)->SetMarkerSize(1.0);
    cc->GetBin(3)->SetMainColor(kGreen);
    cc->GetBin(3)->SetMarkerSize(1.0);
    cc->GetBin(4)->SetMainColor(kYellow);
    cc->GetBin(4)->SetMarkerSize(1.0);
    cc->GetBin(5)->SetMainColor(kRed);
    cc->GetBin(5)->SetMarkerSize(1.0);
    cc->GetBin(6)->SetMainColor(kMagenta);
    cc->GetBin(6)->SetMarkerSize(1.0);
    
    
    // Loop over clusters
    Int_t nentries = clustersTree->GetEntriesFast();
    
    AliTPCClustersRow *clrow = new AliTPCClustersRow();
    clrow->SetClass("AliTPCclusterMI");
    //clrow->SetArray(kMaxCl);
    clustersTree->SetBranchAddress("Segment", &clrow);
    
    for (Int_t i=0; i<nentries; i++) {
        if (!clustersTree->GetEvent(i)) continue;
        
        TClonesArray *cl = clrow->GetArray();
        Int_t ncl = cl->GetEntriesFast();
        
        while (ncl--)
        {
            AliTPCclusterMI* clusterMI = (AliTPCclusterMI*) cl->At(ncl);
            
            AliCluster *c = (AliCluster*) cl->UncheckedAt(ncl);
            Float_t g[3]; //global coordinates
            c->GetGlobalXYZ(g);
            cout<<g[0]<<"\t"<<g[1]<<"\t"<<g[2]<<endl;
            cc->Fill(g[0], g[1], g[2], clusterMI->GetQ());
            clusters->SetNextPoint(g[0], g[1], g[2]);
            AliCluster *atp = new AliCluster(*clusterMI);
            clusters->SetPointId(atp);
            
        }
        cl->Clear();
    }
    
//    delete clrow;
    
    clusters->SetName("TPC Clusters");
    clusters->SetTitle(Form("N=%d", clusters->Size()));
    
//    const TString viz_tag("REC Clusters TPC"); // to be changed
//    clusters->ApplyVizTag(viz_tag, "Clusters");
    
    cc->SetRnrSelf(kTRUE);
    
    gEve->AddElement(cc);
    
    return;
}

void AnalyseEvents(ExRootTreeReader *treeReader, TestPlots *plots)
{
  TClonesArray *branchParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");

  TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack");  // These 3 are used in finding reco jets
  TClonesArray *branchEFlowTower = treeReader->UseBranch("EFlowTower");  // 
  TClonesArray *branchEFlowMuon = treeReader->UseBranch("EFlowMuon");    // 
  TClonesArray *branchJet = treeReader->UseBranch("Jet");

  //  Long64_t allEntries = treeReader->GetEntries();
  Long64_t allEntries = 1000;

  cout << "** Chain contains " << allEntries << " events" << endl;

  GenParticle *particle;
  Electron *electron;
  Photon *photon;
  Muon *muon;

  Track *track;
  Tower *tower;

  Jet *jet;
  TObject *object;

  TLorentzVector momentum;

  Float_t Eem, Ehad;
  Bool_t skip;

  Long64_t entry;

  Int_t i, j, pdgCode;

  // Loop over all events
  for(entry = 0; entry < allEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);

    std::cout << "\n ** NEW ENTRY ** : " << entry << std::endl;

    // Loop over all electrons in event
    for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
    {
      electron = (Electron*) branchElectron->At(i);
      particle = (GenParticle*) electron->Particle.GetObject();

      plots->fElectronDeltaPT->Fill((particle->PT - electron->PT)/particle->PT);
      plots->fElectronDeltaEta->Fill((particle->Eta - electron->Eta)/particle->Eta);
    }

    // Loop over all photons in event
    for(i = 0; i < branchPhoton->GetEntriesFast(); ++i)
    {
      photon = (Photon*) branchPhoton->At(i);

      // skip photons with references to multiple particles
      if(photon->Particles.GetEntriesFast() != 1) continue;

      particle = (GenParticle*) photon->Particles.At(0);

      plots->fPhotonDeltaPT->Fill((particle->PT - photon->PT)/particle->PT);
      plots->fPhotonDeltaEta->Fill((particle->Eta - photon->Eta)/particle->Eta);
      plots->fPhotonDeltaE->Fill((particle->E - photon->E)/particle->E);
    }

    // Loop over all muons in event
    for(i = 0; i < branchMuon->GetEntriesFast(); ++i)
    {
      muon = (Muon*) branchMuon->At(i);
      particle = (GenParticle*) muon->Particle.GetObject();

      plots->fMuonDeltaPT->Fill((particle->PT - muon->PT)/particle->PT);
      plots->fMuonDeltaEta->Fill((particle->Eta - muon->Eta)/particle->Eta);
    }

    // Loop over all tracks in event
    for(i = 0; i < branchEFlowTrack->GetEntriesFast(); ++i)
    {
      track = (Track*) branchEFlowTrack->At(i);
      particle = (GenParticle*) track->Particle.GetObject();

      plots->fTrackDeltaPT->Fill((particle->PT - track->PT)/particle->PT);
      plots->fTrackDeltaEta->Fill((particle->Eta - track->Eta)/particle->Eta);
    }

    // Loop over all towers in event
    for(i = 0; i < branchEFlowTower->GetEntriesFast(); ++i)
    {
      tower = (Tower*) branchEFlowTower->At(i);

      Eem = 0.0;           // Electromagnetic tower energy
      Ehad = 0.0;          // Hadronic tower energy
      skip = kFALSE;
      for(j = 0; j < tower->Particles.GetEntriesFast(); ++j)
      {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                if(hasGen) {
                    genPartArr->Clear();
                    genBr->GetEntry(ientry);
                    genPartBr->GetEntry(ientry);
                    weight*=gen->weight;
                }

                // veto w -> xv decays for signal and w -> ev for bacground samples (needed for inclusive WToLNu sample)
                if (isWrongFlavor && hasGen && fabs(toolbox::flavor(genPartArr, BOSON_ID))==LEPTON_ID) continue;
                else if (isSignal && hasGen && fabs(toolbox::flavor(genPartArr, BOSON_ID))!=LEPTON_ID) continue;

                // check for certified lumi (if applicable)
                baconhep::RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);
                if(hasJSON && !rlrm.hasRunLumi(rl)) continue;

                // trigger requirement
                if (!isEleTrigger(triggerMenu, info->triggerBits, isData)) continue;

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

                //
                // SELECTION PROCEDURE:
                //  (1) Look for 1 good electron matched to trigger
                //  (2) Reject event if another electron is present passing looser cuts
                //
                electronArr->Clear();
                electronBr->GetEntry(ientry);

                Int_t nLooseLep=0;
                const baconhep::TElectron *goodEle=0;
                Bool_t passSel=kFALSE;

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

                    // check ECAL gap
                    if(fabs(ele->scEta)>=ECAL_GAP_LOW && fabs(ele->scEta)<=ECAL_GAP_HIGH) continue;

                    // apply scale and resolution corrections to MC
                    Double_t elescEt_corr = ele->scEt;
                    if(doScaleCorr && snamev[isam].CompareTo("data",TString::kIgnoreCase)!=0)
                        elescEt_corr = gRandom->Gaus(ele->scEt*getEleScaleCorr(ele->scEta,0),getEleResCorr(ele->scEta,0));

                    if(fabs(ele->scEta)   > VETO_ETA) continue;        // loose lepton |eta| cut
                    if(elescEt_corr       < VETO_PT)  continue;        // loose lepton pT cut
                    if(passEleLooseID(ele,info->rhoIso)) nLooseLep++;  // loose lepton selection
                    if(nLooseLep>1) {  // extra lepton veto
                        passSel=kFALSE;
                        break;
                    }

                    if(fabs(ele->scEta)   > ETA_CUT)     continue;  // lepton |eta| cut
                    if(elescEt_corr       < PT_CUT)      continue;  // lepton pT cut
                    if(!passEleID(ele,info->rhoIso))     continue;  // lepton selection
                    if(!isEleTriggerObj(triggerMenu, ele->hltMatchBits, kFALSE, isData)) continue;

                    passSel=kTRUE;
                    goodEle = ele;
                }

                if(passSel) {

                    //******* We have a W candidate! HURRAY! ********
                    nsel+=weight;
                    nselvar+=weight*weight;

void HHToBBGGSelectionMistags(const string inputfile,          // input file
                              const string outputfile,         // output directory
                              Int_t SampleType = 7,
                              Bool_t applyExtrapWeightTo140PU = kFALSE
  ) {
  
  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //============================================================================================================== 
  bool printdebug = false;

  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================  
  Double_t nEvents = 0;
  Double_t nEventsTwoRealPho = 0;
  
  //*****************************************************************************************
  // Set up output ntuple
  //*****************************************************************************************
  TFile *outFile = new TFile(outputfile.c_str(),"RECREATE"); 
  TH1F *NEvents =  new TH1F("NEvents",";;",1,-0.5,0.5);
  //edit
  TH1F *NEventsTwoRealPho =  new TH1F("NEventsTwoRealPho",";;",1,-0.5,0.5);
  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");

  TFile *BJetMistagRateLightJetsFile = TFile::Open("/afs/cern.ch/work/s/sixie/public/releases/analysis/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/BTaggingEfficiency_LightJetsMistagRate.root");
  TFile *BJetMistagRateCharmJetsFile = TFile::Open("/afs/cern.ch/work/s/sixie/public/releases/analysis/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/BTaggingEfficiency_CharmJetsMistagRate.root");
  TFile *PhotonEfficiencyFile = TFile::Open("/afs/cern.ch/work/s/sixie/public/releases/analysis/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/PhotonEfficiency_PromptPhoton.root");
  TH2F *BJetMistagRateLightJetsHist = (TH2F*)BJetMistagRateLightJetsFile->Get("MistagRate_CSVMedium_Pt_Eta");
  TH2F *BJetMistagRateCharmJetsHist = (TH2F*)BJetMistagRateCharmJetsFile->Get("MistagRate_CSVMedium_Pt_Eta");
  TH2F *PhotonEfficiencyHist = (TH2F*)PhotonEfficiencyFile->Get("Efficiency_PtEta");
  TH2F *bjet1FakeRateHist = 0;
  TH2F *bjet2FakeRateHist = 0;
  double bjet1Eff = 0;
  double bjet2Eff = 0;

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

  
  Double_t weight = 1;

  // 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);
    NEventsTwoRealPho->Fill(0);
    NPUMean->Fill(info->nPUMean);

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

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

    genparticleBr->GetEntry(ientry);
    genjetBr->GetEntry(ientry);

    double tmpHT = 0;

    //***********************************************************
    // Find Gen-Level particles
    //***********************************************************
    const cmsana::TGenParticle *genPhoton1 = 0;
    const cmsana::TGenParticle *genPhoton2 = 0;
    for(Int_t i=0; i<genparticleArr->GetEntriesFast(); i++) {
      const cmsana::TGenParticle *p = (cmsana::TGenParticle*)((*genparticleArr)[i]);

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

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

      if (!passCutBasedEleID(tag, ComputeElePFIso04(tag,pfcandidateArr,rhoEleIso,EleEAEra))) continue;

      if(dataType == 0 &&
         !((info->triggerBits & kHLT_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC4_Mass50) && (tag->hltMatchBits & kHLTObject_Ele20_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT)) &&
         !((info->triggerBits & kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30) && (tag->hltMatchBits & kHLTObject_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT)) &&
         !((info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) && (tag->hltMatchBits & kHLTObject_Ele32_CaloIdL_CaloIsoVL)) &&
         !((info->triggerBits & kHLT_Ele32_CaloIdL_CaloIsoVL_SC17) && (tag->hltMatchBits & kHLTObject_Ele32_CaloIdT_CaloIsoT_TrkIdT_TrkIsoT))         
        )
        continue;
      
      if(dataType == 1 &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT)) &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele32_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT)) &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele52_CaloIdVT_TrkIdT)) &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele65_CaloIdVT_TrkIdT)) &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele80_CaloIdVT_TrkIdT)) &&
         !((info->triggerBits & kHLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) && (tag->hltMatchBits & kHLTObject_Ele27_WP80)) &&

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

    genparticleBr->GetEntry(ientry);
    genjetBr->GetEntry(ientry);

    //***********************************************************
    // Find Gen-Level particles
    //***********************************************************
    const cmsana::TGenParticle *genPhoton1 = 0;
    const cmsana::TGenParticle *genPhoton2 = 0;
    const cmsana::TGenParticle* photon1 = 0;
    const cmsana::TGenParticle *genB1 = 0;
    const cmsana::TGenParticle *genB2 = 0;
    const cmsana::TGenJet* genBJet1 = 0;
    const cmsana::TGenJet* genBJet2 = 0;
    const cmsana::TGenJet* bjet1 = 0;
    const cmsana::TGenJet* bjet2 = 0;
    Float_t FakeRate3 = 0;
    Float_t FakeRate4 = 0;
    for(Int_t i=0; i<genparticleArr->GetEntriesFast(); i++) {
      const cmsana::TGenParticle *p = (cmsana::TGenParticle*)((*genparticleArr)[i]);
//       cout << p->pdgid << " " << p->status << " " << p->pt << " " << p->eta << " " << p->phi 
//            << " | " << p->motherPdgID << "\n";

      if ( SampleType == cmsana::HHToBBGGEventTree::HHToBBGG) {
        if (abs(p->pdgid) == 5 && p->motherPdgID == 25 && p->status == 2) {
          if (!genB1) {
            genB1 = p;
          } else {
            if (!genB2) genB2 = p;
          }
        }
      }
      if ( SampleType == cmsana::HHToBBGGEventTree::ttHgg 
           || SampleType == cmsana::HHToBBGGEventTree::ttbar
        ) {
        if (abs(p->pdgid) == 5 && abs(p->motherPdgID) == 6 && p->status == 2) {
          if (!genB1) {
            genB1 = p;
          } else {
            if (!genB2) genB2 = p;
          }
        }
      }
      if ( SampleType == cmsana::HHToBBGGEventTree::ZHgg 
        ) {
        if (abs(p->pdgid) == 5 && p->motherPdgID == 23 && p->status == 2) {
          if (!genB1) {
            genB1 = p;
          } else {
            if (!genB2) genB2 = p;
          }