void IsoMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){
  auto mu_pt = make_auto(new std::vector<float>);
  auto mu_phi = make_auto(new std::vector<float>);
  auto mu_eta = make_auto(new std::vector<float>);
  auto mu_iso = make_auto(new std::vector<float>);
  
  edm::Handle<edm::ValueMap<double> > iso_handle;
  edm::Handle<std::vector<reco::RecoChargedCandidate> > muon_handle;
  iEvent.getByToken(iso_token_, iso_handle);
  iEvent.getByToken(cand_token_, muon_handle);
  
  for(auto it = iso_handle->begin(); it != iso_handle->end(); ++it){
    for(auto ite = it.begin(); ite != it.end(); ++ite){
      mu_iso->push_back(*ite);
    }
  }
  
  for(auto it = muon_handle->begin(); it != muon_handle->end(); ++it){
    mu_pt->push_back(it->pt());
    mu_phi->push_back(it->phi());
    mu_eta->push_back(it->eta());
  }

  iEvent.put(mu_pt, "mupt");
  iEvent.put(mu_phi, "muphi");
  iEvent.put(mu_eta, "mueta");
  iEvent.put(mu_iso, "muiso");
}

void GlobalPositionRcdScan::analyze(const edm::Event& evt, const edm::EventSetup& evtSetup)
{
  lastRun_ = evt.run();
  if (0 == firstRun_) firstRun_ = lastRun_;

  if (watcher_.check(evtSetup)) {
  
    edm::ESHandle<Alignments> globalPositionRcd;
    evtSetup.get<GlobalPositionRcd>().get(globalPositionRcd);

    std::cout << "=====================================================\n" 
	      << "GlobalPositionRcd content starting from run " << evt.run() << ":" << std::endl;
  
    for (std::vector<AlignTransform>::const_iterator i = globalPositionRcd->m_align.begin();
	 i != globalPositionRcd->m_align.end();  ++i) {
      std::cout << "  Component ";
      if (i->rawId() == DetId(DetId::Tracker).rawId()) {
	std::cout << "Tracker";
      }
      else if (i->rawId() == DetId(DetId::Muon).rawId()) {
	std::cout << "Muon   ";
      }
      else if (i->rawId() == DetId(DetId::Ecal).rawId()) {
	std::cout << "Ecal   ";
      }
      else if (i->rawId() == DetId(DetId::Hcal).rawId()) {
	std::cout << "Hcal   ";
      }
      else if (i->rawId() == DetId(DetId::Calo).rawId()) {
	std::cout << "Calo   ";
      }
      else {
	std::cout << "Unknown";
      }
      std::cout << " entry " << i->rawId() 
		<< "\n     translation " << i->translation() << "\n";
      const AlignTransform::Rotation hepRot(i->rotation());
      if (eulerAngles_) {
	std::cout << "     euler angles " << hepRot.eulerAngles() << std::endl;
      }
      if (alignAngles_) {
	const align::RotationType matrix(hepRot.xx(), hepRot.xy(), hepRot.xz(),
					 hepRot.yx(), hepRot.yy(), hepRot.yz(),
					 hepRot.zx(), hepRot.zy(), hepRot.zz());
	const AlgebraicVector angles(align::toAngles(matrix));
	std::cout << "     alpha, beta, gamma (" << angles[0] << ", " << angles[1] << ", "
		  << angles[2] << ')' << std::endl;
      }
      if (matrix_) {
	std::cout << "     rotation matrix " << hepRot << std::endl;
      }
    }
  }
}

 void 
 SleepingProducer::produce(edm::Event& iEvent) {
   //printf("Producer %s\n",label().c_str());
   int sum=0;
   for(std::vector<const Getter*>::iterator it = m_getters.begin(), itEnd=m_getters.end();
       it != itEnd;
       ++it) {
     sum +=iEvent.get(*it);
     //printf("%s got %s with value %i\n",label().c_str(), (*it)->label().c_str(), iEvent.get((*it)));
   }
   wait(iEvent);
   iEvent.put(this,"",static_cast<int>(sum));
 }

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

 bool 
 EventTimesPassFilterBase::filter(const edm::Event& iEvent) {
   //printf("Producer %s\n",label().c_str());
   int sum=0;
   for(std::vector<const Getter*>::iterator it = m_getters.begin(), itEnd=m_getters.end();
       it != itEnd;
       ++it) {
     sum +=iEvent.get(*it);
     //printf("%s got %s with value %i\n",label().c_str(), (*it)->label().c_str(), iEvent.get((*it)));
   }
   unsigned long index = iEvent.index() % m_eventTimes.size();
   
   wait(m_eventTimes[index]);
   return true;
 }

