File indexing completed on 2024-04-06 11:59:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0031 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0032 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0033
0034 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0035 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0036 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0037 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0038 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0039
0040 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0041 #include "DataFormats/TrackReco/interface/Track.h"
0042
0043 #include "CondFormats/L1TObjects/interface/L1GtPrescaleFactors.h"
0044 #include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsAlgoTrigRcd.h"
0045 #include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsTechTrigRcd.h"
0046
0047 #include "DQMServices/Core/interface/DQMStore.h"
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049
0050 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0051 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0052
0053
0054
0055 #include "SimDataFormats/Track/interface/SimTrack.h"
0056 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0057
0058
0059
0060 #include "ValidationHcalIsoTrackAlCaReco.h"
0061 #include <fstream>
0062
0063 #include "TH1F.h"
0064
0065 double ValidationHcalIsoTrackAlCaReco::getDist(double eta1, double phi1, double eta2, double phi2) {
0066 double dphi = fabs(phi1 - phi2);
0067 if (dphi > acos(-1))
0068 dphi = 2 * acos(-1) - dphi;
0069 double dr = sqrt(dphi * dphi + pow(eta1 - eta2, 2));
0070 return dr;
0071 }
0072
0073
0074
0075 double ValidationHcalIsoTrackAlCaReco::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0076 double dR, Rec;
0077 double theta1 = 2 * atan(exp(-eta1));
0078 double theta2 = 2 * atan(exp(-eta2));
0079 if (fabs(eta1) < 1.479)
0080 Rec = 129;
0081 else
0082 Rec = 275;
0083
0084 dR = fabs((Rec / sin(theta1)) * tan(acos(sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) +
0085 cos(theta1) * cos(theta2))));
0086 return dR;
0087 }
0088
0089
0090
0091 std::pair<int, int> ValidationHcalIsoTrackAlCaReco::towerIndex(double eta, double phi) {
0092 int ieta = 0;
0093 int iphi = 0;
0094 for (int i = 1; i < 21; i++) {
0095 if (fabs(eta) < (i * 0.087) && fabs(eta) > (i - 1) * 0.087)
0096 ieta = int(fabs(eta) / eta) * i;
0097 }
0098 if (fabs(eta) > 1.740 && fabs(eta) < 1.830)
0099 ieta = int(fabs(eta) / eta) * 21;
0100 if (fabs(eta) > 1.830 && fabs(eta) < 1.930)
0101 ieta = int(fabs(eta) / eta) * 22;
0102 if (fabs(eta) > 1.930 && fabs(eta) < 2.043)
0103 ieta = int(fabs(eta) / eta) * 23;
0104
0105 double delta = phi + 0.174532925;
0106 if (delta < 0)
0107 delta = delta + 2 * acos(-1);
0108 if (fabs(eta) < 1.740) {
0109 for (int i = 0; i < 72; i++) {
0110 if (delta < (i + 1) * 0.087266462 && delta > i * 0.087266462)
0111 iphi = i;
0112 }
0113 } else {
0114 for (int i = 0; i < 36; i++) {
0115 if (delta < 2 * (i + 1) * 0.087266462 && delta > 2 * i * 0.087266462)
0116 iphi = 2 * i;
0117 }
0118 }
0119
0120 return std::pair<int, int>(ieta, iphi);
0121 }
0122
0123 ValidationHcalIsoTrackAlCaReco::ValidationHcalIsoTrackAlCaReco(const edm::ParameterSet& iConfig)
0124 : folderName_(iConfig.getParameter<std::string>("folderName")),
0125 saveToFile_(iConfig.getParameter<bool>("saveToFile")),
0126 outRootFileName_(iConfig.getParameter<std::string>("outputRootFileName")),
0127 hltFilterTag_(iConfig.getParameter<edm::InputTag>("hltL3FilterLabel")),
0128 recoTrLabel_(iConfig.getParameter<edm::InputTag>("recoTracksLabel")),
0129 pThr_(iConfig.getUntrackedParameter<double>("pThrL3", 0)),
0130 heLow_(iConfig.getUntrackedParameter<double>("lowerHighEnergyCut", 40)),
0131 heUp_(iConfig.getUntrackedParameter<double>("upperHighEnergyCut", 60)),
0132 tok_hlt_(consumes<trigger::TriggerEvent>(iConfig.getParameter<edm::InputTag>("hltTriggerEventLabel"))),
0133 tok_arITr_(consumes<reco::IsolatedPixelTrackCandidateCollection>(
0134 iConfig.getParameter<edm::InputTag>("alcarecoIsoTracksLabel"))),
0135 tok_simTrack_(consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("simTracksTag"))) {
0136 nTotal = 0;
0137 nHLTL3accepts = 0;
0138 }
0139
0140 ValidationHcalIsoTrackAlCaReco::~ValidationHcalIsoTrackAlCaReco() {}
0141
0142 void ValidationHcalIsoTrackAlCaReco::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0143 nTotal++;
0144
0145 const edm::Handle<trigger::TriggerEvent>& trEv = iEvent.getHandle(tok_hlt_);
0146
0147 const edm::Handle<reco::IsolatedPixelTrackCandidateCollection>& recoIsoTracks = iEvent.getHandle(tok_arITr_);
0148
0149 const trigger::TriggerObjectCollection& TOCol(trEv->getObjects());
0150
0151 trigger::Keys KEYS;
0152 const trigger::size_type nFilt(trEv->sizeFilters());
0153 for (trigger::size_type iFilt = 0; iFilt != nFilt; iFilt++) {
0154 if (trEv->filterTag(iFilt) == hltFilterTag_) {
0155 KEYS = trEv->filterKeys(iFilt);
0156 }
0157 }
0158
0159 trigger::size_type nReg = KEYS.size();
0160
0161 std::vector<double> trigEta;
0162 std::vector<double> trigPhi;
0163 bool trig = false;
0164
0165
0166 for (trigger::size_type iReg = 0; iReg < nReg; iReg++) {
0167 const trigger::TriggerObject& TObj(TOCol[KEYS[iReg]]);
0168 if (TObj.p() < pThr_)
0169 continue;
0170 hl3eta->Fill(TObj.eta(), 1);
0171 hl3AbsEta->Fill(fabs(TObj.eta()), 1);
0172 hl3phi->Fill(TObj.phi(), 1);
0173
0174 if (recoIsoTracks->size() > 0) {
0175 double minRecoL3dist = 1000;
0176 reco::IsolatedPixelTrackCandidateCollection::const_iterator mrtr;
0177 for (reco::IsolatedPixelTrackCandidateCollection::const_iterator rtrit = recoIsoTracks->begin();
0178 rtrit != recoIsoTracks->end();
0179 rtrit++) {
0180 double R = getDist(rtrit->eta(), rtrit->phi(), TObj.eta(), TObj.phi());
0181 if (R < minRecoL3dist) {
0182 mrtr = rtrit;
0183 minRecoL3dist = R;
0184 }
0185 }
0186 hOffL3TrackMatch->Fill(minRecoL3dist, 1);
0187 hOffL3TrackPtRat->Fill(TObj.pt() / mrtr->pt(), 1);
0188 }
0189
0190 hl3Pt->Fill(TObj.pt(), 1);
0191 trig = true;
0192 trigEta.push_back(TObj.eta());
0193 trigPhi.push_back(TObj.phi());
0194 }
0195
0196
0197 for (reco::IsolatedPixelTrackCandidateCollection::const_iterator itr = recoIsoTracks->begin();
0198 itr != recoIsoTracks->end();
0199 itr++) {
0200 bool match = false;
0201 for (unsigned int l = 0; l < trigEta.size(); l++) {
0202 if (getDist(itr->eta(), itr->phi(), trigEta[l], trigPhi[l]) < 0.2)
0203 match = true;
0204 }
0205 if (match && trig) {
0206 hOffEtaFP->Fill(itr->eta(), 1);
0207 hOffPhiFP->Fill(itr->phi(), 1);
0208 }
0209
0210 hOffEta->Fill(itr->eta(), 1);
0211 hOffPhi->Fill(itr->phi(), 1);
0212
0213 hOffAbsEta->Fill(fabs(itr->eta()), 1);
0214
0215 hDeposEcalInner->Fill(itr->energyIn(), 1);
0216 hDeposEcalOuter->Fill(itr->energyOut(), 1);
0217
0218 hTracksSumP->Fill(itr->sumPtPxl(), 1);
0219 hTracksMaxP->Fill(itr->maxPtPxl(), 1);
0220
0221 if (fabs(itr->eta()) < 0.5)
0222 hOffP_0005->Fill(itr->p(), 1);
0223 if (fabs(itr->eta()) > 0.5 && fabs(itr->eta()) < 1.0)
0224 hOffP_0510->Fill(itr->p(), 1);
0225 if (fabs(itr->eta()) > 1.0 && fabs(itr->eta()) < 1.5)
0226 hOffP_1015->Fill(itr->p(), 1);
0227 if (fabs(itr->eta()) < 1.5 && fabs(itr->eta()) < 2.0)
0228 hOffP_1520->Fill(itr->p(), 1);
0229
0230 hOffP->Fill(itr->p(), 1);
0231
0232 std::pair<int, int> TI = towerIndex(itr->eta(), itr->phi());
0233 hOccupancyFull->Fill(TI.first, TI.second, 1);
0234 if (itr->p() > heLow_ && itr->p() < heUp_)
0235 hOccupancyHighEn->Fill(TI.first, TI.second, 1);
0236 }
0237
0238
0239
0240 edm::LogVerbatim("HcalIsoTrack") << "\n End / Start ";
0241
0242 const edm::Handle<edm::SimTrackContainer>& simTracks = iEvent.getHandle(tok_simTrack_);
0243
0244 for (reco::IsolatedPixelTrackCandidateCollection::const_iterator bll = recoIsoTracks->begin();
0245 bll != recoIsoTracks->end();
0246 bll++) {
0247 edm::LogVerbatim("HcalIsoTrack") << "ISO Pt " << bll->pt() << " P " << bll->p() << " Eta " << bll->eta()
0248 << " Phi " << bll->phi();
0249
0250 double distanceMin = 1.;
0251 double SimPtMatched = 1.;
0252 double SimPhiMatched = 1.;
0253 double SimEtaMatched = 1.;
0254 double SimDistMatched = 1.;
0255 double SimPMatched = 1.;
0256 double neuen = 0.;
0257 double neuenm = 0.;
0258 int neun = 0;
0259
0260 for (edm::SimTrackContainer::const_iterator tracksCI = simTracks->begin(); tracksCI != simTracks->end();
0261 tracksCI++) {
0262 int partIndex = tracksCI->genpartIndex();
0263 if (tracksCI->momentum().eta() > (bll->eta() - 0.1) && tracksCI->momentum().eta() < (bll->eta() + 0.1) &&
0264 tracksCI->momentum().phi() > (bll->phi() - 0.1) &&
0265 tracksCI->momentum().phi() < (bll->phi() + 0.1)
0266
0267
0268 && tracksCI->momentum().e() > 2. && fabs(tracksCI->charge()) == 1 && partIndex > 0) {
0269 double distance = getDist(tracksCI->momentum().eta(), tracksCI->momentum().phi(), bll->eta(), bll->phi());
0270 double distanceCM = getDistInCM(tracksCI->momentum().eta(), tracksCI->momentum().phi(), bll->eta(), bll->phi());
0271
0272 if (distanceMin > distance) {
0273 distanceMin = distance;
0274 SimPtMatched = tracksCI->momentum().pt();
0275 SimPhiMatched = tracksCI->momentum().phi();
0276 SimEtaMatched = tracksCI->momentum().eta();
0277 SimDistMatched = distance;
0278 SimPMatched = sqrt(tracksCI->momentum().pt() * tracksCI->momentum().pt() +
0279 tracksCI->momentum().pz() * tracksCI->momentum().pz());
0280 }
0281
0282 edm::LogVerbatim("HcalIsoTrack") << " Pt " << tracksCI->momentum().pt() << " Energy "
0283 << tracksCI->momentum().e() << " Eta " << tracksCI->momentum().eta() << " Phi "
0284 << tracksCI->momentum().phi() << " Ind " << partIndex << " Cha "
0285 << tracksCI->charge() << " Dis " << distance << " DCM " << distanceCM;
0286 }
0287
0288 if (tracksCI->momentum().eta() > (bll->eta() - 0.5) && tracksCI->momentum().eta() < (bll->eta() + 0.5) &&
0289 tracksCI->momentum().phi() > (bll->phi() - 0.5) &&
0290 tracksCI->momentum().phi() < (bll->phi() + 0.5)
0291
0292 && tracksCI->charge() == 0 && partIndex > 0) {
0293 double distance = getDist(tracksCI->momentum().eta(), tracksCI->momentum().phi(), bll->eta(), bll->phi());
0294 double distanceCM = getDistInCM(tracksCI->momentum().eta(), tracksCI->momentum().phi(), bll->eta(), bll->phi());
0295
0296 edm::LogVerbatim("HcalIsoTrack") << "NEU Pt " << tracksCI->momentum().pt() << " Energy "
0297 << tracksCI->momentum().e() << " Eta " << tracksCI->momentum().eta() << " Phi "
0298 << tracksCI->momentum().phi() << " Ind " << partIndex << " Cha "
0299 << tracksCI->charge() << " Dis " << distance << " DCM " << distanceCM;
0300
0301 if (distanceCM < 40.) {
0302 neuen = neuen + tracksCI->momentum().e();
0303 neun = neun + 1;
0304 if (neuenm < tracksCI->momentum().e())
0305 neuenm = tracksCI->momentum().e();
0306 }
0307 }
0308 }
0309
0310 hSimNN->Fill(neun, 1);
0311 hSimNE->Fill(neuen, 1);
0312 hSimNM->Fill(neuenm, 1);
0313
0314 if (distanceMin < 0.1) {
0315 hSimPt->Fill(SimPtMatched, 1);
0316 hSimPhi->Fill(SimPhiMatched, 1);
0317 hSimEta->Fill(SimEtaMatched, 1);
0318 hSimAbsEta->Fill(fabs(SimEtaMatched), 1);
0319 hSimDist->Fill(SimDistMatched, 1);
0320 hSimPtRatOff->Fill(SimPtMatched / bll->pt(), 1);
0321 hSimP->Fill(SimPMatched, 1);
0322 hSimN->Fill(1, 1);
0323
0324 edm::LogVerbatim("HcalIsoTrack") << "S Pt " << SimPtMatched;
0325 }
0326
0327 if (distanceMin > 0.1) {
0328 hSimN->Fill(0, 1);
0329 }
0330 }
0331
0332
0333 }
0334
0335 void ValidationHcalIsoTrackAlCaReco::beginJob() {
0336 dbe_ = edm::Service<DQMStore>().operator->();
0337 dbe_->setCurrentFolder(folderName_);
0338
0339 hl3Pt = dbe_->book1D("hl3Pt", "pT of hlt L3 objects", 1000, 0, 1000);
0340 hl3Pt->setAxisTitle("pT(GeV)", 1);
0341 hl3eta = dbe_->book1D("hl3eta", "eta of hlt L3 objects", 16, -2, 2);
0342 hl3eta->setAxisTitle("eta", 1);
0343 hl3AbsEta = dbe_->book1D("hl3AbsEta", "|eta| of hlt L3 objects", 8, 0, 2);
0344 hl3AbsEta->setAxisTitle("eta", 1);
0345 hl3phi = dbe_->book1D("hl3phi", "phi of hlt L3 objects", 16, -3.2, 3.2);
0346 hl3phi->setAxisTitle("phi", 1);
0347
0348 hOffEta = dbe_->book1D("hOffEta", "eta of alcareco objects", 100, -2, 2);
0349 hOffEta->setAxisTitle("eta", 1);
0350 hOffPhi = dbe_->book1D("hOffPhi", "phi of alcareco objects", 100, -3.2, 3.2);
0351 hOffPhi->setAxisTitle("phi", 1);
0352 hOffP = dbe_->book1D("hOffP", "p of alcareco objects", 1000, 0, 1000);
0353 hOffP->setAxisTitle("E(GeV)", 1);
0354 hOffP_0005 = dbe_->book1D("hOffP_0005", "p of alcareco objects, |eta|<0.5", 1000, 0, 1000);
0355 hOffP_0005->setAxisTitle("E(GeV)", 1);
0356 hOffP_0510 = dbe_->book1D("hOffP_0510", "p of alcareco objects, 0.5<|eta|<1.0", 1000, 0, 1000);
0357 hOffP_0510->setAxisTitle("E(GeV)", 1);
0358 hOffP_1015 = dbe_->book1D("hOffP_1015", "p of alcareco objects, 1.0<|eta|<1.5", 1000, 0, 1000);
0359 hOffP_1015->setAxisTitle("E(GeV)", 1);
0360 hOffP_1520 = dbe_->book1D("hOffP_1520", "p of alcareco objects, 1.5<|eta|<2.0", 1000, 0, 1000);
0361 hOffP_1520->setAxisTitle("E(GeV)", 1);
0362 hOffEtaFP = dbe_->book1D("hOffEtaFP", "eta of alcareco objects, FP", 16, -2, 2);
0363 hOffEtaFP->setAxisTitle("eta", 1);
0364 hOffAbsEta = dbe_->book1D("hOffAbsEta", "|eta| of alcareco objects", 8, 0, 2);
0365 hOffAbsEta->setAxisTitle("|eta|", 1);
0366 hOffPhiFP = dbe_->book1D("hOffPhiFP", "phi of alcareco objects, FP", 16, -3.2, 3.2);
0367 hOffPhiFP->setAxisTitle("phi", 1);
0368 hTracksSumP = dbe_->book1D("hTracksSumP", "summary p of tracks in the isolation cone", 100, 0, 20);
0369 hTracksSumP->setAxisTitle("E(GeV)");
0370 hTracksMaxP = dbe_->book1D("hTracksMaxP", "maximum p among tracks in the isolation cone", 100, 0, 20);
0371 hTracksMaxP->setAxisTitle("E(GeV)");
0372
0373 hDeposEcalInner = dbe_->book1D("hDeposEcalInner", "ecal energy deposition in inner cone around track", 1000, 0, 1000);
0374 hDeposEcalInner->setAxisTitle("E(GeV)");
0375 hDeposEcalOuter = dbe_->book1D("hDeposEcalOuter", "ecal energy deposition in outer cone around track", 1000, 0, 1000);
0376 hDeposEcalOuter->setAxisTitle("E(GeV)");
0377
0378 hOccupancyFull =
0379 dbe_->book2D("hOccupancyFull", "number of tracks per tower, full energy range", 48, -25, 25, 73, 0, 73);
0380 hOccupancyFull->setAxisTitle("ieta", 1);
0381 hOccupancyFull->setAxisTitle("iphi", 2);
0382 hOccupancyFull->setOption("colz");
0383 hOccupancyFull->getTH2F()->SetStats(kFALSE);
0384 hOccupancyHighEn =
0385 dbe_->book2D("hOccupancyHighEn", "number of tracks per tower, high energy tracks", 48, -25, 25, 73, 0, 73);
0386 hOccupancyHighEn->setAxisTitle("ieta", 1);
0387 hOccupancyHighEn->setAxisTitle("iphi", 2);
0388 hOccupancyHighEn->setOption("colz");
0389 hOccupancyHighEn->getTH2F()->SetStats(kFALSE);
0390
0391 hOffL3TrackMatch = dbe_->book1D("hOffL3TrackMatch", "Distance from L3 object to offline track", 200, 0, 0.5);
0392 hOffL3TrackMatch->setAxisTitle("R(eta,phi)", 1);
0393 hOffL3TrackPtRat = dbe_->book1D("hOffL3TrackPtRat", "Ratio of pT: L3/offline", 100, 0, 10);
0394 hOffL3TrackPtRat->setAxisTitle("ratio L3/offline", 1);
0395
0396
0397
0398 hSimPt = dbe_->book1D("hSimPt", "pT of matched SimTrack", 1000, 0, 1000);
0399 hSimPt->setAxisTitle("pT(GeV)", 1);
0400
0401 hSimPhi = dbe_->book1D("hSimPhi", "Phi of matched SimTrack", 100, -3.2, 3.2);
0402 hSimPhi->setAxisTitle("phi", 1);
0403
0404 hSimEta = dbe_->book1D("hSimEta", "Eta of matched SimTrack", 100, -2., 2.);
0405 hSimEta->setAxisTitle("eta", 1);
0406
0407 hSimAbsEta = dbe_->book1D("hSimAbsEta", "|eta| of matched SimTrack", 8, 0., 2.);
0408 hSimAbsEta->setAxisTitle("|eta|", 1);
0409
0410 hSimDist = dbe_->book1D("hSimDist", "Distance from matched SimTrack to Offline Track", 200, 0, 0.1);
0411 hSimDist->setAxisTitle("R(eta,phi)", 1);
0412
0413 hSimPtRatOff = dbe_->book1D("hSimPtRatOff", "pT Sim / pT Offline", 100, 0, 10);
0414 hSimPtRatOff->setAxisTitle("pT Sim / pT Offline", 1);
0415
0416 hSimP = dbe_->book1D("hSimP", "p of matched SimTrack", 1000, 0, 1000);
0417 hSimP->setAxisTitle("p(GeV)", 1);
0418
0419 hSimN = dbe_->book1D("hSimN", "Number matched", 2, 0, 2);
0420 hSimN->setAxisTitle("Offline/SimTRack - Matched or Not", 1);
0421
0422 hSimNN = dbe_->book1D("hSimNN", "Number of the neutral particles in cone on ECAL", 100, 0, 100);
0423 hSimNN->setAxisTitle("Number", 1);
0424
0425 hSimNE = dbe_->book1D("hSimNE", "Total energy of the neutral particles in cone on ECAL", 100, 0, 100);
0426 hSimNE->setAxisTitle("Energy", 1);
0427
0428 hSimNM = dbe_->book1D("hSimNM", "Maximum energy of the neutral particles in cone on ECAL", 100, 0, 100);
0429 hSimNM->setAxisTitle("Energy", 1);
0430
0431
0432 }
0433
0434 void ValidationHcalIsoTrackAlCaReco::endJob() {
0435 if (dbe_) {
0436 if (saveToFile_)
0437 dbe_->save(outRootFileName_);
0438 }
0439 }
0440
0441 #include "FWCore/Framework/interface/MakerMacros.h"
0442 DEFINE_FWK_MODULE(ValidationHcalIsoTrackAlCaReco);
