Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-10 01:53:45 

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