File indexing completed on 2022-02-28 01:32:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <memory>
0023
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0032
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0035 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0036 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0037 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0038 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0039 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0040 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0041
0042 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0043 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0044 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0045 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0046 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0047 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0049 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0050 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0051
0052 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
0053
0054 #include "TROOT.h"
0055 #include "TFile.h"
0056 #include "TTree.h"
0057 #include "TH1F.h"
0058
0059 #include <fstream>
0060 #include <map>
0061
0062
0063
0064 class ValidIsoTrkCalib : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0065 public:
0066 explicit ValidIsoTrkCalib(const edm::ParameterSet&);
0067
0068
0069
0070 private:
0071 void beginJob() override;
0072 void analyze(const edm::Event&, const edm::EventSetup&) override;
0073
0074
0075
0076
0077
0078 TrackDetectorAssociator trackAssociator_;
0079 TrackAssociatorParameters parameters_;
0080 double taECALCone_;
0081 double taHCALCone_;
0082
0083 const CaloGeometry* geo;
0084
0085 const bool takeGenTracks_;
0086 const edm::InputTag genhbheLabel_;
0087 const double associationConeSize_;
0088 const bool allowMissingInputs_;
0089 const std::string AxB_;
0090 const double calibrationConeSize_;
0091 const int MinNTrackHitsBarrel;
0092 const int MinNTECHitsEndcap;
0093 const double energyECALmip;
0094 const double energyMinIso;
0095 const double energyMaxIso;
0096 const double maxPNear;
0097
0098 const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0099 const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0100 const edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0101 const edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_track_;
0102 const edm::EDGetTokenT<reco::TrackCollection> tok_track1_;
0103
0104 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_recalibCorrs_;
0105 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0106
0107 edm::InputTag genhoLabel_;
0108 std::vector<edm::InputTag> genecalLabel_;
0109 std::string outputFileName_;
0110
0111 int gen, iso, pix;
0112 float genPt[500], genPhi[500], genEta[500];
0113 float isoPt[500], isoPhi[500], isoEta[500];
0114 float pixPt[500], pixPhi[500], pixEta[500];
0115
0116 int NisoTrk;
0117 float trackPt, trackE, trackEta, trackPhi;
0118 float ptNear;
0119 float ptrack, rvert;
0120
0121
0122 Float_t emEnergy;
0123 Float_t targetE;
0124
0125 TTree *tTree, *fTree;
0126
0127 Float_t xTrkEcal;
0128 Float_t yTrkEcal;
0129 Float_t zTrkEcal;
0130
0131 Float_t xTrkHcal;
0132 Float_t yTrkHcal;
0133 Float_t zTrkHcal;
0134
0135 int Nhits;
0136 float eClustBefore;
0137 float eClustAfter;
0138 float eTrack;
0139 float etaTrack;
0140 float phiTrack;
0141 float eECAL;
0142 int numHits;
0143
0144 float eBeforeDepth1;
0145 float eAfterDepth1;
0146 float eBeforeDepth2;
0147 float eAfterDepth2;
0148 float eCentHitBefore;
0149 float eCentHitAfter;
0150 float CentHitFactor;
0151 int iEta;
0152 int iPhi;
0153 int iEtaTr;
0154 int iPhiTr;
0155 float iDr, delR;
0156 int dietatr;
0157 int diphitr;
0158
0159 float iTime;
0160 float HTime[100];
0161 float e3x3Before;
0162 float e3x3After;
0163 float e5x5Before;
0164 float e5x5After;
0165 int eventNumber;
0166 int runNumber;
0167 float PtNearBy;
0168 float numVH, numVS, numValidTrkHits, numValidTrkStrips;
0169
0170 const HcalRespCorrs* respRecalib;
0171
0172 TH1F* nTracks;
0173 };
0174
0175 ValidIsoTrkCalib::ValidIsoTrkCalib(const edm::ParameterSet& iConfig)
0176 : takeGenTracks_(iConfig.getUntrackedParameter<bool>("takeGenTracks")),
0177 genhbheLabel_(iConfig.getParameter<edm::InputTag>("genHBHE")),
0178 associationConeSize_(iConfig.getParameter<double>("associationConeSize")),
0179 allowMissingInputs_(iConfig.getUntrackedParameter<bool>("allowMissingInputs", true)),
0180 AxB_(iConfig.getParameter<std::string>("AxB")),
0181 calibrationConeSize_(iConfig.getParameter<double>("calibrationConeSize")),
0182 MinNTrackHitsBarrel(iConfig.getParameter<int>("MinNTrackHitsBarrel")),
0183 MinNTECHitsEndcap(iConfig.getParameter<int>("MinNTECHitsEndcap")),
0184 energyECALmip(iConfig.getParameter<double>("energyECALmip")),
0185 energyMinIso(iConfig.getParameter<double>("energyMinIso")),
0186 energyMaxIso(iConfig.getParameter<double>("energyMaxIso")),
0187 maxPNear(iConfig.getParameter<double>("maxPNear")),
0188 tok_genTrack_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("genTracksLabel"))),
0189 tok_hbhe_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"))),
0190 tok_ho_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"))),
0191 tok_track_(consumes<reco::IsolatedPixelTrackCandidateCollection>(
0192 iConfig.getParameter<edm::InputTag>("HcalIsolTrackInput"))),
0193 tok_track1_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trackInput"))),
0194 tok_recalibCorrs_(esConsumes(edm::ESInputTag("", "recalibrate"))),
0195 tok_geom_(esConsumes()) {
0196 usesResource(TFileService::kSharedResource);
0197
0198 edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0199 edm::ConsumesCollector iC = consumesCollector();
0200 parameters_.loadParameters(parameters, iC);
0201 trackAssociator_.useDefaultPropagator();
0202 }
0203
0204
0205 void ValidIsoTrkCalib::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0206 try {
0207 respRecalib = &iSetup.getData(tok_recalibCorrs_);
0208
0209 edm::LogVerbatim("ValidIsoTrkCalib") << " Loaded: OK ";
0210
0211 } catch (const cms::Exception& e) {
0212 edm::LogWarning("ValidIsoTrkCalib") << " Not Found!! ";
0213 }
0214
0215 const edm::Handle<reco::TrackCollection>& generalTracks = iEvent.getHandle(tok_genTrack_);
0216
0217 const edm::Handle<reco::TrackCollection>& isoProdTracks = iEvent.getHandle(tok_track1_);
0218
0219 const edm::Handle<reco::IsolatedPixelTrackCandidateCollection>& isoPixelTracks = iEvent.getHandle(tok_track_);
0220
0221
0222
0223
0224
0225
0226
0227
0228 const edm::Handle<HBHERecHitCollection>& hbhe = iEvent.getHandle(tok_hbhe_);
0229 const HBHERecHitCollection Hithbhe = *(hbhe.product());
0230
0231 geo = &iSetup.getData(tok_geom_);
0232
0233 const HcalGeometry* gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0234
0235
0236
0237 parameters_.useEcal = true;
0238 parameters_.useHcal = true;
0239 parameters_.useCalo = false;
0240 parameters_.useMuon = false;
0241
0242
0243 #ifdef EDM_ML_DEBUG
0244 edm::LogVerbatim("ValidIsoTrkCalib") << "Hello World. TrackCollectionSize: " << isoPixelTracks->size();
0245 #endif
0246 if (isoPixelTracks->empty())
0247 return;
0248
0249 for (reco::TrackCollection::const_iterator trit = isoProdTracks->begin(); trit != isoProdTracks->end(); trit++) {
0250 reco::IsolatedPixelTrackCandidateCollection::const_iterator isoMatched = isoPixelTracks->begin();
0251
0252 bool matched = false;
0253
0254
0255 for (reco::IsolatedPixelTrackCandidateCollection::const_iterator it = isoPixelTracks->begin();
0256 it != isoPixelTracks->end();
0257 it++)
0258
0259 {
0260 if (abs((trit->pt() - it->pt()) / it->pt()) < 0.005 && abs(trit->eta() - it->eta()) < 0.01) {
0261 isoMatched = it;
0262 matched = true;
0263 break;
0264 }
0265 }
0266
0267
0268 if (!matched)
0269 continue;
0270 if (isoMatched->maxPtPxl() > maxPNear)
0271 continue;
0272
0273 ptNear = isoMatched->maxPtPxl();
0274 #ifdef EDM_ML_DEBUG
0275 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.1 isoMatch. ptnear: " << ptNear;
0276 #endif
0277
0278 if (trit->hitPattern().numberOfValidHits() < MinNTrackHitsBarrel)
0279 continue;
0280 if (fabs(trit->eta()) > 1.47 && trit->hitPattern().numberOfValidStripTECHits() < MinNTECHitsEndcap)
0281 continue;
0282
0283 #ifdef EDM_ML_DEBUG
0284 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.2.1 after numofvalidhits HB: "
0285 << trit->hitPattern().numberOfValidHits();
0286 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.2.2 after numofvalidstrips HE: "
0287 << trit->hitPattern().numberOfValidStripTECHits();
0288 #endif
0289 numVH = trit->hitPattern().numberOfValidHits();
0290 numVS = trit->hitPattern().numberOfValidStripTECHits();
0291
0292 trackE = sqrt(trit->px() * trit->px() + trit->py() * trit->py() + trit->pz() * trit->pz() + 0.14 * 0.14);
0293 trackPt = trit->pt();
0294 trackEta = trit->eta();
0295 trackPhi = trit->phi();
0296
0297 emEnergy = isoMatched->energyIn();
0298
0299 #ifdef EDM_ML_DEBUG
0300 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.3. Matched :: pt: " << trit->pt() << " wholeEnergy: " << trackE
0301 << " emEnergy: " << emEnergy << " eta: " << trackEta << " phi: " << trackPhi;
0302 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.4. EM energy in cone: " << emEnergy;
0303 #endif
0304 TrackDetMatchInfo info = trackAssociator_.associate(
0305 iEvent,
0306 iSetup,
0307 trackAssociator_.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *trit),
0308 parameters_);
0309
0310
0311
0312
0313
0314
0315 xTrkEcal = info.trkGlobPosAtEcal.x();
0316 yTrkEcal = info.trkGlobPosAtEcal.y();
0317 zTrkEcal = info.trkGlobPosAtEcal.z();
0318
0319 xTrkHcal = info.trkGlobPosAtHcal.x();
0320 yTrkHcal = info.trkGlobPosAtHcal.y();
0321 zTrkHcal = info.trkGlobPosAtHcal.z();
0322
0323 if (xTrkEcal == 0 && yTrkEcal == 0 && zTrkEcal == 0) {
0324 edm::LogVerbatim("ValidIsoTrkCalib") << "zero point at Ecal";
0325 continue;
0326 }
0327 if (xTrkHcal == 0 && yTrkHcal == 0 && zTrkHcal == 0) {
0328 edm::LogVerbatim("ValidIsoTrkCalib") << "zero point at Hcal";
0329 continue;
0330 }
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340 GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0341 GlobalPoint gPointHcal(xTrkHcal, yTrkHcal, zTrkHcal);
0342
0343 int iphitrue = -10;
0344 int ietatrue = 100;
0345 const HcalDetId tempId = gHcal->getClosestCell(gPointHcal);
0346 ietatrue = tempId.ieta();
0347 iphitrue = tempId.iphi();
0348
0349 MaxHit_struct MaxHit;
0350
0351 MaxHit.hitenergy = -100.;
0352
0353
0354 std::vector<DetId> usedHits;
0355
0356 usedHits.clear();
0357 #ifdef EDM_ML_DEBUG
0358 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 1. Entrance to HBHECollection";
0359 #endif
0360
0361
0362
0363
0364
0365 GlobalPoint gPhot;
0366
0367 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++) {
0368
0369 bool hitIsUsed = false;
0370 for (uint32_t i = 0; i < usedHits.size(); i++) {
0371 if (usedHits[i] == hhit->id())
0372 hitIsUsed = true;
0373 }
0374 if (hitIsUsed)
0375 continue;
0376 usedHits.push_back(hhit->id());
0377
0378
0379
0380 float recal = 1;
0381
0382
0383 GlobalPoint pos = geo->getPosition(hhit->detid());
0384
0385
0386
0387 int iphihitm = (hhit->id()).iphi();
0388 int ietahitm = (hhit->id()).ieta();
0389 int depthhit = (hhit->id()).depth();
0390 float enehit = hhit->energy() * recal;
0391
0392 if (depthhit != 1)
0393 continue;
0394
0395 #ifdef EDM_ML_DEBUG
0396 float dphi = fabs(info.trkGlobPosAtHcal.phi() - pos.phi());
0397 if (dphi > 4. * atan(1.))
0398 dphi = 8. * atan(1.) - dphi;
0399 float deta = fabs(info.trkGlobPosAtHcal.eta() - pos.eta());
0400 float dr = sqrt(dphi * dphi + deta * deta);
0401 #endif
0402
0403
0404 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0405
0406 if (distAtHcal < associationConeSize_) {
0407 for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0408 int iphihitm2 = (hhit2->id()).iphi();
0409 int ietahitm2 = (hhit2->id()).ieta();
0410 int depthhit2 = (hhit2->id()).depth();
0411 float enehit2 = hhit2->energy() * recal;
0412
0413 if (iphihitm == iphihitm2 && ietahitm == ietahitm2 && depthhit != depthhit2)
0414 enehit = enehit + enehit2;
0415 }
0416
0417 #ifdef EDM_ML_DEBUG
0418 edm::LogVerbatim("ValidIsoTrkCalib") << "IN CONE ieta: " << ietahitm << " iphi: " << iphihitm
0419 << " depthhit: " << depthhit << " dr: " << dr << " energy: " << enehit;
0420 #endif
0421
0422
0423 if (enehit > MaxHit.hitenergy) {
0424 MaxHit.hitenergy = enehit;
0425 MaxHit.ietahitm = (hhit->id()).ieta();
0426 MaxHit.iphihitm = (hhit->id()).iphi();
0427 MaxHit.dr = distAtHcal;
0428
0429 MaxHit.depthhit = 1;
0430
0431
0432 }
0433 }
0434 }
0435
0436 usedHits.clear();
0437
0438 #ifdef EDM_ML_DEBUG
0439 edm::LogVerbatim("ValidIsoTrkCalib") << "Hottest ieta: " << MaxHit.ietahitm << " iphi: " << MaxHit.iphihitm
0440 << " dr: " << MaxHit.dr;
0441 edm::LogVerbatim("ValidIsoTrkCalib") << "Track ieta: " << ietatrue << " iphi: " << iphitrue;
0442 edm::LogVerbatim("ValidIsoTrkCalib") << "Point 3. MaxHit :::En " << MaxHit.hitenergy
0443 << " ieta: " << MaxHit.ietahitm << " iphi: " << MaxHit.iphihitm;
0444 #endif
0445
0446 Bool_t passCuts = kFALSE;
0447 if (trackE > energyMinIso && trackE < energyMaxIso && emEnergy < energyECALmip && MaxHit.hitenergy > 0. &&
0448 abs(MaxHit.ietahitm) < 29)
0449 passCuts = kTRUE;
0450
0451 #ifdef EDM_ML_DEBUG
0452 edm::LogVerbatim("ValidIsoTrkCalib") << "Pont 0.1.1. trackE:" << trackE << " emEn: " << emEnergy;
0453 #endif
0454 numHits = 0;
0455
0456 eClustBefore = 0.0;
0457 eClustAfter = 0.0;
0458 eBeforeDepth1 = 0.0;
0459 eAfterDepth1 = 0.0;
0460 eBeforeDepth2 = 0.0;
0461 eAfterDepth2 = 0.0;
0462 CentHitFactor = 0.0;
0463 e3x3After = 0.0;
0464 e3x3Before = 0.0;
0465 e5x5After = 0.0;
0466 e5x5Before = 0.0;
0467
0468 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++) {
0469
0470 bool hitIsUsed = false;
0471 for (uint32_t i = 0; i < usedHits.size(); i++) {
0472 if (usedHits[i] == hhit->id())
0473 hitIsUsed = true;
0474 }
0475 if (hitIsUsed)
0476 continue;
0477 usedHits.push_back(hhit->id());
0478
0479 int DIETA = 100;
0480 if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0481 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0482 }
0483 if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0484 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0485 DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0486 }
0487
0488 int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0489 DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0490
0491 int numbercell = 100;
0492
0493
0494
0495
0496
0497
0498
0499 if (abs(DIETA) <= numbercell &&
0500 (abs(DIPHI) <= numbercell || (abs(MaxHit.ietahitm) >= 20 && abs(DIPHI) <= numbercell + 1))) {
0501 const GlobalPoint pos2 = geo->getPosition(hhit->detid());
0502
0503 if (passCuts && hhit->energy() > 0) {
0504 float factor = 0.0;
0505
0506 factor = respRecalib->getValues(hhit->id())->getValue();
0507
0508 #ifdef EDM_ML_DEBUG
0509 edm::LogVerbatim("ValidIsoTrkCalib") << " calib factors: " << factor;
0510 #endif
0511
0512 if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm)
0513 CentHitFactor = factor;
0514
0515 if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm)
0516 iTime = hhit->time();
0517
0518 if (AxB_ != "3x3" && AxB_ != "5x5" && AxB_ != "Cone")
0519 edm::LogWarning(" AxB ") << " Not supported: " << AxB_;
0520
0521 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0522 !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0523 abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0524 e5x5Before += hhit->energy();
0525 e5x5After += hhit->energy() * factor;
0526 }
0527
0528 if (abs(DIETA) <= 1 && (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0529 !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 &&
0530 abs(DIPHI) > 1)))) {
0531 e3x3Before += hhit->energy();
0532 e3x3After += hhit->energy() * factor;
0533 }
0534
0535 if (AxB_ == "5x5") {
0536 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || (abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4))) {
0537 if (abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 3)
0538 continue;
0539
0540 HTime[numHits] = hhit->time();
0541 numHits++;
0542
0543 eClustBefore += hhit->energy();
0544 eClustAfter += hhit->energy() * factor;
0545
0546 if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0547 eBeforeDepth1 += hhit->energy();
0548 eAfterDepth1 += hhit->energy() * factor;
0549 } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0550 eBeforeDepth2 += hhit->energy();
0551 eAfterDepth2 += hhit->energy() * factor;
0552 }
0553 }
0554 }
0555
0556 if (AxB_ == "3x3") {
0557 if (abs(DIETA) <= 1 && (abs(DIPHI) <= 1 || (abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2))) {
0558 if (abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)
0559 continue;
0560
0561 HTime[numHits] = hhit->time();
0562 numHits++;
0563
0564 eClustBefore += hhit->energy();
0565 eClustAfter += hhit->energy() * factor;
0566
0567 if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0568 eBeforeDepth1 += hhit->energy();
0569 eAfterDepth1 += hhit->energy() * factor;
0570 } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0571 eBeforeDepth2 += hhit->energy();
0572 eAfterDepth2 += hhit->energy() * factor;
0573 }
0574 }
0575 }
0576
0577 if (AxB_ == "Cone" && getDistInPlaneSimple(gPointHcal, pos2) < calibrationConeSize_) {
0578 HTime[numHits] = hhit->time();
0579 numHits++;
0580
0581 eClustBefore += hhit->energy();
0582 eClustAfter += hhit->energy() * factor;
0583
0584 if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0585 eBeforeDepth1 += hhit->energy();
0586 eAfterDepth1 += hhit->energy() * factor;
0587 } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0588 eBeforeDepth2 += hhit->energy();
0589 eAfterDepth2 += hhit->energy() * factor;
0590 }
0591
0592 }
0593
0594 }
0595
0596 }
0597
0598 }
0599
0600 int dieta_M_P = 100;
0601 int diphi_M_P = 100;
0602 if (MaxHit.ietahitm * ietatrue > 0) {
0603 dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0604 }
0605 if (MaxHit.ietahitm * ietatrue < 0) {
0606 dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0607 }
0608 diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0609 diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0610
0611 if (passCuts)
0612
0613 {
0614 eventNumber = iEvent.id().event();
0615 runNumber = iEvent.id().run();
0616
0617 eCentHitBefore = MaxHit.hitenergy;
0618 eCentHitAfter = MaxHit.hitenergy * CentHitFactor;
0619 eECAL = emEnergy;
0620 numValidTrkHits = numVH;
0621 numValidTrkStrips = numVS;
0622 PtNearBy = ptNear;
0623
0624 eTrack = trackE;
0625 phiTrack = trackPhi;
0626 etaTrack = trackEta;
0627
0628 iEta = MaxHit.ietahitm;
0629 iPhi = MaxHit.iphihitm;
0630
0631 iEtaTr = ietatrue;
0632 iPhiTr = iphitrue;
0633 iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0634 delR = MaxHit.dr;
0635 dietatr = dieta_M_P;
0636 diphitr = diphi_M_P;
0637
0638 fTree->Fill();
0639 }
0640
0641 }
0642
0643
0644 #ifdef EDM_ML_DEBUG
0645 edm::LogVerbatim("ValidIsoTrkCalib") << " generalTracks Size: " << generalTracks->size();
0646 #endif
0647 int n = generalTracks->size();
0648 nTracks->Fill(n);
0649
0650 if (takeGenTracks_ && iEvent.id().event() % 10 == 1) {
0651 gen = generalTracks->size();
0652 iso = isoProdTracks->size();
0653 pix = isoPixelTracks->size();
0654
0655 genPt[0] = -33;
0656 genPhi[0] = -33;
0657 genEta[0] = -33;
0658
0659 isoPt[0] = -33;
0660 isoPhi[0] = -33;
0661 isoEta[0] = -33;
0662
0663 pixPt[0] = -33;
0664 pixPhi[0] = -33;
0665 pixEta[0] = -33;
0666
0667 Int_t gencount = 0, isocount = 0, pixcount = 0;
0668 for (reco::TrackCollection::const_iterator gentr = generalTracks->begin(); gentr != generalTracks->end(); gentr++) {
0669 genPt[gencount] = gentr->pt();
0670 genPhi[gencount] = gentr->phi();
0671 genEta[gencount] = gentr->eta();
0672 gencount++;
0673 }
0674
0675 for (reco::TrackCollection::const_iterator isotr = isoProdTracks->begin(); isotr != isoProdTracks->end(); isotr++) {
0676 isoPt[isocount] = isotr->pt();
0677 isoPhi[isocount] = isotr->phi();
0678 isoEta[isocount] = isotr->eta();
0679 isocount++;
0680 }
0681
0682 for (reco::IsolatedPixelTrackCandidateCollection::const_iterator pixtr = isoPixelTracks->begin();
0683 pixtr != isoPixelTracks->end();
0684 pixtr++) {
0685 pixPt[pixcount] = pixtr->pt();
0686 pixPhi[pixcount] = pixtr->phi();
0687 pixEta[pixcount] = pixtr->eta();
0688 pixcount++;
0689 }
0690 }
0691
0692 tTree->Fill();
0693 }
0694
0695
0696 void ValidIsoTrkCalib::beginJob() {
0697 edm::Service<TFileService> fs;
0698
0699 nTracks = fs->make<TH1F>("nTracks", "general;number of general tracks", 11, -0.5, 10.5);
0700
0701 tTree = fs->make<TTree>("tTree", "Tree for gen info");
0702
0703 fTree = fs->make<TTree>("fTree", "Tree for IsoTrack Calibration");
0704
0705 fTree->Branch("eventNumber", &eventNumber, "eventNumber/I");
0706 fTree->Branch("runNumber", &runNumber, "runNumber/I");
0707
0708 fTree->Branch("eClustBefore", &eClustBefore, "eClustBefore/F");
0709 fTree->Branch("eClustAfter", &eClustAfter, "eClustAfter/F");
0710 fTree->Branch("eTrack", &eTrack, "eTrack/F");
0711 fTree->Branch("etaTrack", &etaTrack, "etaTrack/F");
0712 fTree->Branch("phiTrack", &phiTrack, "phiTrack/F");
0713
0714 fTree->Branch("numHits", &numHits, "numHits/I");
0715 fTree->Branch("eECAL", &eECAL, "eECAL/F");
0716 fTree->Branch("PtNearBy", &PtNearBy, "PtNearBy/F");
0717 fTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/F");
0718 fTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/F");
0719
0720 fTree->Branch("eBeforeDepth1", &eBeforeDepth1, "eBeforeDepth1/F");
0721 fTree->Branch("eBeforeDepth2", &eBeforeDepth2, "eBeforeDepth2/F");
0722 fTree->Branch("eAfterDepth1", &eAfterDepth1, "eAfterDepth1/F");
0723 fTree->Branch("eAfterDepth2", &eAfterDepth2, "eAfterDepth2/F");
0724
0725 fTree->Branch("e3x3Before", &e3x3Before, "e3x3Before/F");
0726 fTree->Branch("e3x3After", &e3x3After, "e3x3After/F");
0727 fTree->Branch("e5x5Before", &e5x5Before, "e5x5Before/F");
0728 fTree->Branch("e5x5After", &e5x5After, "e5x5After/F");
0729
0730 fTree->Branch("eCentHitAfter", &eCentHitAfter, "eCentHitAfter/F");
0731 fTree->Branch("eCentHitBefore", &eCentHitBefore, "eCentHitBefore/F");
0732 fTree->Branch("iEta", &iEta, "iEta/I");
0733 fTree->Branch("iPhi", &iPhi, "iPhi/I");
0734
0735 fTree->Branch("iEtaTr", &iEtaTr, "iEtaTr/I");
0736 fTree->Branch("iPhiTr", &iPhiTr, "iPhiTr/I");
0737 fTree->Branch("dietatr", &dietatr, "dietatr/I");
0738 fTree->Branch("diphitr", &diphitr, "diphitr/I");
0739 fTree->Branch("iDr", &iDr, "iDr/F");
0740 fTree->Branch("delR", &delR, "delR/F");
0741
0742 fTree->Branch("iTime", &iTime, "iTime/F");
0743 fTree->Branch("HTime", HTime, "HTime[numHits]/F");
0744
0745 fTree->Branch("xTrkEcal", &xTrkEcal, "xTrkEcal/F");
0746 fTree->Branch("yTrkEcal", &yTrkEcal, "yTrkEcal/F");
0747 fTree->Branch("zTrkEcal", &zTrkEcal, "zTrkEcal/F");
0748 fTree->Branch("xTrkHcal", &xTrkHcal, "xTrkHcal/F");
0749 fTree->Branch("yTrkHcal", &yTrkHcal, "yTrkHcal/F");
0750 fTree->Branch("zTrkHcal", &zTrkHcal, "zTrkHcal/F");
0751
0752 if (takeGenTracks_) {
0753 tTree->Branch("gen", &gen, "gen/I");
0754 tTree->Branch("iso", &iso, "iso/I");
0755 tTree->Branch("pix", &pix, "pix/I");
0756 tTree->Branch("genPt", genPt, "genPt[gen]/F");
0757 tTree->Branch("genPhi", genPhi, "genPhi[gen]/F");
0758 tTree->Branch("genEta", genEta, "genEta[gen]/F");
0759
0760 tTree->Branch("isoPt", isoPt, "isoPt[iso]/F");
0761 tTree->Branch("isoPhi", isoPhi, "isoPhi[iso]/F");
0762 tTree->Branch("isoEta", isoEta, "isoEta[iso]/F");
0763
0764 tTree->Branch("pixPt", pixPt, "pixPt[pix]/F");
0765 tTree->Branch("pixPhi", pixPhi, "pixPhi[pix]/F");
0766 tTree->Branch("pixEta", pixEta, "pixEta[pix]/F");
0767 }
0768 }
0769
0770
0771 DEFINE_FWK_MODULE(ValidIsoTrkCalib);