13 #include <boost/shared_ptr.hpp>
50 #define M_PI 3.14159265358979323846
69 virtual void endJob();
73 int convertPhiToHW(
double iphi,
int steps);
74 int convertEtaToHW(
double ieta,
double minEta,
double maxEta,
int steps);
75 int convertPtToHW(
double ipt,
int maxPt,
double step);
144 produces<BXVector<l1t::EGamma>>();
145 produces<BXVector<l1t::Muon>>();
146 produces<BXVector<l1t::Tau>>();
147 produces<BXVector<l1t::Jet>>();
148 produces<BXVector<l1t::EtSum>>();
154 maxNumMuCands_ = iConfig.
getParameter<
int>(
"maxMuCand");
155 maxNumJetCands_ = iConfig.
getParameter<
int>(
"maxJetCand");
156 maxNumEGCands_ = iConfig.
getParameter<
int>(
"maxEGCand");
157 maxNumTauCands_ = iConfig.
getParameter<
int>(
"maxTauCand");
159 jetEtThreshold_ = iConfig.
getParameter<
double>(
"jetEtThreshold");
160 tauEtThreshold_ = iConfig.
getParameter<
double>(
"tauEtThreshold");
161 egEtThreshold_ = iConfig.
getParameter<
double>(
"egEtThreshold");
162 muEtThreshold_ = iConfig.
getParameter<
double>(
"muEtThreshold");
165 emptyBxTrailer_ = iConfig.
getParameter<
int>(
"emptyBxTrailer");
169 genParticlesToken = consumes <reco::GenParticleCollection> (
std::string(
"genParticles"));
170 genJetsToken = consumes <reco::GenJetCollection> (
std::string(
"ak4GenJets"));
171 genMetToken = consumes <reco::GenMETCollection> (
std::string(
"genMetCalo"));
180 GenToInputProducer::~GenToInputProducer()
197 LogDebug(
"l1t|Global") <<
"GenToInputProducer::produce function called...\n";
200 std::vector<l1t::Muon> muonVec;
201 std::vector<l1t::EGamma> egammaVec;
202 std::vector<l1t::Tau> tauVec;
203 std::vector<l1t::Jet> jetVec;
204 std::vector<l1t::EtSum> etsumVec;
237 std::vector<int> mu_cands_index;
238 std::vector<int> eg_cands_index;
239 std::vector<int> tau_cands_index;
242 if( iEvent.
getByToken(genParticlesToken, genParticles) ){
244 for(
size_t k = 0;
k < genParticles->size();
k++ ){
249 double pt = mcParticle.
pt();
252 if( status!=1 && !(
abs(pdgId)==15 && status==2) )
continue;
254 int absId =
abs(pdgId);
256 if( absId==11 && pt>=egEtThreshold_ ) eg_cands_index.push_back(
k);
257 else if( absId==13 && pt>=muEtThreshold_ ) mu_cands_index.push_back(
k);
258 else if( absId==15 && pt>=tauEtThreshold_ ) tau_cands_index.push_back(
k);
262 LogTrace(
"l1t|Global") <<
">>> GenParticles collection not found!" << std::endl;
268 int numMuCands = int( mu_cands_index.size() );
269 Int_t idxMu[numMuCands];
270 double muPtSorted[numMuCands];
271 for(
int iMu=0; iMu<numMuCands; iMu++ ) muPtSorted[iMu] = genParticles->at(mu_cands_index[iMu]).pt();
273 TMath::Sort(numMuCands,muPtSorted,idxMu);
274 for(
int iMu=0; iMu<numMuCands; iMu++ ){
276 if( iMu>=maxNumMuCands_ )
continue;
278 const reco::Candidate & mcParticle = (*genParticles)[mu_cands_index[idxMu[iMu]]];
283 int qual = gRandom->Integer(16);
284 int iso = gRandom->Integer(4)%2;
291 if( eta>=9999 )
continue;
293 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
295 l1t::Muon mu(*p4, pt, eta, phi, qual, charge, chargeValid, iso, mip, tag);
296 muonVec.push_back(mu);
301 int numEgCands = int( eg_cands_index.size() );
302 Int_t idxEg[numEgCands];
303 double egPtSorted[numEgCands];
304 for(
int iEg=0; iEg<numEgCands; iEg++ ) egPtSorted[iEg] = genParticles->at(eg_cands_index[iEg]).pt();
306 TMath::Sort(numEgCands,egPtSorted,idxEg);
307 for(
int iEg=0; iEg<numEgCands; iEg++ ){
309 if( iEg>=maxNumEGCands_ )
continue;
311 const reco::Candidate & mcParticle = (*genParticles)[eg_cands_index[idxEg[iEg]]];
317 int iso = gRandom->Integer(4)%2;
320 if( eta>=9999 )
continue;
322 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
325 egammaVec.push_back(eg);
331 int numTauCands = int( tau_cands_index.size() );
332 Int_t idxTau[numTauCands];
333 double tauPtSorted[numTauCands];
334 for(
int iTau=0; iTau<numTauCands; iTau++ ) tauPtSorted[iTau] = genParticles->at(tau_cands_index[iTau]).pt();
336 TMath::Sort(numTauCands,tauPtSorted,idxTau);
337 for(
int iTau=0; iTau<numTauCands; iTau++ ){
339 if( iTau>=maxNumTauCands_ )
continue;
341 const reco::Candidate & mcParticle = (*genParticles)[tau_cands_index[idxTau[iTau]]];
347 int iso = gRandom->Integer(4)%2;
350 if( eta>=9999 )
continue;
352 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
355 tauVec.push_back(tau);
361 int maxOtherTaus = 8;
362 int numCurrentEGs = int( egammaVec.size() );
363 int numCurrentTaus = int( tauVec.size() );
365 int numExtraEGs=0, numExtraTaus=0;
375 if( iEvent.
getByToken(genJetsToken, genJets) ){
376 for(reco::GenJetCollection::const_iterator genJet = genJets->begin(); genJet!=genJets->end(); ++genJet ){
379 sumEt += genJet->et();
382 if( genJet->pt()<jetEtThreshold_ )
continue;
385 if( nJet>=maxNumJetCands_ )
continue;
386 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
393 if( eta>=9999 )
continue;
398 jetVec.push_back(jet);
403 if( (numExtraEGs+numCurrentEGs)<maxNumEGCands_ && numExtraEGs<maxOtherEGs ){
408 int EGphi = convertPhiToHW( genJet->phi(),
PhiStepCalo_ );
411 int EGiso = gRandom->Integer(4)%2;
413 l1t::EGamma eg(*p4, EGpt, EGeta, EGphi, EGqual, EGiso);
414 egammaVec.push_back(eg);
417 if( (numExtraTaus+numCurrentTaus)<maxNumTauCands_ && numExtraTaus<maxOtherTaus ){
422 int Tauphi = convertPhiToHW( genJet->phi(),
PhiStepCalo_ );
424 int Tauiso = gRandom->Integer(4)%2;
426 l1t::Tau tau(*p4, Taupt, Taueta, Tauphi, Tauqual, Tauiso);
427 tauVec.push_back(tau);
433 LogTrace(
"l1t|Global") <<
">>> GenJets collection not found!" << std::endl;
443 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
446 l1t::EtSum etmiss(*p4, l1t::EtSum::EtSumType::kMissingEt,pt, 0,phi, 0);
447 etsumVec.push_back(etmiss);
450 pt = convertPtToHW( genMet->front().pt()*0.9,
MaxEt_,
PtStep_ );
451 phi = convertPhiToHW( genMet->front().phi()+ 3.14/5.,
PhiStepCalo_ );
453 l1t::EtSum htmiss(*p4, l1t::EtSum::EtSumType::kMissingHt,pt, 0,phi, 0);
454 etsumVec.push_back(htmiss);
459 LogTrace(
"l1t|Global") <<
">>> GenMet collection not found!" << std::endl;
464 int pt = convertPtToHW( sumEt, 2047, PtStep_ );
465 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *
p4 =
new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
466 l1t::EtSum etTotal(*p4, l1t::EtSum::EtSumType::kTotalEt,pt, 0, 0, 0);
467 etsumVec.push_back(etTotal);
469 pt = convertPtToHW( sumEt*0.9, 2047, PtStep_ );
470 l1t::EtSum htTotal(*p4, l1t::EtSum::EtSumType::kTotalHt,pt, 0, 0, 0);
471 etsumVec.push_back(htTotal);
475 printf(
"Event %i EmptyBxEvt %i emptyBxTrailer %i diff %i \n",eventCnt_,emptyBxEvt_,emptyBxTrailer_,(emptyBxEvt_ - eventCnt_));
478 for(
int iMu=0; iMu<int(muonVec_bxm2.size()); iMu++ ){
479 muons->push_back(-2, muonVec_bxm2[iMu]);
481 for(
int iMu=0; iMu<int(muonVec_bxm1.size()); iMu++ ){
482 muons->push_back(-1, muonVec_bxm1[iMu]);
484 for(
int iMu=0; iMu<int(muonVec_bx0.size()); iMu++ ){
485 muons->push_back(0, muonVec_bx0[iMu]);
487 for(
int iMu=0; iMu<int(muonVec_bxp1.size()); iMu++ ){
488 muons->push_back(1, muonVec_bxp1[iMu]);
490 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
491 for(
int iMu=0; iMu<int(muonVec.size()); iMu++ ){
492 muons->push_back(2, muonVec[iMu]);
500 for(
int iEG=0; iEG<int(egammaVec_bxm2.size()); iEG++ ){
501 egammas->push_back(-2, egammaVec_bxm2[iEG]);
503 for(
int iEG=0; iEG<int(egammaVec_bxm1.size()); iEG++ ){
504 egammas->push_back(-1, egammaVec_bxm1[iEG]);
506 for(
int iEG=0; iEG<int(egammaVec_bx0.size()); iEG++ ){
507 egammas->push_back(0, egammaVec_bx0[iEG]);
509 for(
int iEG=0; iEG<int(egammaVec_bxp1.size()); iEG++ ){
510 egammas->push_back(1, egammaVec_bxp1[iEG]);
512 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
513 for(
int iEG=0; iEG<int(egammaVec.size()); iEG++ ){
514 egammas->push_back(2, egammaVec[iEG]);
522 for(
int iTau=0; iTau<int(tauVec_bxm2.size()); iTau++ ){
523 taus->push_back(-2, tauVec_bxm2[iTau]);
525 for(
int iTau=0; iTau<int(tauVec_bxm1.size()); iTau++ ){
526 taus->push_back(-1, tauVec_bxm1[iTau]);
528 for(
int iTau=0; iTau<int(tauVec_bx0.size()); iTau++ ){
529 taus->push_back(0, tauVec_bx0[iTau]);
531 for(
int iTau=0; iTau<int(tauVec_bxp1.size()); iTau++ ){
532 taus->push_back(1, tauVec_bxp1[iTau]);
534 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
535 for(
int iTau=0; iTau<int(tauVec.size()); iTau++ ){
536 taus->push_back(2, tauVec[iTau]);
544 for(
int iJet=0; iJet<int(jetVec_bxm2.size()); iJet++ ){
545 jets->push_back(-2, jetVec_bxm2[iJet]);
547 for(
int iJet=0; iJet<int(jetVec_bxm1.size()); iJet++ ){
548 jets->push_back(-1, jetVec_bxm1[iJet]);
550 for(
int iJet=0; iJet<int(jetVec_bx0.size()); iJet++ ){
551 jets->push_back(0, jetVec_bx0[iJet]);
553 for(
int iJet=0; iJet<int(jetVec_bxp1.size()); iJet++ ){
554 jets->push_back(1, jetVec_bxp1[iJet]);
556 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
557 for(
int iJet=0; iJet<int(jetVec.size()); iJet++ ){
558 jets->push_back(2, jetVec[iJet]);
566 for(
int iETsum=0; iETsum<int(etsumVec_bxm2.size()); iETsum++ ){
567 etsums->push_back(-2, etsumVec_bxm2[iETsum]);
569 for(
int iETsum=0; iETsum<int(etsumVec_bxm1.size()); iETsum++ ){
570 etsums->push_back(-1, etsumVec_bxm1[iETsum]);
572 for(
int iETsum=0; iETsum<int(etsumVec_bx0.size()); iETsum++ ){
573 etsums->push_back(0, etsumVec_bx0[iETsum]);
575 for(
int iETsum=0; iETsum<int(etsumVec_bxp1.size()); iETsum++ ){
576 etsums->push_back(1, etsumVec_bxp1[iETsum]);
578 if(emptyBxTrailer_<=(emptyBxEvt_ - eventCnt_)) {
579 for(
int iETsum=0; iETsum<int(etsumVec.size()); iETsum++ ){
580 etsums->push_back(2, etsumVec[iETsum]);
595 muonVec_bxm2 = muonVec_bxm1;
596 egammaVec_bxm2 = egammaVec_bxm1;
597 tauVec_bxm2 = tauVec_bxm1;
598 jetVec_bxm2 = jetVec_bxm1;
599 etsumVec_bxm2 = etsumVec_bxm1;
601 muonVec_bxm1 = muonVec_bx0;
602 egammaVec_bxm1 = egammaVec_bx0;
603 tauVec_bxm1 = tauVec_bx0;
604 jetVec_bxm1 = jetVec_bx0;
605 etsumVec_bxm1 = etsumVec_bx0;
607 muonVec_bx0 = muonVec_bxp1;
608 egammaVec_bx0 = egammaVec_bxp1;
609 tauVec_bx0 = tauVec_bxp1;
610 jetVec_bx0 = jetVec_bxp1;
611 etsumVec_bx0 = etsumVec_bxp1;
613 muonVec_bxp1 = muonVec;
614 egammaVec_bxp1 = egammaVec;
615 tauVec_bxp1 = tauVec;
616 jetVec_bxp1 = jetVec;
617 etsumVec_bxp1 = etsumVec;
629 GenToInputProducer::endJob() {
636 LogDebug(
"l1t|Global") <<
"GenToInputProducer::beginRun function called...\n";
641 gRandom =
new TRandom3();
651 int GenToInputProducer::convertPhiToHW(
double iphi,
int steps){
653 double phiMax = 2 *
M_PI;
654 if( iphi < 0 ) iphi += 2*
M_PI;
655 if( iphi > phiMax) iphi -= phiMax;
657 int hwPhi = int( (iphi/phiMax)*steps + 0.00001 );
663 double binWidth = (maxEta -
minEta)/steps;
666 if(ieta < minEta)
return 99999;
667 if(ieta > maxEta)
return 99999;
669 int binNum = (int)(ieta/binWidth);
670 if(ieta<0.) binNum--;
678 int GenToInputProducer::convertPtToHW(
double ipt,
int maxPt,
double step){
680 int hwPt = int( ipt/step + 0.0001 );
682 if( hwPt > maxPt ) hwPt = maxPt;
T getParameter(std::string const &) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
virtual double pt() const =0
transverse momentum
#define DEFINE_FWK_MODULE(type)
virtual int status() const =0
status word
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.)
Geom::Phi< T > phi() const
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity