C++ Event::getByLabel方法代码示例

int NeroAll::analyze(const edm::Event&iEvent)
{
    if(isMc_ <0 )
    {
        isMc_ = ( iEvent.isRealData() ) ? 0 : 1;
        isRealData = ( iEvent.isRealData() ) ? 1 : 0;
    }
    if( isSkim() == 0) 
    {
        //TODO FILL all_
        isRealData = ( iEvent.isRealData() ) ? 1 : 0;	
        runNum = iEvent.id().run();
        lumiNum = iEvent.luminosityBlock(); 	
        eventNum = iEvent.id().event();  
        iEvent.getByLabel(edm::InputTag("generator"), info_handle); // USE TOKEN AND INPUT TAG ?
        iEvent.getByLabel(edm::InputTag("addPileupInfo"), pu_handle);
        mcWeight = info_handle->weight();
        //PU
        puTrueInt = 0;
        for(const auto & pu : *pu_handle)
        {
            //Intime
            if (pu.getBunchCrossing() == 0)
                puTrueInt += pu.getTrueNumInteractions();
            //Out-of-time
        }
    }
    return 0;
}

void UAHiggsTree::GetGenMET(const edm::Event& iEvent , const edm::EventSetup& iSetup,
                            const string GenMETCollection_ , vector<MyGenMET>& METVector )
{
   using namespace std;
   using namespace edm;
   using namespace reco;

   METVector.clear();
  
   Handle<GenMETCollection> met;
   iEvent.getByLabel(GenMETCollection_, met);

   for(reco::GenMETCollection::const_iterator met_iter=met->begin(); 
      met_iter != met->end(); met_iter++)
    {
      
      MyGenMET mymet;
     
      mymet.pt=met_iter->pt();
      mymet.phi=met_iter->phi();
      mymet.eta=met_iter->eta();
      
      
      METVector.push_back(mymet);
         

    } 

}

void ChPartTree::GetL1Trig(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace std;

  // Tests
/*
  if (L1TrigDebug)
  {  
   edm::ESHandle<L1GtTriggerMenu> menuRcd;
   iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
   const L1GtTriggerMenu* menu = menuRcd.product();

   for (CItAlgo algo = menu->gtAlgorithmMap().begin(); 
                algo!=menu->gtAlgorithmMap().end(); 
                ++algo) 
   {
     cout << "Name: " << (algo->second).algoName() 
          << " Alias: " << (algo->second).algoAlias() 
          << " Bit: " << (algo->second).algoBitNumber()
           << " Result: " << l1AlgorithmResult(iEvent, iSetup, (algo->second).algoName() )
          << std::endl;
   }
  }
*/

   edm::Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
   iEvent.getByLabel("gtDigis",L1GTRR);
   for (int i=0 ; i <128 ; i++) 
   {
     if (L1TrigDebug) cout << "PhysicsTriggerWord :" << i << " " << L1GTRR->decisionWord()[i] << endl;
     L1Trig.PhysTrigWord[i] = L1GTRR->decisionWord()[i];
   }
   for (int i=0 ; i <64  ; i++)
   {
    if (L1TrigDebug) cout << "technicalTriggerWord :" << i << " " << L1GTRR->technicalTriggerWord()[i] << endl;
     L1Trig.TechTrigWord[i] = L1GTRR->technicalTriggerWord()[i]; 
   }
}

