13 #include <boost/shared_ptr.hpp>
48 #define M_PI 3.14159265358979323846
67 virtual void endJob();
71 int convertPhiToHW(
double iphi,
int steps);
72 int convertEtaToHW(
double ieta,
double minEta,
double maxEta,
int steps);
73 int convertPtToHW(
double ipt,
int maxPt,
double step);
140 produces<BXVector<l1t::EGamma>>();
141 produces<BXVector<l1t::Muon>>();
142 produces<BXVector<l1t::Tau>>();
143 produces<BXVector<l1t::Jet>>();
144 produces<BXVector<l1t::EtSum>>();
150 maxNumMuCands_ = iConfig.
getParameter<
int>(
"maxMuCand");
151 maxNumJetCands_ = iConfig.
getParameter<
int>(
"maxJetCand");
152 maxNumEGCands_ = iConfig.
getParameter<
int>(
"maxEGCand");
153 maxNumTauCands_ = iConfig.
getParameter<
int>(
"maxTauCand");
155 jetEtThreshold_ = iConfig.
getParameter<
double>(
"jetEtThreshold");
156 tauEtThreshold_ = iConfig.
getParameter<
double>(
"tauEtThreshold");
157 egEtThreshold_ = iConfig.
getParameter<
double>(
"egEtThreshold");
158 muEtThreshold_ = iConfig.
getParameter<
double>(
"muEtThreshold");
161 emptyBxTrailer_ = iConfig.
getParameter<
int>(
"emptyBxTrailer");
165 genParticlesToken = consumes <reco::GenParticleCollection> (
std::string(
"genParticles"));
166 genJetsToken = consumes <reco::GenJetCollection> (
std::string(
"ak5GenJets"));
167 genMetToken = consumes <reco::GenMETCollection> (
std::string(
"genMetCalo"));
176 GenToInputProducer::~GenToInputProducer()
193 LogDebug(
"l1t|Global") <<
"GenToInputProducer::produce function called...\n";
196 std::vector<l1t::Muon> muonVec;
197 std::vector<l1t::EGamma> egammaVec;
198 std::vector<l1t::Tau> tauVec;
199 std::vector<l1t::Jet> jetVec;
200 std::vector<l1t::EtSum> etsumVec;
233 std::vector<int> mu_cands_index;
234 std::vector<int> eg_cands_index;
235 std::vector<int> tau_cands_index;
238 if( iEvent.
getByToken(genParticlesToken, genParticles) ){
240 for(
size_t k = 0;
k < genParticles->size();
k++ ){
245 double pt = mcParticle.
pt();
248 if( status!=1 && !(
abs(pdgId)==15 && status==2) )
continue;
250 int absId =
abs(pdgId);
252 if( absId==11 && pt>=egEtThreshold_ ) eg_cands_index.push_back(
k);
253 else if( absId==13 && pt>=muEtThreshold_ ) mu_cands_index.push_back(
k);
254 else if( absId==15 && pt>=tauEtThreshold_ ) tau_cands_index.push_back(
k);
258 LogTrace(
"l1t|Global") <<
">>> GenParticles collection not found!" << std::endl;
264 int numMuCands = int( mu_cands_index.size() );
265 Int_t idxMu[numMuCands];
266 double muPtSorted[numMuCands];
267 for(
int iMu=0; iMu<numMuCands; iMu++ ) muPtSorted[iMu] = genParticles->at(mu_cands_index[iMu]).pt();
269 TMath::Sort(numMuCands,muPtSorted,idxMu);
270 for(
int iMu=0; iMu<numMuCands; iMu++ ){
272 if( iMu>=maxNumMuCands_ )
continue;
274 const reco::Candidate & mcParticle = (*genParticles)[mu_cands_index[idxMu[iMu]]];
281 int charge = ( mcParticle.
charge()>0 ) ? 1 : 0;
287 if( eta>=9999 )
continue;
289 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
291 l1t::Muon mu(*p4, pt, eta, phi, qual, charge, chargeValid, iso, mip, tag);
292 muonVec.push_back(mu);
297 int numEgCands = int( eg_cands_index.size() );
298 Int_t idxEg[numEgCands];
299 double egPtSorted[numEgCands];
300 for(
int iEg=0; iEg<numEgCands; iEg++ ) egPtSorted[iEg] = genParticles->at(eg_cands_index[iEg]).pt();
302 TMath::Sort(numEgCands,egPtSorted,idxEg);
303 for(
int iEg=0; iEg<numEgCands; iEg++ ){
305 if( iEg>=maxNumEGCands_ )
continue;
307 const reco::Candidate & mcParticle = (*genParticles)[eg_cands_index[idxEg[iEg]]];
316 if( eta>=9999 )
continue;
318 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
321 egammaVec.push_back(eg);
327 int numTauCands = int( tau_cands_index.size() );
328 Int_t idxTau[numTauCands];
329 double tauPtSorted[numTauCands];
330 for(
int iTau=0; iTau<numTauCands; iTau++ ) tauPtSorted[iTau] = genParticles->at(tau_cands_index[iTau]).pt();
332 TMath::Sort(numTauCands,tauPtSorted,idxTau);
333 for(
int iTau=0; iTau<numTauCands; iTau++ ){
335 if( iTau>=maxNumTauCands_ )
continue;
337 const reco::Candidate & mcParticle = (*genParticles)[tau_cands_index[idxTau[iTau]]];
346 if( eta>=9999 )
continue;
348 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
351 tauVec.push_back(tau);
357 int maxOtherTaus = 8;
358 int numCurrentEGs = int( egammaVec.size() );
359 int numCurrentTaus = int( tauVec.size() );
361 int numExtraEGs=0, numExtraTaus=0;
371 if( iEvent.
getByToken(genJetsToken, genJets) ){
372 for(reco::GenJetCollection::const_iterator genJet = genJets->begin(); genJet!=genJets->end(); ++genJet ){
375 sumEt += genJet->et();
378 if( genJet->pt()<jetEtThreshold_ )
continue;
381 if( nJet>=maxNumJetCands_ )
continue;
382 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
389 if( eta>=9999 )
continue;
394 jetVec.push_back(jet);
399 if( (numExtraEGs+numCurrentEGs)<maxNumEGCands_ && numExtraEGs<maxOtherEGs ){
404 int EGphi = convertPhiToHW( genJet->phi(),
PhiStepCalo_ );
409 l1t::EGamma eg(*p4, EGpt, EGeta, EGphi, EGqual, EGiso);
410 egammaVec.push_back(eg);
413 if( (numExtraTaus+numCurrentTaus)<maxNumTauCands_ && numExtraTaus<maxOtherTaus ){
418 int Tauphi = convertPhiToHW( genJet->phi(),
PhiStepCalo_ );
422 l1t::Tau tau(*p4, Taupt, Taueta, Tauphi, Tauqual, Tauiso);
423 tauVec.push_back(tau);
429 LogTrace(
"l1t|Global") <<
">>> GenJets collection not found!" << std::endl;
439 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
442 l1t::EtSum etmiss(*p4, l1t::EtSum::EtSumType::kMissingEt,pt, 0,phi, 0);
443 etsumVec.push_back(etmiss);
446 pt = convertPtToHW( genMet->front().pt()*0.9,
MaxEt_,
PtStep_ );
447 phi = convertPhiToHW( genMet->front().phi()+ 3.14/5.,
PhiStepCalo_ );
449 l1t::EtSum htmiss(*p4, l1t::EtSum::EtSumType::kMissingHt,pt, 0,phi, 0);
450 etsumVec.push_back(htmiss);
455 LogTrace(
"l1t|Global") <<
">>> GenMet collection not found!" << std::endl;
460 int pt = convertPtToHW( sumEt, 2047, PtStep_ );
461 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
462 l1t::EtSum etTotal(*p4, l1t::EtSum::EtSumType::kTotalEt,pt, 0, 0, 0);
463 etsumVec.push_back(etTotal);
465 pt = convertPtToHW( sumEt*0.9, 2047, PtStep_ );
466 l1t::EtSum htTotal(*p4, l1t::EtSum::EtSumType::kTotalHt,pt, 0, 0, 0);
467 etsumVec.push_back(htTotal);
471 printf(
"Event %i EmptyBxEvt %i emptyBxTrailer %i diff %i \n",eventCnt_,emptyBxEvt_,emptyBxTrailer_,(emptyBxEvt_ - eventCnt_));
474 for(
int iMu=0; iMu<int(muonVec_bxm2.size()); iMu++ ){
475 muons->push_back(-2, muonVec_bxm2[iMu]);
477 for(
int iMu=0; iMu<int(muonVec_bxm1.size()); iMu++ ){
478 muons->push_back(-1, muonVec_bxm1[iMu]);
480 for(
int iMu=0; iMu<int(muonVec_bx0.size()); iMu++ ){
481 muons->push_back(0, muonVec_bx0[iMu]);
483 for(
int iMu=0; iMu<int(muonVec_bxp1.size()); iMu++ ){
484 muons->push_back(1, muonVec_bxp1[iMu]);
486 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
487 for(
int iMu=0; iMu<int(muonVec.size()); iMu++ ){
488 muons->push_back(2, muonVec[iMu]);
496 for(
int iEG=0; iEG<int(egammaVec_bxm2.size()); iEG++ ){
497 egammas->push_back(-2, egammaVec_bxm2[iEG]);
499 for(
int iEG=0; iEG<int(egammaVec_bxm1.size()); iEG++ ){
500 egammas->push_back(-1, egammaVec_bxm1[iEG]);
502 for(
int iEG=0; iEG<int(egammaVec_bx0.size()); iEG++ ){
503 egammas->push_back(0, egammaVec_bx0[iEG]);
505 for(
int iEG=0; iEG<int(egammaVec_bxp1.size()); iEG++ ){
506 egammas->push_back(1, egammaVec_bxp1[iEG]);
508 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
509 for(
int iEG=0; iEG<int(egammaVec.size()); iEG++ ){
510 egammas->push_back(2, egammaVec[iEG]);
518 for(
int iTau=0; iTau<int(tauVec_bxm2.size()); iTau++ ){
519 taus->push_back(-2, tauVec_bxm2[iTau]);
521 for(
int iTau=0; iTau<int(tauVec_bxm1.size()); iTau++ ){
522 taus->push_back(-1, tauVec_bxm1[iTau]);
524 for(
int iTau=0; iTau<int(tauVec_bx0.size()); iTau++ ){
525 taus->push_back(0, tauVec_bx0[iTau]);
527 for(
int iTau=0; iTau<int(tauVec_bxp1.size()); iTau++ ){
528 taus->push_back(1, tauVec_bxp1[iTau]);
530 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
531 for(
int iTau=0; iTau<int(tauVec.size()); iTau++ ){
532 taus->push_back(2, tauVec[iTau]);
540 for(
int iJet=0; iJet<int(jetVec_bxm2.size()); iJet++ ){
541 jets->push_back(-2, jetVec_bxm2[iJet]);
543 for(
int iJet=0; iJet<int(jetVec_bxm1.size()); iJet++ ){
544 jets->push_back(-1, jetVec_bxm1[iJet]);
546 for(
int iJet=0; iJet<int(jetVec_bx0.size()); iJet++ ){
547 jets->push_back(0, jetVec_bx0[iJet]);
549 for(
int iJet=0; iJet<int(jetVec_bxp1.size()); iJet++ ){
550 jets->push_back(1, jetVec_bxp1[iJet]);
552 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
553 for(
int iJet=0; iJet<int(jetVec.size()); iJet++ ){
554 jets->push_back(2, jetVec[iJet]);
562 for(
int iETsum=0; iETsum<int(etsumVec_bxm2.size()); iETsum++ ){
563 etsums->push_back(-2, etsumVec_bxm2[iETsum]);
565 for(
int iETsum=0; iETsum<int(etsumVec_bxm1.size()); iETsum++ ){
566 etsums->push_back(-1, etsumVec_bxm1[iETsum]);
568 for(
int iETsum=0; iETsum<int(etsumVec_bx0.size()); iETsum++ ){
569 etsums->push_back(0, etsumVec_bx0[iETsum]);
571 for(
int iETsum=0; iETsum<int(etsumVec_bxp1.size()); iETsum++ ){
572 etsums->push_back(1, etsumVec_bxp1[iETsum]);
574 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
575 for(
int iETsum=0; iETsum<int(etsumVec.size()); iETsum++ ){
576 etsums->push_back(2, etsumVec[iETsum]);
591 muonVec_bxm2 = muonVec_bxm1;
592 egammaVec_bxm2 = egammaVec_bxm1;
593 tauVec_bxm2 = tauVec_bxm1;
594 jetVec_bxm2 = jetVec_bxm1;
595 etsumVec_bxm2 = etsumVec_bxm1;
597 muonVec_bxm1 = muonVec_bx0;
598 egammaVec_bxm1 = egammaVec_bx0;
599 tauVec_bxm1 = tauVec_bx0;
600 jetVec_bxm1 = jetVec_bx0;
601 etsumVec_bxm1 = etsumVec_bx0;
603 muonVec_bx0 = muonVec_bxp1;
604 egammaVec_bx0 = egammaVec_bxp1;
605 tauVec_bx0 = tauVec_bxp1;
606 jetVec_bx0 = jetVec_bxp1;
607 etsumVec_bx0 = etsumVec_bxp1;
609 muonVec_bxp1 = muonVec;
610 egammaVec_bxp1 = egammaVec;
611 tauVec_bxp1 = tauVec;
612 jetVec_bxp1 = jetVec;
613 etsumVec_bxp1 = etsumVec;
625 GenToInputProducer::endJob() {
632 LogDebug(
"l1t|Global") <<
"GenToInputProducer::beginRun function called...\n";
644 int GenToInputProducer::convertPhiToHW(
double iphi,
int steps){
646 double phiMax = 2 *
M_PI;
647 if( iphi < 0 ) iphi += 2*
M_PI;
648 if( iphi > phiMax) iphi -= phiMax;
650 int hwPhi = int( (iphi/phiMax)*steps + 0.00001 );
656 double binWidth = (maxEta -
minEta)/steps;
659 if(ieta < minEta)
return 99999;
660 if(ieta > maxEta)
return 99999;
662 int binNum = (int)(ieta/binWidth);
663 if(ieta<0.) binNum--;
671 int GenToInputProducer::convertPtToHW(
double ipt,
int maxPt,
double step){
673 int hwPt = int( ipt/step + 0.0001 );
675 if( hwPt > maxPt ) hwPt = maxPt;
T getParameter(std::string const &) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
virtual float eta() const =0
momentum pseudorapidity
#define DEFINE_FWK_MODULE(type)
virtual int status() const =0
status word
virtual float phi() const =0
momentum azimuthal angle
virtual float pt() const =0
transverse momentum
void addDefault(ParameterSetDescription const &psetDescription)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Abs< T >::type abs(const T &t)
virtual int charge() const =0
electric charge
virtual int pdgId() const =0
PDG identifier.
TLorentzVector genMet(const HepMC::GenEvent *all, double etamin=-9999., double etamax=9999.)