Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-21 23:31:23

0001 // -*- C++ -*-
0002 //
0003 // Package:    IsolatedParticles
0004 // Class:      StudyCaloResponse
0005 //
0006 /**\class StudyCaloResponse StudyCaloResponse.cc Calibration/IsolatedParticles/plugins/StudyCaloResponse.cc
0007 
0008  Description: Studies single particle response measurements in data/MC
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Sunanda Banerjee
0015 //         Created:  Thu Mar  4 18:52:02 CST 2011
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 
0023 // Root objects
0024 #include "TH1.h"
0025 #include "TH2.h"
0026 #include "TTree.h"
0027 
0028 // user include files
0029 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0030 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0031 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0032 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0033 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0034 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0035 
0036 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0037 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0038 
0039 #include "DataFormats/Common/interface/TriggerResults.h"
0040 #include "DataFormats/DetId/interface/DetId.h"
0041 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0042 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0043 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0044 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0045 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0046 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0047 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0048 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0049 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0050 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0051 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0052 #include "DataFormats/MuonReco/interface/Muon.h"
0053 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0054 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0055 #include "DataFormats/TrackReco/interface/Track.h"
0056 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0057 #include "DataFormats/TrackReco/interface/HitPattern.h"
0058 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0059 #include "DataFormats/VertexReco/interface/Vertex.h"
0060 
0061 #include "FWCore/Framework/interface/Frameworkfwd.h"
0062 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0063 #include "FWCore/Framework/interface/Event.h"
0064 #include "FWCore/Framework/interface/MakerMacros.h"
0065 #include "FWCore/ServiceRegistry/interface/Service.h"
0066 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0067 #include "FWCore/Common/interface/TriggerNames.h"
0068 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0069 
0070 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0071 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0072 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0073 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0074 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0075 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0076 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0077 
0078 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0079 #include "MagneticField/Engine/interface/MagneticField.h"
0080 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0081 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0082 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0083 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0084 
0085 class StudyCaloResponse : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0086 public:
0087   explicit StudyCaloResponse(const edm::ParameterSet&);
0088   ~StudyCaloResponse() override {}
0089 
0090   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0091 
0092 private:
0093   void analyze(edm::Event const&, edm::EventSetup const&) override;
0094   void beginJob() override;
0095   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0096   void endRun(edm::Run const&, edm::EventSetup const&) override;
0097   virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
0098   virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) {}
0099 
0100   void clear();
0101   void fillTrack(int, double, double, double, double);
0102   void fillIsolation(int, double, double, double);
0103   void fillEnergy(int, int, double, double, double, double, double);
0104   std::string truncate_str(const std::string&);
0105   int trackPID(const reco::Track*, const edm::Handle<reco::GenParticleCollection>&);
0106 
0107   // ----------member data ---------------------------
0108   static const int nPBin_ = 15, nEtaBin_ = 4, nPVBin_ = 4;
0109   static const int nGen_ = (nPVBin_ + 12);
0110   HLTConfigProvider hltConfig_;
0111   edm::Service<TFileService> fs_;
0112   const int verbosity_;
0113   const std::vector<std::string> trigNames_, newNames_;
0114   const edm::InputTag labelMuon_, labelGenTrack_;
0115   const std::string theTrackQuality_;
0116   const double minTrackP_, maxTrackEta_;
0117   const double tMinE_, tMaxE_, tMinH_, tMaxH_;
0118   const double eThrEB_, eThrEE_, eThrHB_, eThrHE_;
0119   const bool isItAOD_, vetoTrigger_, doTree_, vetoMuon_, vetoEcal_;
0120   const double cutMuon_, cutEcal_, cutRatio_;
0121   const std::vector<double> puWeights_;
0122   const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
0123   spr::trackSelectionParameters selectionParameters_;
0124   std::vector<std::string> HLTNames_;
0125   bool changed_, firstEvent_;
0126 
0127   edm::EDGetTokenT<LumiDetails> tok_lumi;
0128   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt;
0129   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes;
0130   edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
0131   edm::EDGetTokenT<reco::MuonCollection> tok_Muon_;
0132   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0133   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0134   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0135   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0136   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0137   edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0138 
0139   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0140   edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0141   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_topo_;
0142   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0143   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0144   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0145 
0146   TH1I *h_nHLT, *h_HLTAccept, *h_HLTCorr, *h_numberPV;
0147   TH1I *h_goodPV, *h_goodRun;
0148   TH2I* h_nHLTvsRN;
0149   std::vector<TH1I*> h_HLTAccepts;
0150   TH1D *h_p[nGen_ + 2], *h_pt[nGen_ + 2], *h_counter[8];
0151   TH1D *h_eta[nGen_ + 2], *h_phi[nGen_ + 2], *h_h_pNew[8];
0152   TH1I* h_ntrk[2];
0153   TH1D *h_maxNearP[2], *h_ene1[2], *h_ene2[2], *h_ediff[2];
0154   TH1D* h_energy[nPVBin_ + 8][nPBin_][nEtaBin_][6];
0155   TTree* tree_;
0156   int nRun_, etaBin_[nEtaBin_ + 1], pvBin_[nPVBin_ + 1];
0157   double pBin_[nPBin_ + 1];
0158   int tr_goodPV, tr_goodRun;
0159   double tr_eventWeight;
0160   std::vector<std::string> tr_TrigName;
0161   std::vector<double> tr_TrkPt, tr_TrkP, tr_TrkEta, tr_TrkPhi;
0162   std::vector<double> tr_MaxNearP31X31, tr_MaxNearHcalP7x7;
0163   std::vector<double> tr_H3x3, tr_H5x5, tr_H7x7;
0164   std::vector<double> tr_FE7x7P, tr_FE11x11P, tr_FE15x15P;
0165   std::vector<bool> tr_SE7x7P, tr_SE11x11P, tr_SE15x15P;
0166   std::vector<int> tr_iEta, tr_TrkID;
0167 };
0168 
0169 StudyCaloResponse::StudyCaloResponse(const edm::ParameterSet& iConfig)
0170     : verbosity_(iConfig.getUntrackedParameter<int>("verbosity")),
0171       trigNames_(iConfig.getUntrackedParameter<std::vector<std::string> >("triggers")),
0172       newNames_(iConfig.getUntrackedParameter<std::vector<std::string> >("newNames")),
0173       labelMuon_(iConfig.getUntrackedParameter<edm::InputTag>("labelMuon")),
0174       labelGenTrack_(iConfig.getUntrackedParameter<edm::InputTag>("labelTrack")),
0175       theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality")),
0176       minTrackP_(iConfig.getUntrackedParameter<double>("minTrackP")),
0177       maxTrackEta_(iConfig.getUntrackedParameter<double>("maxTrackEta")),
0178       tMinE_(iConfig.getUntrackedParameter<double>("timeMinCutECAL")),
0179       tMaxE_(iConfig.getUntrackedParameter<double>("timeMaxCutECAL")),
0180       tMinH_(iConfig.getUntrackedParameter<double>("timeMinCutHCAL")),
0181       tMaxH_(iConfig.getUntrackedParameter<double>("timeMaxCutHCAL")),
0182       eThrEB_(iConfig.getUntrackedParameter<double>("thresholdEB")),
0183       eThrEE_(iConfig.getUntrackedParameter<double>("thresholdEE")),
0184       eThrHB_(iConfig.getUntrackedParameter<double>("thresholdHB")),
0185       eThrHE_(iConfig.getUntrackedParameter<double>("thresholdHE")),
0186       isItAOD_(iConfig.getUntrackedParameter<bool>("isItAOD")),
0187       vetoTrigger_(iConfig.getUntrackedParameter<bool>("vetoTrigger")),
0188       doTree_(iConfig.getUntrackedParameter<bool>("doTree")),
0189       vetoMuon_(iConfig.getUntrackedParameter<bool>("vetoMuon")),
0190       vetoEcal_(iConfig.getUntrackedParameter<bool>("vetoEcal")),
0191       cutMuon_(iConfig.getUntrackedParameter<double>("cutMuon")),
0192       cutEcal_(iConfig.getUntrackedParameter<double>("cutEcal")),
0193       cutRatio_(iConfig.getUntrackedParameter<double>("cutRatio")),
0194       puWeights_(iConfig.getUntrackedParameter<std::vector<double> >("puWeights")),
0195       triggerEvent_(edm::InputTag("hltTriggerSummaryAOD", "", "HLT")),
0196       theTriggerResultsLabel_(edm::InputTag("TriggerResults", "", "HLT")),
0197       nRun_(0) {
0198   usesResource(TFileService::kSharedResource);
0199 
0200   reco::TrackBase::TrackQuality trackQuality = reco::TrackBase::qualityByName(theTrackQuality_);
0201   selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt");
0202   selectionParameters_.minQuality = trackQuality;
0203   selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV");
0204   selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV");
0205   selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2");
0206   selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP");
0207   selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("minOuterHit");
0208   selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("minLayerCrossed");
0209   selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("maxInMiss");
0210   selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("maxOutMiss");
0211 
0212   // define tokens for access
0213   tok_lumi = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
0214   tok_trigEvt = consumes<trigger::TriggerEvent>(triggerEvent_);
0215   tok_trigRes = consumes<edm::TriggerResults>(theTriggerResultsLabel_);
0216   tok_Muon_ = consumes<reco::MuonCollection>(labelMuon_);
0217   tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
0218   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0219   tok_parts_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("particleSource"));
0220 
0221   if (isItAOD_) {
0222     tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"));
0223     tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"));
0224     tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("reducedHcalRecHits", "hbhereco"));
0225   } else {
0226     tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0227     tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0228     tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0229   }
0230   tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0231 
0232   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0233   tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0234   tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0235   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0236   tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0237   tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0238 
0239   edm::LogVerbatim("IsoTrack") << "Verbosity " << verbosity_ << " with " << trigNames_.size() << " triggers:";
0240   for (unsigned int k = 0; k < trigNames_.size(); ++k)
0241     edm::LogVerbatim("IsoTrack") << " [" << k << "] " << trigNames_[k];
0242   edm::LogVerbatim("IsoTrack") << "TrackQuality " << theTrackQuality_ << " Minpt " << selectionParameters_.minPt
0243                                << " maxDxy " << selectionParameters_.maxDxyPV << " maxDz "
0244                                << selectionParameters_.maxDzPV << " maxChi2 " << selectionParameters_.maxChi2
0245                                << " maxDp/p " << selectionParameters_.maxDpOverP << " minOuterHit "
0246                                << selectionParameters_.minOuterHit << " minLayerCrossed "
0247                                << selectionParameters_.minLayerCrossed << " maxInMiss "
0248                                << selectionParameters_.maxInMiss << " maxOutMiss " << selectionParameters_.maxOutMiss
0249                                << " minTrackP " << minTrackP_ << " maxTrackEta " << maxTrackEta_ << " tMinE_ " << tMinE_
0250                                << " tMaxE " << tMaxE_ << " tMinH_ " << tMinH_ << " tMaxH_ " << tMaxH_ << " isItAOD "
0251                                << isItAOD_ << " doTree " << doTree_ << " vetoTrigger " << vetoTrigger_ << " vetoMuon "
0252                                << vetoMuon_ << ":" << cutMuon_ << " vetoEcal " << vetoEcal_ << ":" << cutEcal_ << ":"
0253                                << cutRatio_;
0254 
0255   double pBins[nPBin_ + 1] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 11.0, 15.0, 20.0, 25.0, 30.0, 40.0, 60.0, 100.0};
0256   int etaBins[nEtaBin_ + 1] = {1, 7, 13, 17, 23};
0257   int pvBins[nPVBin_ + 1] = {1, 2, 3, 5, 100};
0258   for (int i = 0; i <= nPBin_; ++i)
0259     pBin_[i] = pBins[i];
0260   for (int i = 0; i <= nEtaBin_; ++i)
0261     etaBin_[i] = etaBins[i];
0262   for (int i = 0; i <= nPVBin_; ++i)
0263     pvBin_[i] = pvBins[i];
0264 
0265   firstEvent_ = true;
0266   changed_ = false;
0267 }
0268 
0269 void StudyCaloResponse::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0270   std::vector<std::string> trig;
0271   std::vector<double> weights;
0272   std::vector<std::string> newNames = {"HLT", "PixelTracks_Multiplicity", "HLT_Physics_", "HLT_JetE", "HLT_ZeroBias"};
0273   edm::ParameterSetDescription desc;
0274   desc.add<edm::InputTag>("particleSource", edm::InputTag("genParticles"));
0275   desc.addUntracked<int>("verbosity", 0);
0276   desc.addUntracked<std::vector<std::string> >("triggers", trig);
0277   desc.addUntracked<std::vector<std::string> >("newNames", newNames);
0278   desc.addUntracked<edm::InputTag>("labelMuon", edm::InputTag("muons", "", "RECO"));
0279   desc.addUntracked<edm::InputTag>("labelTrack", edm::InputTag("generalTracks", "", "RECO"));
0280   desc.addUntracked<std::string>("trackQuality", "highPurity");
0281   desc.addUntracked<double>("minTrackPt", 1.0);
0282   desc.addUntracked<double>("maxDxyPV", 0.02);
0283   desc.addUntracked<double>("maxDzPV", 0.02);
0284   desc.addUntracked<double>("maxChi2", 5.0);
0285   desc.addUntracked<double>("maxDpOverP", 0.1);
0286   desc.addUntracked<int>("minOuterHit", 4);
0287   desc.addUntracked<int>("minLayerCrossed", 8);
0288   desc.addUntracked<int>("maxInMiss", 0);
0289   desc.addUntracked<int>("maxOutMiss", 0);
0290   desc.addUntracked<double>("minTrackP", 1.0);
0291   desc.addUntracked<double>("maxTrackEta", 2.6);
0292   desc.addUntracked<double>("timeMinCutECAL", -500.0);
0293   desc.addUntracked<double>("timeMaxCutECAL", 500.0);
0294   desc.addUntracked<double>("timeMinCutHCAL", -500.0);
0295   desc.addUntracked<double>("timeMaxCutHCAL", 500.0);
0296   desc.addUntracked<double>("thresholdEB", 0.030);
0297   desc.addUntracked<double>("thresholdEE", 0.150);
0298   desc.addUntracked<double>("thresholdHB", 0.7);
0299   desc.addUntracked<double>("thresholdHE", 0.8);
0300   desc.addUntracked<bool>("isItAOD", false);
0301   desc.addUntracked<bool>("vetoTrigger", false);
0302   desc.addUntracked<bool>("doTree", false);
0303   desc.addUntracked<bool>("vetoMuon", true);
0304   desc.addUntracked<double>("cutMuon", 0.1);
0305   desc.addUntracked<bool>("vetoEcal", false);
0306   desc.addUntracked<double>("cutEcal", 2.0);
0307   desc.addUntracked<double>("cutRatio", 0.9);
0308   desc.addUntracked<std::vector<double> >("puWeights", weights);
0309   descriptions.add("studyCaloResponse", desc);
0310 }
0311 
0312 void StudyCaloResponse::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0313   clear();
0314   int counter0[1000] = {0};
0315   int counter1[1000] = {0};
0316   int counter2[1000] = {0};
0317   int counter3[1000] = {0};
0318   int counter4[1000] = {0};
0319   int counter5[1000] = {0};
0320   int counter6[1000] = {0};
0321   int counter7[1000] = {0};
0322   if (verbosity_ > 0)
0323     edm::LogVerbatim("IsoTrack") << "Event starts====================================";
0324   int RunNo = iEvent.id().run();
0325   int EvtNo = iEvent.id().event();
0326   int Lumi = iEvent.luminosityBlock();
0327   int Bunch = iEvent.bunchCrossing();
0328 
0329   std::vector<int> newAccept(newNames_.size() + 1, 0);
0330   if (verbosity_ > 0)
0331     edm::LogVerbatim("IsoTrack") << "RunNo " << RunNo << " EvtNo " << EvtNo << " Lumi " << Lumi << " Bunch " << Bunch;
0332 
0333   trigger::TriggerEvent triggerEvent;
0334   edm::Handle<trigger::TriggerEvent> triggerEventHandle;
0335   iEvent.getByToken(tok_trigEvt, triggerEventHandle);
0336 
0337   bool ok(false);
0338   std::string triggerUse("None");
0339   if (!triggerEventHandle.isValid()) {
0340     if (trigNames_.empty()) {
0341       ok = true;
0342     } else {
0343       edm::LogWarning("StudyCaloResponse") << "Error! Can't get the product " << triggerEvent_.label();
0344     }
0345   } else {
0346     triggerEvent = *(triggerEventHandle.product());
0347 
0348     /////////////////////////////TriggerResults
0349     edm::Handle<edm::TriggerResults> triggerResults;
0350     iEvent.getByToken(tok_trigRes, triggerResults);
0351 
0352     if (triggerResults.isValid()) {
0353       h_nHLT->Fill(triggerResults->size());
0354       h_nHLTvsRN->Fill(RunNo, triggerResults->size());
0355 
0356       const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0357       const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
0358       for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0359         int ipos = -1;
0360         std::string newtriggerName = truncate_str(triggerNames_[iHLT]);
0361         for (unsigned int i = 0; i < HLTNames_.size(); ++i) {
0362           if (newtriggerName == HLTNames_[i]) {
0363             ipos = i + 1;
0364             break;
0365           }
0366         }
0367         if (ipos < 0) {
0368           HLTNames_.push_back(newtriggerName);
0369           ipos = (int)(HLTNames_.size());
0370           if (ipos <= h_HLTAccept->GetNbinsX())
0371             h_HLTAccept->GetXaxis()->SetBinLabel(ipos, newtriggerName.c_str());
0372         }
0373         if ((int)(iHLT + 1) > h_HLTAccepts[nRun_]->GetNbinsX()) {
0374           edm::LogVerbatim("IsoTrack") << "Wrong trigger " << RunNo << " Event " << EvtNo << " Hlt " << iHLT;
0375         } else {
0376           if (firstEvent_)
0377             h_HLTAccepts[nRun_]->GetXaxis()->SetBinLabel(iHLT + 1, newtriggerName.c_str());
0378         }
0379         int hlt = triggerResults->accept(iHLT);
0380         if (hlt) {
0381           h_HLTAccepts[nRun_]->Fill(iHLT + 1);
0382           h_HLTAccept->Fill(ipos);
0383         }
0384         if (trigNames_.empty()) {
0385           ok = true;
0386         } else {
0387           for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0388             if (newtriggerName.find(trigNames_[i]) != std::string::npos) {
0389               if (verbosity_ % 10 > 0)
0390                 edm::LogVerbatim("IsoTrack") << newtriggerName;
0391               if (hlt > 0) {
0392                 if (!ok)
0393                   triggerUse = newtriggerName;
0394                 ok = true;
0395                 tr_TrigName.push_back(newtriggerName);
0396               }
0397             }
0398           }
0399           if (vetoTrigger_)
0400             ok = !ok;
0401           for (unsigned int i = 0; i < newNames_.size(); ++i) {
0402             if (newtriggerName.find(newNames_[i]) != std::string::npos) {
0403               if (verbosity_ % 10 > 0)
0404                 edm::LogVerbatim("IsoTrack") << "[" << i << "] " << newNames_[i] << " : " << newtriggerName;
0405               if (hlt > 0)
0406                 newAccept[i] = 1;
0407             }
0408           }
0409         }
0410       }
0411       int iflg(0), indx(1);
0412       for (unsigned int i = 0; i < newNames_.size(); ++i) {
0413         iflg += (indx * newAccept[i]);
0414         indx *= 2;
0415       }
0416       h_HLTCorr->Fill(iflg);
0417     }
0418   }
0419   if ((verbosity_ / 10) % 10 > 0)
0420     edm::LogVerbatim("IsoTrack") << "Trigger check gives " << ok << " with " << triggerUse;
0421 
0422   //Look at the tracks
0423   edm::Handle<reco::TrackCollection> trkCollection;
0424   iEvent.getByToken(tok_genTrack_, trkCollection);
0425 
0426   edm::Handle<reco::MuonCollection> muonEventHandle;
0427   iEvent.getByToken(tok_Muon_, muonEventHandle);
0428 
0429   edm::Handle<reco::VertexCollection> recVtxs;
0430   iEvent.getByToken(tok_recVtx_, recVtxs);
0431 
0432   if ((!trkCollection.isValid()) || (!muonEventHandle.isValid()) || (!recVtxs.isValid())) {
0433     edm::LogWarning("StudyCaloResponse") << "Track collection " << trkCollection.isValid() << " Muon collection "
0434                                          << muonEventHandle.isValid() << " Vertex Collecttion " << recVtxs.isValid();
0435     ok = false;
0436   }
0437 
0438   if (ok) {
0439     h_goodRun->Fill(RunNo);
0440     tr_goodRun = RunNo;
0441     // get handles to calogeometry and calotopology
0442     const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0443     const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0444     const HcalTopology* theHBHETopology = &iSetup.getData(tok_topo_);
0445     const MagneticField* bField = &iSetup.getData(tok_magField_);
0446     const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0447 
0448     int ntrk(0), ngoodPV(0), nPV(-1), nvtxs(0);
0449     nvtxs = (int)(recVtxs->size());
0450     for (int ind = 0; ind < nvtxs; ind++) {
0451       if (!((*recVtxs)[ind].isFake()) && (*recVtxs)[ind].ndof() > 4)
0452         ngoodPV++;
0453     }
0454     for (int i = 0; i < nPVBin_; ++i) {
0455       if (ngoodPV >= pvBin_[i] && ngoodPV < pvBin_[i + 1]) {
0456         nPV = i;
0457         break;
0458       }
0459     }
0460 
0461     tr_eventWeight = 1.0;
0462     edm::Handle<GenEventInfoProduct> genEventInfo;
0463     iEvent.getByToken(tok_ew_, genEventInfo);
0464     if (genEventInfo.isValid())
0465       tr_eventWeight = genEventInfo->weight();
0466 
0467     if ((verbosity_ / 10) % 10 > 0)
0468       edm::LogVerbatim("IsoTrack") << "Number of vertices: " << nvtxs << " Good " << ngoodPV << " Bin " << nPV
0469                                    << " Event weight " << tr_eventWeight;
0470     h_numberPV->Fill(nvtxs, tr_eventWeight);
0471     h_goodPV->Fill(ngoodPV, tr_eventWeight);
0472     tr_goodPV = ngoodPV;
0473 
0474     if (!puWeights_.empty()) {
0475       int npbin = h_goodPV->FindBin(ngoodPV);
0476       if (npbin > 0 && npbin <= (int)(puWeights_.size()))
0477         tr_eventWeight *= puWeights_[npbin - 1];
0478       else
0479         tr_eventWeight = 0;
0480     }
0481 
0482     //=== genParticle information
0483     edm::Handle<reco::GenParticleCollection> genParticles;
0484     iEvent.getByToken(tok_parts_, genParticles);
0485     if (genParticles.isValid()) {
0486       for (const auto& p : (reco::GenParticleCollection)(*genParticles)) {
0487         double pt1 = p.momentum().Rho();
0488         double p1 = p.momentum().R();
0489         double eta1 = p.momentum().Eta();
0490         double phi1 = p.momentum().Phi();
0491         fillTrack(nGen_, pt1, p1, eta1, phi1);
0492         bool match(false);
0493         double phi2(phi1);
0494         if (phi2 < 0)
0495           phi2 += 2.0 * M_PI;
0496         for (const auto& trk : (reco::TrackCollection)(*trkCollection)) {
0497           bool quality = trk.quality(selectionParameters_.minQuality);
0498           if (quality) {
0499             double dEta = trk.eta() - eta1;
0500             double phi0 = trk.phi();
0501             if (phi0 < 0)
0502               phi0 += 2.0 * M_PI;
0503             double dPhi = phi0 - phi2;
0504             if (dPhi > M_PI)
0505               dPhi -= 2. * M_PI;
0506             else if (dPhi < -M_PI)
0507               dPhi += 2. * M_PI;
0508             double dR = sqrt(dEta * dEta + dPhi * dPhi);
0509             if (dR < 0.01) {
0510               match = true;
0511               break;
0512             }
0513           }
0514         }
0515         if (match)
0516           fillTrack(nGen_ + 1, pt1, p1, eta1, phi1);
0517       }
0518     }
0519 
0520     reco::TrackCollection::const_iterator trkItr;
0521     for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); ++trkItr, ++ntrk) {
0522       const reco::Track* pTrack = &(*trkItr);
0523       double pt1 = pTrack->pt();
0524       double p1 = pTrack->p();
0525       double eta1 = pTrack->momentum().eta();
0526       double phi1 = pTrack->momentum().phi();
0527       bool quality = pTrack->quality(selectionParameters_.minQuality);
0528       fillTrack(0, pt1, p1, eta1, phi1);
0529       if (quality)
0530         fillTrack(1, pt1, p1, eta1, phi1);
0531       if (p1 < 1000) {
0532         h_h_pNew[0]->Fill(p1);
0533         ++counter0[(int)(p1)];
0534       }
0535     }
0536     h_ntrk[0]->Fill(ntrk, tr_eventWeight);
0537 
0538     std::vector<spr::propagatedTrackID> trkCaloDets;
0539     spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, ((verbosity_ / 100) % 10 > 0));
0540     std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
0541     for (trkDetItr = trkCaloDets.begin(), ntrk = 0; trkDetItr != trkCaloDets.end(); trkDetItr++, ntrk++) {
0542       const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0543       double pt1 = pTrack->pt();
0544       double p1 = pTrack->p();
0545       double eta1 = pTrack->momentum().eta();
0546       double phi1 = pTrack->momentum().phi();
0547       if ((verbosity_ / 10) % 10 > 0)
0548         edm::LogVerbatim("IsoTrack") << "track: p " << p1 << " pt " << pt1 << " eta " << eta1 << " phi " << phi1
0549                                      << " okEcal " << trkDetItr->okECAL;
0550       fillTrack(2, pt1, p1, eta1, phi1);
0551 
0552       bool vetoMuon(false);
0553       double chiGlobal(0), dr(0);
0554       bool goodGlob(false);
0555       if (vetoMuon_) {
0556         if (muonEventHandle.isValid()) {
0557           for (reco::MuonCollection::const_iterator recMuon = muonEventHandle->begin();
0558                recMuon != muonEventHandle->end();
0559                ++recMuon) {
0560             if (((recMuon->isPFMuon()) && (recMuon->isGlobalMuon() || recMuon->isTrackerMuon())) &&
0561                 (recMuon->innerTrack()->validFraction() > 0.49) && (recMuon->innerTrack().isNonnull())) {
0562               chiGlobal = ((recMuon->globalTrack().isNonnull()) ? recMuon->globalTrack()->normalizedChi2() : 999);
0563               goodGlob = (recMuon->isGlobalMuon() && chiGlobal < 3 &&
0564                           recMuon->combinedQuality().chi2LocalPosition < 12 && recMuon->combinedQuality().trkKink < 20);
0565               if (muon::segmentCompatibility(*recMuon) > (goodGlob ? 0.303 : 0.451)) {
0566                 const reco::Track* pTrack0 = (recMuon->innerTrack()).get();
0567                 dr = reco::deltaR(pTrack0->eta(), pTrack0->phi(), pTrack->eta(), pTrack->phi());
0568                 if (dr < cutMuon_) {
0569                   vetoMuon = true;
0570                   break;
0571                 }
0572               }
0573             }
0574           }
0575         }
0576       }
0577       if ((verbosity_ / 10) % 10 > 0)
0578         edm::LogVerbatim("IsoTrack") << "vetoMuon: " << vetoMuon_ << ":" << vetoMuon << " chi:good:dr " << chiGlobal
0579                                      << ":" << goodGlob << ":" << dr;
0580       if (pt1 > minTrackP_ && std::abs(eta1) < maxTrackEta_ && trkDetItr->okECAL && (!vetoMuon)) {
0581         fillTrack(3, pt1, p1, eta1, phi1);
0582         double maxNearP31x31 =
0583             spr::chargeIsolationEcal(ntrk, trkCaloDets, geo, caloTopology, 15, 15, ((verbosity_ / 1000) % 10 > 0));
0584 
0585         const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0586 
0587         edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0588         edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0589         iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0590         iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0591         // get ECal Tranverse Profile
0592         std::pair<double, bool> e7x7P, e11x11P, e15x15P;
0593         const DetId isoCell = trkDetItr->detIdECAL;
0594         e7x7P = spr::eECALmatrix(isoCell,
0595                                  barrelRecHitsHandle,
0596                                  endcapRecHitsHandle,
0597                                  *theEcalChStatus,
0598                                  geo,
0599                                  caloTopology,
0600                                  sevlv,
0601                                  3,
0602                                  3,
0603                                  eThrEB_,
0604                                  eThrEE_,
0605                                  tMinE_,
0606                                  tMaxE_,
0607                                  ((verbosity_ / 10000) % 10 > 0));
0608         e11x11P = spr::eECALmatrix(isoCell,
0609                                    barrelRecHitsHandle,
0610                                    endcapRecHitsHandle,
0611                                    *theEcalChStatus,
0612                                    geo,
0613                                    caloTopology,
0614                                    sevlv,
0615                                    5,
0616                                    5,
0617                                    eThrEB_,
0618                                    eThrEE_,
0619                                    tMinE_,
0620                                    tMaxE_,
0621                                    ((verbosity_ / 10000) % 10 > 0));
0622         e15x15P = spr::eECALmatrix(isoCell,
0623                                    barrelRecHitsHandle,
0624                                    endcapRecHitsHandle,
0625                                    *theEcalChStatus,
0626                                    geo,
0627                                    caloTopology,
0628                                    sevlv,
0629                                    7,
0630                                    7,
0631                                    eThrEB_,
0632                                    eThrEE_,
0633                                    tMinE_,
0634                                    tMaxE_,
0635                                    ((verbosity_ / 10000) % 10 > 0));
0636 
0637         double maxNearHcalP7x7 =
0638             spr::chargeIsolationHcal(ntrk, trkCaloDets, theHBHETopology, 3, 3, ((verbosity_ / 1000) % 10 > 0));
0639         int ieta(0);
0640         double h3x3(0), h5x5(0), h7x7(0);
0641         fillIsolation(0, maxNearP31x31, e11x11P.first, e15x15P.first);
0642         if ((verbosity_ / 10) % 10 > 0)
0643           edm::LogVerbatim("IsoTrack") << "Accepted Tracks reaching Ecal maxNearP31x31 " << maxNearP31x31 << " e11x11P "
0644                                        << e11x11P.first << " e15x15P " << e15x15P.first << " okHCAL "
0645                                        << trkDetItr->okHCAL;
0646 
0647         int trackID = trackPID(pTrack, genParticles);
0648         if (trkDetItr->okHCAL) {
0649           edm::Handle<HBHERecHitCollection> hbhe;
0650           iEvent.getByToken(tok_hbhe_, hbhe);
0651           const DetId ClosestCell(trkDetItr->detIdHCAL);
0652           ieta = ((HcalDetId)(ClosestCell)).ietaAbs();
0653           h3x3 = spr::eHCALmatrix(theHBHETopology,
0654                                   ClosestCell,
0655                                   hbhe,
0656                                   1,
0657                                   1,
0658                                   false,
0659                                   true,
0660                                   eThrHB_,
0661                                   eThrHE_,
0662                                   -100.0,
0663                                   -100.0,
0664                                   tMinH_,
0665                                   tMaxH_,
0666                                   ((verbosity_ / 10000) % 10 > 0));
0667           h5x5 = spr::eHCALmatrix(theHBHETopology,
0668                                   ClosestCell,
0669                                   hbhe,
0670                                   2,
0671                                   2,
0672                                   false,
0673                                   true,
0674                                   eThrHB_,
0675                                   eThrHE_,
0676                                   -100.0,
0677                                   -100.0,
0678                                   tMinH_,
0679                                   tMaxH_,
0680                                   ((verbosity_ / 10000) % 10 > 0));
0681           h7x7 = spr::eHCALmatrix(theHBHETopology,
0682                                   ClosestCell,
0683                                   hbhe,
0684                                   3,
0685                                   3,
0686                                   false,
0687                                   true,
0688                                   eThrHB_,
0689                                   eThrHE_,
0690                                   -100.0,
0691                                   -100.0,
0692                                   tMinH_,
0693                                   tMaxH_,
0694                                   ((verbosity_ / 10000) % 10 > 0));
0695           fillIsolation(1, maxNearHcalP7x7, h5x5, h7x7);
0696           double eByh = ((e11x11P.second) ? (e11x11P.first / std::max(h3x3, 0.001)) : 0.0);
0697           bool notAnElec = ((vetoEcal_ && e11x11P.second) ? ((e11x11P.first < cutEcal_) || (eByh < cutRatio_)) : true);
0698           if ((verbosity_ / 10) % 10 > 0)
0699             edm::LogVerbatim("IsoTrack") << "Tracks Reaching Hcal maxNearHcalP7x7/h5x5/h7x7 " << maxNearHcalP7x7 << "/"
0700                                          << h5x5 << "/" << h7x7 << " eByh " << eByh << " notAnElec " << notAnElec;
0701           tr_TrkPt.push_back(pt1);
0702           tr_TrkP.push_back(p1);
0703           tr_TrkEta.push_back(eta1);
0704           tr_TrkPhi.push_back(phi1);
0705           tr_TrkID.push_back(trackID);
0706           tr_MaxNearP31X31.push_back(maxNearP31x31);
0707           tr_MaxNearHcalP7x7.push_back(maxNearHcalP7x7);
0708           tr_FE7x7P.push_back(e7x7P.first);
0709           tr_FE11x11P.push_back(e11x11P.first);
0710           tr_FE15x15P.push_back(e15x15P.first);
0711           tr_SE7x7P.push_back(e7x7P.second);
0712           tr_SE11x11P.push_back(e11x11P.second);
0713           tr_SE15x15P.push_back(e15x15P.second);
0714           tr_iEta.push_back(ieta);
0715           tr_H3x3.push_back(h3x3);
0716           tr_H5x5.push_back(h5x5);
0717           tr_H7x7.push_back(h7x7);
0718 
0719           if (maxNearP31x31 < 0 && notAnElec) {
0720             fillTrack(4, pt1, p1, eta1, phi1);
0721             fillEnergy(0, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0722             if (maxNearHcalP7x7 < 0) {
0723               fillTrack(5, pt1, p1, eta1, phi1);
0724               fillEnergy(1, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0725               if ((e11x11P.second) && (e15x15P.second) && (e15x15P.first - e11x11P.first) < 2.0) {
0726                 fillTrack(6, pt1, p1, eta1, phi1);
0727                 fillEnergy(2, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0728                 if (h7x7 - h5x5 < 2.0) {
0729                   fillTrack(7, pt1, p1, eta1, phi1);
0730                   fillEnergy(3, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0731                   if (nPV >= 0) {
0732                     fillTrack(nPV + 8, pt1, p1, eta1, phi1);
0733                     fillEnergy(nPV + 4, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0734                   }
0735                   if (trackID > 0) {
0736                     fillTrack(nPVBin_ + trackID + 7, pt1, p1, eta1, phi1);
0737                     fillEnergy(nPVBin_ + trackID + 3, ieta, p1, e7x7P.first, h3x3, e11x11P.first, h5x5);
0738                   }
0739                   if (p1 < 1000) {
0740                     h_h_pNew[7]->Fill(p1);
0741                     ++counter7[(int)(p1)];
0742                   }
0743                 }
0744                 if (p1 < 1000) {
0745                   h_h_pNew[6]->Fill(p1);
0746                   ++counter6[(int)(p1)];
0747                 }
0748               }
0749               if (p1 < 1000) {
0750                 h_h_pNew[5]->Fill(p1);
0751                 ++counter5[(int)(p1)];
0752               }
0753             }
0754             if (p1 < 1000) {
0755               h_h_pNew[4]->Fill(p1);
0756               ++counter4[(int)(p1)];
0757             }
0758           }
0759           if (p1 < 1000) {
0760             h_h_pNew[3]->Fill(p1);
0761             ++counter3[(int)(p1)];
0762           }
0763         }
0764         if (p1 < 1000) {
0765           h_h_pNew[2]->Fill(p1);
0766           ++counter2[(int)(p1)];
0767         }
0768       }
0769       if (p1 < 1000) {
0770         h_h_pNew[1]->Fill(p1);
0771         ++counter1[(int)(p1)];
0772       }
0773     }
0774     h_ntrk[1]->Fill(ntrk, tr_eventWeight);
0775     if ((!tr_TrkPt.empty()) && doTree_)
0776       tree_->Fill();
0777     for (int i = 0; i < 1000; ++i) {
0778       if (counter0[i])
0779         h_counter[0]->Fill(i, counter0[i]);
0780       if (counter1[i])
0781         h_counter[1]->Fill(i, counter1[i]);
0782       if (counter2[i])
0783         h_counter[2]->Fill(i, counter2[i]);
0784       if (counter3[i])
0785         h_counter[3]->Fill(i, counter3[i]);
0786       if (counter4[i])
0787         h_counter[4]->Fill(i, counter4[i]);
0788       if (counter5[i])
0789         h_counter[5]->Fill(i, counter5[i]);
0790       if (counter6[i])
0791         h_counter[6]->Fill(i, counter6[i]);
0792       if (counter7[i])
0793         h_counter[7]->Fill(i, counter7[i]);
0794     }
0795   }
0796   firstEvent_ = false;
0797 }
0798 
0799 void StudyCaloResponse::beginJob() {
0800   // Book histograms
0801   h_nHLT = fs_->make<TH1I>("h_nHLT", "size of trigger Names", 1000, 0, 1000);
0802   h_HLTAccept = fs_->make<TH1I>("h_HLTAccept", "HLT Accepts for all runs", 500, 0, 500);
0803   for (int i = 1; i <= 500; ++i)
0804     h_HLTAccept->GetXaxis()->SetBinLabel(i, " ");
0805   h_nHLTvsRN = fs_->make<TH2I>("h_nHLTvsRN", "size of trigger Names vs RunNo", 2168, 190949, 193116, 100, 400, 500);
0806   h_HLTCorr = fs_->make<TH1I>("h_HLTCorr", "Correlation among different paths", 100, 0, 100);
0807   h_numberPV = fs_->make<TH1I>("h_numberPV", "Number of Primary Vertex", 100, 0, 100);
0808   h_goodPV = fs_->make<TH1I>("h_goodPV", "Number of good Primary Vertex", 100, 0, 100);
0809   h_goodRun = fs_->make<TH1I>("h_goodRun", "Number of accepted events for Run", 4000, 190000, 1940000);
0810   char hname[60], htit[200];
0811   std::string CollectionNames[2] = {"Reco", "Propagated"};
0812   for (unsigned int i = 0; i < 2; i++) {
0813     sprintf(hname, "h_nTrk_%s", CollectionNames[i].c_str());
0814     sprintf(htit, "Number of %s tracks", CollectionNames[i].c_str());
0815     h_ntrk[i] = fs_->make<TH1I>(hname, htit, 500, 0, 500);
0816   }
0817   std::string TrkNames[8] = {
0818       "All", "Quality", "NoIso", "okEcal", "EcalCharIso", "HcalCharIso", "EcalNeutIso", "HcalNeutIso"};
0819   std::string particle[4] = {"Electron", "Pion", "Kaon", "Proton"};
0820   for (unsigned int i = 0; i <= nGen_ + 1; i++) {
0821     if (i < 8) {
0822       sprintf(hname, "h_pt_%s", TrkNames[i].c_str());
0823       sprintf(htit, "p_{T} of %s tracks", TrkNames[i].c_str());
0824     } else if (i < 8 + nPVBin_) {
0825       sprintf(hname, "h_pt_%s_%d", TrkNames[7].c_str(), i - 8);
0826       sprintf(htit, "p_{T} of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i - 8], pvBin_[i - 7] - 1);
0827     } else if (i >= nGen_) {
0828       sprintf(hname, "h_pt_%s_%d", TrkNames[0].c_str(), i - nGen_);
0829       sprintf(htit, "p_{T} of %s Generator tracks", TrkNames[0].c_str());
0830     } else {
0831       sprintf(hname, "h_pt_%s_%s", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0832       sprintf(htit, "p_{T} of %s tracks (%s)", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0833     }
0834     h_pt[i] = fs_->make<TH1D>(hname, htit, 400, 0, 200.0);
0835     h_pt[i]->Sumw2();
0836 
0837     if (i < 8) {
0838       sprintf(hname, "h_p_%s", TrkNames[i].c_str());
0839       sprintf(htit, "Momentum of %s tracks", TrkNames[i].c_str());
0840     } else if (i < 8 + nPVBin_) {
0841       sprintf(hname, "h_p_%s_%d", TrkNames[7].c_str(), i - 8);
0842       sprintf(htit, "Momentum of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i - 8], pvBin_[i - 7] - 1);
0843     } else if (i >= nGen_) {
0844       sprintf(hname, "h_p_%s_%d", TrkNames[0].c_str(), i - nGen_);
0845       sprintf(htit, "Momentum of %s Generator tracks", TrkNames[0].c_str());
0846     } else {
0847       sprintf(hname, "h_p_%s_%s", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0848       sprintf(htit, "Momentum of %s tracks (%s)", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0849     }
0850     h_p[i] = fs_->make<TH1D>(hname, htit, 400, 0, 200.0);
0851     h_p[i]->Sumw2();
0852 
0853     if (i < 8) {
0854       sprintf(hname, "h_eta_%s", TrkNames[i].c_str());
0855       sprintf(htit, "Eta of %s tracks", TrkNames[i].c_str());
0856     } else if (i < 8 + nPVBin_) {
0857       sprintf(hname, "h_eta_%s_%d", TrkNames[7].c_str(), i - 8);
0858       sprintf(htit, "Eta of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i - 8], pvBin_[i - 7] - 1);
0859     } else if (i >= nGen_) {
0860       sprintf(hname, "h_eta_%s_%d", TrkNames[0].c_str(), i - nGen_);
0861       sprintf(htit, "Eta of %s Generator tracks", TrkNames[0].c_str());
0862     } else {
0863       sprintf(hname, "h_eta_%s_%s", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0864       sprintf(htit, "Eta of %s tracks (%s)", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0865     }
0866     h_eta[i] = fs_->make<TH1D>(hname, htit, 60, -3.0, 3.0);
0867     h_eta[i]->Sumw2();
0868 
0869     if (i < 8) {
0870       sprintf(hname, "h_phi_%s", TrkNames[i].c_str());
0871       sprintf(htit, "Phi of %s tracks", TrkNames[i].c_str());
0872     } else if (i < 8 + nPVBin_) {
0873       sprintf(hname, "h_phi_%s_%d", TrkNames[7].c_str(), i - 8);
0874       sprintf(htit, "Phi of %s tracks (PV=%d:%d)", TrkNames[7].c_str(), pvBin_[i - 8], pvBin_[i - 7] - 1);
0875     } else if (i >= nGen_) {
0876       sprintf(hname, "h_phi_%s_%d", TrkNames[0].c_str(), i - nGen_);
0877       sprintf(htit, "Phi of %s Generator tracks", TrkNames[0].c_str());
0878     } else {
0879       sprintf(hname, "h_phi_%s_%s", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0880       sprintf(htit, "Phi of %s tracks (%s)", TrkNames[7].c_str(), particle[i - 8 - nPVBin_].c_str());
0881     }
0882     h_phi[i] = fs_->make<TH1D>(hname, htit, 100, -3.15, 3.15);
0883     h_phi[i]->Sumw2();
0884   }
0885   std::string IsolationNames[2] = {"Ecal", "Hcal"};
0886   for (unsigned int i = 0; i < 2; i++) {
0887     sprintf(hname, "h_maxNearP_%s", IsolationNames[i].c_str());
0888     sprintf(htit, "Energy in ChargeIso region for %s", IsolationNames[i].c_str());
0889     h_maxNearP[i] = fs_->make<TH1D>(hname, htit, 120, -1.5, 10.5);
0890     h_maxNearP[i]->Sumw2();
0891 
0892     sprintf(hname, "h_ene1_%s", IsolationNames[i].c_str());
0893     sprintf(htit, "Energy in smaller cone for %s", IsolationNames[i].c_str());
0894     h_ene1[i] = fs_->make<TH1D>(hname, htit, 400, 0.0, 200.0);
0895     h_ene1[i]->Sumw2();
0896 
0897     sprintf(hname, "h_ene2_%s", IsolationNames[i].c_str());
0898     sprintf(htit, "Energy in bigger cone for %s", IsolationNames[i].c_str());
0899     h_ene2[i] = fs_->make<TH1D>(hname, htit, 400, 0.0, 200.0);
0900     h_ene2[i]->Sumw2();
0901 
0902     sprintf(hname, "h_ediff_%s", IsolationNames[i].c_str());
0903     sprintf(htit, "Energy in NeutralIso region for %s", IsolationNames[i].c_str());
0904     h_ediff[i] = fs_->make<TH1D>(hname, htit, 100, -0.5, 19.5);
0905     h_ediff[i]->Sumw2();
0906   }
0907   std::string energyNames[6] = {
0908       "E_{7x7}", "H_{3x3}", "(E_{7x7}+H_{3x3})", "E_{11x11}", "H_{5x5}", "{E_{11x11}+H_{5x5})"};
0909   for (int i = 0; i < 4 + nPVBin_ + 4; ++i) {
0910     for (int ip = 0; ip < nPBin_; ++ip) {
0911       for (int ie = 0; ie < nEtaBin_; ++ie) {
0912         for (int j = 0; j < 6; ++j) {
0913           sprintf(hname, "h_energy_%d_%d_%d_%d", i, ip, ie, j);
0914           if (i < 4) {
0915             sprintf(htit,
0916                     "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d) for tracks with %s",
0917                     energyNames[j].c_str(),
0918                     pBin_[ip],
0919                     pBin_[ip + 1],
0920                     etaBin_[ie],
0921                     (etaBin_[ie + 1] - 1),
0922                     TrkNames[i + 4].c_str());
0923           } else if (i < 4 + nPVBin_) {
0924             sprintf(htit,
0925                     "%s/p (p=%4.1f:%4.1f, i#eta=%d:%d, PV=%d:%d) for tracks with %s",
0926                     energyNames[j].c_str(),
0927                     pBin_[ip],
0928                     pBin_[ip + 1],
0929                     etaBin_[ie],
0930                     (etaBin_[ie + 1] - 1),
0931                     pvBin_[i - 4],
0932                     (pvBin_[i - 3] - 1),
0933                     TrkNames[7].c_str());
0934           } else {
0935             sprintf(htit,
0936                     "%s/p (p=%4.1f:%4.1f, i#eta=%d:%d %s) for tracks with %s",
0937                     energyNames[j].c_str(),
0938                     pBin_[ip],
0939                     pBin_[ip + 1],
0940                     etaBin_[ie],
0941                     (etaBin_[ie + 1] - 1),
0942                     particle[i - 4 - nPVBin_].c_str(),
0943                     TrkNames[7].c_str());
0944           }
0945           h_energy[i][ip][ie][j] = fs_->make<TH1D>(hname, htit, 5000, -0.1, 49.9);
0946           h_energy[i][ip][ie][j]->Sumw2();
0947         }
0948       }
0949     }
0950   }
0951 
0952   for (int i = 0; i < 8; ++i) {
0953     sprintf(hname, "counter%d", i);
0954     sprintf(htit, "Counter with cut %d", i);
0955     h_counter[i] = fs_->make<TH1D>(hname, htit, 1000, 0, 1000);
0956     sprintf(hname, "h_pTNew%d", i);
0957     sprintf(htit, "Track momentum with cut %d", i);
0958     h_h_pNew[i] = fs_->make<TH1D>(hname, htit, 1000, 0, 1000);
0959   }
0960 
0961   // Now the tree
0962   if (doTree_) {
0963     tree_ = fs_->make<TTree>("testTree", "new HLT Tree");
0964     tree_->Branch("tr_goodRun", &tr_goodRun, "tr_goodRun/I");
0965     tree_->Branch("tr_goodPV", &tr_goodPV, "tr_goodPV/I");
0966     tree_->Branch("tr_eventWeight", &tr_eventWeight, "tr_eventWeight/D");
0967     tree_->Branch("tr_tr_TrigName", &tr_TrigName);
0968     tree_->Branch("tr_TrkPt", &tr_TrkPt);
0969     tree_->Branch("tr_TrkP", &tr_TrkP);
0970     tree_->Branch("tr_TrkEta", &tr_TrkEta);
0971     tree_->Branch("tr_TrkPhi", &tr_TrkPhi);
0972     tree_->Branch("tr_TrkID", &tr_TrkID);
0973     tree_->Branch("tr_MaxNearP31X31", &tr_MaxNearP31X31);
0974     tree_->Branch("tr_MaxNearHcalP7x7", &tr_MaxNearHcalP7x7);
0975     tree_->Branch("tr_FE7x7P", &tr_FE7x7P);
0976     tree_->Branch("tr_FE11x11P", &tr_FE11x11P);
0977     tree_->Branch("tr_FE15x15P", &tr_FE15x15P);
0978     tree_->Branch("tr_SE7x7P", &tr_SE7x7P);
0979     tree_->Branch("tr_SE11x11P", &tr_SE11x11P);
0980     tree_->Branch("tr_SE15x15P", &tr_SE15x15P);
0981     tree_->Branch("tr_H3x3", &tr_H3x3);
0982     tree_->Branch("tr_H5x5", &tr_H5x5);
0983     tree_->Branch("tr_H7x7", &tr_H7x7);
0984     tree_->Branch("tr_iEta", &tr_iEta);
0985   }
0986 }
0987 
0988 // ------------ method called when starting to processes a run  ------------
0989 void StudyCaloResponse::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0990   char hname[100], htit[400];
0991   edm::LogVerbatim("IsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
0992                                << hltConfig_.init(iRun, iSetup, "HLT", changed_);
0993   sprintf(hname, "h_HLTAccepts_%i", iRun.run());
0994   sprintf(htit, "HLT Accepts for Run No %i", iRun.run());
0995   TH1I* hnew = fs_->make<TH1I>(hname, htit, 500, 0, 500);
0996   for (int i = 1; i <= 500; ++i)
0997     hnew->GetXaxis()->SetBinLabel(i, " ");
0998   h_HLTAccepts.push_back(hnew);
0999   edm::LogVerbatim("IsoTrack") << "beginRun " << iRun.run();
1000   firstEvent_ = true;
1001   changed_ = false;
1002 }
1003 
1004 // ------------ method called when ending the processing of a run  ------------
1005 void StudyCaloResponse::endRun(edm::Run const& iRun, edm::EventSetup const&) {
1006   ++nRun_;
1007   edm::LogVerbatim("IsoTrack") << "endRun[" << nRun_ << "] " << iRun.run();
1008 }
1009 
1010 void StudyCaloResponse::clear() {
1011   tr_TrigName.clear();
1012   tr_TrkPt.clear();
1013   tr_TrkP.clear();
1014   tr_TrkEta.clear();
1015   tr_TrkPhi.clear();
1016   tr_TrkID.clear();
1017   tr_MaxNearP31X31.clear();
1018   tr_MaxNearHcalP7x7.clear();
1019   tr_FE7x7P.clear();
1020   tr_FE11x11P.clear();
1021   tr_FE15x15P.clear();
1022   tr_SE7x7P.clear();
1023   tr_SE11x11P.clear();
1024   tr_SE15x15P.clear();
1025   tr_H3x3.clear();
1026   tr_H5x5.clear();
1027   tr_H7x7.clear();
1028   tr_iEta.clear();
1029 }
1030 
1031 void StudyCaloResponse::fillTrack(int i, double pt, double p, double eta, double phi) {
1032   h_pt[i]->Fill(pt, tr_eventWeight);
1033   h_p[i]->Fill(p, tr_eventWeight);
1034   h_eta[i]->Fill(eta, tr_eventWeight);
1035   h_phi[i]->Fill(phi, tr_eventWeight);
1036 }
1037 
1038 void StudyCaloResponse::fillIsolation(int i, double emaxnearP, double eneutIso1, double eneutIso2) {
1039   h_maxNearP[i]->Fill(emaxnearP, tr_eventWeight);
1040   h_ene1[i]->Fill(eneutIso1, tr_eventWeight);
1041   h_ene2[i]->Fill(eneutIso2, tr_eventWeight);
1042   h_ediff[i]->Fill(eneutIso2 - eneutIso1, tr_eventWeight);
1043 }
1044 
1045 void StudyCaloResponse::fillEnergy(
1046     int flag, int ieta, double p, double enEcal1, double enHcal1, double enEcal2, double enHcal2) {
1047   int ip(-1), ie(-1);
1048   for (int i = 0; i < nPBin_; ++i) {
1049     if (p >= pBin_[i] && p < pBin_[i + 1]) {
1050       ip = i;
1051       break;
1052     }
1053   }
1054   for (int i = 0; i < nEtaBin_; ++i) {
1055     if (ieta >= etaBin_[i] && ieta < etaBin_[i + 1]) {
1056       ie = i;
1057       break;
1058     }
1059   }
1060   if (ip >= 0 && ie >= 0 && enEcal1 > 0.02 && enHcal1 > 0.1) {
1061     h_energy[flag][ip][ie][0]->Fill(enEcal1 / p, tr_eventWeight);
1062     h_energy[flag][ip][ie][1]->Fill(enHcal1 / p, tr_eventWeight);
1063     h_energy[flag][ip][ie][2]->Fill((enEcal1 + enHcal1) / p, tr_eventWeight);
1064     h_energy[flag][ip][ie][3]->Fill(enEcal2 / p, tr_eventWeight);
1065     h_energy[flag][ip][ie][4]->Fill(enHcal2 / p, tr_eventWeight);
1066     h_energy[flag][ip][ie][5]->Fill((enEcal2 + enHcal2) / p, tr_eventWeight);
1067   }
1068 }
1069 
1070 std::string StudyCaloResponse::truncate_str(const std::string& str) {
1071   std::string truncated_str(str);
1072   int length = str.length();
1073   for (int i = 0; i < length - 2; i++) {
1074     if (str[i] == '_' && str[i + 1] == 'v' && isdigit(str.at(i + 2))) {
1075       int z = i + 1;
1076       truncated_str = str.substr(0, z);
1077     }
1078   }
1079   return (truncated_str);
1080 }
1081 
1082 int StudyCaloResponse::trackPID(const reco::Track* pTrack,
1083                                 const edm::Handle<reco::GenParticleCollection>& genParticles) {
1084   int id(0);
1085   if (genParticles.isValid()) {
1086     unsigned int indx;
1087     reco::GenParticleCollection::const_iterator p;
1088     double mindR(999.9);
1089     for (p = genParticles->begin(), indx = 0; p != genParticles->end(); ++p, ++indx) {
1090       int pdgId = std::abs(p->pdgId());
1091       int idx = (pdgId == 11) ? 1 : ((pdgId == 211) ? 2 : ((pdgId == 321) ? 3 : ((pdgId == 2212) ? 4 : 0)));
1092       if (idx > 0) {
1093         double dEta = pTrack->eta() - p->momentum().Eta();
1094         double phi1 = pTrack->phi();
1095         double phi2 = p->momentum().Phi();
1096         if (phi1 < 0)
1097           phi1 += 2.0 * M_PI;
1098         if (phi2 < 0)
1099           phi2 += 2.0 * M_PI;
1100         double dPhi = phi1 - phi2;
1101         if (dPhi > M_PI)
1102           dPhi -= 2. * M_PI;
1103         else if (dPhi < -M_PI)
1104           dPhi += 2. * M_PI;
1105         double dR = sqrt(dEta * dEta + dPhi * dPhi);
1106         if (dR < mindR) {
1107           mindR = dR;
1108           id = idx;
1109         }
1110       }
1111     }
1112   }
1113   return id;
1114 }
1115 
1116 DEFINE_FWK_MODULE(StudyCaloResponse);