File indexing completed on 2022-02-28 01:32:06
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009
0010 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0011 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0012
0013 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0014 #include "DataFormats/HcalRecHit/interface/HcalSourcePositionData.h"
0015 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0017
0018 #include "DataFormats/HcalDetId/interface/HcalElectronicsId.h"
0019 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0022
0023 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0024 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0025 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0026 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0027 #include "CondFormats/HcalObjects/interface/HcalPFCorrs.h"
0028 #include "CondFormats/DataRecord/interface/HcalPFCorrsRcd.h"
0029 #include "MagneticField/Engine/interface/MagneticField.h"
0030 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0031 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0032 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0033 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0034 #include "DataFormats/GeometrySurface/interface/Plane.h"
0035
0036 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
0037
0038 #include <iostream>
0039 #include "TProfile.h"
0040 #include "TTree.h"
0041
0042
0043
0044 class HcalCorrPFCalculation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0045 public:
0046 HcalCorrPFCalculation(edm::ParameterSet const& conf);
0047 void analyze(edm::Event const& ev, edm::EventSetup const& c) override;
0048 void beginJob() override;
0049
0050 private:
0051 double RecalibFactor(HcalDetId id);
0052
0053 const bool Respcorr_;
0054 const bool PFcorr_;
0055 const bool Conecorr_;
0056
0057
0058 const double clusterConeSize_, associationConeSize_, ecalCone_, neutralIsolationCone_, trackIsolationCone_;
0059 float eECAL, eECAL09cm, eECAL40cm;
0060
0061
0062 float eCentHit, eTrack, eParticle;
0063 float CentHitFactor;
0064 int iEta, iPhi;
0065
0066 float iDr, delR;
0067 float e3x3, e5x5;
0068 float phiTrack, etaTrack;
0069 float phiGPoint, etaGPoint;
0070
0071 Bool_t doHF;
0072 Bool_t AddRecalib;
0073 int nevtot;
0074
0075 const HcalRespCorrs* respRecalib;
0076 const HcalPFCorrs* pfRecalib;
0077
0078 SteppingHelixPropagator* stepPropF;
0079 MagneticField* theMagField;
0080
0081 TrackDetectorAssociator trackAssociator_;
0082 TrackAssociatorParameters parameters_;
0083
0084 const CaloGeometry* geo;
0085 const HcalGeometry* gHcal;
0086
0087 Float_t xTrkHcal, yTrkHcal, zTrkHcal;
0088 Float_t xTrkEcal, yTrkEcal, zTrkEcal;
0089 Float_t xAtHcal, yAtHcal, zAtHcal;
0090
0091 float eEcalCone, eHcalCone, eHcalConeNoise;
0092 int UsedCells, UsedCellsNoise;
0093 float phiParticle, etaParticle;
0094
0095 Int_t numValidTrkHits[50], numValidTrkStrips[50], numLayers[50];
0096
0097 Bool_t trkQual[50];
0098
0099 TProfile *nCells, *nCellsNoise, *enHcal, *enHcalNoise;
0100
0101 TTree *pfTree, *tracksTree;
0102
0103 edm::Service<TFileService> fs;
0104 Int_t nTracks;
0105 Float_t genEta, genPhi, trackEta[50], trackPhi[50], trackP[50], delRmc[50];
0106
0107 const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0108 const edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0109 const edm::EDGetTokenT<HFRecHitCollection> tok_hf_;
0110
0111 const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0112 const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0113 const edm::EDGetTokenT<reco::TrackCollection> tok_tracks_;
0114 const edm::EDGetTokenT<edm::HepMCProduct> tok_gen_;
0115
0116 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0117 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0118 const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0119 const edm::ESGetToken<HcalPFCorrs, HcalPFCorrsRcd> tok_pfcorr_;
0120 };
0121
0122 HcalCorrPFCalculation::HcalCorrPFCalculation(edm::ParameterSet const& iConfig)
0123 : Respcorr_(iConfig.getUntrackedParameter<bool>("RespcorrAdd", false)),
0124 PFcorr_(iConfig.getUntrackedParameter<bool>("PFcorrAdd", false)),
0125 Conecorr_(iConfig.getUntrackedParameter<bool>("ConeCorrAdd", true)),
0126 clusterConeSize_(iConfig.getParameter<double>("clusterConeSize")),
0127 associationConeSize_(iConfig.getParameter<double>("associationConeSize")),
0128 ecalCone_(iConfig.getParameter<double>("ecalCone")),
0129 neutralIsolationCone_(iConfig.getParameter<double>("neutralIsolationCone")),
0130 trackIsolationCone_(iConfig.getParameter<double>("trackIsolationCone")),
0131 tok_hbhe_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheRecHitCollectionTag"))),
0132 tok_ho_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoRecHitCollectionTag"))),
0133 tok_hf_(consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfRecHitCollectionTag"))),
0134 tok_EE_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"))),
0135 tok_EB_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"))),
0136 tok_tracks_(consumes<reco::TrackCollection>(edm::InputTag("generalTracks"))),
0137 tok_gen_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0138 tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0139 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0140 tok_resp_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0141 tok_pfcorr_(esConsumes<HcalPFCorrs, HcalPFCorrsRcd>()) {
0142 usesResource(TFileService::kSharedResource);
0143
0144
0145
0146
0147
0148
0149
0150
0151 edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0152 edm::ConsumesCollector iC = consumesCollector();
0153 parameters_.loadParameters(parameters, iC);
0154 trackAssociator_.useDefaultPropagator();
0155
0156
0157 }
0158
0159 double HcalCorrPFCalculation::RecalibFactor(HcalDetId id) {
0160 Float_t resprecal = 1.;
0161 Float_t pfrecal = 1.;
0162 if (AddRecalib) {
0163 if (Respcorr_)
0164 resprecal = respRecalib->getValues(id)->getValue();
0165 if (PFcorr_)
0166 pfrecal = pfRecalib->getValues(id)->getValue();
0167 }
0168 Float_t factor = resprecal * pfrecal;
0169 return factor;
0170 }
0171
0172 void HcalCorrPFCalculation::analyze(edm::Event const& ev, edm::EventSetup const& c) {
0173 AddRecalib = kFALSE;
0174
0175 try {
0176 respRecalib = &c.getData(tok_resp_);
0177 pfRecalib = &c.getData(tok_pfcorr_);
0178
0179 AddRecalib = kTRUE;
0180 #ifdef EDM_ML_DEBUG
0181 edm::LogVerbatim("CalibConstants") << " OK ";
0182 #endif
0183
0184 } catch (const cms::Exception& e) {
0185 edm::LogWarning("CalibConstants") << " Not Found!! ";
0186 }
0187
0188 const HBHERecHitCollection& Hithbhe = ev.get(tok_hbhe_);
0189 const HFRecHitCollection& Hithf = ev.get(tok_hf_);
0190
0191
0192 const edm::Handle<EERecHitCollection>& ecalEE = ev.getHandle(tok_EE_);
0193 const EERecHitCollection HitecalEE = *(ecalEE.product());
0194
0195 const edm::Handle<EBRecHitCollection>& ecalEB = ev.getHandle(tok_EB_);
0196 const EBRecHitCollection HitecalEB = *(ecalEB.product());
0197
0198
0199 auto tmpEcalRecHitCollection = std::make_unique<EcalRecHitCollection>();
0200 for (EcalRecHitCollection::const_iterator recHit = (*ecalEB).begin(); recHit != (*ecalEB).end(); ++recHit) {
0201 tmpEcalRecHitCollection->push_back(*recHit);
0202 }
0203 for (EcalRecHitCollection::const_iterator recHit = (*ecalEE).begin(); recHit != (*ecalEE).end(); ++recHit) {
0204 tmpEcalRecHitCollection->push_back(*recHit);
0205 }
0206 const EcalRecHitCollection Hitecal = *tmpEcalRecHitCollection;
0207
0208 const edm::Handle<reco::TrackCollection>& generalTracks = ev.getHandle(tok_tracks_);
0209
0210 geo = &c.getData(tok_geom_);
0211
0212 gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0213
0214 parameters_.useEcal = true;
0215 parameters_.useHcal = true;
0216 parameters_.useCalo = false;
0217 parameters_.useMuon = false;
0218 parameters_.dREcal = 0.5;
0219 parameters_.dRHcal = 0.6;
0220
0221
0222
0223
0224 const MagneticField* bField = &c.getData(tok_bFieldH_);
0225 stepPropF = new SteppingHelixPropagator(bField, alongMomentum);
0226 stepPropF->setMaterialMode(false);
0227 stepPropF->applyRadX0Correction(true);
0228
0229
0230
0231
0232
0233
0234
0235
0236 double mom_MC = 50.;
0237
0238
0239
0240
0241 const edm::Handle<edm::HepMCProduct>& evtMC = ev.getHandle(tok_gen_);
0242 if (!evtMC.isValid()) {
0243 edm::LogVerbatim("HcalCalib") << "no HepMCProduct found";
0244 } else {
0245
0246 #ifdef EDM_ML_DEBUG
0247 edm::LogVerbatim("HcalCalib") << "*** source HepMCProduct found";
0248 #endif
0249 }
0250
0251
0252 double maxPt = -99999.;
0253 int npart = 0;
0254
0255 GlobalPoint initpos(0, 0, 0);
0256
0257 HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evtMC->GetEvent()));
0258 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0259 phiParticle = (*p)->momentum().phi();
0260 etaParticle = (*p)->momentum().eta();
0261 double pt = (*p)->momentum().perp();
0262 mom_MC = (*p)->momentum().rho();
0263 if (pt > maxPt) {
0264 npart++;
0265 maxPt = pt;
0266 }
0267 GlobalVector mom((*p)->momentum().x(), (*p)->momentum().y(), (*p)->momentum().z());
0268 int charge = -1;
0269
0270 if (abs((*p)->pdg_id()) == 211)
0271 charge = (*p)->pdg_id() / abs((*p)->pdg_id());
0272 else
0273 continue;
0274
0275
0276
0277 const FreeTrajectoryState* freetrajectorystate_ = new FreeTrajectoryState(initpos, mom, charge, &(*theMagField));
0278
0279 TrackDetMatchInfo info = trackAssociator_.associate(ev, c, *freetrajectorystate_, parameters_);
0280
0281 GlobalPoint barrelMC(0, 0, 0), endcapMC(0, 0, 0), forwardMC1(0, 0, 0), forwardMC2(0, 0, 0);
0282
0283 GlobalPoint gPointHcal(0., 0., 0.);
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296 if (fabs(etaParticle) < 1.392) {
0297 Cylinder* cylinder = new Cylinder(181.1, Surface::PositionType(0, 0, 0), Surface::RotationType());
0298
0299 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*cylinder));
0300 if (steppingHelixstateinfo_.isValid()) {
0301 barrelMC = steppingHelixstateinfo_.freeState()->position();
0302 }
0303 }
0304
0305 doHF = kFALSE;
0306 if (fabs(etaParticle) >= 1.392 && fabs(etaParticle) < 5.191) {
0307 Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0308
0309 Surface::PositionType pos(0., 0., 400.5 * fabs(etaParticle) / etaParticle);
0310 PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
0311
0312 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane));
0313 if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) < 3.0)
0314 endcapMC = steppingHelixstateinfo_.freeState()->position();
0315 if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) > 3.0)
0316 doHF = kTRUE;
0317 }
0318 if (doHF) {
0319 if (abs(etaParticle) > 5.191)
0320 continue;
0321
0322 if (abs(etaParticle) > 2.99) {
0323 Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0324
0325 Surface::PositionType pos1(0., 0., 1125 * fabs(etaParticle) / etaParticle);
0326
0327 Surface::PositionType pos2(0., 0., 1137 * fabs(etaParticle) / etaParticle);
0328 PlaneBuilder::ReturnType aPlane1 = PlaneBuilder().plane(pos1, rot);
0329 PlaneBuilder::ReturnType aPlane2 = PlaneBuilder().plane(pos2, rot);
0330
0331 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane1));
0332 if (steppingHelixstateinfo_.isValid())
0333 forwardMC1 = steppingHelixstateinfo_.freeState()->position();
0334
0335 steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
0336 if (steppingHelixstateinfo_.isValid())
0337 forwardMC2 = steppingHelixstateinfo_.freeState()->position();
0338 }
0339 }
0340
0341
0342
0343 Int_t iphitrue = -10;
0344 Int_t ietatrue = 100;
0345 HcalDetId tempId;
0346
0347 if (abs(etaParticle) < 1.392) {
0348 gPointHcal = barrelMC;
0349 tempId = gHcal->getClosestCell(gPointHcal);
0350 }
0351 if (abs(etaParticle) >= 1.392 && abs(etaParticle) < 3.0) {
0352 gPointHcal = endcapMC;
0353 tempId = gHcal->getClosestCell(gPointHcal);
0354 }
0355 if (abs(etaParticle) >= 3.0 && abs(etaParticle) < 5.191) {
0356 gPointHcal = forwardMC1;
0357 tempId = gHcal->getClosestCell(gPointHcal);
0358
0359 }
0360
0361 tempId = gHcal->getClosestCell(gPointHcal);
0362
0363 ietatrue = tempId.ieta();
0364 iphitrue = tempId.iphi();
0365
0366 etaGPoint = gPointHcal.eta();
0367 phiGPoint = gPointHcal.phi();
0368
0369
0370
0371
0372
0373
0374 if (gPointHcal.x() == 0 && gPointHcal.y() == 0 && gPointHcal.z() == 0) {
0375 #ifdef EDM_ML_DEBUG
0376 edm::LogVerbatim("HcalCalib") << "gPointHcal is Zero!";
0377 #endif
0378 continue;
0379 }
0380
0381 float etahcal = gPointHcal.eta();
0382 if (abs(etahcal) > 5.192)
0383 continue;
0384 #ifdef EDM_ML_DEBUG
0385 if (std::abs(etahcal) > 3.0 && std::abs(etahcal) < 5.191) {
0386 edm::LogVerbatim("HcalCalib") << gPointHcal.x() << " " << gPointHcal.y() << " " << gPointHcal.z() << " "
0387 << gPointHcal.eta() << " " << gPointHcal.phi() << " " << ietatrue << " "
0388 << iphitrue;
0389 if (ietatrue == 100 || iphitrue == -10) {
0390 edm::LogVerbatim("HcalCalib") << "ietatrue: " << ietatrue << " iphitrue: " << iphitrue
0391 << " etahcal: " << etahcal << " phihcal: " << gPointHcal.phi();
0392 }
0393 }
0394 #endif
0395
0396
0397
0398
0399 xTrkEcal = info.trkGlobPosAtEcal.x();
0400 yTrkEcal = info.trkGlobPosAtEcal.y();
0401 zTrkEcal = info.trkGlobPosAtEcal.z();
0402
0403 GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0404
0405 eECAL = ecalEnergyInCone(gPointEcal, ecalCone_, Hitecal, geo);
0406 eECAL09cm = ecalEnergyInCone(gPointEcal, 9, Hitecal, geo);
0407 eECAL40cm = ecalEnergyInCone(gPointEcal, 40, Hitecal, geo);
0408
0409 eEcalCone = eECAL;
0410
0411
0412
0413
0414
0415
0416 MaxHit_struct MaxHit;
0417
0418 MaxHit.hitenergy = -100.;
0419 Float_t recal = 1.0;
0420
0421
0422 eHcalCone = 0.;
0423 eHcalConeNoise = 0.;
0424 UsedCells = 0;
0425 UsedCellsNoise = 0;
0426 e3x3 = 0.;
0427 e5x5 = 0.;
0428
0429 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0430
0431 {
0432 recal = RecalibFactor(hhit->detid());
0433 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0434
0435 int iphihit = (hhit->id()).iphi();
0436 int ietahit = (hhit->id()).ieta();
0437 int depthhit = (hhit->id()).depth();
0438 float enehit = hhit->energy() * recal;
0439
0440 if (depthhit != 1)
0441 continue;
0442
0443 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0444
0445 if (distAtHcal < clusterConeSize_) {
0446 for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++)
0447
0448 {
0449 int iphihit2 = (hhit2->id()).iphi();
0450 int ietahit2 = (hhit2->id()).ieta();
0451 int depthhit2 = (hhit2->id()).depth();
0452 float enehit2 = hhit2->energy() * recal;
0453
0454 if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0455 enehit = enehit + enehit2;
0456 }
0457
0458 if (enehit > MaxHit.hitenergy) {
0459 MaxHit.hitenergy = enehit;
0460 MaxHit.ietahitm = (hhit->id()).ieta();
0461 MaxHit.iphihitm = (hhit->id()).iphi();
0462 MaxHit.dr = distAtHcal;
0463
0464 MaxHit.depthhit = 1;
0465 }
0466 }
0467 }
0468
0469 for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0470 recal = RecalibFactor(hhit->detid());
0471
0472 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0473
0474 int iphihit = (hhit->id()).iphi();
0475 int ietahit = (hhit->id()).ieta();
0476 int depthhit = (hhit->id()).depth();
0477 float enehit = hhit->energy() * recal;
0478
0479 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0480
0481 if (distAtHcal < associationConeSize_) {
0482 for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0483 int iphihit2 = (hhit2->id()).iphi();
0484 int ietahit2 = (hhit2->id()).ieta();
0485 int depthhit2 = (hhit2->id()).depth();
0486 float enehit2 = hhit2->energy() * recal;
0487
0488 if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0489 enehit = enehit + enehit2;
0490 }
0491
0492 if (enehit > MaxHit.hitenergy) {
0493 MaxHit.hitenergy = enehit;
0494 MaxHit.ietahitm = (hhit->id()).ieta();
0495 MaxHit.iphihitm = (hhit->id()).iphi();
0496 MaxHit.dr = distAtHcal;
0497 MaxHit.depthhit = 1;
0498 }
0499 }
0500 }
0501
0502
0503
0504
0505
0506 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0507
0508 {
0509 recal = RecalibFactor(hhit->detid());
0510 #ifdef EDM_ML_DEBUG
0511 edm::LogVerbatim("HcalCalib") << "recal: " << recal;
0512 #endif
0513 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0514
0515 int iphihit = (hhit->id()).iphi();
0516 int ietahit = (hhit->id()).ieta();
0517 int depthhit = (hhit->id()).depth();
0518 float enehit = hhit->energy() * recal;
0519
0520
0521
0522
0523 int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0524 int ietahitNoise = -ietahit;
0525 int depthhitNoise = depthhit;
0526
0527 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0528 if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.) {
0529 eHcalCone += enehit;
0530 UsedCells++;
0531
0532 int DIETA = 100;
0533 if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0534 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0535 }
0536 if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0537 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0538 DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0539 }
0540 int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0541 DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0542
0543 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0544 !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0545 abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0546 e5x5 += hhit->energy();
0547 }
0548 if (abs(DIETA) <= 1 &&
0549 (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0550 !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0551 e3x3 += hhit->energy();
0552 }
0553 #ifdef EDM_ML_DEBUG
0554 edm::LogVerbatim("HcalCalib") << "track: ieta " << ietahit << " iphi: " << iphihit << " depth: " << depthhit
0555 << " energydepos: " << enehit;
0556 #endif
0557 for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0558 recal = RecalibFactor(hhit2->detid());
0559 int iphihit2 = (hhit2->id()).iphi();
0560 int ietahit2 = (hhit2->id()).ieta();
0561 int depthhit2 = (hhit2->id()).depth();
0562 float enehit2 = hhit2->energy() * recal;
0563
0564 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.) {
0565 eHcalConeNoise += hhit2->energy() * recal;
0566 UsedCellsNoise++;
0567 #ifdef EDM_ML_DEBUG
0568 edm::LogVerbatim("HcalCalib") << "Noise: ieta " << ietahit2 << " iphi: " << iphihit2
0569 << " depth: " << depthhit2 << " energydepos: " << enehit2;
0570 #endif
0571 }
0572 }
0573 }
0574 }
0575
0576 for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0577 recal = RecalibFactor(hhit->detid());
0578
0579 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0580
0581
0582
0583 int iphihit = (hhit->id()).iphi();
0584 int ietahit = (hhit->id()).ieta();
0585 int depthhit = (hhit->id()).depth();
0586 float enehit = hhit->energy() * recal;
0587
0588
0589 int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0590 int ietahitNoise = -ietahit;
0591 int depthhitNoise = depthhit;
0592
0593 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0594
0595 if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
0596
0597 {
0598 eHcalCone += enehit;
0599 UsedCells++;
0600
0601 int DIETA = 100;
0602 if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0603 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0604 }
0605 if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0606 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0607 DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0608 }
0609 int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0610 DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0611
0612 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0613 !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0614 abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0615 e5x5 += hhit->energy();
0616 }
0617 if (abs(DIETA) <= 1 &&
0618 (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0619 !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0620 e3x3 += hhit->energy();
0621 }
0622
0623 for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0624 recal = RecalibFactor(hhit2->detid());
0625
0626 int iphihit2 = (hhit2->id()).iphi();
0627 int ietahit2 = (hhit2->id()).ieta();
0628 int depthhit2 = (hhit2->id()).depth();
0629 float enehit2 = hhit2->energy() * recal;
0630
0631 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.01) {
0632 eHcalConeNoise += hhit2->energy() * recal;
0633 UsedCellsNoise++;
0634 }
0635 }
0636 }
0637 }
0638
0639
0640
0641
0642
0643 nTracks = 0;
0644
0645 delRmc[0] = 5;
0646
0647 float delR_track_particle = 100;
0648
0649 for (reco::TrackCollection::const_iterator track1 = generalTracks->begin(); track1 != generalTracks->end();
0650 track1++) {
0651 delR_track_particle = deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
0652
0653 trackEta[nTracks] = track1->eta();
0654 trackPhi[nTracks] = track1->phi();
0655 trackP[nTracks] = sqrt(track1->px() * track1->px() + track1->py() * track1->py() + track1->pz() * track1->pz());
0656
0657 delRmc[nTracks] = delR_track_particle;
0658 numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
0659 numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
0660 numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement();
0661 trkQual[nTracks] = track1->quality(reco::TrackBase::highPurity);
0662
0663 nTracks++;
0664 }
0665
0666
0667
0668 int dieta_M_P = 100;
0669 int diphi_M_P = 100;
0670 if (MaxHit.ietahitm * ietatrue > 0) {
0671 dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0672 }
0673 if (MaxHit.ietahitm * ietatrue < 0) {
0674 dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0675 }
0676 diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0677
0678 diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0679 iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0680
0681 #ifdef EDM_ML_DEBUG
0682 if (iDr > 15) {
0683 edm::LogVerbatim("HcalCalib") << "diphi: " << diphi_M_P << " dieta: " << dieta_M_P << " iDr: " << iDr
0684 << " ietatrue:" << ietatrue << " iphitrue:" << iphitrue;
0685 edm::LogVerbatim("HcalCalib") << "M ieta: " << MaxHit.ietahitm << " M iphi: " << MaxHit.iphihitm;
0686 }
0687 #endif
0688
0689 Bool_t passCuts = kFALSE;
0690 passCuts = kTRUE;
0691
0692
0693
0694 if (passCuts) {
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704 iEta = ietatrue;
0705 iPhi = iphitrue;
0706
0707
0708
0709 delR = MaxHit.dr;
0710 eCentHit = MaxHit.hitenergy;
0711
0712 eParticle = mom_MC;
0713
0714
0715
0716
0717 pfTree->Fill();
0718 }
0719
0720 }
0721 }
0722
0723 void HcalCorrPFCalculation::beginJob() {
0724 pfTree = fs->make<TTree>("pfTree", "Tree for pf info");
0725
0726 pfTree->Branch("nTracks", &nTracks, "nTracks/I");
0727 pfTree->Branch("trackEta", trackEta, "trackEta[nTracks]/F");
0728 pfTree->Branch("trackPhi", trackPhi, "trackPhi[nTracks]/F");
0729 pfTree->Branch("trackP", trackP, "trackP[nTracks]/F");
0730
0731 pfTree->Branch("delRmc", delRmc, "delRmc[nTracks]/F");
0732 pfTree->Branch("numValidTrkHits", numValidTrkHits, "numValidTrkHits[nTracks]/I");
0733 pfTree->Branch("numValidTrkStrips", numValidTrkStrips, "numValidTrkStrips[nTracks]/I");
0734 pfTree->Branch("numLayers", numLayers, "numLayers[nTracks]/I");
0735 pfTree->Branch("trkQual", trkQual, "trkQual[nTracks]/O");
0736
0737 pfTree->Branch("eEcalCone", &eEcalCone, "eEcalCone/F");
0738 pfTree->Branch("eHcalCone", &eHcalCone, "eHcalCone/F");
0739 pfTree->Branch("eHcalConeNoise", &eHcalConeNoise, "eHcalConeNoise/F");
0740
0741 pfTree->Branch("UsedCellsNoise", &UsedCellsNoise, "UsedCellsNoise/I");
0742 pfTree->Branch("UsedCells", &UsedCells, "UsedCells/I");
0743
0744 pfTree->Branch("eCentHit", &eCentHit, "eCentHit/F");
0745
0746 pfTree->Branch("eParticle", &eParticle, "eParticle/F");
0747 pfTree->Branch("etaParticle", &etaParticle, "etaParticle/F");
0748 pfTree->Branch("phiParticle", &phiParticle, "phiParticle/F");
0749
0750 pfTree->Branch("etaGPoint", &etaGPoint, "etaGPoint/F");
0751 pfTree->Branch("phiGPoint", &phiGPoint, "phiGPoint/F");
0752
0753 pfTree->Branch("xAtHcal", &xAtHcal, "xAtHcal/F");
0754 pfTree->Branch("yAtHcal", &yAtHcal, "yAtHcal/F");
0755 pfTree->Branch("zAtHcal", &zAtHcal, "zAtHcal/F");
0756
0757 pfTree->Branch("eECAL09cm", &eECAL09cm, "eECAL09cm/F");
0758 pfTree->Branch("eECAL40cm", &eECAL40cm, "eECAL40cm/F");
0759 pfTree->Branch("eECAL", &eECAL, "eECAL/F");
0760
0761 pfTree->Branch("e3x3 ", &e3x3, "e3x3/F");
0762 pfTree->Branch("e5x5", &e5x5, "e5x5/F");
0763
0764 pfTree->Branch("iDr", &iDr, "iDr/F");
0765 pfTree->Branch("delR", &delR, "delR/F");
0766
0767 pfTree->Branch("iEta", &iEta, "iEta/I");
0768 pfTree->Branch("iPhi", &iPhi, "iPhi/I");
0769
0770
0771
0772
0773
0774 }
0775
0776 #include "FWCore/Framework/interface/MakerMacros.h"
0777
0778
0779 DEFINE_FWK_MODULE(HcalCorrPFCalculation);