HcalHits/src/HcalSimHitsValidation.cc

0001 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0002 #include "Validation/HcalHits/interface/HcalSimHitsValidation.h"
0003
0004 HcalSimHitsValidation::HcalSimHitsValidation(edm::ParameterSet const &conf)
0005     : g4Label_(conf.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0006       hcalHits_(conf.getUntrackedParameter<std::string>("HcalHitCollection", "HcalHits")),
0007       ebHits_(conf.getUntrackedParameter<std::string>("EBHitCollection", "EcalHitsEB")),
0008       eeHits_(conf.getUntrackedParameter<std::string>("EEHitCollection", "EcalHitsEE")),
0009       hf1_(conf.getParameter<double>("hf1")),
0010       hf2_(conf.getParameter<double>("hf2")),
0011       outputFile_(conf.getUntrackedParameter<std::string>("outputFile", "myfile.root")),
0012       testNumber_(conf.getUntrackedParameter<bool>("TestNumber", false)),
0013       auxPlots_(conf.getUntrackedParameter<bool>("auxiliaryPlots", false)),
0014       tok_evt_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0015       tok_hcal_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
0016       tok_ecalEB_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, ebHits_))),
0017       tok_ecalEE_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, eeHits_))),
0018       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0019       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0020   // DQM ROOT output
0021
0022   // register for data access
0023
0024   // import sampling factors
0025
0026   if (!outputFile_.empty()) {
0027     edm::LogVerbatim("OutputInfo") << " Hcal SimHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
0028   } else {
0029     edm::LogVerbatim("OutputInfo") << " Hcal SimHit Task histograms will NOT be saved";
0030   }
0031
0032   nevtot = 0;
0033 }
0034
0035 void HcalSimHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &run, edm::EventSetup const &es) {
0036   hcons_ = &es.getData(tok_HRNDC_);
0037   maxDepthHB_ = hcons_->getMaxDepth(0);
0038   maxDepthHE_ = hcons_->getMaxDepth(1);
0039   maxDepthHF_ = hcons_->getMaxDepth(2);
0040   maxDepthHO_ = hcons_->getMaxDepth(3);
0041
0042   // Get Phi segmentation from geometry, use the max phi number so that all iphi
0043   // values are included.
0044
0045   //int NphiMax = hcons_->getNPhi(0);
0046
0047   //NphiMax = (hcons_->getNPhi(1) > NphiMax ? hcons_->getNPhi(1) : NphiMax);
0048   //NphiMax = (hcons_->getNPhi(2) > NphiMax ? hcons_->getNPhi(2) : NphiMax);
0049   //NphiMax = (hcons_->getNPhi(3) > NphiMax ? hcons_->getNPhi(3) : NphiMax);
0050
0051   // Center the iphi bins on the integers
0052   // float iphi_min = 0.5;
0053   // float iphi_max = NphiMax + 0.5;
0054   // int iphi_bins = (int) (iphi_max - iphi_min);
0055
0056   int iEtaHBMax = hcons_->getEtaRange(0).second;
0057   int iEtaHEMax = std::max(hcons_->getEtaRange(1).second, 1);
0058   int iEtaHFMax = hcons_->getEtaRange(2).second;
0059   int iEtaHOMax = hcons_->getEtaRange(3).second;
0060
0061   // Retain classic behavior, all plots have same ieta range.
0062   // Comment out    code to allow each subdetector to have its on range
0063
0064   int iEtaMax = (iEtaHBMax > iEtaHEMax ? iEtaHBMax : iEtaHEMax);
0065   iEtaMax = (iEtaMax > iEtaHFMax ? iEtaMax : iEtaHFMax);
0066   iEtaMax = (iEtaMax > iEtaHOMax ? iEtaMax : iEtaHOMax);
0067
0068   iEtaHBMax = iEtaMax;
0069   iEtaHEMax = iEtaMax;
0070   iEtaHFMax = iEtaMax;
0071   iEtaHOMax = iEtaMax;
0072
0073   // Give an empty bin around the subdet ieta range to make it clear that all
0074   // ieta rings have been included
0075   float ieta_min_HB = -iEtaHBMax - 1.5;
0076   float ieta_max_HB = iEtaHBMax + 1.5;
0077   int ieta_bins_HB = (int)(ieta_max_HB - ieta_min_HB);
0078
0079   float ieta_min_HE = -iEtaHEMax - 1.5;
0080   float ieta_max_HE = iEtaHEMax + 1.5;
0081   int ieta_bins_HE = (int)(ieta_max_HE - ieta_min_HE);
0082
0083   float ieta_min_HF = -iEtaHFMax - 1.5;
0084   float ieta_max_HF = iEtaHFMax + 1.5;
0085   int ieta_bins_HF = (int)(ieta_max_HF - ieta_min_HF);
0086
0087   float ieta_min_HO = -iEtaHOMax - 1.5;
0088   float ieta_max_HO = iEtaHOMax + 1.5;
0089   int ieta_bins_HO = (int)(ieta_max_HO - ieta_min_HO);
0090
0091   Char_t histo[200];
0092
0093   ib.setCurrentFolder("HcalHitsV/HcalSimHitTask");
0094
0095   if (auxPlots_) {
0096     // General counters
0097     for (int depth = 0; depth <= maxDepthHB_; depth++) {
0098       if (depth == 0) {
0099         sprintf(histo, "N_HB");
0100       } else {
0101         sprintf(histo, "N_HB%d", depth);
0102       }
0103
0104       Nhb.push_back(ib.book1D(histo, histo, 2600, 0., 2600.));
0105     }
0106     for (int depth = 0; depth <= maxDepthHE_; depth++) {
0107       if (depth == 0) {
0108         sprintf(histo, "N_HE");
0109       } else {
0110         sprintf(histo, "N_HE%d", depth);
0111       }
0112
0113       Nhe.push_back(ib.book1D(histo, histo, 2600, 0., 2600.));
0114     }
0115
0116     sprintf(histo, "N_HO");
0117     Nho = ib.book1D(histo, histo, 2200, 0., 2200.);
0118
0119     for (int depth = 0; depth <= maxDepthHF_; depth++) {
0120       if (depth == 0) {
0121         sprintf(histo, "N_HF");
0122       } else {
0123         sprintf(histo, "N_HF%d", depth);
0124       }
0125
0126       Nhf.push_back(ib.book1D(histo, histo, 1800, 0., 1800.));
0127     }
0128
0129     // Mean energy vs iEta TProfiles
0130     for (int depth = 0; depth <= maxDepthHB_; depth++) {
0131       if (depth == 0) {
0132         sprintf(histo, "emean_vs_ieta_HB");
0133       } else {
0134         sprintf(histo, "emean_vs_ieta_HB%d", depth);
0135       }
0136
0137       emean_vs_ieta_HB.push_back(
0138           ib.bookProfile(histo, histo, ieta_bins_HB, ieta_min_HB, ieta_max_HB, -10., 2000., " "));
0139     }
0140     for (int depth = 0; depth <= maxDepthHE_; depth++) {
0141       if (depth == 0) {
0142         sprintf(histo, "emean_vs_ieta_HE");
0143       } else {
0144         sprintf(histo, "emean_vs_ieta_HE%d", depth);
0145       }
0146
0147       emean_vs_ieta_HE.push_back(
0148           ib.bookProfile(histo, histo, ieta_bins_HE, ieta_min_HE, ieta_max_HE, -10., 2000., " "));
0149     }
0150
0151     sprintf(histo, "emean_vs_ieta_HO");
0152     emean_vs_ieta_HO = ib.bookProfile(histo, histo, ieta_bins_HO, ieta_min_HO, ieta_max_HO, -10., 2000., " ");
0153
0154     for (int depth = 0; depth <= maxDepthHF_; depth++) {
0155       if (depth == 0) {
0156         sprintf(histo, "emean_vs_ieta_HF");
0157       } else {
0158         sprintf(histo, "emean_vs_ieta_HF%d", depth);
0159       }
0160
0161       emean_vs_ieta_HF.push_back(
0162           ib.bookProfile(histo, histo, ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10., 2000., " "));
0163     }
0164
0165     // Occupancy vs. iEta TH1Fs
0166     for (int depth = 0; depth <= maxDepthHB_; depth++) {
0167       if (depth == 0) {
0168         sprintf(histo, "occupancy_vs_ieta_HB");
0169       } else {
0170         sprintf(histo, "occupancy_vs_ieta_HB%d", depth);
0171       }
0172
0173       occupancy_vs_ieta_HB.push_back(ib.book1D(histo, histo, ieta_bins_HB, ieta_min_HB, ieta_max_HB));
0174     }
0175     for (int depth = 0; depth <= maxDepthHE_; depth++) {
0176       if (depth == 0) {
0177         sprintf(histo, "occupancy_vs_ieta_HE");
0178       } else {
0179         sprintf(histo, "occupancy_vs_ieta_HE%d", depth);
0180       }
0181
0182       occupancy_vs_ieta_HE.push_back(ib.book1D(histo, histo, ieta_bins_HE, ieta_min_HE, ieta_max_HE));
0183     }
0184
0185     sprintf(histo, "occupancy_vs_ieta_HO");
0186     occupancy_vs_ieta_HO = ib.book1D(histo, histo, ieta_bins_HO, ieta_min_HO, ieta_max_HO);
0187
0188     for (int depth = 0; depth <= maxDepthHF_; depth++) {
0189       if (depth == 0) {
0190         sprintf(histo, "occupancy_vs_ieta_HF");
0191       } else {
0192         sprintf(histo, "occupancy_vs_ieta_HF%d", depth);
0193       }
0194
0195       occupancy_vs_ieta_HF.push_back(ib.book1D(histo, histo, ieta_bins_HF, ieta_min_HF, ieta_max_HF));
0196     }
0197
0198     // Energy spectra
0199     for (int depth = 0; depth <= maxDepthHB_; depth++) {
0200       if (depth == 0) {
0201         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HB");
0202       } else {
0203         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HB%d", depth);
0204       }
0205
0206       meSimHitsEnergyHB.push_back(ib.book1D(histo, histo, 510, -0.1, 5.));
0207     }
0208     for (int depth = 0; depth <= maxDepthHE_; depth++) {
0209       if (depth == 0) {
0210         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HE");
0211       } else {
0212         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HE%d", depth);
0213       }
0214
0215       meSimHitsEnergyHE.push_back(ib.book1D(histo, histo, 510, -0.1, 5.));
0216     }
0217
0218     sprintf(histo, "HcalSimHitTask_energy_of_simhits_HO");
0219     meSimHitsEnergyHO = ib.book1D(histo, histo, 510, -0.1, 5.);
0220
0221     for (int depth = 0; depth <= maxDepthHF_; depth++) {
0222       if (depth == 0) {
0223         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HF");
0224       } else {
0225         sprintf(histo, "HcalSimHitTask_energy_of_simhits_HF%d", depth);
0226       }
0227
0228       meSimHitsEnergyHF.push_back(ib.book1D(histo, histo, 1010, -5., 500.));
0229     }
0230
0231   }  // auxPlots_
0232
0233   // Energy in Cone
0234   sprintf(histo, "HcalSimHitTask_En_simhits_cone_profile_vs_ieta_all_depths");
0235   meEnConeEtaProfile = ib.bookProfile(histo, histo, ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10., 200., " ");
0236
0237   sprintf(histo, "HcalSimHitTask_En_simhits_cone_profile_vs_ieta_all_depths_E");
0238   meEnConeEtaProfile_E = ib.bookProfile(histo, histo, ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10., 200., " ");
0239
0240   sprintf(histo, "HcalSimHitTask_En_simhits_cone_profile_vs_ieta_all_depths_EH");
0241   meEnConeEtaProfile_EH = ib.bookProfile(histo, histo, ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10., 200., " ");
0242 }
0243
0244 void HcalSimHitsValidation::endJob() {
0245   if (auxPlots_) {
0246     for (int i = 1; i <= occupancy_vs_ieta_HB[0]->getNbinsX(); i++) {
0247       int ieta = i - 43;  // -41 -1, 1 41
0248
0249       float phi_factor;
0250
0251       if (std::abs(ieta) <= 20)
0252         phi_factor = 72.;
0253       else if (std::abs(ieta) < 40)
0254         phi_factor = 36.;
0255       else
0256         phi_factor = 18.;
0257
0258       float cnorm;
0259
0260       // Occupancy vs. iEta TH1Fs
0261       for (int depth = 0; depth <= maxDepthHB_; depth++) {
0262         cnorm = occupancy_vs_ieta_HB[depth]->getBinContent(i) / (phi_factor * nevtot);
0263         occupancy_vs_ieta_HB[depth]->setBinContent(i, cnorm);
0264       }
0265       for (int depth = 0; depth <= maxDepthHE_; depth++) {
0266         cnorm = occupancy_vs_ieta_HE[depth]->getBinContent(i) / (phi_factor * nevtot);
0267         occupancy_vs_ieta_HE[depth]->setBinContent(i, cnorm);
0268       }
0269
0270       cnorm = occupancy_vs_ieta_HO->getBinContent(i) / (phi_factor * nevtot);
0271       occupancy_vs_ieta_HO->setBinContent(i, cnorm);
0272
0273       for (int depth = 0; depth <= maxDepthHF_; depth++) {
0274         cnorm = occupancy_vs_ieta_HF[depth]->getBinContent(i) / (phi_factor * nevtot);
0275         occupancy_vs_ieta_HF[depth]->setBinContent(i, cnorm);
0276       }
0277     }
0278   }
0279
0280   // let's see if this breaks anything
0281   // if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);
0282 }
0283
0284 void HcalSimHitsValidation::analyze(edm::Event const &ev, edm::EventSetup const &c) {
0285   //===========================================================================
0286   // Getting SimHits
0287   //===========================================================================
0288
0289   double phi_MC = -999.;  // phi of initial particle from HepMC
0290   double eta_MC = -999.;  // eta of initial particle from HepMC
0291
0292   const edm::Handle<edm::HepMCProduct> &evtMC = ev.getHandle(tok_evt_);  // generator in late 310_preX
0293   if (!evtMC.isValid()) {
0294     edm::LogVerbatim("OutputInfo") << "no HepMCProduct found";
0295   }
0296
0297   // MC particle with highest pt is taken as a direction reference
0298   double maxPt = -99999.;
0299
0300   const HepMC::GenEvent *myGenEvent = evtMC->GetEvent();
0301   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0302        ++p) {
0303     double phip = (*p)->momentum().phi();
0304     double etap = (*p)->momentum().eta();
0305     double pt = (*p)->momentum().perp();
0306     if (pt > maxPt) {
0307       maxPt = pt;
0308       phi_MC = phip;
0309       eta_MC = etap;
0310     }
0311   }
0312
0313   double partR = 0.3;
0314
0315   // Hcal SimHits
0316
0317   // Approximate calibration constants
0318   const float calib_HB = 120.;
0319   const float calib_HE = 190.;
0320   const float calib_HF1 = hf1_;  // 1.0/0.383;
0321   const float calib_HF2 = hf2_;  // 1.0/0.368;
0322
0323   const edm::Handle<edm::PCaloHitContainer> &hcalHits = ev.getHandle(tok_hcal_);
0324   const auto SimHitResult = hcalHits.product();
0325
0326   float eta_diff;
0327   float etaMax = 9999;
0328   int ietaMax = 0;
0329
0330   double HcalCone = 0;
0331
0332   geometry_ = &c.getData(tok_geom_);
0333
0334   for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin(); SimHits != SimHitResult->end();
0335        ++SimHits) {
0336     HcalDetId cell;
0337     if (testNumber_)
0338       cell = HcalHitRelabeller::relabel(SimHits->id(), hcons_);
0339     else
0340       cell = HcalDetId(SimHits->id());
0341
0342     auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0343     double etaS = cellGeometry->getPosition().eta();
0344     double phiS = cellGeometry->getPosition().phi();
0345     double en = SimHits->energy();
0346
0347     int sub = cell.subdet();
0348     int depth = cell.depth();
0349     double ieta = cell.ieta();
0350
0351     // Energy in Cone
0352     double r = dR(eta_MC, phi_MC, etaS, phiS);
0353
0354     if (r < partR) {
0355       eta_diff = std::abs(eta_MC - etaS);
0356       if (eta_diff < etaMax) {
0357         etaMax = eta_diff;
0358         ietaMax = cell.ieta();
0359       }
0360       // Approximation of calibration
0361       if (sub == 1)
0362         HcalCone += en * calib_HB;
0363       else if (sub == 2)
0364         HcalCone += en * calib_HE;
0365       else if (sub == 4 && (depth == 1 || depth == 3))
0366         HcalCone += en * calib_HF1;
0367       else if (sub == 4 && (depth == 2 || depth == 4))
0368         HcalCone += en * calib_HF2;
0369     }
0370
0371     if (auxPlots_) {
0372       // HB
0373       if (sub == 1) {
0374         meSimHitsEnergyHB[0]->Fill(en);
0375         meSimHitsEnergyHB[depth]->Fill(en);
0376
0377         emean_vs_ieta_HB[0]->Fill(double(ieta), en);
0378         emean_vs_ieta_HB[depth]->Fill(double(ieta), en);
0379
0380         occupancy_vs_ieta_HB[0]->Fill(double(ieta));
0381         occupancy_vs_ieta_HB[depth]->Fill(double(ieta));
0382       }
0383       // HE
0384       if (sub == 2 && maxDepthHE_ > 0) {
0385         meSimHitsEnergyHE[0]->Fill(en);
0386         meSimHitsEnergyHE[depth]->Fill(en);
0387
0388         emean_vs_ieta_HE[0]->Fill(double(ieta), en);
0389         emean_vs_ieta_HE[depth]->Fill(double(ieta), en);
0390
0391         occupancy_vs_ieta_HE[0]->Fill(double(ieta));
0392         occupancy_vs_ieta_HE[depth]->Fill(double(ieta));
0393       }
0394       // HO
0395       if (sub == 3) {
0396         meSimHitsEnergyHO->Fill(en);
0397
0398         emean_vs_ieta_HO->Fill(double(ieta), en);
0399
0400         occupancy_vs_ieta_HO->Fill(double(ieta));
0401       }
0402       // HF
0403       if (sub == 4) {
0404         meSimHitsEnergyHF[0]->Fill(en);
0405         meSimHitsEnergyHF[depth]->Fill(en);
0406
0407         emean_vs_ieta_HF[0]->Fill(double(ieta), en);
0408         emean_vs_ieta_HF[depth]->Fill(double(ieta), en);
0409
0410         occupancy_vs_ieta_HF[0]->Fill(double(ieta));
0411         occupancy_vs_ieta_HF[depth]->Fill(double(ieta));
0412       }
0413
0414     }  // auxPlots_
0415
0416   }  // Loop over SimHits
0417
0418   // Ecal EB SimHits
0419   double EcalCone = 0;
0420
0421   if (!ebHits_.empty()) {
0422     const auto &SimHitResultEB = &ev.get(tok_ecalEB_);
0423
0424     for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResultEB->begin(); SimHits != SimHitResultEB->end();
0425          ++SimHits) {
0426       EBDetId EBid = EBDetId(SimHits->id());
0427
0428       auto cellGeometry = geometry_->getSubdetectorGeometry(EBid)->getGeometry(EBid);
0429       double etaS = cellGeometry->getPosition().eta();
0430       double phiS = cellGeometry->getPosition().phi();
0431       double en = SimHits->energy();
0432
0433       double r = dR(eta_MC, phi_MC, etaS, phiS);
0434
0435       if (r < partR)
0436         EcalCone += en;
0437     }
0438   }  // ebHits_
0439
0440   // Ecal EE SimHits
0441   if (!eeHits_.empty()) {
0442     const auto &SimHitResultEE = &ev.get(tok_ecalEE_);
0443
0444     for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResultEE->begin(); SimHits != SimHitResultEE->end();
0445          ++SimHits) {
0446       EEDetId EEid = EEDetId(SimHits->id());
0447
0448       auto cellGeometry = geometry_->getSubdetectorGeometry(EEid)->getGeometry(EEid);
0449       double etaS = cellGeometry->getPosition().eta();
0450       double phiS = cellGeometry->getPosition().phi();
0451       double en = SimHits->energy();
0452
0453       double r = dR(eta_MC, phi_MC, etaS, phiS);
0454
0455       if (r < partR)
0456         EcalCone += en;
0457     }
0458   }  // eeHits_
0459
0460   if (ietaMax != 0) {  // If ietaMax == 0, there were no good HCAL SimHits
0461     meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);
0462     meEnConeEtaProfile_E->Fill(double(ietaMax), EcalCone);
0463     meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + EcalCone);
0464   }
0465
0466   nevtot++;
0467 }
0468
0469 double HcalSimHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0470   double PI = 3.1415926535898;
0471   double deltaphi = phi1 - phi2;
0472   if (phi2 > phi1) {
0473     deltaphi = phi2 - phi1;
0474   }
0475   if (deltaphi > PI) {
0476     deltaphi = 2. * PI - deltaphi;
0477   }
0478   double deltaeta = eta2 - eta1;
0479   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0480   return tmp;
0481 }
0482
0483 double HcalSimHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
0484   // weighted mean value of phi1 and phi2
0485
0486   double tmp;
0487   double PI = 3.1415926535898;
0488   double a1 = phi1;
0489   double a2 = phi2;
0490
0491   if (a1 > 0.5 * PI && a2 < 0.)
0492     a2 += 2 * PI;
0493   if (a2 > 0.5 * PI && a1 < 0.)
0494     a1 += 2 * PI;
0495   tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
0496   if (tmp > PI)
0497     tmp -= 2. * PI;
0498
0499   return tmp;
0500 }
0501
0502 double HcalSimHitsValidation::dPhiWsign(double phi1, double phi2) {
0503   // clockwise      phi2 w.r.t phi1 means "+" phi distance
0504   // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
0505
0506   double PI = 3.1415926535898;
0507   double a1 = phi1;
0508   double a2 = phi2;
0509   double tmp = a2 - a1;
0510   if (a1 * a2 < 0.) {
0511     if (a1 > 0.5 * PI)
0512       tmp += 2. * PI;
0513     if (a2 > 0.5 * PI)
0514       tmp -= 2. * PI;
0515   }
0516   return tmp;
0517 }
0518
0519 DEFINE_FWK_MODULE(HcalSimHitsValidation);