Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-12 09:07:16

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