File indexing completed on 2023-01-13 23:40:10
0001 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0002 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0003 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0004 #include "Validation/HcalRecHits/interface/HcalRecHitsValidation.h"
0005
0006
0007
0008 HcalRecHitsValidation::HcalRecHitsValidation(edm::ParameterSet const &conf)
0009 : topFolderName_(conf.getParameter<std::string>("TopFolderName")),
0010 outputFile_(conf.getUntrackedParameter<std::string>("outputFile", "myfile.root")),
0011 hcalselector_(conf.getUntrackedParameter<std::string>("hcalselector", "all")),
0012 ecalselector_(conf.getUntrackedParameter<std::string>("ecalselector", "yes")),
0013 sign_(conf.getUntrackedParameter<std::string>("sign", "*")),
0014 mc_(conf.getUntrackedParameter<std::string>("mc", "yes")),
0015 testNumber_(conf.getParameter<bool>("TestNumber")),
0016 EBRecHitCollectionLabel_(conf.getParameter<edm::InputTag>("EBRecHitCollectionLabel")),
0017 EERecHitCollectionLabel_(conf.getParameter<edm::InputTag>("EERecHitCollectionLabel")),
0018 tok_evt_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0019 tok_EB_(consumes<EBRecHitCollection>(EBRecHitCollectionLabel_)),
0020 tok_EE_(consumes<EERecHitCollection>(EERecHitCollectionLabel_)),
0021 tok_hh_(consumes<edm::PCaloHitContainer>(conf.getUntrackedParameter<edm::InputTag>("SimHitCollectionLabel"))),
0022 tok_hbhe_(consumes<HBHERecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HBHERecHitCollectionLabel"))),
0023 tok_hf_(consumes<HFRecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HFRecHitCollectionLabel"))),
0024 tok_ho_(consumes<HORecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HORecHitCollectionLabel"))),
0025 tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()),
0026 tok_Geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0027
0028 if (!outputFile_.empty()) {
0029 edm::LogVerbatim("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
0030 } else {
0031 edm::LogVerbatim("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
0032 }
0033
0034 nevtot = 0;
0035
0036
0037
0038
0039
0040 subdet_ = 5;
0041 if (hcalselector_ == "noise")
0042 subdet_ = 0;
0043 if (hcalselector_ == "HB")
0044 subdet_ = 1;
0045 if (hcalselector_ == "HE")
0046 subdet_ = 2;
0047 if (hcalselector_ == "HO")
0048 subdet_ = 3;
0049 if (hcalselector_ == "HF")
0050 subdet_ = 4;
0051 if (hcalselector_ == "all")
0052 subdet_ = 5;
0053 if (hcalselector_ == "ZS")
0054 subdet_ = 6;
0055
0056 iz = 1;
0057 if (sign_ == "-")
0058 iz = -1;
0059 if (sign_ == "*")
0060 iz = 0;
0061
0062 imc = 1;
0063 if (mc_ == "no")
0064 imc = 0;
0065 }
0066
0067 void HcalRecHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &run, edm::EventSetup const &es) {
0068 Char_t histo[200];
0069
0070 ib.setCurrentFolder(topFolderName_);
0071
0072
0073
0074
0075 if (subdet_ != 0 && imc != 0) {
0076 sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths");
0077 meEnConeEtaProfile = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0078
0079 sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_E");
0080 meEnConeEtaProfile_E = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0081
0082 sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_EH");
0083 meEnConeEtaProfile_EH = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0084 }
0085
0086
0087 if (subdet_ == 1 || subdet_ == 5) {
0088 sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HB");
0089 meRecHitsM2Chi2HB = ib.book1D(histo, histo, 120, -2., 10.);
0090
0091 sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HB");
0092 meLog10Chi2profileHB = ib.bookProfile(histo, histo, 300, -5., 295., -2., 9.9, " ");
0093
0094 sprintf(histo, "HcalRecHitTask_energy_of_rechits_HB");
0095 meRecHitsEnergyHB = ib.book1D(histo, histo, 2010, -10., 2000.);
0096
0097 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HB");
0098 meTEprofileHB = ib.bookProfile(histo, histo, 150, -5., 295., -48., 92., " ");
0099
0100 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HB");
0101 meTEprofileHB_Low = ib.bookProfile(histo, histo, 150, -5., 295., -48., 92., " ");
0102
0103 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HB");
0104 meTEprofileHB_High = ib.bookProfile(histo, histo, 150, -5., 295., 48., 92., " ");
0105
0106 if (imc != 0) {
0107 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HB");
0108 meRecHitSimHitHB = ib.book2D(histo, histo, 120, 0., 1.2, 300, 0., 150.);
0109 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HB");
0110 meRecHitSimHitProfileHB = ib.bookProfile(histo, histo, 120, 0., 1.2, 0., 500., " ");
0111 }
0112 }
0113
0114
0115 if (subdet_ == 2 || subdet_ == 5) {
0116 sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HE");
0117 meRecHitsM2Chi2HE = ib.book1D(histo, histo, 120, -2., 10.);
0118
0119 sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HE");
0120 meLog10Chi2profileHE = ib.bookProfile(histo, histo, 1000, -5., 995., -2., 9.9, " ");
0121
0122 sprintf(histo, "HcalRecHitTask_energy_of_rechits_HE");
0123 meRecHitsEnergyHE = ib.book1D(histo, histo, 2010, -10., 2000.);
0124
0125 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HE");
0126 meTEprofileHE_Low = ib.bookProfile(histo, histo, 80, -5., 75., -48., 92., " ");
0127
0128 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HE");
0129 meTEprofileHE = ib.bookProfile(histo, histo, 200, -5., 2995., -48., 92., " ");
0130
0131 if (imc != 0) {
0132 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HE");
0133 meRecHitSimHitHE = ib.book2D(histo, histo, 120, 0., 0.6, 300, 0., 150.);
0134 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HE");
0135 meRecHitSimHitProfileHE = ib.bookProfile(histo, histo, 120, 0., 0.6, 0., 500., " ");
0136 }
0137 }
0138
0139
0140 if (subdet_ == 3 || subdet_ == 5) {
0141 sprintf(histo, "HcalRecHitTask_energy_of_rechits_HO");
0142 meRecHitsEnergyHO = ib.book1D(histo, histo, 2010, -10., 2000.);
0143
0144 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HO");
0145 meTEprofileHO = ib.bookProfile(histo, histo, 60, -5., 55., -48., 92., " ");
0146
0147 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HO");
0148 meTEprofileHO_High = ib.bookProfile(histo, histo, 100, -5., 995., -48., 92., " ");
0149
0150 if (imc != 0) {
0151 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HO");
0152 meRecHitSimHitHO = ib.book2D(histo, histo, 150, 0., 1.5, 350, 0., 350.);
0153 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HO");
0154 meRecHitSimHitProfileHO = ib.bookProfile(histo, histo, 150, 0., 1.5, 0., 500., " ");
0155 }
0156 }
0157
0158
0159 if (subdet_ == 4 || subdet_ == 5) {
0160 sprintf(histo, "HcalRecHitTask_energy_of_rechits_HF");
0161 meRecHitsEnergyHF = ib.book1D(histo, histo, 2010, -10., 2000.);
0162
0163 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HF");
0164 meTEprofileHF_Low = ib.bookProfile(histo, histo, 100, -5., 195., -48., 92., " ");
0165
0166 sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HF");
0167 meTEprofileHF = ib.bookProfile(histo, histo, 200, -5., 995., -48., 92., " ");
0168
0169 if (imc != 0) {
0170 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HF");
0171 meRecHitSimHitHF = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0172 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFL");
0173 meRecHitSimHitHFL = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0174 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFS");
0175 meRecHitSimHitHFS = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0176 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HF");
0177 meRecHitSimHitProfileHF = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0178 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFL");
0179 meRecHitSimHitProfileHFL = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0180 sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFS");
0181 meRecHitSimHitProfileHFS = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0182 }
0183 }
0184 }
0185
0186 void HcalRecHitsValidation::analyze(edm::Event const &ev, edm::EventSetup const &c) {
0187 using namespace edm;
0188
0189 const HcalDDDRecConstants *hcons = &c.getData(tok_HRNDC_);
0190
0191
0192
0193
0194
0195 double eHcalConeHB = 0.;
0196 double eHcalConeHE = 0.;
0197 double eHcalConeHO = 0.;
0198 double eHcalConeHF = 0.;
0199 double eHcalConeHFL = 0.;
0200 double eHcalConeHFS = 0.;
0201
0202
0203 int nrechits = 0;
0204 int nrechitsCone = 0;
0205 int nrechitsThresh = 0;
0206
0207
0208 double eEcalCone = 0.;
0209 int numrechitsEcal = 0;
0210
0211
0212 double phi_MC = -999999.;
0213 double eta_MC = -999999.;
0214
0215
0216 double partR = 0.3;
0217
0218 if (imc != 0) {
0219 const edm::Handle<edm::HepMCProduct> &evtMC = ev.getHandle(tok_evt_);
0220 if (!evtMC.isValid()) {
0221 edm::LogVerbatim("HcalRecHitsValidation") << "no HepMCProduct found";
0222 } else {
0223 #ifdef EDM_ML_DEBUG
0224 edm::LogVerbatim("HcalRecHitsValidation") << "*** source HepMCProduct found";
0225 #endif
0226 }
0227
0228
0229 double maxPt = -99999.;
0230 int npart = 0;
0231 const HepMC::GenEvent *myGenEvent = evtMC->GetEvent();
0232 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0233 ++p) {
0234 double phip = (*p)->momentum().phi();
0235 double etap = (*p)->momentum().eta();
0236
0237
0238 double pt = (*p)->momentum().perp();
0239 if (pt > maxPt) {
0240 npart++;
0241 maxPt = pt;
0242 phi_MC = phip;
0243 eta_MC = etap;
0244 }
0245 }
0246 #ifdef EDM_ML_DEBUG
0247 edm::LogVerbatim("HcalRecHitsValidation") << "*** Max pT = " << maxPt;
0248 #endif
0249 }
0250
0251 #ifdef EDM_ML_DEBUG
0252 edm::LogVerbatim("HcalRecHitsValidation") << "*** 2";
0253 #endif
0254
0255 geometry_ = &c.getData(tok_Geom_);
0256
0257
0258 fillRecHitsTmp(subdet_, ev);
0259
0260 #ifdef EDM_ML_DEBUG
0261 edm::LogVerbatim("HcalRecHitsValidation") << "*** 3";
0262 #endif
0263
0264
0265
0266
0267
0268
0269 if (ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
0270 const edm::Handle<EBRecHitCollection> &rhitEB = ev.getHandle(tok_EB_);
0271
0272 EcalRecHitCollection::const_iterator RecHit;
0273 EcalRecHitCollection::const_iterator RecHitEnd;
0274
0275 if (rhitEB.isValid()) {
0276 RecHit = rhitEB.product()->begin();
0277 RecHitEnd = rhitEB.product()->end();
0278
0279 for (; RecHit != RecHitEnd; ++RecHit) {
0280 EBDetId EBid = EBDetId(RecHit->id());
0281
0282 auto cellGeometry = geometry_->getSubdetectorGeometry(EBid)->getGeometry(EBid);
0283 double eta = cellGeometry->getPosition().eta();
0284 double phi = cellGeometry->getPosition().phi();
0285 double en = RecHit->energy();
0286
0287 double r = dR(eta_MC, phi_MC, eta, phi);
0288 if (r < partR) {
0289 eEcalCone += en;
0290 numrechitsEcal++;
0291 }
0292 }
0293 }
0294
0295 const edm::Handle<EERecHitCollection> &rhitEE = ev.getHandle(tok_EE_);
0296
0297 if (rhitEE.isValid()) {
0298 RecHit = rhitEE.product()->begin();
0299 RecHitEnd = rhitEE.product()->end();
0300
0301 for (; RecHit != RecHitEnd; ++RecHit) {
0302 EEDetId EEid = EEDetId(RecHit->id());
0303
0304 auto cellGeometry = geometry_->getSubdetectorGeometry(EEid)->getGeometry(EEid);
0305 double eta = cellGeometry->getPosition().eta();
0306 double phi = cellGeometry->getPosition().phi();
0307 double en = RecHit->energy();
0308
0309 double r = dR(eta_MC, phi_MC, eta, phi);
0310 if (r < partR) {
0311 eEcalCone += en;
0312 numrechitsEcal++;
0313 }
0314 }
0315 }
0316 }
0317
0318 #ifdef EDM_ML_DEBUG
0319 edm::LogVerbatim("HcalRecHitsValidation") << "*** 4";
0320 #endif
0321
0322
0323
0324
0325
0326 if ((subdet_ != 6) && (subdet_ != 0)) {
0327 #ifdef EDM_ML_DEBUG
0328 edm::LogVerbatim("HcalRecHitsValidation") << "*** 6";
0329 #endif
0330
0331 double HcalCone = 0.;
0332
0333 int ietaMax = 9999;
0334 double etaMax = 9999.;
0335
0336
0337
0338 for (unsigned int i = 0; i < cen.size(); i++) {
0339 int sub = csub[i];
0340 int depth = cdepth[i];
0341 double eta = ceta[i];
0342 double phi = cphi[i];
0343 double en = cen[i];
0344 double t = ctime[i];
0345 int ieta = cieta[i];
0346 double chi2 = cchi2[i];
0347
0348 double chi2_log10 = 9.99;
0349 if (chi2 > 0.)
0350 chi2_log10 = log10(chi2);
0351
0352 nrechits++;
0353 if (en > 1.)
0354 nrechitsThresh++;
0355
0356 double r = dR(eta_MC, phi_MC, eta, phi);
0357 if (r < partR) {
0358 if (sub == 1)
0359 eHcalConeHB += en;
0360 if (sub == 2)
0361 eHcalConeHE += en;
0362 if (sub == 3)
0363 eHcalConeHO += en;
0364 if (sub == 4) {
0365 eHcalConeHF += en;
0366 if (depth == 1)
0367 eHcalConeHFL += en;
0368 else
0369 eHcalConeHFS += en;
0370 }
0371 nrechitsCone++;
0372
0373 HcalCone += en;
0374
0375 float eta_diff = fabs(eta_MC - eta);
0376 if (eta_diff < etaMax) {
0377 etaMax = eta_diff;
0378 ietaMax = ieta;
0379 }
0380 }
0381
0382
0383
0384 if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0385 meRecHitsM2Chi2HB->Fill(chi2_log10);
0386 meLog10Chi2profileHB->Fill(en, chi2_log10);
0387
0388 meRecHitsEnergyHB->Fill(en);
0389
0390 meTEprofileHB_Low->Fill(en, t);
0391 meTEprofileHB->Fill(en, t);
0392 meTEprofileHB_High->Fill(en, t);
0393 }
0394 if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0395 meRecHitsM2Chi2HE->Fill(chi2_log10);
0396 meLog10Chi2profileHE->Fill(en, chi2_log10);
0397
0398 meRecHitsEnergyHE->Fill(en);
0399
0400 meTEprofileHE_Low->Fill(en, t);
0401 meTEprofileHE->Fill(en, t);
0402 }
0403 if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0404 meRecHitsEnergyHF->Fill(en);
0405
0406 meTEprofileHF_Low->Fill(en, t);
0407 meTEprofileHF->Fill(en, t);
0408 }
0409 if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0410 meRecHitsEnergyHO->Fill(en);
0411
0412 meTEprofileHO->Fill(en, t);
0413 meTEprofileHO_High->Fill(en, t);
0414 }
0415 }
0416
0417 if (imc != 0) {
0418 meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);
0419 meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0420 meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0421 }
0422
0423 #ifdef EDM_ML_DEBUG
0424 edm::LogVerbatim("HcalRecHitsValidation") << "*** 7";
0425 #endif
0426 }
0427
0428
0429 if (subdet_ > 0 && subdet_ < 6 && imc != 0) {
0430
0431 const edm::Handle<PCaloHitContainer> &hcalHits = ev.getHandle(tok_hh_);
0432 if (hcalHits.isValid()) {
0433 const PCaloHitContainer *SimHitResult = hcalHits.product();
0434
0435 double enSimHitsHB = 0.;
0436 double enSimHitsHE = 0.;
0437 double enSimHitsHO = 0.;
0438 double enSimHitsHF = 0.;
0439 double enSimHitsHFL = 0.;
0440 double enSimHitsHFS = 0.;
0441
0442
0443 for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin(); SimHits != SimHitResult->end();
0444 ++SimHits) {
0445 int sub, depth;
0446 HcalDetId cell;
0447
0448 if (testNumber_)
0449 cell = HcalHitRelabeller::relabel(SimHits->id(), hcons);
0450 else
0451 cell = HcalDetId(SimHits->id());
0452
0453 sub = cell.subdet();
0454 depth = cell.depth();
0455
0456 if (sub != subdet_ && subdet_ != 5)
0457 continue;
0458
0459
0460
0461 const HcalGeometry *cellGeometry =
0462 dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0463 double etaS = cellGeometry->getPosition(cell).eta();
0464 double phiS = cellGeometry->getPosition(cell).phi();
0465 double en = SimHits->energy();
0466
0467 double r = dR(eta_MC, phi_MC, etaS, phiS);
0468
0469 if (r < partR) {
0470
0471 if (sub == static_cast<int>(HcalBarrel))
0472 enSimHitsHB += en;
0473 if (sub == static_cast<int>(HcalEndcap))
0474 enSimHitsHE += en;
0475 if (sub == static_cast<int>(HcalOuter))
0476 enSimHitsHO += en;
0477 if (sub == static_cast<int>(HcalForward)) {
0478 enSimHitsHF += en;
0479 if (depth == 1)
0480 enSimHitsHFL += en;
0481 else
0482 enSimHitsHFS += en;
0483 }
0484 }
0485 }
0486
0487
0488
0489 if (subdet_ == 4 || subdet_ == 5) {
0490 meRecHitSimHitHF->Fill(enSimHitsHF, eHcalConeHF);
0491 meRecHitSimHitProfileHF->Fill(enSimHitsHF, eHcalConeHF);
0492
0493 meRecHitSimHitHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0494 meRecHitSimHitProfileHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0495 meRecHitSimHitHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0496 meRecHitSimHitProfileHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0497 }
0498 if (subdet_ == 1 || subdet_ == 5) {
0499 meRecHitSimHitHB->Fill(enSimHitsHB, eHcalConeHB);
0500 meRecHitSimHitProfileHB->Fill(enSimHitsHB, eHcalConeHB);
0501 }
0502 if (subdet_ == 2 || subdet_ == 5) {
0503 meRecHitSimHitHE->Fill(enSimHitsHE, eHcalConeHE);
0504 meRecHitSimHitProfileHE->Fill(enSimHitsHE, eHcalConeHE);
0505 }
0506 if (subdet_ == 3 || subdet_ == 5) {
0507 meRecHitSimHitHO->Fill(enSimHitsHO, eHcalConeHO);
0508 meRecHitSimHitProfileHO->Fill(enSimHitsHO, eHcalConeHO);
0509 }
0510 }
0511 }
0512
0513 nevtot++;
0514 }
0515
0516
0517 void HcalRecHitsValidation::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0518 using namespace edm;
0519
0520
0521 csub.clear();
0522 cen.clear();
0523 ceta.clear();
0524 cphi.clear();
0525 ctime.clear();
0526 cieta.clear();
0527 ciphi.clear();
0528 cdepth.clear();
0529 cz.clear();
0530 cchi2.clear();
0531
0532 if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0533
0534 const edm::Handle<HBHERecHitCollection> &hbhecoll = ev.getHandle(tok_hbhe_);
0535 if (hbhecoll.isValid()) {
0536 for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
0537 HcalDetId cell(j->id());
0538 const HcalGeometry *cellGeometry =
0539 dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0540 double eta = cellGeometry->getPosition(cell).eta();
0541 double phi = cellGeometry->getPosition(cell).phi();
0542 double zc = cellGeometry->getPosition(cell).z();
0543 int sub = cell.subdet();
0544 int depth = cell.depth();
0545 int inteta = cell.ieta();
0546 int intphi = cell.iphi() - 1;
0547 double en = j->energy();
0548 double t = j->time();
0549 double chi2 = j->chi2();
0550
0551 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0552 csub.push_back(sub);
0553 cen.push_back(en);
0554 ceta.push_back(eta);
0555 cphi.push_back(phi);
0556 ctime.push_back(t);
0557 cieta.push_back(inteta);
0558 ciphi.push_back(intphi);
0559 cdepth.push_back(depth);
0560 cz.push_back(zc);
0561 cchi2.push_back(chi2);
0562 }
0563 }
0564 }
0565 }
0566
0567 if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0568
0569 const edm::Handle<HFRecHitCollection> &hfcoll = ev.getHandle(tok_hf_);
0570 if (hfcoll.isValid()) {
0571 for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
0572 HcalDetId cell(j->id());
0573 auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0574
0575 double eta = cellGeometry->getPosition().eta();
0576 double phi = cellGeometry->getPosition().phi();
0577 double zc = cellGeometry->getPosition().z();
0578 int sub = cell.subdet();
0579 int depth = cell.depth();
0580 int inteta = cell.ieta();
0581 int intphi = cell.iphi() - 1;
0582 double en = j->energy();
0583 double t = j->time();
0584
0585 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0586 csub.push_back(sub);
0587 cen.push_back(en);
0588 ceta.push_back(eta);
0589 cphi.push_back(phi);
0590 ctime.push_back(t);
0591 cieta.push_back(inteta);
0592 ciphi.push_back(intphi);
0593 cdepth.push_back(depth);
0594 cz.push_back(zc);
0595 cchi2.push_back(0.);
0596 }
0597 }
0598 }
0599 }
0600
0601
0602 if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0603 const edm::Handle<HORecHitCollection> &hocoll = ev.getHandle(tok_ho_);
0604 if (hocoll.isValid()) {
0605 for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
0606 HcalDetId cell(j->id());
0607 auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0608
0609 double eta = cellGeometry->getPosition().eta();
0610 double phi = cellGeometry->getPosition().phi();
0611 double zc = cellGeometry->getPosition().z();
0612 int sub = cell.subdet();
0613 int depth = cell.depth();
0614 int inteta = cell.ieta();
0615 int intphi = cell.iphi() - 1;
0616 double t = j->time();
0617 double en = j->energy();
0618
0619 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0620 csub.push_back(sub);
0621 cen.push_back(en);
0622 ceta.push_back(eta);
0623 cphi.push_back(phi);
0624 ctime.push_back(t);
0625 cieta.push_back(inteta);
0626 ciphi.push_back(intphi);
0627 cdepth.push_back(depth);
0628 cz.push_back(zc);
0629 cchi2.push_back(0.);
0630 }
0631 }
0632 }
0633 }
0634 }
0635
0636 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0637 double PI = 3.1415926535898;
0638 double deltaphi = phi1 - phi2;
0639 if (phi2 > phi1) {
0640 deltaphi = phi2 - phi1;
0641 }
0642 if (deltaphi > PI) {
0643 deltaphi = 2. * PI - deltaphi;
0644 }
0645 double deltaeta = eta2 - eta1;
0646 double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0647 return tmp;
0648 }
0649
0650 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
0651
0652
0653 double tmp;
0654 double PI = 3.1415926535898;
0655 double a1 = phi1;
0656 double a2 = phi2;
0657
0658 if (a1 > 0.5 * PI && a2 < 0.)
0659 a2 += 2 * PI;
0660 if (a2 > 0.5 * PI && a1 < 0.)
0661 a1 += 2 * PI;
0662 tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
0663 if (tmp > PI)
0664 tmp -= 2. * PI;
0665
0666 return tmp;
0667 }
0668
0669 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
0670
0671
0672
0673 double PI = 3.1415926535898;
0674 double a1 = phi1;
0675 double a2 = phi2;
0676 double tmp = a2 - a1;
0677 if (a1 * a2 < 0.) {
0678 if (a1 > 0.5 * PI)
0679 tmp += 2. * PI;
0680 if (a2 > 0.5 * PI)
0681 tmp -= 2. * PI;
0682 }
0683 return tmp;
0684 }
0685
0686 DEFINE_FWK_MODULE(HcalRecHitsValidation);