int NeroPF::analyze(const edm::Event& iEvent)
{
	
    	iEvent.getByToken( token , handle);
        if ( not handle.isValid() ) cout<<"[NeroPF]::[analyze]::[ERROR] handle is not valid"<<endl;
	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);
         

    } 

}

EventInfo::EventInfo(const edm::Event& iEvent, const edm::Handle<reco::BeamSpot>& beamSpot, const edm::Handle<HcalNoiseSummary>& hcalNoiseSummary, const edm::Handle< std::vector<PileupSummaryInfo> >& puSummaryInfo, bool firstVertexIsGood_, float rho_):firstVertexIsGood(firstVertexIsGood_),rho(rho_)
{
  
  evt         = iEvent.id().event();
  run         = iEvent.id().run();
  lumiBlock   = iEvent.luminosityBlock();
  
  BSx = beamSpot->x0();
  BSy = beamSpot->y0();
  BSz = beamSpot->z0();
  
 
  if( hcalNoiseSummary.isValid() ){
    hcalnoiseLoose = hcalNoiseSummary->passLooseNoiseFilter();
    hcalnoiseTight = hcalNoiseSummary->passTightNoiseFilter();
  }
  
  if( puSummaryInfo.isValid() ){
    for(std::vector<PileupSummaryInfo>::const_iterator PVI = puSummaryInfo->begin(); PVI != puSummaryInfo->end(); ++PVI) {

      int BX = PVI->getBunchCrossing();

      sumNVtx  += float(PVI->getPU_NumInteractions());
      sumTrueNVtx += float(PVI->getTrueNumInteractions());

      if( BX==0 ){
	      numGenPV = PVI->getPU_NumInteractions();
	      numTruePV = PVI->getTrueNumInteractions();
      }

      if(BX == -1) { 
	      nm1 = PVI->getPU_NumInteractions();
	      nm1_true = PVI->getTrueNumInteractions();
      }
      else if(BX == 0) { 
	      n0 = PVI->getPU_NumInteractions();
	      n0_true = PVI->getTrueNumInteractions();
      }
      else if(BX == 1) { 
	      np1 = PVI->getPU_NumInteractions();
	      np1_true = PVI->getTrueNumInteractions();
      }
    }
  }
    
  
}

 bool 
 SleepingPassFilter::filter(const edm::Event& iEvent) {
   //printf("Producer %s\n",label().c_str());
   int sum=0;
   for(std::vector<const Getter*>::iterator it = m_getters.begin(), itEnd=m_getters.end();
       it != itEnd;
       ++it) {
     sum +=iEvent.get(*it);
     //printf("%s got %s with value %i\n",label().c_str(), (*it)->label().c_str(), iEvent.get((*it)));
   }
   wait(iEvent);
   return true;
 }

int NeroPuppiJets::analyze(const edm::Event& iEvent, const edm::EventSetup &iSetup){

    if ( mOnlyMc  ) return 0;

    // maybe handle should be taken before
    iEvent.getByToken(token, handle);

    if ( not handle.isValid() ) cout<<"[NeroPuppiJets]::[analyze]::[ERROR] handle is not valid"<<endl;

    for (const pat::Jet& j : *handle)
    {

        if (j.pt() < mMinPt ) continue;
        if ( fabs(j.eta()) > mMinEta ) continue;
        if ( !JetId(j,mMinId) ) continue;


        float charge =  0.;
        float charge_den = 0.;

        for( size_t idx =0; idx < j.numberOfDaughters() ; ++idx)
        {
            pat::PackedCandidate *cand = ( pat::PackedCandidate* ) j.daughter(idx) ; 
            if (cand->charge() !=0 ) {  
                charge     += cand->charge() * cand->puppiWeight() * ( j.px()*cand->px() + j.py()*cand->py() + j.pz()*cand->pz()  ) ;
                charge_den +=                  cand->puppiWeight() * ( j.px()*cand->px() + j.py()*cand->py() + j.pz()*cand->pz()  ) ;
            }
        }

        if ( charge_den == 0 ) { charge=0.0 ; charge_den =1.0;}  //  guard, if no jet id

        // Fill output object	
        new ( (*p4)[p4->GetEntriesFast()]) TLorentzVector(j.px(), j.py(), j.pz(), j.energy());
        rawPt  -> push_back (j.pt()*j.jecFactor("Uncorrected"));
        bDiscr -> push_back( j.bDiscriminator("pfCombinedInclusiveSecondaryVertexV2BJetTags") );
        unsigned bits=0;
        bits |=  (1 * JetBaseline);
        bits |= JetId(j,"monojet") * mjId;
        bits |= JetId(j,"monojetloose") * mjIdLoose;
        bits |= JetId(j,"monojet2015") * mjId2015;
        bits |= JetId(j,"loose") * JetLoose;
        bits |= JetId(j,"tight") * JetTight;

        selBits -> push_back( bits);
        Q -> push_back (charge/charge_den);
    }

    if ( int(rawPt -> size()) < mMinNjets ) return 1;

    return 0;
}

