Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:16 

0001 // system include files
0002 #include <memory>
0003 
0004 // Root objects
0005 #include "TROOT.h"
0006 #include "TSystem.h"
0007 #include "TFile.h"
0008 #include "TProfile.h"
0009 #include "TDirectory.h"
0010 #include "TTree.h"
0011 #include "TLorentzVector.h"
0012 #include "TInterpreter.h"
0013 
0014 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0015 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0016 #include "Calibration/IsolatedParticles/interface/eCone.h"
0017 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0018 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0019 
0020 //L1 trigger Menus etc
0021 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0022 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenuFwd.h"
0023 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0024 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0025 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0026 
0027 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0028 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0029 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0030 
0031 //Tracks
0032 #include "DataFormats/TrackReco/interface/Track.h"
0033 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0034 #include "DataFormats/TrackReco/interface/HitPattern.h"
0035 #include "DataFormats/TrackReco/interface/TrackBase.h"
0036 // Vertices
0037 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0040 // Jets
0041 #include "DataFormats/JetReco/interface/GenJet.h"
0042 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0043 #include "DataFormats/JetReco/interface/PFJet.h"
0044 
0045 //Triggers
0046 #include "DataFormats/Common/interface/TriggerResults.h"
0047 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0048 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0049 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0050 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0051 #include "DataFormats/Math/interface/LorentzVector.h"
0052 
0053 #include "FWCore/Framework/interface/Frameworkfwd.h"
0054 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0055 #include "FWCore/Framework/interface/Event.h"
0056 #include "FWCore/Framework/interface/MakerMacros.h"
0057 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0058 #include "FWCore/ServiceRegistry/interface/Service.h"
0059 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0060 #include "FWCore/Common/interface/TriggerNames.h"
0061 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0062 
0063 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0064 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0065 
0066 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0067 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0068 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0069 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0070 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0071 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0072 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0073 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0074 
0075 #include "MagneticField/Engine/interface/MagneticField.h"
0076 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0077 
0078 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0079 
0080 class IsoTrackCalib : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0081 public:
0082   explicit IsoTrackCalib(const edm::ParameterSet&);
0083   ~IsoTrackCalib() override;
0084 
0085   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0086 
0087 private:
0088   void beginJob() override;
0089   void analyze(const edm::Event&, const edm::EventSetup&) override;
0090   void endJob() override {}
0091   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0092   void endRun(edm::Run const&, edm::EventSetup const&) override;
0093 
0094   double dEta(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0095   double dPhi(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0096   double dR(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0097   double deltaR(double eta1, double eta2, double phi1, double phi2);
0098 
0099   edm::Service<TFileService> fs_;
0100   HLTConfigProvider hltConfig_;
0101   L1GtUtils m_l1GtUtils;
0102   const L1GtTriggerMenu* m_l1GtMenu;
0103   const int verbosity_;
0104   const std::vector<std::string> l1Names_;
0105   spr::trackSelectionParameters selectionParameters_;
0106   const std::string theTrackQuality_;
0107   const double a_coneR_, a_charIsoR_, a_mipR_;
0108   std::vector<bool>* t_l1bits;
0109 
0110   const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0111   const edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0112   const edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0113   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0114   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0115   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0116   const edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0117   const edm::EDGetTokenT<reco::GenJetCollection> tok_jets_;
0118   const edm::EDGetTokenT<reco::PFJetCollection> tok_pfjets_;
0119   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extTauJet_;
0120   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extCenJet_;
0121   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extFwdJet_;
0122 
0123   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0124   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0125 
0126   TTree* tree;
0127   int t_Run, t_Event, t_ieta;
0128   double t_EventWeight, t_l1pt, t_l1eta, t_l1phi;
0129   double t_l3pt, t_l3eta, t_l3phi, t_p, t_mindR1;
0130   double t_mindR2, t_eMipDR, t_eHcal, t_hmaxNearP;
0131 
0132   bool t_selectTk, t_qltyFlag, t_qltyMissFlag, t_qltyPVFlag;
0133   std::vector<unsigned int>* t_DetIds;
0134   std::vector<double>*t_HitEnergies, pbin;
0135   TProfile *h_RecHit_iEta, *h_RecHit_num;
0136   TH1I *h_iEta, *h_tkEta0[5], *h_tkEta1[5], *h_tkEta2[5];
0137   TH1I *h_tkEta3[5], *h_tkEta4[5], *h_tkEta5[5];
0138   TH1F *h_Rechit_E, *h_jetp;
0139   TH1F* h_jetpt[4];
0140   TH1I *h_tketa0[6], *h_tketa1[6], *h_tketa2[6];
0141   TH1I *h_tketa3[6], *h_tketa4[6], *h_tketa5[6];
0142   std::map<std::pair<unsigned int, std::string>, int> l1AlgoMap_;
0143 };
0144 
0145 static const bool useL1GtTriggerMenuLite(true);
0146 IsoTrackCalib::IsoTrackCalib(const edm::ParameterSet& iConfig)
0147     : m_l1GtUtils(iConfig, consumesCollector(), useL1GtTriggerMenuLite, *this, L1GtUtils::UseEventSetupIn::Event),
0148       verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0149       l1Names_(iConfig.getUntrackedParameter<std::vector<std::string> >("L1Seed")),
0150       theTrackQuality_(iConfig.getUntrackedParameter<std::string>("TrackQuality", "highPurity")),
0151       a_coneR_(iConfig.getUntrackedParameter<double>("ConeRadius", 34.98)),
0152       a_charIsoR_(a_coneR_ + 28.9),
0153       a_mipR_(iConfig.getUntrackedParameter<double>("ConeRadiusMIP", 14.0)),
0154       tok_genTrack_(consumes<reco::TrackCollection>(edm::InputTag("generalTracks"))),
0155       tok_recVtx_(consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"))),
0156       tok_bs_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
0157       tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generatorSmeared"))),
0158       tok_jets_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("JetSource"))),
0159       tok_pfjets_(consumes<reco::PFJetCollection>(edm::InputTag("ak5PFJets"))) {
0160   usesResource(TFileService::kSharedResource);
0161 
0162   //now do whatever initialization is needed
0163   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0164   selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
0165   selectionParameters_.minQuality = trackQuality_;
0166   selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
0167   selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("MaxDzPV", 5.0);
0168   selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("MaxChi2", 5.0);
0169   selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("MaxDpOverP", 0.1);
0170   selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
0171   selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
0172   selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
0173   selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
0174   bool isItAOD = iConfig.getUntrackedParameter<bool>("IsItAOD", false);
0175   edm::InputTag L1extraTauJetSource_ = iConfig.getParameter<edm::InputTag>("L1extraTauJetSource");
0176   edm::InputTag L1extraCenJetSource_ = iConfig.getParameter<edm::InputTag>("L1extraCenJetSource");
0177   edm::InputTag L1extraFwdJetSource_ = iConfig.getParameter<edm::InputTag>("L1extraFwdJetSource");
0178 
0179   // define tokens for access
0180   if (isItAOD) {
0181     tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"));
0182     tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"));
0183     tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("reducedHcalRecHits", "hbhereco"));
0184   } else {
0185     tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0186     tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0187     tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0188   }
0189   tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
0190   tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
0191   tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
0192   if (verbosity_ >= 0) {
0193     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0194                                  << "\t minPt " << selectionParameters_.minPt << "\t theTrackQuality "
0195                                  << theTrackQuality_ << "\t minQuality " << selectionParameters_.minQuality
0196                                  << "\t maxDxyPV " << selectionParameters_.maxDxyPV << "\t maxDzPV "
0197                                  << selectionParameters_.maxDzPV << "\t maxChi2 " << selectionParameters_.maxChi2
0198                                  << "\t maxDpOverP " << selectionParameters_.maxDpOverP << "\t minOuterHit "
0199                                  << selectionParameters_.minOuterHit << "\t minLayerCrossed "
0200                                  << selectionParameters_.minLayerCrossed << "\t maxInMiss "
0201                                  << selectionParameters_.maxInMiss << "\t maxOutMiss "
0202                                  << selectionParameters_.maxOutMiss << "\t a_coneR " << a_coneR_ << "\t a_charIsoR "
0203                                  << a_charIsoR_ << "\t a_mipR " << a_mipR_ << "\t isItAOD " << isItAOD;
0204     edm::LogVerbatim("IsoTrack") << l1Names_.size() << " triggers to be studied";
0205     for (unsigned int k = 0; k < l1Names_.size(); ++k)
0206       edm::LogVerbatim("IsoTrack") << "[" << k << "]: " << l1Names_[k];
0207   }
0208 
0209   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0210   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0211 }
0212 
0213 IsoTrackCalib::~IsoTrackCalib() {
0214   // do anything here that needs to be done at desctruction time
0215   // (e.g. close files, deallocate resources etc.)
0216 }
0217 
0218 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0219 void IsoTrackCalib::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0220   std::vector<std::string> seeds = {
0221       "L1_SingleJet36", "L1_SingleJet52", "L1_SingleJet68", "L1_SingleJet92", "L1_SingleJet128"};
0222   edm::ParameterSetDescription desc;
0223   desc.addUntracked<int>("Verbosity", 0);
0224   desc.addUntracked<std::vector<std::string> >("L1Seed", seeds);
0225   desc.addUntracked<std::string>("TrackQuality", "highPurity");
0226   desc.addUntracked<double>("MinTrackPt", 10.0);
0227   desc.addUntracked<double>("MaxDxyPV", 0.02);
0228   desc.addUntracked<double>("MaxDzPV", 0.02);
0229   desc.addUntracked<double>("MaxChi2", 5.0);
0230   desc.addUntracked<double>("MaxDpOverP", 0.1);
0231   desc.addUntracked<int>("MinOuterHit", 4);
0232   desc.addUntracked<int>("MinLayerCrossed", 8);
0233   desc.addUntracked<int>("MaxInMiss", 0);
0234   desc.addUntracked<int>("MaxOutMiss", 0);
0235   desc.addUntracked<double>("ConeRadius", 34.98);
0236   desc.addUntracked<double>("ConeRadiusMIP", 14.0);
0237   desc.addUntracked<bool>("IsItAOD", false);
0238   descriptions.add("isoTrackCalib", desc);
0239 }
0240 
0241 void IsoTrackCalib::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0242   t_Run = iEvent.id().run();
0243   t_Event = iEvent.id().event();
0244   if (verbosity_ % 10 > 0)
0245     edm::LogVerbatim("IsoTrack") << "Run " << t_Run << " Event " << t_Event << " Luminosity "
0246                                  << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing()
0247                                  << " starts ==========";
0248 
0249   //Get magnetic field and geometry
0250   const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0251   const MagneticField* bField = &iSetup.getData(tok_magField_);
0252 
0253   //Get track collection
0254   edm::Handle<reco::TrackCollection> trkCollection;
0255   iEvent.getByToken(tok_genTrack_, trkCollection);
0256 
0257   //event weight for FLAT sample
0258   t_EventWeight = 1.0;
0259   edm::Handle<GenEventInfoProduct> genEventInfo;
0260   iEvent.getByToken(tok_ew_, genEventInfo);
0261   if (genEventInfo.isValid())
0262     t_EventWeight = genEventInfo->weight();
0263 
0264   // genJet information
0265   edm::Handle<reco::GenJetCollection> genJets;
0266   iEvent.getByToken(tok_jets_, genJets);
0267   if (genJets.isValid()) {
0268     for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
0269       const reco::GenJet& genJet = (*genJets)[iGenJet];
0270       double genJetPt = genJet.pt();
0271       double genJetEta = genJet.eta();
0272       h_jetpt[0]->Fill(genJetPt);
0273       h_jetpt[1]->Fill(genJetPt, t_EventWeight);
0274       if (genJetEta > -2.5 && genJetEta < 2.5) {
0275         h_jetpt[2]->Fill(genJetPt);
0276         h_jetpt[3]->Fill(genJetPt, t_EventWeight);
0277       }
0278       break;
0279     }
0280   }
0281 
0282   //pf jets
0283   edm::Handle<reco::PFJetCollection> pfJets;
0284   iEvent.getByToken(tok_pfjets_, pfJets);
0285 
0286   //Define the best vertex and the beamspot
0287   edm::Handle<reco::VertexCollection> recVtxs;
0288   iEvent.getByToken(tok_recVtx_, recVtxs);
0289   edm::Handle<reco::BeamSpot> beamSpotH;
0290   iEvent.getByToken(tok_bs_, beamSpotH);
0291   math::XYZPoint leadPV(0, 0, 0);
0292   if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
0293     leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
0294   } else if (beamSpotH.isValid()) {
0295     leadPV = beamSpotH->position();
0296   }
0297   if (verbosity_ > 10) {
0298     if ((verbosity_ % 100) / 10 > 2)
0299       edm::LogVerbatim("IsoTrack") << "Primary Vertex " << leadPV;
0300     if (beamSpotH.isValid())
0301       edm::LogVerbatim("IsoTrack") << " Beam Spot " << beamSpotH->position();
0302   }
0303 
0304   // RecHits
0305   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0306   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0307   iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0308   iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0309   edm::Handle<HBHERecHitCollection> hbhe;
0310   iEvent.getByToken(tok_hbhe_, hbhe);
0311   HBHERecHitCollection::const_iterator rhitItr;
0312 
0313   for (rhitItr = hbhe->begin(); rhitItr != hbhe->end(); rhitItr++) {
0314     double rec_energy = rhitItr->energy();
0315     int rec_ieta = rhitItr->id().ieta();
0316     int rec_depth = rhitItr->id().depth();
0317     int rec_zside = rhitItr->id().zside();
0318     double num1_1 = rec_zside * (rec_ieta + 0.2 * (rec_depth - 1));
0319     if (verbosity_ % 10 > 0)
0320       edm::LogVerbatim("IsoTrack") << "detid/rechit/ieta/zside/depth/num "
0321                                    << " = " << rhitItr->id() << "/" << rec_energy << "/" << rec_ieta << "/" << rec_zside
0322                                    << "/" << rec_depth << "/" << num1_1;
0323     h_iEta->Fill(rec_ieta);
0324     h_Rechit_E->Fill(rec_energy);
0325     h_RecHit_iEta->Fill(rec_ieta, rec_energy);
0326     h_RecHit_num->Fill(num1_1, rec_energy);
0327   }
0328 
0329   //Propagate tracks to calorimeter surface)
0330   std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0331   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, ((verbosity_ / 100) % 10 > 2));
0332   std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0333   for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++) {
0334     if (trkDetItr->okHCAL) {
0335       HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0336       int tk_ieta = detId.ieta();
0337       const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0338       double tk_p = pTrack->p();
0339       h_tketa0[0]->Fill(tk_ieta);
0340       for (unsigned int k = 1; k < pbin.size(); ++k) {
0341         if (tk_p >= pbin[k - 1] && tk_p < pbin[k]) {
0342           h_tketa0[k]->Fill(tk_ieta);
0343           break;
0344         }
0345       }
0346     }
0347   }
0348   //////////////////////////////L1 Trigger Results//////////////////////////////////////////////////
0349   t_l1bits->clear();
0350   for (unsigned int i = 0; i < l1Names_.size(); ++i)
0351     t_l1bits->push_back(false);
0352   bool useL1EventSetup = true;
0353   bool useL1GtTriggerMenuLite = true;
0354 
0355   m_l1GtUtils.getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
0356   int iErrorCode = -1;
0357   int l1ConfCode = -1;
0358   const bool l1Conf = m_l1GtUtils.availableL1Configuration(iErrorCode, l1ConfCode);
0359   if (!l1Conf) {
0360     edm::LogVerbatim("IsoTrack") << "\nL1 configuration code:" << l1ConfCode
0361                                  << "\nNo valid L1 trigger configuration available."
0362                                  << "\nSee text above for error code interpretation"
0363                                  << "\nNo return here, in order to test each method"
0364                                  << ", protected against configuration error.";
0365   }
0366 
0367   const AlgorithmMap& algorithmMap = m_l1GtMenu->gtAlgorithmMap();
0368   const std::string& menuName = m_l1GtMenu->gtTriggerMenuName();
0369   if (verbosity_ % 10 > 0)
0370     edm::LogVerbatim("IsoTrack") << "menuName " << menuName << std::endl;
0371 
0372   std::vector<int> algbits;
0373   for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0374     std::string algName = itAlgo->first;
0375     int algBitNumber = (itAlgo->second).algoBitNumber();
0376     bool decision = m_l1GtUtils.decision(iEvent, itAlgo->first, iErrorCode);
0377 
0378     bool l1ok(false);
0379     if (verbosity_ % 10 > 0)
0380       edm::LogVerbatim("IsoTrack") << algName << "  " << algBitNumber << "  " << decision;
0381     for (unsigned int i = 0; i < l1Names_.size(); ++i) {
0382       if (algName.find(l1Names_[i]) != std::string::npos) {
0383         if (verbosity_ % 10 > 0)
0384           edm::LogVerbatim("IsoTrack") << "match found"
0385                                        << " " << algName << "  " << decision;
0386         t_l1bits->at(i) = (decision > 0);
0387         if (decision > 0)
0388           l1ok = true;
0389       }
0390     }
0391     if (verbosity_ % 10 > 0)
0392       edm::LogVerbatim("IsoTrack") << "l1 ok =" << l1ok;
0393 
0394     if (l1ok) {
0395       edm::Handle<l1extra::L1JetParticleCollection> l1TauHandle;
0396       iEvent.getByToken(tok_L1extTauJet_, l1TauHandle);
0397       l1extra::L1JetParticleCollection::const_iterator itr;
0398       double ptTriggered = -10;
0399       double etaTriggered = -100;
0400       double phiTriggered = -100;
0401 
0402       for (itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr) {
0403         if (itr->pt() > ptTriggered) {
0404           ptTriggered = itr->pt();
0405           etaTriggered = itr->eta();
0406           phiTriggered = itr->phi();
0407         }
0408         if (verbosity_ % 10 > 0)
0409           edm::LogVerbatim("IsoTrack") << "tauJ pt " << itr->pt() << "  eta/phi " << itr->eta() << " " << itr->phi();
0410       }
0411       edm::Handle<l1extra::L1JetParticleCollection> l1CenJetHandle;
0412       iEvent.getByToken(tok_L1extCenJet_, l1CenJetHandle);
0413       for (itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr) {
0414         if (itr->pt() > ptTriggered) {
0415           ptTriggered = itr->pt();
0416           etaTriggered = itr->eta();
0417           phiTriggered = itr->phi();
0418         }
0419         if (verbosity_ % 10 > 0)
0420           edm::LogVerbatim("IsoTrack") << "cenJ pt     " << itr->pt() << "  eta/phi " << itr->eta() << " "
0421                                        << itr->phi();
0422       }
0423       edm::Handle<l1extra::L1JetParticleCollection> l1FwdJetHandle;
0424       iEvent.getByToken(tok_L1extFwdJet_, l1FwdJetHandle);
0425       for (itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr) {
0426         if (itr->pt() > ptTriggered) {
0427           ptTriggered = itr->pt();
0428           etaTriggered = itr->eta();
0429           phiTriggered = itr->phi();
0430         }
0431         if (verbosity_ % 10 > 0)
0432           edm::LogVerbatim("IsoTrack") << "forJ pt     " << itr->pt() << " eta/phi " << itr->eta() << " " << itr->phi();
0433       }
0434       if (verbosity_ % 10 > 0)
0435         edm::LogVerbatim("IsoTrack") << "jets pt/eta/phi = " << ptTriggered << "/" << etaTriggered << "/"
0436                                      << phiTriggered;
0437       //////////////////////loop over tracks////////////////////////////////////////
0438       unsigned int nTracks(0);
0439       for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0440            trkDetItr++, nTracks++) {
0441         const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0442         math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0443 
0444         t_mindR1 = deltaR(etaTriggered, v4.eta(), phiTriggered, v4.phi());
0445         t_mindR2 = -999;
0446         if (verbosity_ % 10 > 0)
0447           edm::LogVerbatim("IsoTrack") << "This track : " << nTracks << " (pt/eta/phi/p) :" << pTrack->pt() << "/"
0448                                        << pTrack->eta() << "/" << pTrack->phi() << "/" << pTrack->p();
0449 
0450         if (verbosity_ % 10 > 0)
0451           edm::LogVerbatim("IsoTrack") << "dr values are = " << t_mindR1;
0452 
0453         t_l1pt = ptTriggered;
0454         t_l1eta = etaTriggered;
0455         t_l1phi = phiTriggered;
0456         t_l3pt = -999;
0457         t_l3eta = -999;
0458         t_l3phi = -999;
0459 
0460         //Selection of good track
0461         t_selectTk = spr::goodTrack(pTrack, leadPV, selectionParameters_, ((verbosity_ / 100) % 10 > 2));
0462         spr::trackSelectionParameters oneCutParameters = selectionParameters_;
0463         oneCutParameters.maxDxyPV = 10;
0464         oneCutParameters.maxDzPV = 100;
0465         oneCutParameters.maxInMiss = 2;
0466         oneCutParameters.maxOutMiss = 2;
0467         bool qltyFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, ((verbosity_ / 100) % 10 > 2));
0468         oneCutParameters = selectionParameters_;
0469         oneCutParameters.maxDxyPV = 10;
0470         oneCutParameters.maxDzPV = 100;
0471         t_qltyMissFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, ((verbosity_ / 100) % 10 > 2));
0472         oneCutParameters = selectionParameters_;
0473         oneCutParameters.maxInMiss = 2;
0474         oneCutParameters.maxOutMiss = 2;
0475         t_qltyPVFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, ((verbosity_ / 100) % 10 > 2));
0476         t_ieta = 0;
0477         if (trkDetItr->okHCAL) {
0478           HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0479           t_ieta = detId.ieta();
0480         }
0481         if (verbosity_ % 10 > 0)
0482           edm::LogVerbatim("IsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "/"
0483                                        << trkDetItr->okHCAL;
0484         t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
0485         t_p = pTrack->p();
0486         h_tketa1[0]->Fill(t_ieta);
0487         for (unsigned int k = 1; k < pbin.size(); ++k) {
0488           if (t_p >= pbin[k - 1] && t_p < pbin[k]) {
0489             h_tketa1[k]->Fill(t_ieta);
0490             break;
0491           }
0492         }
0493         if (t_qltyFlag) {
0494           h_tketa2[0]->Fill(t_ieta);
0495           for (unsigned int k = 1; k < pbin.size(); ++k) {
0496             if (t_p >= pbin[k - 1] && t_p < pbin[k]) {
0497               h_tketa2[k]->Fill(t_ieta);
0498               break;
0499             }
0500           }
0501           int nRH_eMipDR(0), nNearTRKs(0), nRecHits(-999);
0502           t_eMipDR = spr::eCone_ecal(geo,
0503                                      barrelRecHitsHandle,
0504                                      endcapRecHitsHandle,
0505                                      trkDetItr->pointHCAL,
0506                                      trkDetItr->pointECAL,
0507                                      a_mipR_,
0508                                      trkDetItr->directionECAL,
0509                                      nRH_eMipDR);
0510           t_DetIds->clear();
0511           t_HitEnergies->clear();
0512           std::vector<DetId> ids;
0513           t_eHcal = spr::eCone_hcal(geo,
0514                                     hbhe,
0515                                     trkDetItr->pointHCAL,
0516                                     trkDetItr->pointECAL,
0517                                     a_coneR_,
0518                                     trkDetItr->directionHCAL,
0519                                     nRecHits,
0520                                     ids,
0521                                     *t_HitEnergies);
0522           for (unsigned int k = 0; k < ids.size(); ++k) {
0523             t_DetIds->push_back(ids[k].rawId());
0524           }
0525           t_hmaxNearP = spr::chargeIsolationCone(
0526               nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 2));
0527           if (t_hmaxNearP < 2) {
0528             h_tketa3[0]->Fill(t_ieta);
0529             for (unsigned int k = 1; k < pbin.size(); ++k) {
0530               if (t_p >= pbin[k - 1] && t_p < pbin[k]) {
0531                 h_tketa3[k]->Fill(t_ieta);
0532                 break;
0533               }
0534             }
0535             if (t_eMipDR < 1) {
0536               h_tketa4[0]->Fill(t_ieta);
0537               for (unsigned int k = 1; k < pbin.size(); ++k) {
0538                 if (t_p >= pbin[k - 1] && t_p < pbin[k]) {
0539                   h_tketa4[k]->Fill(t_ieta);
0540                   break;
0541                 }
0542               }
0543               if (t_mindR1 > 1) {
0544                 h_tketa5[0]->Fill(t_ieta);
0545                 for (unsigned int k = 1; k < pbin.size(); ++k) {
0546                   if (t_p >= pbin[k - 1] && t_p < pbin[k]) {
0547                     h_tketa5[k]->Fill(t_ieta);
0548                     break;
0549                   }
0550                 }
0551               }
0552             }
0553           }
0554           if (verbosity_ % 10 > 0) {
0555             edm::LogVerbatim("IsoTrack") << "This track : " << nTracks << " (pt/eta/phi/p) :" << pTrack->pt() << "/"
0556                                          << pTrack->eta() << "/" << pTrack->phi() << "/" << t_p;
0557             edm::LogVerbatim("IsoTrack") << "e_MIP " << t_eMipDR << " Chg Isolation " << t_hmaxNearP << " eHcal"
0558                                          << t_eHcal << " ieta " << t_ieta << " Quality " << t_qltyMissFlag << ":"
0559                                          << t_qltyPVFlag << ":" << t_selectTk;
0560             for (unsigned int lll = 0; lll < t_DetIds->size(); lll++) {
0561               edm::LogVerbatim("IsoTrack") << "det id is = " << t_DetIds->at(lll) << "  "
0562                                            << " hit enery is  = " << t_HitEnergies->at(lll);
0563             }
0564           }
0565           if (t_p > 20.0 && t_eMipDR < 2.0 && t_hmaxNearP < 10.0) {
0566             tree->Fill();
0567           }
0568         }
0569       }
0570     }
0571   }
0572 }
0573 
0574 void IsoTrackCalib::beginJob() {
0575   h_RecHit_iEta = fs_->make<TProfile>("rechit_ieta", "Rec hit vs. ieta", 60, -30, 30, 0, 1000);
0576   h_RecHit_num = fs_->make<TProfile>("rechit_num", "Rec hit vs. num", 100, 0, 20, 0, 1000);
0577   h_iEta = fs_->make<TH1I>("iEta", "iEta", 60, -30, 30);
0578   h_Rechit_E = fs_->make<TH1F>("Rechit_E", "Rechit_E", 100, 0, 1000);
0579 
0580   double prange[5] = {20, 30, 40, 60, 100};
0581   for (int k = 0; k < 5; ++k)
0582     pbin.push_back(prange[k]);
0583   std::string type[6] = {"All", "Trigger OK", "Tree Selected", "Charge Isolation", "MIP Cut", "L1 Cut"};
0584   for (unsigned int k = 0; k < pbin.size(); ++k) {
0585     char name[20], namp[20], title[100];
0586     if (k == 0)
0587       sprintf(namp, "all momentum");
0588     else
0589       sprintf(namp, "p = %4.0f:%4.0f GeV", pbin[k - 1], pbin[k]);
0590     sprintf(name, "TrackEta0%d", k);
0591     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[0].c_str());
0592     h_tketa0[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0593     sprintf(name, "TrackEta1%d", k);
0594     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[1].c_str());
0595     h_tketa1[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0596     sprintf(name, "TrackEta2%d", k);
0597     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[2].c_str());
0598     h_tketa2[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0599     sprintf(name, "TrackEta3%d", k);
0600     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[3].c_str());
0601     h_tketa3[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0602     sprintf(name, "TrackEta4%d", k);
0603     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[4].c_str());
0604     h_tketa4[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0605     sprintf(name, "TrackEta5%d", k);
0606     sprintf(title, "Track #eta for tracks with %s (%s)", namp, type[5].c_str());
0607     h_tketa5[k] = fs_->make<TH1I>(name, title, 60, -30, 30);
0608   }
0609   h_jetpt[0] = fs_->make<TH1F>("Jetpt0", "Jet p_T (All)", 500, 0., 2500.);
0610   h_jetpt[1] = fs_->make<TH1F>("Jetpt1", "Jet p_T (All Weighted)", 500, 0., 2500.);
0611   h_jetpt[2] = fs_->make<TH1F>("Jetpt2", "Jet p_T (|#eta| < 2.5)", 500, 0., 2500.);
0612   h_jetpt[3] = fs_->make<TH1F>("Jetpt3", "Jet p_T (|#eta| < 2.5 Weighted)", 500, 0., 2500.);
0613 
0614   tree = fs_->make<TTree>("CalibTree", "CalibTree");
0615 
0616   tree->Branch("t_Run", &t_Run, "t_Run/I");
0617   tree->Branch("t_Event", &t_Event, "t_Event/I");
0618   tree->Branch("t_ieta", &t_ieta, "t_ieta/I");
0619   tree->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
0620   tree->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
0621   tree->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
0622   tree->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
0623   tree->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
0624   tree->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
0625   tree->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
0626   tree->Branch("t_p", &t_p, "t_p/D");
0627   tree->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
0628   tree->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
0629   tree->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
0630   tree->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
0631   tree->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
0632   tree->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
0633   tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
0634   tree->Branch("t_qltyMissFlag", &t_qltyMissFlag, "t_qltyMissFlag/O");
0635   tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O)");
0636 
0637   t_DetIds = new std::vector<unsigned int>();
0638   t_HitEnergies = new std::vector<double>();
0639   t_l1bits = new std::vector<bool>();
0640   tree->Branch("t_DetIds", "std::vector<unsigned int>", &t_DetIds);
0641   tree->Branch("t_HitEnergies", "std::vector<double>", &t_HitEnergies);
0642   tree->Branch("t_l1bits", "std::vector<bool>", &t_l1bits);
0643 }
0644 
0645 // ------------ method called when starting to processes a run  ------------
0646 void IsoTrackCalib::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0647   bool changed = false;
0648   bool ok = hltConfig_.init(iRun, iSetup, "HLT", changed);
0649   edm::LogVerbatim("IsoTrack") << "Run " << iRun.run() << " hltconfig.init " << ok;
0650 
0651   int iErrorCode = -1;
0652   m_l1GtMenu = m_l1GtUtils.ptrL1TriggerMenuEventSetup(iErrorCode);
0653   const AlgorithmMap& algorithmMap = m_l1GtMenu->gtAlgorithmMap();
0654   const std::string& menuName = m_l1GtMenu->gtTriggerMenuName();
0655 
0656   if (verbosity_ % 10 > 0)
0657     edm::LogVerbatim("IsoTrack") << "menuName " << menuName;
0658   for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0659     std::string algName = itAlgo->first;
0660     int algBitNumber = (itAlgo->second).algoBitNumber();
0661     l1AlgoMap_.insert(std::pair<std::pair<unsigned int, std::string>, int>(
0662         std::pair<unsigned int, std::string>(algBitNumber, algName), 0));
0663   }
0664   std::map<std::pair<unsigned int, std::string>, int>::iterator itr;
0665   for (itr = l1AlgoMap_.begin(); itr != l1AlgoMap_.end(); itr++) {
0666     if (verbosity_ % 10 > 0)
0667       edm::LogVerbatim("IsoTrack") << " ********** " << (itr->first).first << " " << (itr->first).second << " "
0668                                    << itr->second;
0669   }
0670 }
0671 
0672 // ------------ method called when ending the processing of a run  ------------
0673 void IsoTrackCalib::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0674   edm::LogVerbatim("IsoTrack") << "endRun " << iRun.run() << std::endl;
0675 }
0676 
0677 double IsoTrackCalib::dEta(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
0678   return (vec1.eta() - vec2.eta());
0679 }
0680 
0681 double IsoTrackCalib::dPhi(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
0682   double phi1 = vec1.phi();
0683   if (phi1 < 0)
0684     phi1 += 2.0 * M_PI;
0685   double phi2 = vec2.phi();
0686   if (phi2 < 0)
0687     phi2 += 2.0 * M_PI;
0688   double dphi = phi1 - phi2;
0689   if (dphi > M_PI)
0690     dphi -= 2. * M_PI;
0691   else if (dphi < -M_PI)
0692     dphi += 2. * M_PI;
0693   return dphi;
0694 }
0695 
0696 double IsoTrackCalib::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
0697   double deta = dEta(vec1, vec2);
0698   double dphi = dPhi(vec1, vec2);
0699   return std::sqrt(deta * deta + dphi * dphi);
0700 }
0701 
0702 double IsoTrackCalib::deltaR(double eta1, double eta2, double phi1, double phi2) {
0703   double deta = eta1 - eta2;
0704   if (phi1 < 0)
0705     phi1 += 2.0 * M_PI;
0706   if (phi2 < 0)
0707     phi2 += 2.0 * M_PI;
0708   double dphi = phi1 - phi2;
0709   if (dphi > M_PI)
0710     dphi -= 2. * M_PI;
0711   else if (dphi < -M_PI)
0712     dphi += 2. * M_PI;
0713   return std::sqrt(deta * deta + dphi * dphi);
0714 }
0715 
0716 //define this as a plug-in
0717 DEFINE_FWK_MODULE(IsoTrackCalib);
[ 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!