File indexing completed on 2024-04-06 11:59: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
0254 GlobalPoint initpos(0, 0, 0);
0255
0256 HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evtMC->GetEvent()));
0257 for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0258 phiParticle = (*p)->momentum().phi();
0259 etaParticle = (*p)->momentum().eta();
0260 double pt = (*p)->momentum().perp();
0261 mom_MC = (*p)->momentum().rho();
0262 if (pt > maxPt) {
0263 maxPt = pt;
0264 }
0265 GlobalVector mom((*p)->momentum().x(), (*p)->momentum().y(), (*p)->momentum().z());
0266 int charge = -1;
0267
0268 if (abs((*p)->pdg_id()) == 211)
0269 charge = (*p)->pdg_id() / abs((*p)->pdg_id());
0270 else
0271 continue;
0272
0273
0274
0275 const FreeTrajectoryState* freetrajectorystate_ = new FreeTrajectoryState(initpos, mom, charge, &(*theMagField));
0276
0277 TrackDetMatchInfo info = trackAssociator_.associate(ev, c, *freetrajectorystate_, parameters_);
0278
0279 GlobalPoint barrelMC(0, 0, 0), endcapMC(0, 0, 0), forwardMC1(0, 0, 0), forwardMC2(0, 0, 0);
0280
0281 GlobalPoint gPointHcal(0., 0., 0.);
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294 if (fabs(etaParticle) < 1.392) {
0295 Cylinder* cylinder = new Cylinder(181.1, Surface::PositionType(0, 0, 0), Surface::RotationType());
0296
0297 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*cylinder));
0298 if (steppingHelixstateinfo_.isValid()) {
0299 barrelMC = steppingHelixstateinfo_.freeState()->position();
0300 }
0301 }
0302
0303 doHF = kFALSE;
0304 if (fabs(etaParticle) >= 1.392 && fabs(etaParticle) < 5.191) {
0305 Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0306
0307 Surface::PositionType pos(0., 0., 400.5 * fabs(etaParticle) / etaParticle);
0308 PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
0309
0310 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane));
0311 if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) < 3.0)
0312 endcapMC = steppingHelixstateinfo_.freeState()->position();
0313 if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) > 3.0)
0314 doHF = kTRUE;
0315 }
0316 if (doHF) {
0317 if (abs(etaParticle) > 5.191)
0318 continue;
0319
0320 if (abs(etaParticle) > 2.99) {
0321 Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0322
0323 Surface::PositionType pos1(0., 0., 1125 * fabs(etaParticle) / etaParticle);
0324
0325 Surface::PositionType pos2(0., 0., 1137 * fabs(etaParticle) / etaParticle);
0326 PlaneBuilder::ReturnType aPlane1 = PlaneBuilder().plane(pos1, rot);
0327 PlaneBuilder::ReturnType aPlane2 = PlaneBuilder().plane(pos2, rot);
0328
0329 TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane1));
0330 if (steppingHelixstateinfo_.isValid())
0331 forwardMC1 = steppingHelixstateinfo_.freeState()->position();
0332
0333 steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
0334 if (steppingHelixstateinfo_.isValid())
0335 forwardMC2 = steppingHelixstateinfo_.freeState()->position();
0336 }
0337 }
0338
0339
0340
0341 Int_t iphitrue = -10;
0342 Int_t ietatrue = 100;
0343 HcalDetId tempId;
0344
0345 if (abs(etaParticle) < 1.392) {
0346 gPointHcal = barrelMC;
0347 tempId = gHcal->getClosestCell(gPointHcal);
0348 }
0349 if (abs(etaParticle) >= 1.392 && abs(etaParticle) < 3.0) {
0350 gPointHcal = endcapMC;
0351 tempId = gHcal->getClosestCell(gPointHcal);
0352 }
0353 if (abs(etaParticle) >= 3.0 && abs(etaParticle) < 5.191) {
0354 gPointHcal = forwardMC1;
0355 tempId = gHcal->getClosestCell(gPointHcal);
0356
0357 }
0358
0359 tempId = gHcal->getClosestCell(gPointHcal);
0360
0361 ietatrue = tempId.ieta();
0362 iphitrue = tempId.iphi();
0363
0364 etaGPoint = gPointHcal.eta();
0365 phiGPoint = gPointHcal.phi();
0366
0367
0368
0369
0370
0371
0372 if (gPointHcal.x() == 0 && gPointHcal.y() == 0 && gPointHcal.z() == 0) {
0373 #ifdef EDM_ML_DEBUG
0374 edm::LogVerbatim("HcalCalib") << "gPointHcal is Zero!";
0375 #endif
0376 continue;
0377 }
0378
0379 float etahcal = gPointHcal.eta();
0380 if (abs(etahcal) > 5.192)
0381 continue;
0382 #ifdef EDM_ML_DEBUG
0383 if (std::abs(etahcal) > 3.0 && std::abs(etahcal) < 5.191) {
0384 edm::LogVerbatim("HcalCalib") << gPointHcal.x() << " " << gPointHcal.y() << " " << gPointHcal.z() << " "
0385 << gPointHcal.eta() << " " << gPointHcal.phi() << " " << ietatrue << " "
0386 << iphitrue;
0387 if (ietatrue == 100 || iphitrue == -10) {
0388 edm::LogVerbatim("HcalCalib") << "ietatrue: " << ietatrue << " iphitrue: " << iphitrue
0389 << " etahcal: " << etahcal << " phihcal: " << gPointHcal.phi();
0390 }
0391 }
0392 #endif
0393
0394
0395
0396
0397 xTrkEcal = info.trkGlobPosAtEcal.x();
0398 yTrkEcal = info.trkGlobPosAtEcal.y();
0399 zTrkEcal = info.trkGlobPosAtEcal.z();
0400
0401 GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0402
0403 eECAL = ecalEnergyInCone(gPointEcal, ecalCone_, Hitecal, geo);
0404 eECAL09cm = ecalEnergyInCone(gPointEcal, 9, Hitecal, geo);
0405 eECAL40cm = ecalEnergyInCone(gPointEcal, 40, Hitecal, geo);
0406
0407 eEcalCone = eECAL;
0408
0409
0410
0411
0412
0413
0414 MaxHit_struct MaxHit;
0415
0416 MaxHit.hitenergy = -100.;
0417 Float_t recal = 1.0;
0418
0419
0420 eHcalCone = 0.;
0421 eHcalConeNoise = 0.;
0422 UsedCells = 0;
0423 UsedCellsNoise = 0;
0424 e3x3 = 0.;
0425 e5x5 = 0.;
0426
0427 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0428
0429 {
0430 recal = RecalibFactor(hhit->detid());
0431 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0432
0433 int iphihit = (hhit->id()).iphi();
0434 int ietahit = (hhit->id()).ieta();
0435 int depthhit = (hhit->id()).depth();
0436 float enehit = hhit->energy() * recal;
0437
0438 if (depthhit != 1)
0439 continue;
0440
0441 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0442
0443 if (distAtHcal < clusterConeSize_) {
0444 for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++)
0445
0446 {
0447 int iphihit2 = (hhit2->id()).iphi();
0448 int ietahit2 = (hhit2->id()).ieta();
0449 int depthhit2 = (hhit2->id()).depth();
0450 float enehit2 = hhit2->energy() * recal;
0451
0452 if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0453 enehit = enehit + enehit2;
0454 }
0455
0456 if (enehit > MaxHit.hitenergy) {
0457 MaxHit.hitenergy = enehit;
0458 MaxHit.ietahitm = (hhit->id()).ieta();
0459 MaxHit.iphihitm = (hhit->id()).iphi();
0460 MaxHit.dr = distAtHcal;
0461
0462 MaxHit.depthhit = 1;
0463 }
0464 }
0465 }
0466
0467 for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0468 recal = RecalibFactor(hhit->detid());
0469
0470 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0471
0472 int iphihit = (hhit->id()).iphi();
0473 int ietahit = (hhit->id()).ieta();
0474 int depthhit = (hhit->id()).depth();
0475 float enehit = hhit->energy() * recal;
0476
0477 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0478
0479 if (distAtHcal < associationConeSize_) {
0480 for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0481 int iphihit2 = (hhit2->id()).iphi();
0482 int ietahit2 = (hhit2->id()).ieta();
0483 int depthhit2 = (hhit2->id()).depth();
0484 float enehit2 = hhit2->energy() * recal;
0485
0486 if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0487 enehit = enehit + enehit2;
0488 }
0489
0490 if (enehit > MaxHit.hitenergy) {
0491 MaxHit.hitenergy = enehit;
0492 MaxHit.ietahitm = (hhit->id()).ieta();
0493 MaxHit.iphihitm = (hhit->id()).iphi();
0494 MaxHit.dr = distAtHcal;
0495 MaxHit.depthhit = 1;
0496 }
0497 }
0498 }
0499
0500
0501
0502
0503
0504 for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0505
0506 {
0507 recal = RecalibFactor(hhit->detid());
0508 #ifdef EDM_ML_DEBUG
0509 edm::LogVerbatim("HcalCalib") << "recal: " << recal;
0510 #endif
0511 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0512
0513 int iphihit = (hhit->id()).iphi();
0514 int ietahit = (hhit->id()).ieta();
0515 int depthhit = (hhit->id()).depth();
0516 float enehit = hhit->energy() * recal;
0517
0518
0519
0520
0521 int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0522 int ietahitNoise = -ietahit;
0523 int depthhitNoise = depthhit;
0524
0525 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0526 if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.) {
0527 eHcalCone += enehit;
0528 UsedCells++;
0529
0530 int DIETA = 100;
0531 if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0532 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0533 }
0534 if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0535 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0536 DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0537 }
0538 int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0539 DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0540
0541 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0542 !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0543 abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0544 e5x5 += hhit->energy();
0545 }
0546 if (abs(DIETA) <= 1 &&
0547 (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0548 !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0549 e3x3 += hhit->energy();
0550 }
0551 #ifdef EDM_ML_DEBUG
0552 edm::LogVerbatim("HcalCalib") << "track: ieta " << ietahit << " iphi: " << iphihit << " depth: " << depthhit
0553 << " energydepos: " << enehit;
0554 #endif
0555 for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0556 recal = RecalibFactor(hhit2->detid());
0557 int iphihit2 = (hhit2->id()).iphi();
0558 int ietahit2 = (hhit2->id()).ieta();
0559 int depthhit2 = (hhit2->id()).depth();
0560 float enehit2 = hhit2->energy() * recal;
0561
0562 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.) {
0563 eHcalConeNoise += hhit2->energy() * recal;
0564 UsedCellsNoise++;
0565 #ifdef EDM_ML_DEBUG
0566 edm::LogVerbatim("HcalCalib") << "Noise: ieta " << ietahit2 << " iphi: " << iphihit2
0567 << " depth: " << depthhit2 << " energydepos: " << enehit2;
0568 #endif
0569 }
0570 }
0571 }
0572 }
0573
0574 for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0575 recal = RecalibFactor(hhit->detid());
0576
0577 GlobalPoint pos = gHcal->getPosition(hhit->detid());
0578
0579
0580
0581 int iphihit = (hhit->id()).iphi();
0582 int ietahit = (hhit->id()).ieta();
0583 int depthhit = (hhit->id()).depth();
0584 float enehit = hhit->energy() * recal;
0585
0586
0587 int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0588 int ietahitNoise = -ietahit;
0589 int depthhitNoise = depthhit;
0590
0591 double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0592
0593 if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
0594
0595 {
0596 eHcalCone += enehit;
0597 UsedCells++;
0598
0599 int DIETA = 100;
0600 if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0601 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0602 }
0603 if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0604 DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0605 DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0606 }
0607 int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0608 DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0609
0610 if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0611 !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0612 abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0613 e5x5 += hhit->energy();
0614 }
0615 if (abs(DIETA) <= 1 &&
0616 (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0617 !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0618 e3x3 += hhit->energy();
0619 }
0620
0621 for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0622 recal = RecalibFactor(hhit2->detid());
0623
0624 int iphihit2 = (hhit2->id()).iphi();
0625 int ietahit2 = (hhit2->id()).ieta();
0626 int depthhit2 = (hhit2->id()).depth();
0627 float enehit2 = hhit2->energy() * recal;
0628
0629 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.01) {
0630 eHcalConeNoise += hhit2->energy() * recal;
0631 UsedCellsNoise++;
0632 }
0633 }
0634 }
0635 }
0636
0637
0638
0639
0640
0641 nTracks = 0;
0642
0643 delRmc[0] = 5;
0644
0645 float delR_track_particle = 100;
0646
0647 for (reco::TrackCollection::const_iterator track1 = generalTracks->begin(); track1 != generalTracks->end();
0648 track1++) {
0649 delR_track_particle = deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
0650
0651 trackEta[nTracks] = track1->eta();
0652 trackPhi[nTracks] = track1->phi();
0653 trackP[nTracks] = sqrt(track1->px() * track1->px() + track1->py() * track1->py() + track1->pz() * track1->pz());
0654
0655 delRmc[nTracks] = delR_track_particle;
0656 numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
0657 numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
0658 numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement();
0659 trkQual[nTracks] = track1->quality(reco::TrackBase::highPurity);
0660
0661 nTracks++;
0662 }
0663
0664
0665
0666 int dieta_M_P = 100;
0667 int diphi_M_P = 100;
0668 if (MaxHit.ietahitm * ietatrue > 0) {
0669 dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0670 }
0671 if (MaxHit.ietahitm * ietatrue < 0) {
0672 dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0673 }
0674 diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0675
0676 diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0677 iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0678
0679 #ifdef EDM_ML_DEBUG
0680 if (iDr > 15) {
0681 edm::LogVerbatim("HcalCalib") << "diphi: " << diphi_M_P << " dieta: " << dieta_M_P << " iDr: " << iDr
0682 << " ietatrue:" << ietatrue << " iphitrue:" << iphitrue;
0683 edm::LogVerbatim("HcalCalib") << "M ieta: " << MaxHit.ietahitm << " M iphi: " << MaxHit.iphihitm;
0684 }
0685 #endif
0686
0687 Bool_t passCuts = kFALSE;
0688 passCuts = kTRUE;
0689
0690
0691
0692 if (passCuts) {
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702 iEta = ietatrue;
0703 iPhi = iphitrue;
0704
0705
0706
0707 delR = MaxHit.dr;
0708 eCentHit = MaxHit.hitenergy;
0709
0710 eParticle = mom_MC;
0711
0712
0713
0714
0715 pfTree->Fill();
0716 }
0717
0718 }
0719 }
0720
0721 void HcalCorrPFCalculation::beginJob() {
0722 pfTree = fs->make<TTree>("pfTree", "Tree for pf info");
0723
0724 pfTree->Branch("nTracks", &nTracks, "nTracks/I");
0725 pfTree->Branch("trackEta", trackEta, "trackEta[nTracks]/F");
0726 pfTree->Branch("trackPhi", trackPhi, "trackPhi[nTracks]/F");
0727 pfTree->Branch("trackP", trackP, "trackP[nTracks]/F");
0728
0729 pfTree->Branch("delRmc", delRmc, "delRmc[nTracks]/F");
0730 pfTree->Branch("numValidTrkHits", numValidTrkHits, "numValidTrkHits[nTracks]/I");
0731 pfTree->Branch("numValidTrkStrips", numValidTrkStrips, "numValidTrkStrips[nTracks]/I");
0732 pfTree->Branch("numLayers", numLayers, "numLayers[nTracks]/I");
0733 pfTree->Branch("trkQual", trkQual, "trkQual[nTracks]/O");
0734
0735 pfTree->Branch("eEcalCone", &eEcalCone, "eEcalCone/F");
0736 pfTree->Branch("eHcalCone", &eHcalCone, "eHcalCone/F");
0737 pfTree->Branch("eHcalConeNoise", &eHcalConeNoise, "eHcalConeNoise/F");
0738
0739 pfTree->Branch("UsedCellsNoise", &UsedCellsNoise, "UsedCellsNoise/I");
0740 pfTree->Branch("UsedCells", &UsedCells, "UsedCells/I");
0741
0742 pfTree->Branch("eCentHit", &eCentHit, "eCentHit/F");
0743
0744 pfTree->Branch("eParticle", &eParticle, "eParticle/F");
0745 pfTree->Branch("etaParticle", &etaParticle, "etaParticle/F");
0746 pfTree->Branch("phiParticle", &phiParticle, "phiParticle/F");
0747
0748 pfTree->Branch("etaGPoint", &etaGPoint, "etaGPoint/F");
0749 pfTree->Branch("phiGPoint", &phiGPoint, "phiGPoint/F");
0750
0751 pfTree->Branch("xAtHcal", &xAtHcal, "xAtHcal/F");
0752 pfTree->Branch("yAtHcal", &yAtHcal, "yAtHcal/F");
0753 pfTree->Branch("zAtHcal", &zAtHcal, "zAtHcal/F");
0754
0755 pfTree->Branch("eECAL09cm", &eECAL09cm, "eECAL09cm/F");
0756 pfTree->Branch("eECAL40cm", &eECAL40cm, "eECAL40cm/F");
0757 pfTree->Branch("eECAL", &eECAL, "eECAL/F");
0758
0759 pfTree->Branch("e3x3 ", &e3x3, "e3x3/F");
0760 pfTree->Branch("e5x5", &e5x5, "e5x5/F");
0761
0762 pfTree->Branch("iDr", &iDr, "iDr/F");
0763 pfTree->Branch("delR", &delR, "delR/F");
0764
0765 pfTree->Branch("iEta", &iEta, "iEta/I");
0766 pfTree->Branch("iPhi", &iPhi, "iPhi/I");
0767
0768
0769
0770
0771
0772 }
0773
0774 #include "FWCore/Framework/interface/MakerMacros.h"
0775
0776
0777 DEFINE_FWK_MODULE(HcalCorrPFCalculation);