File indexing completed on 2024-04-06 11:59:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <string>
0022
0023
0024 #include "TH1.h"
0025 #include "TH2.h"
0026 #include "TTree.h"
0027
0028
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
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
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
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
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
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
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
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
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
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
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
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);