int NeroTaus::analyze(const edm::Event & iEvent)
{
    if ( mOnlyMc  ) return 0;

    iEvent.getByToken(token, handle);
    for (const pat::Tau &tau : *handle) {

        if (tau.pt() < 18 ) continue;	
        if (tau.pt() < mMinPt ) continue;	
        
        /// miniaod taus = decayModeFindingNewDMs
        if ( mMinId != "" and !(tau.tauID(mMinId)) ) continue; // minimum requirement to be saved.
        if ( mMaxIso >=0 and tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") >= mMaxIso ) continue;

        if ( fabs(tau.eta()) > mMinEta ) continue;
        
        // ------------ END SELECTION 

         float phoIso = 0.; for(auto cand : tau.isolationGammaCands() ) phoIso += cand->pt();//tau.isolationPFGammaCandsEtSum() ;
         float chIso  = 0.; for(auto cand : tau.isolationChargedHadrCands() ) chIso += cand->pt();//tau.isolationPFChargedHadrCandsPtSum();
         float nhIso  = 0.; for(auto cand : tau.isolationNeutrHadrCands() ) nhIso += cand->pt(); // PF Cands not exists in miniAOD
         float totIso = phoIso + chIso + nhIso;
        
        //FILL
        new ( (*p4)[p4->GetEntriesFast()]) TLorentzVector(tau.px(), tau.py(), tau.pz(), tau.energy());
        id -> push_back( tau.tauID("decayModeFinding"));
        Q -> push_back( tau.charge() );
        M -> push_back( tau.mass() );
        iso -> push_back( totIso ) ; 

        if (IsExtend() ){
            chargedIsoPtSum  -> push_back( tau.tauID("chargedIsoPtSum") );
            neutralIsoPtSum  -> push_back( tau.tauID("neutralIsoPtSum") );
            isoDeltaBetaCorr -> push_back( tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));

            againstEleLoose  -> push_back( tau.tauID("againstElectronLooseMVA5") );
            againstEleMedium -> push_back( tau.tauID("againstElectronMediumMVA5") );
            
            againstMuLoose   -> push_back( tau.tauID("againstMuonLoose3"));
            againstMuTight   -> push_back( tau.tauID("againstMuonTight3"));
        }


    }
    if( int(id->size()) < mMinNtaus) return 1;
    return 0;
}

int NeroFatJets::analyze(const edm::Event& iEvent){

	if ( mOnlyMc  ) return 0;

	// maybe handle should be taken before
	iEvent.getByToken(token, handle);

	int ijetRef = -1;
    int nsubjet = 0;
	for (const pat::Jet& j : *handle)
		{
		ijetRef++;
		if (j.pt() < 100 ) continue;

		// JET ID
		if ( !NeroJets::JetId(j,"loose") ) continue;
		
		// GET  ValueMaps
		
		// Fill output object	
		//p4 -> AddLast(new TLorentzVector(j.px(), j.py(), j.pz(), j.energy())  );
		new ( (*p4)[p4->GetEntriesFast()]) TLorentzVector(j.px(), j.py(), j.pz(), j.energy());

		rawPt -> push_back (j.pt()*j.jecFactor("Uncorrected"));
        	flavour -> push_back( j.partonFlavour() );
		tau1 -> push_back(j.userFloat("NjettinessAK8:tau1"));
		tau2 -> push_back(j.userFloat("NjettinessAK8:tau2"));
		tau3 -> push_back(j.userFloat("NjettinessAK8:tau3"));

		trimmedMass ->push_back(j.userFloat("ak8PFJetsCHSTrimmedMass"));
		prunedMass  ->push_back(j.userFloat("ak8PFJetsCHSPrunedMass"));
		filteredMass->push_back(j.userFloat("ak8PFJetsCHSFilteredMass"));
		softdropMass->push_back(j.userFloat("ak8PFJetsCHSSoftDropMass"));
        ak8jet_hasSubjet->push_back(j.hasSubjets("SoftDrop"));

        auto Subjets = j.subjets("SoftDrop");
        for ( auto const & i : Subjets ) {
            new ( (*ak8_subjet)[nsubjet]) TLorentzVector(i->px(), i->py(), i->pz(), i->energy());
            ak8subjet_btag->push_back(i->bDiscriminator("pfCombinedInclusiveSecondaryVertexV2BJetTags"));
            nsubjet++;
        }
        
		}
	return 0;
}

