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