File indexing completed on 2022-03-24 05:59:48
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 eHcal = 0.;
0196 double eHcalCone = 0.;
0197 double eHcalConeHB = 0.;
0198 double eHcalConeHE = 0.;
0199 double eHcalConeHO = 0.;
0200 double eHcalConeHF = 0.;
0201 double eHcalConeHFL = 0.;
0202 double eHcalConeHFS = 0.;
0203
0204
0205 int nrechits = 0;
0206 int nrechitsCone = 0;
0207 int nrechitsThresh = 0;
0208
0209
0210 double eEcal = 0.;
0211 double eEcalB = 0.;
0212 double eEcalE = 0.;
0213 double eEcalCone = 0.;
0214 int numrechitsEcal = 0;
0215
0216
0217 double phi_MC = -999999.;
0218 double eta_MC = -999999.;
0219
0220
0221 double partR = 0.3;
0222
0223 if (imc != 0) {
0224 const edm::Handle<edm::HepMCProduct> &evtMC = ev.getHandle(tok_evt_);
0225 if (!evtMC.isValid()) {
0226 edm::LogVerbatim("HcalRecHitsValidation") << "no HepMCProduct found";
0227 } else {
0228 #ifdef EDM_ML_DEBUG
0229 edm::LogVerbatim("HcalRecHitsValidation") << "*** source HepMCProduct found";
0230 #endif
0231 }
0232
0233
0234 double maxPt = -99999.;
0235 int npart = 0;
0236 const HepMC::GenEvent *myGenEvent = evtMC->GetEvent();
0237 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0238 ++p) {
0239 double phip = (*p)->momentum().phi();
0240 double etap = (*p)->momentum().eta();
0241
0242
0243 double pt = (*p)->momentum().perp();
0244 if (pt > maxPt) {
0245 npart++;
0246 maxPt = pt;
0247 phi_MC = phip;
0248 eta_MC = etap;
0249 }
0250 }
0251 #ifdef EDM_ML_DEBUG
0252 edm::LogVerbatim("HcalRecHitsValidation") << "*** Max pT = " << maxPt;
0253 #endif
0254 }
0255
0256 #ifdef EDM_ML_DEBUG
0257 edm::LogVerbatim("HcalRecHitsValidation") << "*** 2";
0258 #endif
0259
0260 geometry_ = &c.getData(tok_Geom_);
0261
0262
0263 fillRecHitsTmp(subdet_, ev);
0264
0265 #ifdef EDM_ML_DEBUG
0266 edm::LogVerbatim("HcalRecHitsValidation") << "*** 3";
0267 #endif
0268
0269
0270
0271
0272
0273
0274 if (ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
0275 const edm::Handle<EBRecHitCollection> &rhitEB = ev.getHandle(tok_EB_);
0276
0277 EcalRecHitCollection::const_iterator RecHit;
0278 EcalRecHitCollection::const_iterator RecHitEnd;
0279
0280 if (rhitEB.isValid()) {
0281 RecHit = rhitEB.product()->begin();
0282 RecHitEnd = rhitEB.product()->end();
0283
0284 for (; RecHit != RecHitEnd; ++RecHit) {
0285 EBDetId EBid = EBDetId(RecHit->id());
0286
0287 auto cellGeometry = geometry_->getSubdetectorGeometry(EBid)->getGeometry(EBid);
0288 double eta = cellGeometry->getPosition().eta();
0289 double phi = cellGeometry->getPosition().phi();
0290 double en = RecHit->energy();
0291 eEcal += en;
0292 eEcalB += en;
0293
0294 double r = dR(eta_MC, phi_MC, eta, phi);
0295 if (r < partR) {
0296 eEcalCone += en;
0297 numrechitsEcal++;
0298 }
0299 }
0300 }
0301
0302 const edm::Handle<EERecHitCollection> &rhitEE = ev.getHandle(tok_EE_);
0303
0304 if (rhitEE.isValid()) {
0305 RecHit = rhitEE.product()->begin();
0306 RecHitEnd = rhitEE.product()->end();
0307
0308 for (; RecHit != RecHitEnd; ++RecHit) {
0309 EEDetId EEid = EEDetId(RecHit->id());
0310
0311 auto cellGeometry = geometry_->getSubdetectorGeometry(EEid)->getGeometry(EEid);
0312 double eta = cellGeometry->getPosition().eta();
0313 double phi = cellGeometry->getPosition().phi();
0314 double en = RecHit->energy();
0315 eEcal += en;
0316 eEcalE += en;
0317
0318 double r = dR(eta_MC, phi_MC, eta, phi);
0319 if (r < partR) {
0320 eEcalCone += en;
0321 numrechitsEcal++;
0322 }
0323 }
0324 }
0325 }
0326
0327 #ifdef EDM_ML_DEBUG
0328 edm::LogVerbatim("HcalRecHitsValidation") << "*** 4";
0329 #endif
0330
0331
0332
0333
0334
0335 if ((subdet_ != 6) && (subdet_ != 0)) {
0336 #ifdef EDM_ML_DEBUG
0337 edm::LogVerbatim("HcalRecHitsValidation") << "*** 6";
0338 #endif
0339
0340 double HcalCone = 0.;
0341
0342 int ietaMax = 9999;
0343 double etaMax = 9999.;
0344
0345
0346
0347 for (unsigned int i = 0; i < cen.size(); i++) {
0348 int sub = csub[i];
0349 int depth = cdepth[i];
0350 double eta = ceta[i];
0351 double phi = cphi[i];
0352 double en = cen[i];
0353 double t = ctime[i];
0354 int ieta = cieta[i];
0355 double chi2 = cchi2[i];
0356
0357 double chi2_log10 = 9.99;
0358 if (chi2 > 0.)
0359 chi2_log10 = log10(chi2);
0360
0361 nrechits++;
0362 eHcal += en;
0363 if (en > 1.)
0364 nrechitsThresh++;
0365
0366 double r = dR(eta_MC, phi_MC, eta, phi);
0367 if (r < partR) {
0368 if (sub == 1)
0369 eHcalConeHB += en;
0370 if (sub == 2)
0371 eHcalConeHE += en;
0372 if (sub == 3)
0373 eHcalConeHO += en;
0374 if (sub == 4) {
0375 eHcalConeHF += en;
0376 if (depth == 1)
0377 eHcalConeHFL += en;
0378 else
0379 eHcalConeHFS += en;
0380 }
0381 eHcalCone += en;
0382 nrechitsCone++;
0383
0384 HcalCone += en;
0385
0386
0387 float eta_diff = fabs(eta_MC - eta);
0388 if (eta_diff < etaMax) {
0389 etaMax = eta_diff;
0390 ietaMax = ieta;
0391 }
0392 }
0393
0394
0395
0396 if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0397 meRecHitsM2Chi2HB->Fill(chi2_log10);
0398 meLog10Chi2profileHB->Fill(en, chi2_log10);
0399
0400 meRecHitsEnergyHB->Fill(en);
0401
0402 meTEprofileHB_Low->Fill(en, t);
0403 meTEprofileHB->Fill(en, t);
0404 meTEprofileHB_High->Fill(en, t);
0405 }
0406 if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0407 meRecHitsM2Chi2HE->Fill(chi2_log10);
0408 meLog10Chi2profileHE->Fill(en, chi2_log10);
0409
0410 meRecHitsEnergyHE->Fill(en);
0411
0412 meTEprofileHE_Low->Fill(en, t);
0413 meTEprofileHE->Fill(en, t);
0414 }
0415 if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0416 meRecHitsEnergyHF->Fill(en);
0417
0418 meTEprofileHF_Low->Fill(en, t);
0419 meTEprofileHF->Fill(en, t);
0420 }
0421 if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0422 meRecHitsEnergyHO->Fill(en);
0423
0424 meTEprofileHO->Fill(en, t);
0425 meTEprofileHO_High->Fill(en, t);
0426 }
0427 }
0428
0429 if (imc != 0) {
0430 meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);
0431 meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0432 meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0433 }
0434
0435 #ifdef EDM_ML_DEBUG
0436 edm::LogVerbatim("HcalRecHitsValidation") << "*** 7";
0437 #endif
0438 }
0439
0440
0441 if (subdet_ > 0 && subdet_ < 6 && imc != 0) {
0442
0443 const edm::Handle<PCaloHitContainer> &hcalHits = ev.getHandle(tok_hh_);
0444 if (hcalHits.isValid()) {
0445 const PCaloHitContainer *SimHitResult = hcalHits.product();
0446
0447 double enSimHits = 0.;
0448 double enSimHitsHB = 0.;
0449 double enSimHitsHE = 0.;
0450 double enSimHitsHO = 0.;
0451 double enSimHitsHF = 0.;
0452 double enSimHitsHFL = 0.;
0453 double enSimHitsHFS = 0.;
0454
0455
0456 for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin(); SimHits != SimHitResult->end();
0457 ++SimHits) {
0458 int sub, depth;
0459 HcalDetId cell;
0460
0461 if (testNumber_)
0462 cell = HcalHitRelabeller::relabel(SimHits->id(), hcons);
0463 else
0464 cell = HcalDetId(SimHits->id());
0465
0466 sub = cell.subdet();
0467 depth = cell.depth();
0468
0469 if (sub != subdet_ && subdet_ != 5)
0470 continue;
0471
0472
0473
0474 const HcalGeometry *cellGeometry =
0475 dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0476 double etaS = cellGeometry->getPosition(cell).eta();
0477 double phiS = cellGeometry->getPosition(cell).phi();
0478 double en = SimHits->energy();
0479
0480 double r = dR(eta_MC, phi_MC, etaS, phiS);
0481
0482 if (r < partR) {
0483
0484 enSimHits += en;
0485 if (sub == static_cast<int>(HcalBarrel))
0486 enSimHitsHB += en;
0487 if (sub == static_cast<int>(HcalEndcap))
0488 enSimHitsHE += en;
0489 if (sub == static_cast<int>(HcalOuter))
0490 enSimHitsHO += en;
0491 if (sub == static_cast<int>(HcalForward)) {
0492 enSimHitsHF += en;
0493 if (depth == 1)
0494 enSimHitsHFL += en;
0495 else
0496 enSimHitsHFS += en;
0497 }
0498 }
0499 }
0500
0501
0502
0503 if (subdet_ == 4 || subdet_ == 5) {
0504 meRecHitSimHitHF->Fill(enSimHitsHF, eHcalConeHF);
0505 meRecHitSimHitProfileHF->Fill(enSimHitsHF, eHcalConeHF);
0506
0507 meRecHitSimHitHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0508 meRecHitSimHitProfileHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0509 meRecHitSimHitHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0510 meRecHitSimHitProfileHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0511 }
0512 if (subdet_ == 1 || subdet_ == 5) {
0513 meRecHitSimHitHB->Fill(enSimHitsHB, eHcalConeHB);
0514 meRecHitSimHitProfileHB->Fill(enSimHitsHB, eHcalConeHB);
0515 }
0516 if (subdet_ == 2 || subdet_ == 5) {
0517 meRecHitSimHitHE->Fill(enSimHitsHE, eHcalConeHE);
0518 meRecHitSimHitProfileHE->Fill(enSimHitsHE, eHcalConeHE);
0519 }
0520 if (subdet_ == 3 || subdet_ == 5) {
0521 meRecHitSimHitHO->Fill(enSimHitsHO, eHcalConeHO);
0522 meRecHitSimHitProfileHO->Fill(enSimHitsHO, eHcalConeHO);
0523 }
0524 }
0525 }
0526
0527 nevtot++;
0528 }
0529
0530
0531 void HcalRecHitsValidation::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0532 using namespace edm;
0533
0534
0535 csub.clear();
0536 cen.clear();
0537 ceta.clear();
0538 cphi.clear();
0539 ctime.clear();
0540 cieta.clear();
0541 ciphi.clear();
0542 cdepth.clear();
0543 cz.clear();
0544 cchi2.clear();
0545
0546 if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0547
0548 const edm::Handle<HBHERecHitCollection> &hbhecoll = ev.getHandle(tok_hbhe_);
0549 if (hbhecoll.isValid()) {
0550 for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
0551 HcalDetId cell(j->id());
0552 const HcalGeometry *cellGeometry =
0553 dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0554 double eta = cellGeometry->getPosition(cell).eta();
0555 double phi = cellGeometry->getPosition(cell).phi();
0556 double zc = cellGeometry->getPosition(cell).z();
0557 int sub = cell.subdet();
0558 int depth = cell.depth();
0559 int inteta = cell.ieta();
0560 int intphi = cell.iphi() - 1;
0561 double en = j->energy();
0562 double t = j->time();
0563 double chi2 = j->chi2();
0564
0565 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0566 csub.push_back(sub);
0567 cen.push_back(en);
0568 ceta.push_back(eta);
0569 cphi.push_back(phi);
0570 ctime.push_back(t);
0571 cieta.push_back(inteta);
0572 ciphi.push_back(intphi);
0573 cdepth.push_back(depth);
0574 cz.push_back(zc);
0575 cchi2.push_back(chi2);
0576 }
0577 }
0578 }
0579 }
0580
0581 if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0582
0583 const edm::Handle<HFRecHitCollection> &hfcoll = ev.getHandle(tok_hf_);
0584 if (hfcoll.isValid()) {
0585 for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
0586 HcalDetId cell(j->id());
0587 auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0588
0589 double eta = cellGeometry->getPosition().eta();
0590 double phi = cellGeometry->getPosition().phi();
0591 double zc = cellGeometry->getPosition().z();
0592 int sub = cell.subdet();
0593 int depth = cell.depth();
0594 int inteta = cell.ieta();
0595 int intphi = cell.iphi() - 1;
0596 double en = j->energy();
0597 double t = j->time();
0598
0599 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0600 csub.push_back(sub);
0601 cen.push_back(en);
0602 ceta.push_back(eta);
0603 cphi.push_back(phi);
0604 ctime.push_back(t);
0605 cieta.push_back(inteta);
0606 ciphi.push_back(intphi);
0607 cdepth.push_back(depth);
0608 cz.push_back(zc);
0609 cchi2.push_back(0.);
0610 }
0611 }
0612 }
0613 }
0614
0615
0616 if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0617 const edm::Handle<HORecHitCollection> &hocoll = ev.getHandle(tok_ho_);
0618 if (hocoll.isValid()) {
0619 for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
0620 HcalDetId cell(j->id());
0621 auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0622
0623 double eta = cellGeometry->getPosition().eta();
0624 double phi = cellGeometry->getPosition().phi();
0625 double zc = cellGeometry->getPosition().z();
0626 int sub = cell.subdet();
0627 int depth = cell.depth();
0628 int inteta = cell.ieta();
0629 int intphi = cell.iphi() - 1;
0630 double t = j->time();
0631 double en = j->energy();
0632
0633 if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0634 csub.push_back(sub);
0635 cen.push_back(en);
0636 ceta.push_back(eta);
0637 cphi.push_back(phi);
0638 ctime.push_back(t);
0639 cieta.push_back(inteta);
0640 ciphi.push_back(intphi);
0641 cdepth.push_back(depth);
0642 cz.push_back(zc);
0643 cchi2.push_back(0.);
0644 }
0645 }
0646 }
0647 }
0648 }
0649
0650 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0651 double PI = 3.1415926535898;
0652 double deltaphi = phi1 - phi2;
0653 if (phi2 > phi1) {
0654 deltaphi = phi2 - phi1;
0655 }
0656 if (deltaphi > PI) {
0657 deltaphi = 2. * PI - deltaphi;
0658 }
0659 double deltaeta = eta2 - eta1;
0660 double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0661 return tmp;
0662 }
0663
0664 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
0665
0666
0667 double tmp;
0668 double PI = 3.1415926535898;
0669 double a1 = phi1;
0670 double a2 = phi2;
0671
0672 if (a1 > 0.5 * PI && a2 < 0.)
0673 a2 += 2 * PI;
0674 if (a2 > 0.5 * PI && a1 < 0.)
0675 a1 += 2 * PI;
0676 tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
0677 if (tmp > PI)
0678 tmp -= 2. * PI;
0679
0680 return tmp;
0681 }
0682
0683 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
0684
0685
0686
0687 double PI = 3.1415926535898;
0688 double a1 = phi1;
0689 double a2 = phi2;
0690 double tmp = a2 - a1;
0691 if (a1 * a2 < 0.) {
0692 if (a1 > 0.5 * PI)
0693 tmp += 2. * PI;
0694 if (a2 > 0.5 * PI)
0695 tmp -= 2. * PI;
0696 }
0697 return tmp;
0698 }
0699
0700 DEFINE_FWK_MODULE(HcalRecHitsValidation);