int NeroPhotons::analyze(const edm::Event& iEvent){

    if ( mOnlyMc  ) return 0;

    // maybe handle should be taken before
    iEvent.getByToken(token, handle);

    // ID and ISO
    iEvent.getByToken(loose_id_token,loose_id);
    iEvent.getByToken(medium_id_token,medium_id);
    iEvent.getByToken(tight_id_token,tight_id);
    iEvent.getByToken(iso_ch_token, iso_ch);
    iEvent.getByToken(iso_nh_token, iso_nh);
    iEvent.getByToken(iso_pho_token, iso_pho);

    int iPho = -1;	
    for (auto &pho : *handle)
    {
        ++iPho;

        if (pho.pt() <15 or pho.chargedHadronIso()/pho.pt() > 0.3) continue;		
        if (fabs(pho.eta()) > mMinEta ) continue;
        if (pho.pt() < mMinPt) continue;

        edm::RefToBase<pat::Photon> ref ( edm::Ref< pat::PhotonCollection >(handle, iPho) ) ;
        float chIso =  (*iso_ch) [ref];
        float nhIso =  (*iso_nh) [ref];
        float phIso =  (*iso_pho)[ref];	
        float totIso = chIso + nhIso + phIso;

        bool isPassLoose  = (*loose_id)[ref];	
        bool isPassMedium = (*medium_id)[ref];	
        bool isPassTight  = (*tight_id)[ref];	

        if (not isPassLoose) continue;
        if (mMaxIso >=0 and totIso > mMaxIso) continue;

        //FILL
        new ( (*p4)[p4->GetEntriesFast()]) TLorentzVector(pho.px(),pho.py(),pho.pz(),pho.energy());
        iso->push_back(totIso);	
        sieie -> push_back( pho. full5x5_sigmaIetaIeta() );


        tightid->push_back(isPassTight);
        mediumid->push_back(isPassMedium);
    }

    return 0;
}

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