void UAHiggsTree::GetGenPart(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace std;
   using namespace edm;
   using namespace reco;

   // Clear GenPart List
   // (evt->GenPart).clear();
   GenPart.clear();
   GenElec.clear();
   GenMu.clear();
   GenNu.clear();

   // Handle to access PDG data from iSetup
   ESHandle <ParticleDataTable> pdt;
   iSetup.getData( pdt );

   // Handle to access GenParticleCollection
   Handle<GenParticleCollection> genParticles;
   iEvent.getByLabel(genPartColl_, genParticles);
   if (!genParticles.isValid()) {
     cerr << "[UAHiggsTree::GetGenPart] Error: non valid GenParticleCollection " << endl;
     return;
   }

   // List for Daugther/Mother Search
   vector<const reco::Candidate *> cands;
   vector<const Candidate *>::const_iterator found = cands.begin();
   for(GenParticleCollection::const_iterator p = genParticles->begin();
     p != genParticles->end(); ++ p) {
     cands.push_back(&*p);
   }

   // Loop on generated particle and store if status=3
   if ( GenPartDebug )
     cout << "GenPart # : " << genParticles->size() << endl;
   
   for(GenParticleCollection::const_iterator p = genParticles->begin();
     p != genParticles->end(); ++ p) {
     int st = p->status();

       MyGenPart genpart;

      if ( GenPartDebug )
      {
       if(( fabs(p->pdgId())==11 || fabs(p->pdgId())== 13) && st==1) 
         cout    << p-genParticles->begin()
                 << " " << st  
                 << " " << p->pdgId()
                 << " " << (pdt->particle(p->pdgId()))->name()
                 << " " << p->pt()
                 << " " << p->eta()
                 << " " << p->phi()
                 << " " << p->charge()
                 << endl ;

       if(( fabs(p->pdgId())==11 || fabs(p->pdgId())== 13) && st==3)
         cout    << p-genParticles->begin()
                 << " " << st
                 << " " << p->pdgId()
                 << " " << (pdt->particle(p->pdgId()))->name()
                 << " " << p->pt()
                 << " " << p->eta()
                 << " " << p->phi()
                 << " " << p->charge()
                 << endl ;
       }

       // Four vector

       genpart.pt   = p->pt();
       genpart.eta  = p->eta();
       genpart.phi  = p->phi();
       genpart.e    = p->energy();
       genpart.px   = p->px();
       genpart.py   = p->py();
       genpart.pz   = p->pz();
       genpart.m    = p->mass();

       //genpart.v.SetPxPyPzE( p->px() , p->py() , p->pz() , p->energy() );
       genpart.Part.v.SetPxPyPzE( p->px() , p->py() , p->pz() , p->energy() );

       // Extra properties

       genpart.Part.charge  = p->charge();
       genpart.charge  = p->charge();
       genpart.pdgId   = p->pdgId();
       genpart.name    = (pdt->particle(p->pdgId()))->name();
       genpart.status  = p->status();

       // Mother Daughter relations

       int nMo = p->numberOfMothers();
       int nDa = p->numberOfDaughters();
       int iMo1 = -1 , iMo2 = -1 ; 
       int iDa1 = -1 , iDa2 = -1 ; 

       found = find(cands.begin(), cands.end(), p->mother(0));
       if(found != cands.end()) iMo1 = found - cands.begin() ;
 
//.........这里部分代码省略.........

// ------------ method called for each event  ------------
void
ThreePointCorrelatorEtaTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;
  using namespace std;

  edm::Handle<reco::VertexCollection> vertices;
  iEvent.getByLabel(vertexSrc_,vertices);
  double bestvz=-999.9, bestvx=-999.9, bestvy=-999.9;
  double bestvzError=-999.9, bestvxError=-999.9, bestvyError=-999.9;
  const reco::Vertex & vtx = (*vertices)[0];
  bestvz = vtx.z(); bestvx = vtx.x(); bestvy = vtx.y();
  bestvzError = vtx.zError(); bestvxError = vtx.xError(); bestvyError = vtx.yError();
  
  //first selection; vertices
  if(bestvz < vzLow_ || bestvz > vzHigh_ ) return;
  
  Handle<CaloTowerCollection> towers;
  iEvent.getByLabel(towerSrc_, towers);

  Handle<reco::TrackCollection> tracks;
  iEvent.getByLabel(trackSrc_, tracks);
  
  if( useCentrality_ ){

    CentralityProvider * centProvider = 0;
    if (!centProvider) centProvider = new CentralityProvider(iSetup);
    centProvider->newEvent(iEvent,iSetup);
    int hiBin = centProvider->getBin();
    cbinHist->Fill( hiBin );
    if( hiBin < Nmin_ || hiBin >= Nmax_ ) return;

  }

  const int NetaBins = etaBins_.size() - 1 ;

// initialize Qcos and Qsin

  double QcosTRK = 0.;
  double QsinTRK = 0.;
  double QcountsTrk = 0;

  double Q1[NetaBins][2][2];
  double Q1_count[NetaBins][2];

  double Q2[NetaBins][2][2];
  double Q2_count[NetaBins][2];

  for(int i = 0; i < NetaBins; i++){
    for(int j = 0; j < 2; j++){
      Q1_count[i][j] = 0.0;
      Q2_count[i][j] = 0.0;
      for(int k = 0; k < 2; k++){
        Q1[i][j][k] = 0.0;
        Q2[i][j][k] = 0.0;
      }
    }
  }

  int nTracks = 0;
  for(unsigned it = 0; it < tracks->size(); it++){

     const reco::Track & trk = (*tracks)[it];

     math::XYZPoint bestvtx(bestvx,bestvy,bestvz);
        
        double dzvtx = trk.dz(bestvtx);
        double dxyvtx = trk.dxy(bestvtx);
        double dzerror = sqrt(trk.dzError()*trk.dzError()+bestvzError*bestvzError);
        double dxyerror = sqrt(trk.d0Error()*trk.d0Error()+bestvxError*bestvyError);
        double nhits = trk.numberOfValidHits();
        double chi2n = trk.normalizedChi2();
        double nlayers = trk.hitPattern().trackerLayersWithMeasurement();
        chi2n = chi2n/nlayers;

        double weight = 1.0;

        if(!trk.quality(reco::TrackBase::highPurity)) continue;
        if(fabs(trk.ptError())/trk.pt() > offlineptErr_ ) continue;
        if(fabs(dzvtx/dzerror) > offlineDCA_) continue;
        if(fabs(dxyvtx/dxyerror) > offlineDCA_) continue;
        if(fabs(trk.eta()) < 2.4 && trk.pt() > 0.4 ){nTracks++;}// NtrkOffline        
        if(fabs(trk.eta()) > 2.4 || trk.pt() < ptLow_ || trk.pt() > ptHigh_) continue;
        if(chi2n > offlineChi2_) continue;
        if(nhits < offlinenhits_) continue;
        if( messAcceptance_ ) { if( trk.phi() < holeRight_ && trk.phi() > holeLeft_ ) continue;}
        if( doEffCorrection_ ){ weight = 1.0/effTable->GetBinContent( effTable->FindBin(trk.eta(), trk.pt()) );}
        else{ weight = 1.0; }
       
        trkPhi->Fill( trk.phi() );//make sure if messAcceptance is on or off

        QcosTRK += weight*cos( 2*trk.phi() );
        QsinTRK += weight*sin( 2*trk.phi() );
        QcountsTrk += weight;

        for(int eta = 0; eta < NetaBins; eta++){
          if( trk.eta() > etaBins_[eta] && trk.eta() < etaBins_[eta+1] ){

            if( trk.charge() == 1){
//.........这里部分代码省略.........

// ------------ method called to skim the data  ------------
bool MCBarrelEndcapFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
    using namespace edm;
    bool accepted = false;
    Handle<HepMCProduct> evt;
    iEvent.getByLabel(label_, evt);

    vector<HepMC::GenParticle*> barrel_passed;
    vector<HepMC::GenParticle*> endcap_passed;


    const HepMC::GenEvent * myGenEvent = evt->GetEvent();

    if(verbose) cout << "Start Event\n";

    for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
            p != myGenEvent->particles_end(); ++p ) {


        if(verbose) cout << abs((*p)->pdg_id()) << ' ' 
            << (*p)->momentum().perp() << ' ' 
            << (*p)->momentum().eta() << ' ' 
            << (*p)->status() << endl;
        // check for barrel conditions
        if ((abs((*p)->pdg_id()) == abs(barrelID) || barrelID == 0) &&
                (*p)->momentum().perp() > ptMinBarrel && (*p)->momentum().eta() > etaMinBarrel 
                && (*p)->momentum().eta() < etaMaxBarrel && ((*p)->status() == statusBarrel || statusBarrel == 0))
        { 
            if(verbose) cout << "Found Barrel\n";
            // passed barrel conditions ...
            // ... now check pair-conditions with endcap passed particles
            //
            unsigned int i=0;
            double invmass =0.;
            while(!accepted && i<endcap_passed.size()) {
                invmass = ( math::PtEtaPhiMLorentzVector((*p)->momentum()) + math::PtEtaPhiMLorentzVector(endcap_passed[i]->momentum())).mass();
                if(verbose) cout << "Invariant Mass is" << invmass << endl;
                if(invmass > minInvMass && invmass < maxInvMass) {
                    accepted = true;
                }
                i++;
            }    
            // if we found a matching pair quit the loop
            if(accepted) break;
            else{
                barrel_passed.push_back(*p);   // else remember the particle to have passed barrel conditions
            }
        }

        // check for endcap conditions

        if ((abs((*p)->pdg_id()) == abs(endcapID) || endcapID == 0) && 
                (*p)->momentum().perp() > ptMinEndcap && (*p)->momentum().eta() > etaMinEndcap
                && (*p)->momentum().eta() < etaMaxEndcap && ((*p)->status() == statusEndcap || statusEndcap == 0)) { 
            if(verbose) cout << "Found Endcap\n";
            // passed endcap conditions ...
            // ... now check pair-conditions with barrel passed particles vector
            unsigned int i=0;
            double invmass =0.;
            while(!accepted && i<barrel_passed.size()) {
                if((*p) != barrel_passed[i]) {
                    if(verbose) cout << "Checking to match Barrel\n";
                    invmass = ( math::PtEtaPhiMLorentzVector((*p)->momentum()) + math::PtEtaPhiMLorentzVector(barrel_passed[i]->momentum())).mass();
                    if(verbose) cout << "Invariant Mass is" << invmass << endl;
                    if(invmass > minInvMass && invmass < maxInvMass) {
                        accepted = true;
                    }
                    i++;
                }
            }    
            // if we found a matching pair quit the loop
            if(accepted) break;
            else {
                endcap_passed.push_back(*p);   // else remember the particle to have passed endcap conditions
            }
        }
    }

    if(verbose)
    {
        cout << "This event ";
        if (accepted){ cout <<  "passed "; } else {cout << "failed ";}
        cout << "the filter\n";
    }

    return accepted;

}

void UAHiggsTree::GetGenKin(const edm::Event& iEvent)
{
   using namespace std;
   using namespace edm;
   using namespace reco;

/* 
   // MC Process Id and PtHat for AOD

   Handle< int > genProcessID;
   iEvent.getByLabel( "genEventProcID", genProcessID );
   int processId = *genProcessID;
   cout<<"Process ID : "<<processId<<endl;

*/

   //K Factor For Signal
   edm::Handle<double> KFactor;
   try {
    iEvent.getByLabel("KFactorProducer",KFactor);
    double kfactor = *KFactor;
    GenKin.kfactor=kfactor;
   }
   catch (...){;}

   /*
   // MC Process Id and PtHat for RECO
   Handle<HepMCProduct> hepMCHandle;
   iEvent.getByLabel(hepMCColl_, hepMCHandle ) ;
   const HepMC::GenEvent* GenEvt = hepMCHandle->GetEvent() ;
   int processId = GenEvt->signal_process_id();
   if (GenKinDebug) cout<<"Process ID2: "<<processId<<endl; 
   double ptHat  = GenEvt->event_scale();
   if (GenKinDebug) cout<<"PtHat MC : "<<ptHat<<endl; 

   GenKin.MCProcId  = processId ; 
   GenKin.Scale     = ptHat ;

   // PDF Info -> x,Q, parton id's -> see HepMC manual
   const HepMC::PdfInfo* PdfInfo = GenEvt->pdf_info();
   double x1 = PdfInfo->x1();
   double x2 = PdfInfo->x2();
   double Q  = PdfInfo->scalePDF(); 

   int id1 = PdfInfo->id1();
   int id2 = PdfInfo->id2();

   GenKin.x1 = x1 ;
   GenKin.x2 = x2 ;
   GenKin.Q  = Q  ;
   GenKin.Part1Id = id1 ; 
   GenKin.Part2Id = id2 ; */

   // Gen MET: From GenMETCollection  

   Handle<GenMETCollection> genmet;
//   iEvent.getByLabel(GenMetColl_ , genmet);
   iEvent.getByLabel("genMetTrue" , genmet);

   if (GenKinDebug) 
   { 
     cout << "GenMET et()  = " << (genmet->front()).et() << endl ;
     cout << "GenMET px()  = " << (genmet->front()).px() << endl ;
     cout << "GenMET py()  = " << (genmet->front()).py() << endl ;
     cout << "GenMET phi() = " << (genmet->front()).phi() << endl;
   }

   GenKin.Met    = (genmet->front()).et() ; 
   GenKin.MetX   = (genmet->front()).px() ; 
   GenKin.MetY   = (genmet->front()).py() ; 
   GenKin.MetPhi = (genmet->front()).phi() ; 

   // Gen MET: From GenPart Neutrinos ( no SUSY !!! )

   // ... Handle to access GenParticleCollection
   Handle<GenParticleCollection> genParticles;
   iEvent.getByLabel(genPartColl_, genParticles);
   if (!genParticles.isValid()) {
     cerr << "[UAHiggsTree::GetGenKin] Error: non valid GenParticleCollection " << endl;
     return;
   }
 
   double met1Px = 0 ;
   double met1Py = 0 ;
   double met1Pz = 0 ;
   double met1E  = 0 ;

   double met3Px = 0 ;
   double met3Py = 0 ;
   double met3Pz = 0 ;
   double met3E  = 0 ;

   // ... Loop on stable (final) particles
   for(GenParticleCollection::const_iterator p = genParticles->begin(); p != genParticles->end(); ++ p) 
   {
     int st  = p->status();
     int pid = p->pdgId();
     if ( st == 1 && ( abs(pid)==12 || abs(pid)==14 || abs(pid)==16 ) )
     {
       met1Px += p->px();
//.........这里部分代码省略.........

void UAHiggsTree::GetRecoMuon(const edm::Event& iEvent, const edm::EventSetup& iSetup,
                                                       const string globalMuonCollection_, vector<MyMuon>& MuonVector )
{
   using namespace std;
   using namespace edm;
   using namespace reco;

   // BORIS 3D IP stuff: Declaration
   ESHandle<TransientTrackBuilder> theTTBuilder;
   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theTTBuilder);


   MuonVector.clear();

   Handle<MuonCollection> muonsHandle;
   try {
     iEvent.getByLabel(globalMuonCollection_,muonsHandle);
   } catch ( cms::Exception& ex ) {
    printf("Error! can't get globalMuon collection\n");
   }
   const MuonCollection & muons = *(muonsHandle.product());

   for (reco::MuonCollection::const_iterator iMuon = muons.begin(); 
        iMuon != muons.end(); iMuon++) {

     MyMuon muon;     
     
     // Global Muon
   
     muon.pt  = iMuon->pt(); 
     muon.px  = iMuon->px(); 
     muon.py  = iMuon->py(); 
     muon.pz  = iMuon->pz(); 
     muon.e   = sqrt(iMuon->momentum().mag2()+MASS_MU*MASS_MU); 
     muon.eta = iMuon->eta(); 
     muon.phi = iMuon->phi(); 
     
     
     muon.Part.v.SetPxPyPzE(iMuon->px(),
                            iMuon->py(),
                            iMuon->pz(), 
                            sqrt(iMuon->momentum().mag2()+MASS_MU*MASS_MU));
     muon.Part.charge = iMuon->charge();

     // Global Properties

     muon.nChambers             = iMuon -> numberOfChambers();
     muon.nChambersMatched      = iMuon -> numberOfMatches();
     
     muon.AllGlobalMuons                          =  muon::isGoodMuon(*iMuon,muon::AllGlobalMuons);
     muon.AllStandAloneMuons                      =  muon::isGoodMuon(*iMuon,muon::AllStandAloneMuons);
     muon.AllTrackerMuons                         =  muon::isGoodMuon(*iMuon,muon::AllTrackerMuons);
     muon.TrackerMuonArbitrated                   =  muon::isGoodMuon(*iMuon,muon::TrackerMuonArbitrated);
     muon.AllArbitrated                           =  muon::isGoodMuon(*iMuon,muon::AllArbitrated);
     muon.GlobalMuonPromptTight                   =  muon::isGoodMuon(*iMuon,muon::GlobalMuonPromptTight);
     muon.TMLastStationLoose                      =  muon::isGoodMuon(*iMuon,muon::TMLastStationLoose);
     muon.TMLastStationTight                      =  muon::isGoodMuon(*iMuon,muon::TMLastStationTight);
     muon.TM2DCompatibilityLoose                  =  muon::isGoodMuon(*iMuon,muon::TM2DCompatibilityLoose);
     muon.TM2DCompatibilityTight                  =  muon::isGoodMuon(*iMuon,muon::TM2DCompatibilityTight);
     muon.TMOneStationLoose                       =  muon::isGoodMuon(*iMuon,muon::TMOneStationLoose);
     muon.TMOneStationTight                       =  muon::isGoodMuon(*iMuon,muon::TMOneStationTight);
     muon.TMLastStationOptimizedLowPtLoose        =  muon::isGoodMuon(*iMuon,muon::TMLastStationOptimizedLowPtLoose);
     muon.TMLastStationOptimizedLowPtTight        =  muon::isGoodMuon(*iMuon,muon::TMLastStationOptimizedLowPtTight);
     muon.GMTkChiCompatibility                    =  muon::isGoodMuon(*iMuon,muon::GMTkChiCompatibility);
     muon.GMStaChiCompatibility                   =  muon::isGoodMuon(*iMuon,muon::GMStaChiCompatibility);
     muon.GMTkKinkTight                           =  muon::isGoodMuon(*iMuon,muon::GMTkKinkTight);
     muon.TMLastStationAngLoose                   =  muon::isGoodMuon(*iMuon,muon::TMLastStationAngLoose);
     muon.TMLastStationAngTight                   =  muon::isGoodMuon(*iMuon,muon::TMLastStationAngTight);
     muon.TMOneStationAngLoose                    =  muon::isGoodMuon(*iMuon,muon::TMOneStationAngLoose);
     muon.TMOneStationAngTight                    =  muon::isGoodMuon(*iMuon,muon::TMOneStationAngTight);
     muon.TMLastStationOptimizedBarrelLowPtLoose  =  muon::isGoodMuon(*iMuon,muon::TMLastStationOptimizedBarrelLowPtLoose);
     muon.TMLastStationOptimizedBarrelLowPtTight  =  muon::isGoodMuon(*iMuon,muon::TMLastStationOptimizedBarrelLowPtTight);
     
 
     muon.isoR03sumPt   = iMuon->isolationR03().sumPt   ;
     muon.isoR03emEt    = iMuon->isolationR03().emEt    ;
     muon.isoR03hadEt   = iMuon->isolationR03().hadEt   ;
     muon.isoR03hoEt    = iMuon->isolationR03().hoEt    ;
     muon.isoR03nTracks = iMuon->isolationR03().nTracks ;
     muon.isoR03nJets   = iMuon->isolationR03().nJets   ;

     muon.isoR05sumPt   = iMuon->isolationR05().sumPt   ;
     muon.isoR05emEt    = iMuon->isolationR05().emEt    ;
     muon.isoR05hadEt   = iMuon->isolationR05().hadEt   ;
     muon.isoR05hoEt    = iMuon->isolationR05().hoEt    ;
     muon.isoR05nTracks = iMuon->isolationR05().nTracks ;
     muon.isoR05nJets   = iMuon->isolationR05().nJets   ;

     muon.calEnergyEm   = iMuon->calEnergy().em   ; 
     muon.calEnergyHad  = iMuon->calEnergy().had  ;
     muon.calEnergyHo   = iMuon->calEnergy().ho   ;
     muon.calEnergyEmS9 = iMuon->calEnergy().emS9 ;
     muon.calEnergyHadS9= iMuon->calEnergy().hadS9;
     muon.calEnergyHoS9 = iMuon->calEnergy().hoS9 ;

     muon.IsGlobalMuon      = iMuon->isGlobalMuon()     ;
     muon.IsTrackerMuon     = iMuon->isTrackerMuon()    ;
     muon.IsStandaloneMuon  = iMuon->isStandAloneMuon() ;
     muon.IsCaloMuon        = iMuon->isCaloMuon()       ;

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

void ChPartTree::GetFwdGap(const edm::Event& iEvent , const edm::EventSetup& iSetup )
{
   using namespace std;


  int nTowersHF_plus = 0;
  int nTowersHF_minus = 0;
  int nTowersHE_plus = 0;
  int nTowersHE_minus = 0;
  int nTowersHB_plus = 0;
  int nTowersHB_minus = 0;
  int nTowersEE_plus = 0;
  int nTowersEE_minus = 0;
  int nTowersEB_plus = 0;
  int nTowersEB_minus = 0; 
  //Sum(E)
  double sumEHF_plus = 0.;
  double sumEHF_minus = 0.;
  double sumEHE_plus = 0.;
  double sumEHE_minus = 0.;
  double sumEHB_plus = 0.;
  double sumEHB_minus = 0.;
  double sumEEE_plus = 0.;
  double sumEEE_minus = 0.;
  double sumEEB_plus = 0.;
  double sumEEB_minus = 0.;
  // Sum(ET)
  double sumETHF_plus = 0.;
  double sumETHF_minus = 0.;
  double sumETHE_plus = 0.;
  double sumETHE_minus = 0.;
  double sumETHB_plus = 0.;
  double sumETHB_minus = 0.;
  double sumETEE_plus = 0.;
  double sumETEE_minus = 0.;
  double sumETEB_plus = 0.;
  double sumETEB_minus = 0.;

  FwdGap.Reset();

  // Calo tower collection from event
  edm::Handle<CaloTowerCollection> towerCollectionH;
  iEvent.getByLabel(caloTowerTag_,towerCollectionH);
  const CaloTowerCollection& towerCollection = *towerCollectionH;

  // Loop over calo towers
  CaloTowerCollection::const_iterator calotower = towerCollection.begin();
  CaloTowerCollection::const_iterator calotowers_end = towerCollection.end();
  for(; calotower != calotowers_end; ++calotower) {
     bool hasHCAL = false;
     bool hasHF = false;
     bool hasHE = false;
     bool hasHB = false;
     bool hasHO = false;
     bool hasECAL = false;
     bool hasEE = false;
     bool hasEB = false;     
     for(size_t iconst = 0; iconst < calotower->constituentsSize(); iconst++){
        DetId detId = calotower->constituent(iconst);
        if(detId.det()==DetId::Hcal){
           hasHCAL = true;
           HcalDetId hcalDetId(detId);
           if(hcalDetId.subdet()==HcalForward) hasHF = true;
           else if(hcalDetId.subdet()==HcalEndcap) hasHE = true;
           else if(hcalDetId.subdet()==HcalBarrel) hasHB = true;
           else if(hcalDetId.subdet()==HcalOuter) hasHO = true;  
        } else if(detId.det()==DetId::Ecal){
           hasECAL = true;
           EcalSubdetector ecalSubDet = (EcalSubdetector)detId.subdetId();
           if(ecalSubDet == EcalEndcap) hasEE = true;
           else if(ecalSubDet == EcalBarrel) hasEB = true;
        }
     }

     int zside = calotower->zside();
     double caloTowerEnergy = calotower->energy();
     // FIXME
     //double caloTowerET = calotower->et(primVtx.position());
     //double caloTowerET = calotower->et(primVtx.z());
     double caloTowerET = calotower->et();

     // HCAL: Towers made of at least one component from HB,HE,HF
     if( hasHF && !hasHE ){
        if( caloTowerEnergy >= energyThresholdHF_ ){
           if(zside >= 0){
              ++nTowersHF_plus;
              sumEHF_plus += caloTowerEnergy; 
              sumETHF_plus += caloTowerET;
           } else{
              ++nTowersHF_minus;
              sumEHF_minus += caloTowerEnergy;
              sumETHF_minus += caloTowerET;
           } 
        }
     } else if( hasHE && !hasHF && !hasHB ){
        if( caloTowerEnergy >= energyThresholdHE_ ){
           if(zside >= 0){
              ++nTowersHE_plus;
              sumEHE_plus += caloTowerEnergy;
              sumETHE_plus += caloTowerET;
//.........这里部分代码省略.........

void UAHiggsTree::GetGenJet(const edm::Event& iEvent , const edm::EventSetup& iSetup,
                            const string GenJetCollection_ , vector<MyGenJet>& JetVector )
{
   using namespace std;
   using namespace edm;
   using namespace reco;

   // Clear

   JetVector.clear();
   

   // Handle to access PDG data from iSetup
   ESHandle <ParticleDataTable> pdt;
   iSetup.getData( pdt );


   // get gen jet collection
   Handle<GenJetCollection> genjets;
   iEvent.getByLabel(GenJetCollection_, genjets);
   
   for(GenJetCollection::const_iterator genjet=genjets->begin();genjet!=genjets->end();genjet++){ 
     
      MyGenJet Jet;
      
      Jet.et  = genjet->et();
      Jet.pt  = genjet->pt();
      Jet.eta = genjet->eta();
      Jet.phi = genjet->phi();
      Jet.e   = genjet->energy();
      Jet.px  = genjet->px();
      Jet.py  = genjet->py();
      Jet.pz  = genjet->pz();
      Jet.npart = genjet->nConstituents();

      Jet.JetPart.clear();

      if (GenJetDebug) cout << genjet->nConstituents() << endl;    
      std::vector <const reco::GenParticle*> mcparts = genjet->getGenConstituents();
      for (unsigned i = 0; i < mcparts.size (); i++) {
        const reco::GenParticle* p = mcparts[i]; 
        if (GenJetDebug) 
          cout    << i
                  << " " <<  p->pdgId()
                  << " " << (pdt->particle(p->pdgId()))->name()
                  << " " << p->pt()
                  << " " << p->eta()
                  << " " << p->phi()
                  << " " << p->charge()
                  << endl ;

         MyGenPart genpart;

         // Four vector

         genpart.pt   = p->pt();
         genpart.eta  = p->eta();
         genpart.phi  = p->phi();
         genpart.e    = p->energy();
         genpart.px   = p->px();
         genpart.py   = p->py();
         genpart.pz   = p->pz();
         genpart.m    = p->mass();
         genpart.Part.v.SetPxPyPzE( p->px() , p->py() , p->pz() , p->energy() );

         // Extra properties

         genpart.Part.charge  = p->charge();
         genpart.charge  = p->charge();
         genpart.pdgId   = p->pdgId();
         genpart.name    = (pdt->particle(p->pdgId()))->name();
         genpart.status  = p->status();

         Jet.JetPart.push_back(genpart);

      }
      vector<int> i;

      JetVector.push_back(Jet);

   }

}