Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-14 02:53:34 

0001 ////
0002 /// \class l1t::GenToInputProducer
0003 ///
0004 /// Description: Create Input Collections for the GT from MC gen particles.  Allows testing of emulation.
0005 ///
0006 ///
0007 /// \author: D. Puigh OSU
0008 ///
0009 ///  Modeled after FakeInputProducer.cc
0010 ///
0011 /// \new features:  R. Cavanaugh
0012 ///          - added unconstrained pT and impact parameter for displaced muons in barrel
0013 ///          - added LLP DISP (quality) bit for displaced jets in calorimeter
0014 
0015 // system include files
0016 #include <memory>
0017 
0018 // user include files
0019 
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/one/EDProducer.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/EDGetToken.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 
0031 //#include <vector>
0032 #include "DataFormats/L1Trigger/interface/BXVector.h"
0033 
0034 #include "DataFormats/L1Trigger/interface/EGamma.h"
0035 #include "DataFormats/L1Trigger/interface/Muon.h"
0036 #include "DataFormats/L1Trigger/interface/MuonShower.h"
0037 #include "DataFormats/L1Trigger/interface/Tau.h"
0038 #include "DataFormats/L1Trigger/interface/Jet.h"
0039 #include "DataFormats/L1Trigger/interface/EtSum.h"
0040 #include "DataFormats/L1TGlobal/interface/GlobalExtBlk.h"
0041 
0042 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0043 #include "DataFormats/JetReco/interface/GenJet.h"
0044 #include "DataFormats/METReco/interface/GenMETCollection.h"
0045 #include "DataFormats/METReco/interface/GenMET.h"
0046 
0047 #include "TMath.h"
0048 #include "TRandom3.h"
0049 #include <cstdlib>
0050 
0051 using namespace std;
0052 using namespace edm;
0053 
0054 #ifndef M_PI
0055 #define M_PI 3.14159265358979323846
0056 #endif
0057 
0058 namespace l1t {
0059 
0060   //
0061   // class declaration
0062   //
0063 
0064   class GenToInputProducer : public one::EDProducer<edm::one::WatchRuns> {
0065   public:
0066     explicit GenToInputProducer(const ParameterSet&);
0067     ~GenToInputProducer() override;
0068 
0069     static void fillDescriptions(ConfigurationDescriptions& descriptions);
0070 
0071   private:
0072     void produce(Event&, EventSetup const&) override;
0073     void beginJob() override;
0074     void endJob() override;
0075     void beginRun(Run const& iR, EventSetup const& iE) override;
0076     void endRun(Run const& iR, EventSetup const& iE) override;
0077 
0078     int convertPhiToHW(double iphi, int steps);
0079     int convertEtaToHW(double ieta, double minEta, double maxEta, int steps);
0080     int convertPtToHW(double ipt, int maxPt, double step);
0081 
0082     // ----------member data ---------------------------
0083     unsigned long long m_paramsCacheId;  // Cache-ID from current parameters, to check if needs to be updated.
0084     //std::shared_ptr<const CaloParams> m_dbpars; // Database parameters for the trigger, to be updated as needed.
0085     //std::shared_ptr<const FirmwareVersion> m_fwv;
0086     //std::shared_ptr<FirmwareVersion> m_fwv; //not const during testing.
0087 
0088     TRandom3* gRandom;
0089     TRandom3* mushowerRandom;
0090 
0091     // BX parameters
0092     int bxFirst_;
0093     int bxLast_;
0094 
0095     int maxNumMuCands_;
0096     int maxNumMuShowerCands_;
0097     int maxNumJetCands_;
0098     int maxNumEGCands_;
0099     int maxNumTauCands_;
0100 
0101     double jetEtThreshold_;
0102     double tauEtThreshold_;
0103     double egEtThreshold_;
0104     double muEtThreshold_;
0105 
0106     // Control how to end the job
0107     int emptyBxTrailer_;
0108     int emptyBxEvt_;
0109     int eventCnt_;
0110 
0111     // Tokens
0112     edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken;
0113     edm::EDGetTokenT<reco::GenJetCollection> genJetsToken;
0114     edm::EDGetTokenT<reco::GenMETCollection> genMetToken;
0115 
0116     std::vector<l1t::Muon> muonVec_bxm2;
0117     std::vector<l1t::Muon> muonVec_bxm1;
0118     std::vector<l1t::Muon> muonVec_bx0;
0119     std::vector<l1t::Muon> muonVec_bxp1;
0120 
0121     std::vector<l1t::MuonShower> muonShowerVec_bxm2;
0122     std::vector<l1t::MuonShower> muonShowerVec_bxm1;
0123     std::vector<l1t::MuonShower> muonShowerVec_bx0;
0124     std::vector<l1t::MuonShower> muonShowerVec_bxp1;
0125 
0126     std::vector<l1t::EGamma> egammaVec_bxm2;
0127     std::vector<l1t::EGamma> egammaVec_bxm1;
0128     std::vector<l1t::EGamma> egammaVec_bx0;
0129     std::vector<l1t::EGamma> egammaVec_bxp1;
0130 
0131     std::vector<l1t::Tau> tauVec_bxm2;
0132     std::vector<l1t::Tau> tauVec_bxm1;
0133     std::vector<l1t::Tau> tauVec_bx0;
0134     std::vector<l1t::Tau> tauVec_bxp1;
0135 
0136     std::vector<l1t::Jet> jetVec_bxm2;
0137     std::vector<l1t::Jet> jetVec_bxm1;
0138     std::vector<l1t::Jet> jetVec_bx0;
0139     std::vector<l1t::Jet> jetVec_bxp1;
0140 
0141     std::vector<l1t::EtSum> etsumVec_bxm2;
0142     std::vector<l1t::EtSum> etsumVec_bxm1;
0143     std::vector<l1t::EtSum> etsumVec_bx0;
0144     std::vector<l1t::EtSum> etsumVec_bxp1;
0145 
0146     GlobalExtBlk extCond_bxm2;
0147     GlobalExtBlk extCond_bxm1;
0148     GlobalExtBlk extCond_bx0;
0149     GlobalExtBlk extCond_bxp1;
0150   };
0151 
0152   //
0153   // constructors and destructor
0154   //
0155   GenToInputProducer::GenToInputProducer(const ParameterSet& iConfig) {
0156     // register what you produce
0157     produces<BXVector<l1t::EGamma>>();
0158     produces<BXVector<l1t::Muon>>();
0159     produces<BXVector<l1t::MuonShower>>();
0160     produces<BXVector<l1t::Tau>>();
0161     produces<BXVector<l1t::Jet>>();
0162     produces<BXVector<l1t::EtSum>>();
0163     produces<GlobalExtBlkBxCollection>();
0164 
0165     // Setup parameters
0166     bxFirst_ = iConfig.getParameter<int>("bxFirst");
0167     bxLast_ = iConfig.getParameter<int>("bxLast");
0168 
0169     maxNumMuCands_ = iConfig.getParameter<int>("maxMuCand");
0170     maxNumMuShowerCands_ = iConfig.getParameter<int>("maxMuShowerCand");
0171     maxNumJetCands_ = iConfig.getParameter<int>("maxJetCand");
0172     maxNumEGCands_ = iConfig.getParameter<int>("maxEGCand");
0173     maxNumTauCands_ = iConfig.getParameter<int>("maxTauCand");
0174 
0175     jetEtThreshold_ = iConfig.getParameter<double>("jetEtThreshold");
0176     tauEtThreshold_ = iConfig.getParameter<double>("tauEtThreshold");
0177     egEtThreshold_ = iConfig.getParameter<double>("egEtThreshold");
0178     muEtThreshold_ = iConfig.getParameter<double>("muEtThreshold");
0179 
0180     emptyBxTrailer_ = iConfig.getParameter<int>("emptyBxTrailer");
0181     emptyBxEvt_ = iConfig.getParameter<int>("emptyBxEvt");
0182 
0183     genParticlesToken = consumes<reco::GenParticleCollection>(std::string("genParticles"));
0184     genJetsToken = consumes<reco::GenJetCollection>(std::string("ak4GenJets"));
0185     genMetToken = consumes<reco::GenMETCollection>(std::string("genMetCalo"));
0186 
0187     // set cache id to zero, will be set at first beginRun:
0188     m_paramsCacheId = 0;
0189     eventCnt_ = 0;
0190   }
0191 
0192   GenToInputProducer::~GenToInputProducer() {}
0193 
0194   //
0195   // member functions
0196   //
0197 
0198   // ------------ method called to produce the data ------------
0199   void GenToInputProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0200     eventCnt_++;
0201 
0202     LogDebug("GtGenToInputProducer") << "GenToInputProducer::produce function called...\n";
0203 
0204     // Setup vectors
0205     std::vector<l1t::Muon> muonVec;
0206     std::vector<l1t::MuonShower> muonShowerVec;
0207     std::vector<l1t::EGamma> egammaVec;
0208     std::vector<l1t::Tau> tauVec;
0209     std::vector<l1t::Jet> jetVec;
0210     std::vector<l1t::EtSum> etsumVec;
0211     GlobalExtBlk extCond_bx;
0212 
0213     // Set the range of BX....TO DO...move to Params or determine from param set.
0214     int bxFirst = bxFirst_;
0215     int bxLast = bxLast_;
0216 
0217     // Default values objects
0218     double MaxLepPt_ = 255;
0219     double MaxJetPt_ = 1023;
0220     double MaxEt_ = 2047;
0221 
0222     double MaxCaloEta_ = 5.0;
0223     double MaxMuonEta_ = 2.45;
0224 
0225     double PhiStepCalo_ = 144;
0226     double PhiStepMuon_ = 576;
0227 
0228     // eta scale
0229     double EtaStepCalo_ = 230;
0230     double EtaStepMuon_ = 450;
0231 
0232     // Et scale (in GeV)
0233     double PtStep_ = 0.5;
0234 
0235     //outputs
0236     std::unique_ptr<l1t::EGammaBxCollection> egammas(new l1t::EGammaBxCollection(0, bxFirst, bxLast));
0237     std::unique_ptr<l1t::MuonBxCollection> muons(new l1t::MuonBxCollection(0, bxFirst, bxLast));
0238     std::unique_ptr<l1t::MuonShowerBxCollection> muonShowers(new l1t::MuonShowerBxCollection(0, bxFirst, bxLast));
0239     std::unique_ptr<l1t::TauBxCollection> taus(new l1t::TauBxCollection(0, bxFirst, bxLast));
0240     std::unique_ptr<l1t::JetBxCollection> jets(new l1t::JetBxCollection(0, bxFirst, bxLast));
0241     std::unique_ptr<l1t::EtSumBxCollection> etsums(new l1t::EtSumBxCollection(0, bxFirst, bxLast));
0242     std::unique_ptr<GlobalExtBlkBxCollection> extCond(new GlobalExtBlkBxCollection(0, bxFirst, bxLast));
0243 
0244     std::vector<int> mu_cands_index;
0245     std::vector<int> eg_cands_index;
0246     std::vector<int> tau_cands_index;
0247     edm::Handle<reco::GenParticleCollection> genParticles;
0248     // Make sure that you can get genParticles
0249     if (iEvent.getByToken(genParticlesToken, genParticles)) {
0250       for (size_t k = 0; k < genParticles->size(); k++) {
0251         const reco::Candidate& mcParticle = (*genParticles)[k];
0252 
0253         int status = mcParticle.status();
0254         int pdgId = mcParticle.pdgId();
0255         double pt = mcParticle.pt();
0256 
0257         // Only use status 1 particles  (Tau's need to be allowed through..take status 2 taus)
0258         if (status != 1 && !(abs(pdgId) == 15 && status == 2))
0259           continue;
0260 
0261         int absId = abs(pdgId);
0262 
0263         if (absId == 11 && pt >= egEtThreshold_)
0264           eg_cands_index.push_back(k);
0265         else if (absId == 13 && pt >= muEtThreshold_)
0266           mu_cands_index.push_back(k);
0267         else if (absId == 15 && pt >= tauEtThreshold_)
0268           tau_cands_index.push_back(k);
0269       }
0270     } else {
0271       LogTrace("GtGenToInputProducer") << ">>> GenParticles collection not found!" << std::endl;
0272     }
0273 
0274     // Muon Shower Collection
0275 
0276     bool mus0 = (bool)mushowerRandom->Integer(2);  // should be [0,1] = 1 bit;
0277     bool mus1 = (bool)mushowerRandom->Integer(2);
0278     bool mus2 = (bool)mushowerRandom->Integer(2);
0279     bool mus2Loose = (bool)mushowerRandom->Integer(2);
0280     bool musoot0 = (bool)mushowerRandom->Integer(2);
0281     bool musoot1 = (bool)mushowerRandom->Integer(2);
0282     bool musoot2Loose = (bool)mushowerRandom->Integer(2);
0283     //fill a vector of MuonShower objs with only one obj per BX, not one obj per muon obj
0284     cout << "GenToInputProducer MuonShower = (MUS0, MUS1, MUS2, MUSOOT0, MUSOOT1) = (" << mus0 << "," << mus1 << ","
0285          << mus2 << "," << musoot0 << "," << musoot1 << ")" << endl;
0286     l1t::MuonShower muShower(mus0, musoot0, mus2Loose, musoot2Loose, mus1, musoot1, mus2);
0287     muShower.setMusOutOfTime0(musoot0);
0288     muShower.setMusOutOfTime1(musoot1);
0289     muonShowerVec.push_back(muShower);
0290 
0291     // Muon Collection
0292     int numMuCands = int(mu_cands_index.size());
0293     Int_t idxMu[numMuCands];
0294     double muPtSorted[numMuCands];
0295     for (int iMu = 0; iMu < numMuCands; iMu++)
0296       muPtSorted[iMu] = genParticles->at(mu_cands_index[iMu]).pt();
0297 
0298     TMath::Sort(numMuCands, muPtSorted, idxMu);
0299     for (int iMu = 0; iMu < numMuCands; iMu++) {
0300       if (iMu >= maxNumMuCands_)
0301         continue;
0302 
0303       const reco::Candidate& mcParticle = (*genParticles)[mu_cands_index[idxMu[iMu]]];
0304 
0305       int pt = convertPtToHW(mcParticle.pt(), MaxLepPt_, PtStep_);
0306       int eta = convertEtaToHW(mcParticle.eta(), -MaxMuonEta_, MaxMuonEta_, EtaStepMuon_);
0307       int phi = convertPhiToHW(mcParticle.phi(), PhiStepMuon_);
0308       int qual = gRandom->Integer(16);    //4;
0309       int iso = gRandom->Integer(4) % 2;  //1;
0310       int charge = (mcParticle.charge() < 0) ? 1 : 0;
0311       int chargeValid = 1;
0312       int tfMuIdx = 0;
0313       int tag = 1;
0314       bool debug = false;
0315       int isoSum = 0;
0316       int dPhi = 0;
0317       int dEta = 0;
0318       int rank = 0;
0319       int hwEtaAtVtx = eta;
0320       int hwPhiAtVtx = phi;
0321       double etaAtVtx = 0.0;
0322       double phiAtVtx = 0.0;
0323       int hwPtUnconstrained =
0324           convertPtToHW(mcParticle.pt(), MaxLepPt_, PtStep_) / 2;  // word is 8 bits wide so divide 9 bit word by 2
0325       double ptUnconstrained = 0.0;
0326       int dXY = gRandom->Integer(4);  // should be [0,3] = 2 bits
0327 
0328       // Eta outside of acceptance
0329       if (eta >= 9999)
0330         continue;
0331 
0332       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>* p4 =
0333           new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>();
0334 
0335       l1t::Muon mu(*p4,
0336                    pt,
0337                    eta,
0338                    phi,
0339                    qual,
0340                    charge,
0341                    chargeValid,
0342                    iso,
0343                    tfMuIdx,
0344                    tag,
0345                    debug,
0346                    isoSum,
0347                    dPhi,
0348                    dEta,
0349                    rank,
0350                    hwEtaAtVtx,
0351                    hwPhiAtVtx,
0352                    etaAtVtx,
0353                    phiAtVtx,
0354                    hwPtUnconstrained,
0355                    ptUnconstrained,
0356                    dXY);  // modified to conform to latest Muon.h interface
0357       muonVec.push_back(mu);
0358     }
0359 
0360     // EG Collection
0361     int numEgCands = int(eg_cands_index.size());
0362     Int_t idxEg[numEgCands];
0363     double egPtSorted[numEgCands];
0364     for (int iEg = 0; iEg < numEgCands; iEg++)
0365       egPtSorted[iEg] = genParticles->at(eg_cands_index[iEg]).pt();
0366 
0367     TMath::Sort(numEgCands, egPtSorted, idxEg);
0368     for (int iEg = 0; iEg < numEgCands; iEg++) {
0369       if (iEg >= maxNumEGCands_)
0370         continue;
0371 
0372       const reco::Candidate& mcParticle = (*genParticles)[eg_cands_index[idxEg[iEg]]];
0373 
0374       int pt = convertPtToHW(mcParticle.pt(), MaxLepPt_, PtStep_);
0375       int eta = convertEtaToHW(mcParticle.eta(), -MaxCaloEta_, MaxCaloEta_, EtaStepCalo_);
0376       int phi = convertPhiToHW(mcParticle.phi(), PhiStepCalo_);
0377       int qual = gRandom->Integer(2);  // modified for LLP Jets
0378       int iso = gRandom->Integer(4) % 2;
0379 
0380       // Eta outside of acceptance
0381       if (eta >= 9999)
0382         continue;
0383 
0384       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>* p4 =
0385           new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>();
0386 
0387       l1t::EGamma eg(*p4, pt, eta, phi, qual, iso);
0388       egammaVec.push_back(eg);
0389     }
0390 
0391     // Tau Collection
0392     int numTauCands = int(tau_cands_index.size());
0393     Int_t idxTau[numTauCands];
0394     double tauPtSorted[numTauCands];
0395     for (int iTau = 0; iTau < numTauCands; iTau++)
0396       tauPtSorted[iTau] = genParticles->at(tau_cands_index[iTau]).pt();
0397 
0398     TMath::Sort(numTauCands, tauPtSorted, idxTau);
0399     for (int iTau = 0; iTau < numTauCands; iTau++) {
0400       if (iTau >= maxNumTauCands_)
0401         continue;
0402 
0403       const reco::Candidate& mcParticle = (*genParticles)[tau_cands_index[idxTau[iTau]]];
0404 
0405       int pt = convertPtToHW(mcParticle.pt(), MaxLepPt_, PtStep_);
0406       int eta = convertEtaToHW(mcParticle.eta(), -MaxCaloEta_, MaxCaloEta_, EtaStepCalo_);
0407       int phi = convertPhiToHW(mcParticle.phi(), PhiStepCalo_);
0408       int qual = gRandom->Integer(2);  // modified for LLP Jets
0409       int iso = gRandom->Integer(4) % 2;
0410 
0411       // Eta outside of acceptance
0412       if (eta >= 9999)
0413         continue;
0414 
0415       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>* p4 =
0416           new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>();
0417 
0418       l1t::Tau tau(*p4, pt, eta, phi, qual, iso);
0419       tauVec.push_back(tau);
0420     }
0421 
0422     // Temporary hack to increase number of EGs and taus
0423     int maxOtherEGs = 4;
0424     int maxOtherTaus = 8;
0425     int numCurrentEGs = int(egammaVec.size());
0426     int numCurrentTaus = int(tauVec.size());
0427 
0428     int numExtraEGs = 0, numExtraTaus = 0;
0429     // end hack
0430 
0431     // Use to sum the energy of the objects in the event for ETT and HTT
0432     // sum all jets
0433     double sumEt = 0;
0434 
0435     int nJet = 0;
0436     edm::Handle<reco::GenJetCollection> genJets;
0437     // Make sure that you can get genJets
0438     if (iEvent.getByToken(genJetsToken, genJets)) {  // Jet Collection
0439       for (reco::GenJetCollection::const_iterator genJet = genJets->begin(); genJet != genJets->end(); ++genJet) {
0440         //Keep running sum of total Et
0441         sumEt += genJet->et();
0442 
0443         // Apply pt and eta cut?
0444         if (genJet->pt() < jetEtThreshold_)
0445           continue;
0446 
0447         //
0448         if (nJet >= maxNumJetCands_)
0449           continue;
0450         ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>* p4 =
0451             new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>();
0452 
0453         int pt = convertPtToHW(genJet->et(), MaxJetPt_, PtStep_);
0454         int eta = convertEtaToHW(genJet->eta(), -MaxCaloEta_, MaxCaloEta_, EtaStepCalo_);
0455         int phi = convertPhiToHW(genJet->phi(), PhiStepCalo_);
0456 
0457         // Eta outside of acceptance
0458         if (eta >= 9999)
0459           continue;
0460 
0461         int qual = gRandom->Integer(2);  // modified for LLP Jets
0462 
0463         l1t::Jet jet(*p4, pt, eta, phi, qual);
0464         jetVec.push_back(jet);
0465 
0466         nJet++;
0467 
0468         // Temporary hack to increase number of EGs and taus
0469         if ((numExtraEGs + numCurrentEGs) < maxNumEGCands_ && numExtraEGs < maxOtherEGs) {
0470           numExtraEGs++;
0471 
0472           int EGpt = convertPtToHW(genJet->et(), MaxLepPt_, PtStep_);
0473           int EGeta = convertEtaToHW(genJet->eta(), -MaxCaloEta_, MaxCaloEta_, EtaStepCalo_);
0474           int EGphi = convertPhiToHW(genJet->phi(), PhiStepCalo_);
0475 
0476           int EGqual = gRandom->Integer(2);  // modified for LLP Jets
0477           int EGiso = gRandom->Integer(4) % 2;
0478 
0479           l1t::EGamma eg(*p4, EGpt, EGeta, EGphi, EGqual, EGiso);
0480           egammaVec.push_back(eg);
0481         }
0482 
0483         if ((numExtraTaus + numCurrentTaus) < maxNumTauCands_ && numExtraTaus < maxOtherTaus) {
0484           numExtraTaus++;
0485 
0486           int Taupt = convertPtToHW(genJet->et(), MaxLepPt_, PtStep_);
0487           int Taueta = convertEtaToHW(genJet->eta(), -MaxCaloEta_, MaxCaloEta_, EtaStepCalo_);
0488           int Tauphi = convertPhiToHW(genJet->phi(), PhiStepCalo_);
0489           int Tauqual = gRandom->Integer(2);  // modified for LLP Jets
0490           int Tauiso = gRandom->Integer(4) % 2;
0491 
0492           l1t::Tau tau(*p4, Taupt, Taueta, Tauphi, Tauqual, Tauiso);
0493           tauVec.push_back(tau);
0494         }
0495         // end hack
0496       }
0497     } else {
0498       LogTrace("GtGenToInputProducer") << ">>> GenJets collection not found!" << std::endl;
0499     }
0500 
0501     // Put the total Et into EtSums  (Make HTT slightly smaller to tell them apart....not supposed to be realistic)
0502     int pt = convertPtToHW(sumEt, 2047, PtStep_);
0503     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>* p4 =
0504         new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>();
0505     l1t::EtSum etTotal(*p4, l1t::EtSum::EtSumType::kTotalEt, pt, 0, 0, 0);
0506 
0507     // Scale down ETTem as an estimate
0508     pt = convertPtToHW(sumEt * 0.6, 2047, PtStep_);
0509     l1t::EtSum etEmTotal(*p4, l1t::EtSum::EtSumType::kTotalEtEm, pt, 0, 0, 0);
0510 
0511     //ccla Generate uniform distribution of tower counts
0512     int nTowers = 4095 * gRandom->Rndm();
0513     l1t::EtSum towerCounts(*p4, l1t::EtSum::EtSumType::kTowerCount, nTowers, 0, 0, 0);
0514 
0515     //ccla Generate uniform distributions of AsymEt, AsymHt, AsymEtHF, AsymHtF
0516     int nAsymEt = 255 * gRandom->Rndm();
0517     l1t::EtSum AsymEt(*p4, l1t::EtSum::EtSumType::kAsymEt, nAsymEt, 0, 0, 0);
0518     int nAsymHt = 255 * gRandom->Rndm();
0519     l1t::EtSum AsymHt(*p4, l1t::EtSum::EtSumType::kAsymHt, nAsymHt, 0, 0, 0);
0520     int nAsymEtHF = 255 * gRandom->Rndm();
0521     l1t::EtSum AsymEtHF(*p4, l1t::EtSum::EtSumType::kAsymEtHF, nAsymEtHF, 0, 0, 0);
0522     int nAsymHtHF = 255 * gRandom->Rndm();
0523     l1t::EtSum AsymHtHF(*p4, l1t::EtSum::EtSumType::kAsymHtHF, nAsymHtHF, 0, 0, 0);
0524 
0525     pt = convertPtToHW(sumEt * 0.9, 2047, PtStep_);
0526     l1t::EtSum htTotal(*p4, l1t::EtSum::EtSumType::kTotalHt, pt, 0, 0, 0);
0527 
0528     // Add EtSums for testing the MinBias Trigger (use some random numbers)
0529     int hfP0val = gRandom->Poisson(4.);
0530     if (hfP0val > 15)
0531       hfP0val = 15;
0532     l1t::EtSum hfP0(*p4, l1t::EtSum::EtSumType::kMinBiasHFP0, hfP0val, 0, 0, 0);
0533 
0534     int hfM0val = gRandom->Poisson(4.);
0535     if (hfM0val > 15)
0536       hfM0val = 15;
0537     l1t::EtSum hfM0(*p4, l1t::EtSum::EtSumType::kMinBiasHFM0, hfM0val, 0, 0, 0);
0538 
0539     int hfP1val = gRandom->Poisson(4.);
0540     if (hfP1val > 15)
0541       hfP1val = 15;
0542     l1t::EtSum hfP1(*p4, l1t::EtSum::EtSumType::kMinBiasHFP1, hfP1val, 0, 0, 0);
0543 
0544     int hfM1val = gRandom->Poisson(4.);
0545     if (hfM1val > 15)
0546       hfM1val = 15;
0547     l1t::EtSum hfM1(*p4, l1t::EtSum::EtSumType::kMinBiasHFM1, hfM1val, 0, 0, 0);
0548 
0549     // Do same for Centrality
0550     int cent30val(0), cent74val(0);
0551     int centa = gRandom->Poisson(2.);
0552     int centb = gRandom->Poisson(2.);
0553     if (centa >= centb) {
0554       cent30val = centa;
0555       cent74val = centb;
0556     } else {
0557       cent30val = centb;
0558       cent74val = centa;
0559     }
0560 
0561     if (cent30val > 15)
0562       cent30val = 15;
0563     if (cent74val > 15)
0564       cent74val = 15;
0565 
0566     int shift = 4;
0567     int centralval = 0;
0568     centralval |= cent30val & 0xF;
0569     centralval |= (cent74val & 0xF) << shift;
0570 
0571     l1t::EtSum centrality(*p4, l1t::EtSum::EtSumType::kCentrality, centralval, 0, 0, 0);
0572 
0573     int mpt = 0;
0574     int mphi = 0;
0575     int mptHf = 0;
0576     int mphiHf = 0;
0577     int mhpt = 0;
0578     int mhphi = 0;
0579     int mhptHf = 0;
0580     int mhphiHf = 0;
0581 
0582     edm::Handle<reco::GenMETCollection> genMet;
0583     // Make sure that you can get genMET
0584     if (iEvent.getByToken(genMetToken, genMet)) {
0585       mpt = convertPtToHW(genMet->front().pt(), MaxEt_, PtStep_);
0586       mphi = convertPhiToHW(genMet->front().phi(), PhiStepCalo_);
0587 
0588       // Make Missing Et with HF slightly largeer and rotated (These are all fake inputs anyway...not supposed to be realistic)
0589       mptHf = convertPtToHW(genMet->front().pt() * 1.1, MaxEt_, PtStep_);
0590       mphiHf = convertPhiToHW(genMet->front().phi() + 3.14 / 7., PhiStepCalo_);
0591 
0592       // Make Missing Ht slightly smaller and rotated (These are all fake inputs anyway...not supposed to be realistic)
0593       mhpt = convertPtToHW(genMet->front().pt() * 0.9, MaxEt_, PtStep_);
0594       mhphi = convertPhiToHW(genMet->front().phi() + 3.14 / 5., PhiStepCalo_);
0595 
0596       // Ditto with Hissing Ht with HF
0597       mhptHf = convertPtToHW(genMet->front().pt() * 0.95, MaxEt_, PtStep_);
0598       mhphiHf = convertPhiToHW(genMet->front().phi() + 3.14 / 6., PhiStepCalo_);
0599     } else {
0600       LogTrace("GtGenToInputProducer") << ">>> GenMet collection not found!" << std::endl;
0601     }
0602 
0603     // Missing Et and missing htt
0604     l1t::EtSum etmiss(*p4, l1t::EtSum::EtSumType::kMissingEt, mpt, 0, mphi, 0);
0605     l1t::EtSum etmissHF(*p4, l1t::EtSum::EtSumType::kMissingEtHF, mptHf, 0, mphiHf, 0);
0606     l1t::EtSum htmiss(*p4, l1t::EtSum::EtSumType::kMissingHt, mhpt, 0, mhphi, 0);
0607     l1t::EtSum htmissHF(*p4, l1t::EtSum::EtSumType::kMissingHtHF, mhptHf, 0, mhphiHf, 0);
0608 
0609     // Fill the EtSums in the Correct order
0610     etsumVec.push_back(etTotal);
0611     etsumVec.push_back(etEmTotal);
0612     etsumVec.push_back(hfP0);  // Frame0
0613 
0614     etsumVec.push_back(htTotal);
0615     etsumVec.push_back(towerCounts);
0616     etsumVec.push_back(hfM0);  //Frame1
0617 
0618     etsumVec.push_back(etmiss);
0619     etsumVec.push_back(AsymEt);
0620     etsumVec.push_back(hfP1);  //Frame2
0621 
0622     etsumVec.push_back(htmiss);
0623     etsumVec.push_back(AsymHt);
0624     etsumVec.push_back(hfM1);  //Frame3
0625 
0626     etsumVec.push_back(etmissHF);
0627     etsumVec.push_back(AsymEtHF);  // Frame4
0628 
0629     etsumVec.push_back(htmissHF);
0630     etsumVec.push_back(AsymHtHF);
0631     etsumVec.push_back(centrality);  // Frame5
0632 
0633     // Fill in some external conditions for testing
0634     if ((iEvent.id().event()) % 2 == 0) {
0635       for (int i = 0; i < 255; i = i + 2)
0636         extCond_bx.setExternalDecision(i, true);
0637     } else {
0638       for (int i = 1; i < 255; i = i + 2)
0639         extCond_bx.setExternalDecision(i, true);
0640     }
0641 
0642     // Insert all the bx into the L1 Collections
0643     //printf("Event %i  EmptyBxEvt %i emptyBxTrailer %i diff %i \n",eventCnt_,emptyBxEvt_,emptyBxTrailer_,(emptyBxEvt_ - eventCnt_));
0644 
0645     // Fill Muons
0646     for (int iMu = 0; iMu < int(muonVec_bxm2.size()); iMu++) {
0647       muons->push_back(-2, muonVec_bxm2[iMu]);
0648     }
0649     for (int iMu = 0; iMu < int(muonVec_bxm1.size()); iMu++) {
0650       muons->push_back(-1, muonVec_bxm1[iMu]);
0651     }
0652     for (int iMu = 0; iMu < int(muonVec_bx0.size()); iMu++) {
0653       muons->push_back(0, muonVec_bx0[iMu]);
0654     }
0655     for (int iMu = 0; iMu < int(muonVec_bxp1.size()); iMu++) {
0656       muons->push_back(1, muonVec_bxp1[iMu]);
0657     }
0658     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0659       for (int iMu = 0; iMu < int(muonVec.size()); iMu++) {
0660         muons->push_back(2, muonVec[iMu]);
0661       }
0662     } else {
0663       // this event is part of empty trailer...clear out data
0664       muonVec.clear();
0665     }
0666 
0667     // Fill MuonShowers
0668     for (int iMuShower = 0; iMuShower < int(muonShowerVec_bxm2.size()); iMuShower++) {
0669       muonShowers->push_back(-2, muonShowerVec_bxm2[iMuShower]);
0670     }
0671     for (int iMuShower = 0; iMuShower < int(muonShowerVec_bxm1.size()); iMuShower++) {
0672       muonShowers->push_back(-1, muonShowerVec_bxm1[iMuShower]);
0673     }
0674     for (int iMuShower = 0; iMuShower < int(muonShowerVec_bx0.size()); iMuShower++) {
0675       muonShowers->push_back(0, muonShowerVec_bx0[iMuShower]);
0676     }
0677     for (int iMuShower = 0; iMuShower < int(muonShowerVec_bxp1.size()); iMuShower++) {
0678       muonShowers->push_back(1, muonShowerVec_bxp1[iMuShower]);
0679     }
0680     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0681       for (int iMuShower = 0; iMuShower < int(muonShowerVec.size()); iMuShower++) {
0682         muonShowers->push_back(2, muonShowerVec[iMuShower]);
0683       }
0684     } else {
0685       // this event is part of empty trailer...clear out data
0686       muonShowerVec.clear();
0687     }
0688 
0689     // Fill Egammas
0690     for (int iEG = 0; iEG < int(egammaVec_bxm2.size()); iEG++) {
0691       egammas->push_back(-2, egammaVec_bxm2[iEG]);
0692     }
0693     for (int iEG = 0; iEG < int(egammaVec_bxm1.size()); iEG++) {
0694       egammas->push_back(-1, egammaVec_bxm1[iEG]);
0695     }
0696     for (int iEG = 0; iEG < int(egammaVec_bx0.size()); iEG++) {
0697       egammas->push_back(0, egammaVec_bx0[iEG]);
0698     }
0699     for (int iEG = 0; iEG < int(egammaVec_bxp1.size()); iEG++) {
0700       egammas->push_back(1, egammaVec_bxp1[iEG]);
0701     }
0702     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0703       for (int iEG = 0; iEG < int(egammaVec.size()); iEG++) {
0704         egammas->push_back(2, egammaVec[iEG]);
0705       }
0706     } else {
0707       // this event is part of empty trailer...clear out data
0708       egammaVec.clear();
0709     }
0710 
0711     // Fill Taus
0712     for (int iTau = 0; iTau < int(tauVec_bxm2.size()); iTau++) {
0713       taus->push_back(-2, tauVec_bxm2[iTau]);
0714     }
0715     for (int iTau = 0; iTau < int(tauVec_bxm1.size()); iTau++) {
0716       taus->push_back(-1, tauVec_bxm1[iTau]);
0717     }
0718     for (int iTau = 0; iTau < int(tauVec_bx0.size()); iTau++) {
0719       taus->push_back(0, tauVec_bx0[iTau]);
0720     }
0721     for (int iTau = 0; iTau < int(tauVec_bxp1.size()); iTau++) {
0722       taus->push_back(1, tauVec_bxp1[iTau]);
0723     }
0724     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0725       for (int iTau = 0; iTau < int(tauVec.size()); iTau++) {
0726         taus->push_back(2, tauVec[iTau]);
0727       }
0728     } else {
0729       // this event is part of empty trailer...clear out data
0730       tauVec.clear();
0731     }
0732 
0733     // Fill Jets
0734     for (int iJet = 0; iJet < int(jetVec_bxm2.size()); iJet++) {
0735       jets->push_back(-2, jetVec_bxm2[iJet]);
0736     }
0737     for (int iJet = 0; iJet < int(jetVec_bxm1.size()); iJet++) {
0738       jets->push_back(-1, jetVec_bxm1[iJet]);
0739     }
0740     for (int iJet = 0; iJet < int(jetVec_bx0.size()); iJet++) {
0741       jets->push_back(0, jetVec_bx0[iJet]);
0742     }
0743     for (int iJet = 0; iJet < int(jetVec_bxp1.size()); iJet++) {
0744       jets->push_back(1, jetVec_bxp1[iJet]);
0745     }
0746     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0747       for (int iJet = 0; iJet < int(jetVec.size()); iJet++) {
0748         jets->push_back(2, jetVec[iJet]);
0749       }
0750     } else {
0751       // this event is part of empty trailer...clear out data
0752       jetVec.clear();
0753     }
0754 
0755     // Fill Etsums
0756     for (int iETsum = 0; iETsum < int(etsumVec_bxm2.size()); iETsum++) {
0757       etsums->push_back(-2, etsumVec_bxm2[iETsum]);
0758     }
0759     for (int iETsum = 0; iETsum < int(etsumVec_bxm1.size()); iETsum++) {
0760       etsums->push_back(-1, etsumVec_bxm1[iETsum]);
0761     }
0762     for (int iETsum = 0; iETsum < int(etsumVec_bx0.size()); iETsum++) {
0763       etsums->push_back(0, etsumVec_bx0[iETsum]);
0764     }
0765     for (int iETsum = 0; iETsum < int(etsumVec_bxp1.size()); iETsum++) {
0766       etsums->push_back(1, etsumVec_bxp1[iETsum]);
0767     }
0768     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0769       for (int iETsum = 0; iETsum < int(etsumVec.size()); iETsum++) {
0770         etsums->push_back(2, etsumVec[iETsum]);
0771       }
0772     } else {
0773       // this event is part of empty trailer...clear out data
0774       etsumVec.clear();
0775     }
0776 
0777     // Fill Externals
0778     extCond->push_back(-2, extCond_bxm2);
0779     extCond->push_back(-1, extCond_bxm1);
0780     extCond->push_back(0, extCond_bx0);
0781     extCond->push_back(1, extCond_bxp1);
0782     if (emptyBxTrailer_ <= (emptyBxEvt_ - eventCnt_)) {
0783       extCond->push_back(2, extCond_bx);
0784     } else {
0785       // this event is part of the empty trailer...clear out data
0786       extCond_bx.reset();
0787     }
0788 
0789     iEvent.put(std::move(egammas));
0790     iEvent.put(std::move(muons));
0791     iEvent.put(std::move(muonShowers));
0792     iEvent.put(std::move(taus));
0793     iEvent.put(std::move(jets));
0794     iEvent.put(std::move(etsums));
0795     iEvent.put(std::move(extCond));
0796 
0797     // Now shift the bx data by one to prepare for next event.
0798     muonVec_bxm2 = muonVec_bxm1;
0799     muonShowerVec_bxm2 = muonShowerVec_bxm1;
0800     egammaVec_bxm2 = egammaVec_bxm1;
0801     tauVec_bxm2 = tauVec_bxm1;
0802     jetVec_bxm2 = jetVec_bxm1;
0803     etsumVec_bxm2 = etsumVec_bxm1;
0804     extCond_bxm2 = extCond_bxm1;
0805 
0806     muonVec_bxm1 = muonVec_bx0;
0807     muonShowerVec_bxm1 = muonShowerVec_bx0;
0808     egammaVec_bxm1 = egammaVec_bx0;
0809     tauVec_bxm1 = tauVec_bx0;
0810     jetVec_bxm1 = jetVec_bx0;
0811     etsumVec_bxm1 = etsumVec_bx0;
0812     extCond_bxm1 = extCond_bx0;
0813 
0814     muonVec_bx0 = muonVec_bxp1;
0815     muonShowerVec_bx0 = muonShowerVec_bxp1;
0816     egammaVec_bx0 = egammaVec_bxp1;
0817     tauVec_bx0 = tauVec_bxp1;
0818     jetVec_bx0 = jetVec_bxp1;
0819     etsumVec_bx0 = etsumVec_bxp1;
0820     extCond_bx0 = extCond_bxp1;
0821 
0822     muonVec_bxp1 = muonVec;
0823     muonShowerVec_bxp1 = muonShowerVec;
0824     egammaVec_bxp1 = egammaVec;
0825     tauVec_bxp1 = tauVec;
0826     jetVec_bxp1 = jetVec;
0827     etsumVec_bxp1 = etsumVec;
0828     extCond_bxp1 = extCond_bx;
0829   }
0830 
0831   // ------------ method called once each job just before starting event loop ------------
0832   void GenToInputProducer::beginJob() {}
0833 
0834   // ------------ method called once each job just after ending the event loop ------------
0835   void GenToInputProducer::endJob() {}
0836 
0837   // ------------ method called when starting to processes a run ------------
0838 
0839   void GenToInputProducer::beginRun(Run const& iR, EventSetup const& iE) {
0840     LogDebug("GtGenToInputProducer") << "GenToInputProducer::beginRun function called...\n";
0841 
0842     srand(0);
0843 
0844     gRandom = new TRandom3();
0845     mushowerRandom = new TRandom3();
0846   }
0847 
0848   // ------------ method called when ending the processing of a run ------------
0849   void GenToInputProducer::endRun(Run const& iR, EventSetup const& iE) {}
0850 
0851   // ------------ methods to convert from physical to HW values ------------
0852   int GenToInputProducer::convertPhiToHW(double iphi, int steps) {
0853     double phiMax = 2 * M_PI;
0854     if (iphi < 0)
0855       iphi += 2 * M_PI;
0856     if (iphi > phiMax)
0857       iphi -= phiMax;
0858 
0859     int hwPhi = int((iphi / phiMax) * steps + 0.00001);
0860     return hwPhi;
0861   }
0862 
0863   int GenToInputProducer::convertEtaToHW(double ieta, double minEta, double maxEta, int steps) {
0864     double binWidth = (maxEta - minEta) / steps;
0865 
0866     //if we are outside the limits, set error
0867     if (ieta < minEta)
0868       return 99999;  //ieta = minEta+binWidth/2.;
0869     if (ieta > maxEta)
0870       return 99999;  //ieta = maxEta-binWidth/2.;
0871 
0872     int binNum = (int)(ieta / binWidth);
0873     if (ieta < 0.)
0874       binNum--;
0875 
0876     //   unsigned int hwEta = binNum & bitMask;
0877     //   Remove masking for BXVectors...only assume in raw data
0878 
0879     return binNum;
0880   }
0881 
0882   int GenToInputProducer::convertPtToHW(double ipt, int maxPt, double step) {
0883     int hwPt = int(ipt / step + 0.0001);
0884     // if above max Pt, set to largest value
0885     if (hwPt > maxPt)
0886       hwPt = maxPt;
0887 
0888     return hwPt;
0889   }
0890 
0891   // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
0892   void GenToInputProducer::fillDescriptions(ConfigurationDescriptions& descriptions) {
0893     //The following says we do not know what parameters are allowed so do no validation
0894     // Please change this to state exactly what you do use, even if it is no parameters
0895     ParameterSetDescription desc;
0896     desc.setUnknown();
0897     descriptions.addDefault(desc);
0898   }
0899 
0900 }  // namespace l1t
0901 
0902 //define this as a plug-in
0903 DEFINE_FWK_MODULE(l1t::GenToInputProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!