int NeroTrigger::analyze(const edm::Event& iEvent){

	if ( mOnlyMc  ) return 0;

	iEvent.getByToken(token,handle);
	iEvent.getByToken(prescale_token,prescale_handle);

	 const pat::helper::TriggerMatchHelper matchHelper; 

	 triggerFired   ->resize( triggerNames->size(),0);
	 triggerPrescale->resize( triggerNames->size(),0);

	 const edm::TriggerNames &names = iEvent.triggerNames(*handle);

	 if ( handle.isValid() and not handle.failedToGet() ) {
		int ntrigs = (int)names.size();
		for (int i = 0; i < ntrigs; i++) {
			string name = names.triggerName(i);
			for(unsigned int j=0;j< triggerNames->size() ;++j) // destination loop
				{
				if (name.find( (*triggerNames)[j]) != string::npos)
					(*triggerFired)[j] = 1;
					(*triggerPrescale)[j] = prescale_handle -> getPrescaleForIndex(i) ;

				} // my trigger end
		} // trigger loop
	 	} // handle is valid

	 // ---- TRIGGER MATCHING ---
	 if (leps_ !=NULL) triggerLeps -> resize(leps_ -> p4 -> GetEntries()  ,0);
	 if (jets_ !=NULL) triggerJets -> resize(jets_ -> p4 -> GetEntries()  ,0);
	 if (taus_ !=NULL) triggerTaus -> resize(taus_ -> p4 -> GetEntries()  ,0);
	 if (photons_ !=NULL) triggerPhotons -> resize(photons_ -> p4 -> GetEntries()  ,0);
	

	 iEvent.getByToken(object_token,object_handle);
	 for (pat::TriggerObjectStandAlone obj : *object_handle) { 
		obj.unpackPathNames(names);

		bool isPhoton  =false;
		bool isElectron=false;
		bool isMuon    =false;
		bool isTau     =false;
		//bool isMet     =false;
		bool isJet     =false;
		bool isBJet    =false;

        	for (unsigned h = 0; h < obj.filterIds().size(); ++h)  // can be more than one
		{
			// HLT: DataFormats/HLTReco/interface/TriggerTypeDefs.h
			if      ( obj.filterIds()[h] == trigger::TriggerPhoton ) isPhoton=true;
			else if ( obj.filterIds()[h] == trigger::TriggerElectron) isElectron=true;
			else if ( obj.filterIds()[h] == trigger::TriggerMuon ) isMuon=true;
			else if ( obj.filterIds()[h] == trigger::TriggerTau ) isTau=true;
			else if ( obj.filterIds()[h] == trigger::TriggerJet ) isJet=true;
			else if ( obj.filterIds()[h] == trigger::TriggerBJet ) isBJet=true;
			//else if ( obj.filterIds()[h] == trigger::TriggerMET ) isMet=true;
		}
		std::vector<std::string> pathNamesAll  = obj.pathNames(false);
        std::vector<std::string> pathNamesLast = obj.pathNames(true);
            
		// match : if more than one ? best
		if(VERBOSE)cout<<"[NeroTrigger]::[analize] Matching"<<endl;

		int muonIndex = -1 ; if( isMuon ) muonIndex = match(leps_,obj,13);

		int eleIndex=-1 ; if (isElectron) eleIndex = match(leps_,obj,11);

		int jetIndex=-1; if (isJet or isBJet) jetIndex = match(jets_,obj);

		int tauIndex=-1; if(isTau) tauIndex = match(taus_,obj);

		int photonIndex=-1; if (isPhoton) photonIndex = match(photons_,obj);
		//
        // Print all trigger paths, for each one record also if the object is associated to a 'l3' filter (always true for the
        // definition used in the PAT trigger producer) and if it's associated to the last filter of a successfull path (which
        // means that this object did cause this trigger to succeed; however, it doesn't work on some multi-object triggers)
		if(VERBOSE)cout<<"[NeroTrigger]::[analize] Lopping on path names and filling"<<endl;
        	for (unsigned h = 0, n = pathNamesAll.size(); h < n; ++h) {
			if(VERBOSE)cout<<"[NeroTrigger]::[analize] pathNames h="<<h<<endl;
           		bool isBoth = obj.hasPathName( pathNamesAll[h], true, true ); 
           		bool isL3   = obj.hasPathName( pathNamesAll[h], false, true ); 
           		bool isLF   = obj.hasPathName( pathNamesAll[h], true, false ); 
           		bool isNone = obj.hasPathName( pathNamesAll[h], false, false ); 
            		//std::cout << "   " << pathNamesAll[h];
            		if (isNone && !isBoth && !isL3 && !isLF) continue;
			for ( unsigned i =0  ; i<triggerNames->size() ;++i)
				{
				if(VERBOSE)cout<<"[NeroTrigger]::[analize] triggerNames i="<<i<<endl;
				if (  pathNamesAll[h].find((*triggerNames)[i] ) !=string::npos )
					{
					if(VERBOSE) cout <<"Considering SUCCESS Path: eleIndex"<< eleIndex <<"/"<<triggerLeps->size()<<" :"<<pathNamesAll[h]<<" and Match string "<<(*triggerNames)[i]<<endl;
					if(VERBOSE)cout<<"[NeroTrigger]::[analize] ----- Mu "<<muonIndex<<endl;
					if (muonIndex >=0) (*triggerLeps)[muonIndex] |= 1<<i;
					if(VERBOSE)cout<<"[NeroTrigger]::[analize] ----- E "<<eleIndex<<endl;
					if (eleIndex >=0) (*triggerLeps)[eleIndex] |= 1<<i;
					if(VERBOSE)cout<<"[NeroTrigger]::[analize] ----- Tau "<<tauIndex <<endl;
					if (tauIndex >=0) (*triggerTaus)[tauIndex] |= 1<<i;
					if(VERBOSE)cout<<"[NeroTrigger]::[analize] ----- Jet "<< jetIndex<<endl;
					if (jetIndex >=0) (*triggerJets)[jetIndex] |= 1<<i;
//.........这里部分代码省略.........

int NeroPhotons::analyze(const edm::Event& iEvent,const edm::EventSetup &iSetup){

    if ( mOnlyMc  ) return 0;

    // maybe handle should be taken before
    iEvent.getByToken(token, handle);

    // ID and ISO
    iEvent.getByToken(loose_id_token,loose_id);
    iEvent.getByToken(medium_id_token,medium_id);
    iEvent.getByToken(tight_id_token,tight_id);
    iEvent.getByToken(iso_ch_token, iso_ch);
    iEvent.getByToken(iso_nh_token, iso_nh);
    iEvent.getByToken(iso_pho_token, iso_pho);

    if ( not handle.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] handle is not valid"<<endl;
    if ( not loose_id.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] loose_id is not valid"<<endl;
    if ( not medium_id.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] medium_id is not valid"<<endl;
    if ( not tight_id.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] tight_id is not valid"<<endl;
    if ( not iso_ch.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] iso_ch is not valid"<<endl;
    if ( not iso_nh.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] iso_nh is not valid"<<endl;
    if ( not iso_pho.isValid() ) cout<<"[NeroPhotons]::[analyze]::[ERROR] iso_pho is not valid"<<endl;

    int iPho = -1;	
    for (auto &pho : *handle)
    {
        ++iPho;
        #ifdef VERBOSE
        if (VERBOSE>0) cout<<"[NeroPhotons]::[analyze]::[DEBUG] analyzing photon"<<iPho<<" pt="<<pho.pt() <<" pz"<<pho.pz() <<endl;
        #endif

        //if ( not pho.passElectronVeto ()  ) continue;

        // r9()>0.8 , chargedHadronIso()<20, chargedHadronIso()<0.3*pt()
        if (pho.pt() <15 or pho.chargedHadronIso()/pho.pt() > 0.3) continue; // 10 -- 14  GeV photons are saved if chargedHadronIso()<10
        if (fabs(pho.eta()) > mMaxEta ) continue;
        if (pho.pt() < mMinPt) continue;

        #ifdef VERBOSE
        if (VERBOSE>1) cout<<"[NeroPhotons]::[analize]::[DEBUG2] photonInfo:" <<endl
            <<" \t pho.chargedHadronIso()/pho.pt() (0.3) "<<pho.chargedHadronIso()/pho.pt() <<endl
            <<" \t chargedHadronIso() (20) "<<pho.chargedHadronIso()<<endl
            <<" \t r9 (0.8) "<<pho.r9()<<endl
            <<" \t SC is non null? "<< pho.superCluster().isNonnull()<<endl
            <<endl;
        #endif

        edm::RefToBase<pat::Photon> ref ( edm::Ref< pat::PhotonCollection >(handle, iPho) ) ;
        float _chIso_ =  (*iso_ch) [ref];
        float _nhIso_ =  (*iso_nh) [ref];
        float _phIso_ =  (*iso_pho)[ref];	
        float _puIso_ =  pho.puChargedHadronIso() ;  // the other are eff area corrected, no need for this correction
        float totIso = _chIso_ + _nhIso_ + _phIso_;

        bool isPassLoose  = (*loose_id)[ref];	
        bool isPassMedium = (*medium_id)[ref];	
        bool isPassTight  = (*tight_id)[ref];	
        bool isPassVLoose50 = cutBasedPhotonId( pho, "loose_50ns", false, false); // no pho iso , no sieie
        bool isPassVLoose25 = cutBasedPhotonId( pho, "loose_25ns", false, false); // no pho iso , no sieie

        //if (not isPassVLoose) continue;
        if (mMaxIso >=0 and totIso > mMaxIso) continue;

        unsigned bits=0;

        bits |= isPassTight * PhoTight;
        bits |= isPassMedium * PhoMedium;
        bits |= isPassLoose * PhoLoose;
        bits |= isPassVLoose50 * PhoVLoose50;
        bits |= isPassVLoose25 * PhoVLoose25;

        if (not bits) continue; // even if there is some misalignment ntuples will not be corrupted

        bits |= pho.passElectronVeto() * PhoElectronVeto;

        // RC -- with FPR
        //
        float _chIsoRC_ = 0;
        float _nhIsoRC_ = 0;
        float _phIsoRC_ = 0;
        float _puIsoRC_ = 0;// not fill for the moment in the FPR TODO
        
        if (  pho.chargedHadronIso()< 20 )
        {
                                        //<<" \t r9 (0.8) "<<pho.r9()<<endl
        #ifdef VERBOSE
            if (VERBOSE >0 ) cout <<"[NeroPhotons]::[analyze]::[DEBUG] FPR START"<<endl;
        #endif

        fpr -> Config(iSetup);
        fpr -> SetHandles(
                pf  -> handle,
                handle,
                jets-> handle,
                leps->mu_handle,
                leps->el_handle
                );

        PFIsolation_struct FPR_out = fpr -> PFIsolation(pho.superCluster(), edm::Ptr<reco::Vertex>(vtx->handle,vtx->firstGoodVertexIdx) );
        _chIsoRC_ = FPR_out.chargediso_primvtx_rcone;
//.........这里部分代码省略.........

void AlignmentRcdScan::inspectRecord(const std::string & rcdname,
				     const edm::Event& evt, 
				     const edm::ESHandle<Alignments> & alignments)
{
  std::cout << rcdname << " content starting from run " << evt.run();
  
  if (verbose_==false) {
    std::cout << std::endl;
    return;
  }

  std::cout << " with " << alignments->m_align.size() << " entries" << std::endl;
  
  if (refAlignments_) {

    std::cout << "  Compared to previous record:" << std::endl;
    
    double meanX = 0;
    double rmsX = 0;
    double meanY = 0;
    double rmsY = 0;
    double meanZ = 0;
    double rmsZ = 0;
    double meanR = 0;
    double rmsR = 0;
    double dPhi;
    double meanPhi = 0;
    double rmsPhi = 0;
    
    std::vector<AlignTransform>::const_iterator iref = refAlignments_->m_align.begin();
    for (std::vector<AlignTransform>::const_iterator i = alignments->m_align.begin();
	 i != alignments->m_align.end();
	 ++i, ++iref) {
      
      meanX += i->translation().x() - iref->translation().x();
      rmsX += pow(i->translation().x() - iref->translation().x(), 2);
 
      meanY += i->translation().y() - iref->translation().y();
      rmsY += pow(i->translation().y() - iref->translation().y(), 2);
      
      meanZ += i->translation().z() - iref->translation().z();
      rmsZ += pow(i->translation().z() - iref->translation().z(), 2);
      
      meanR += i->translation().perp() - iref->translation().perp();
      rmsR += pow(i->translation().perp() - iref->translation().perp(), 2);

      dPhi = i->translation().phi() - iref->translation().phi();
      if (dPhi>TMath::Pi()) dPhi -= 2.0*TMath::Pi();
      if (dPhi<-TMath::Pi()) dPhi += 2.0*TMath::Pi();

      meanPhi += dPhi;
      rmsPhi += dPhi*dPhi;
    }

    meanX /= alignments->m_align.size();
    rmsX /= alignments->m_align.size();
    meanY /= alignments->m_align.size();
    rmsY /= alignments->m_align.size();
    meanZ /= alignments->m_align.size();
    rmsZ /= alignments->m_align.size();
    meanR /= alignments->m_align.size();
    rmsR /= alignments->m_align.size();
    meanPhi /= alignments->m_align.size();
    rmsPhi /= alignments->m_align.size();
    
    std::cout << "    mean X shift:   " 
	      << std::setw(12) << std::scientific << std::setprecision(3) << meanX
	      << " (RMS = " << sqrt(rmsX) << ")" << std::endl;
    std::cout << "    mean Y shift:   " 
	      << std::setw(12) << std::scientific << std::setprecision(3) << meanY
	      << " (RMS = " << sqrt(rmsY) << ")" << std::endl;
    std::cout << "    mean Z shift:   " 
	      << std::setw(12) << std::scientific << std::setprecision(3) << meanZ
	      << " (RMS = " << sqrt(rmsZ) << ")" << std::endl;
    std::cout << "    mean R shift:   " 
	      << std::setw(12) << std::scientific << std::setprecision(3) << meanR
	      << " (RMS = " << sqrt(rmsR) << ")" << std::endl;
    std::cout << "    mean Phi shift: " 
	      << std::setw(12) << std::scientific << std::setprecision(3) << meanPhi
	      << " (RMS = " << sqrt(rmsPhi) << ")" << std::endl;
    
    delete refAlignments_;
  }

  refAlignments_ = new Alignments(*alignments);

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

void
TreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   // get gen particles before photons for the truth match
   edm::Handle<edm::View<reco::GenParticle>> prunedGenParticles;
   iEvent.getByToken(prunedGenToken_,prunedGenParticles);

   edm::Handle<edm::ValueMap<bool>> loose_id_dec;
   edm::Handle<edm::ValueMap<bool>> medium_id_dec;
   edm::Handle<edm::ValueMap<bool>> tight_id_dec;
   edm::Handle<edm::ValueMap<float>> mva_value;
   edm::Handle<edm::ValueMap<vid::CutFlowResult>> loose_id_cutflow;
   iEvent.getByToken(photonLooseIdMapToken_, loose_id_dec);
   iEvent.getByToken(photonMediumIdMapToken_, medium_id_dec);
   iEvent.getByToken(photonTightIdMapToken_, tight_id_dec);
   iEvent.getByToken(photonMvaValuesMapToken_, mva_value);
   iEvent.getByToken(phoLooseIdFullInfoMapToken_, loose_id_cutflow);

   // photon collection
   edm::Handle<edm::View<pat::Photon>> photonColl;
   iEvent.getByToken(photonCollectionToken_, photonColl);

   for (edm::View<pat::Photon>::const_iterator pho = photonColl->begin(); pho != photonColl->end(); pho++) {
      // some basic selections
      if (pho->pt() < 15 || pho->hadTowOverEm() > 0.15) continue;
      const edm::Ptr<pat::Photon> phoPtr(photonColl, pho - photonColl->begin());

      pt = pho->pt();
      eta = pho->eta();
      sigmaIetaIeta=pho->full5x5_sigmaIetaIeta();
      hOverE=pho->hadTowOverEm() ;
      hasPixelSeed=pho->hasPixelSeed() ;
      passElectronVeto= pho->passElectronVeto() ;
      r9  = pho->r9();

      vid::CutFlowResult cutFlow = (*loose_id_cutflow)[phoPtr];
      cIso = cutFlow.getValueCutUpon(4);
      nIso = cutFlow.getValueCutUpon(5);
      pIso = cutFlow.getValueCutUpon(6);

      mvaValue = (*mva_value)    [phoPtr];
      isLoose  = (*loose_id_dec) [phoPtr];
      isMedium = (*medium_id_dec)[phoPtr];
      isTight  = (*tight_id_dec) [phoPtr];

      // mc matching
      bool foundPhoton = false;
      bool vetoPhoton = false;
      for (auto& genP: *prunedGenParticles) {
         if (fabs(genP.pdgId()) == 22
            && genP.statusFlags().fromHardProcess()
            && genP.status() == 1
            && ROOT::Math::VectorUtil::DeltaR(pho->p4(), genP.p4())<0.1) {
            foundPhoton = true;
         } else if (fabs(genP.pdgId()) == 22
            && ROOT::Math::VectorUtil::DeltaR(pho->p4(), genP.p4())<0.4) {
            vetoPhoton = true;
         }
      }
      if (foundPhoton) {
         isTrue = true;
         eventTree_->Fill();
      }
      if (!vetoPhoton) {
         isTrue = false;
         eventTree_->Fill();
      }

   }

   // Get the electron ID data from the event stream
   edm::Handle<edm::ValueMap<bool>> veto_id_decisions;
   edm::Handle<edm::ValueMap<bool>> loose_id_decisions;
   edm::Handle<edm::ValueMap<bool>> medium_id_decisions;
   edm::Handle<edm::ValueMap<bool>> tight_id_decisions;
   iEvent.getByToken(electronVetoIdMapToken_, veto_id_decisions);
   iEvent.getByToken(electronLooseIdMapToken_, loose_id_decisions);
   iEvent.getByToken(electronMediumIdMapToken_, medium_id_decisions);
   iEvent.getByToken(electronTightIdMapToken_, tight_id_decisions);

   edm::Handle<edm::View<pat::Electron>> electronColl;
   iEvent.getByToken(electronCollectionToken_, electronColl);

   for (edm::View<pat::Electron>::const_iterator el = electronColl->begin();el != electronColl->end(); el++) {
      const edm::Ptr<pat::Electron> elPtr (electronColl, el - electronColl->begin() );
      bool isLoose = (*loose_id_decisions) [elPtr];
      if (isLoose && false ) std::cout << el->pt() << std::endl;
   }
}

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;