Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:12:43 

0001 #include "DQMOffline/Hcal/interface/HcalRecHitsAnalyzer.h"
0002 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0003 
0004 #include "FWCore/Utilities/interface/ESInputTag.h"
0005 #include "FWCore/Utilities/interface/Transition.h"
0006 
0007 HcalRecHitsAnalyzer::HcalRecHitsAnalyzer(edm::ParameterSet const &conf)
0008     : topFolderName_(conf.getParameter<std::string>("TopFolderName")),
0009       hcalDDDRecConstantsToken_{esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()},
0010       caloGeometryRunToken_{esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()},
0011       caloGeometryEventToken_{esConsumes<CaloGeometry, CaloGeometryRecord>()},
0012       hcalTopologyToken_{esConsumes<HcalTopology, HcalRecNumberingRecord>()},
0013       hcalChannelQualityToken_{esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"))},
0014       hcalSeverityLevelComputerToken_{esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>()} {
0015   // DQM ROOT output
0016   outputFile_ = conf.getUntrackedParameter<std::string>("outputFile", "myfile.root");
0017 
0018   if (!outputFile_.empty()) {
0019     edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
0020   } else {
0021     edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
0022   }
0023 
0024   nevtot = 0;
0025 
0026   hcalselector_ = conf.getUntrackedParameter<std::string>("hcalselector", "all");
0027   ecalselector_ = conf.getUntrackedParameter<std::string>("ecalselector", "yes");
0028   eventype_ = conf.getUntrackedParameter<std::string>("eventype", "single");
0029   sign_ = conf.getUntrackedParameter<std::string>("sign", "*");
0030   // useAllHistos_ = conf.getUntrackedParameter<bool>("useAllHistos", false);
0031 
0032   // HEP17 configuration
0033   hep17_ = conf.getUntrackedParameter<bool>("hep17");
0034 
0035   // Collections
0036   tok_hbhe_ = consumes<HBHERecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HBHERecHitCollectionLabel"));
0037   tok_hf_ = consumes<HFRecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HFRecHitCollectionLabel"));
0038   tok_ho_ = consumes<HORecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HORecHitCollectionLabel"));
0039   edm::InputTag EBRecHitCollectionLabel = conf.getParameter<edm::InputTag>("EBRecHitCollectionLabel");
0040   tok_EB_ = consumes<EBRecHitCollection>(EBRecHitCollectionLabel);
0041   edm::InputTag EERecHitCollectionLabel = conf.getParameter<edm::InputTag>("EERecHitCollectionLabel");
0042   tok_EE_ = consumes<EERecHitCollection>(EERecHitCollectionLabel);
0043 
0044   subdet_ = 5;
0045   if (hcalselector_ == "noise")
0046     subdet_ = 0;
0047   if (hcalselector_ == "HB")
0048     subdet_ = 1;
0049   if (hcalselector_ == "HE")
0050     subdet_ = 2;
0051   if (hcalselector_ == "HO")
0052     subdet_ = 3;
0053   if (hcalselector_ == "HF")
0054     subdet_ = 4;
0055   if (hcalselector_ == "all")
0056     subdet_ = 5;
0057   if (hcalselector_ == "ZS")
0058     subdet_ = 6;
0059 
0060   etype_ = 1;
0061   if (eventype_ == "multi")
0062     etype_ = 2;
0063 
0064   iz = 1;
0065   if (sign_ == "-")
0066     iz = -1;
0067   if (sign_ == "*")
0068     iz = 0;
0069 
0070   imc = 0;
0071 }
0072 
0073 void HcalRecHitsAnalyzer::dqmBeginRun(const edm::Run &run, const edm::EventSetup &iSetup) {
0074   HcalDDDRecConstants const &hcons = iSetup.getData(hcalDDDRecConstantsToken_);
0075   maxDepthHB_ = hcons.getMaxDepth(0);
0076   maxDepthHE_ = hcons.getMaxDepth(1);
0077   maxDepthHF_ = std::max(hcons.getMaxDepth(2), 1);
0078   maxDepthHO_ = hcons.getMaxDepth(3);
0079 
0080   CaloGeometry const &geo = iSetup.getData(caloGeometryRunToken_);
0081 
0082   const HcalGeometry *gHB = static_cast<const HcalGeometry *>(geo.getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0083   const HcalGeometry *gHE = static_cast<const HcalGeometry *>(geo.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0084   const HcalGeometry *gHO = static_cast<const HcalGeometry *>(geo.getSubdetectorGeometry(DetId::Hcal, HcalOuter));
0085   const HcalGeometry *gHF = static_cast<const HcalGeometry *>(geo.getSubdetectorGeometry(DetId::Hcal, HcalForward));
0086 
0087   nChannels_[1] = gHB->getHxSize(1);
0088   nChannels_[2] = std::max(int(gHE->getHxSize(2)), 1);
0089   nChannels_[3] = gHO->getHxSize(3);
0090   nChannels_[4] = gHF->getHxSize(4);
0091 
0092   nChannels_[0] = nChannels_[1] + nChannels_[2] + nChannels_[3] + nChannels_[4];
0093 
0094   // std::cout << "Channels HB:" << nChannels_[1] << " HE:" << nChannels_[2] <<
0095   // " HO:" << nChannels_[3] << " HF:" << nChannels_[4] << std::endl;
0096 
0097   // We hardcode the HF depths because in the dual readout configuration,
0098   // rechits are not defined for depths 3&4
0099   maxDepthHF_ = (maxDepthHF_ > 2 ? 2 : maxDepthHF_);  // We retain the dynamic possibility
0100                                                       // that HF might have 0 or 1 depths
0101 
0102   maxDepthAll_ = (maxDepthHB_ + maxDepthHO_ > maxDepthHE_ ? maxDepthHB_ + maxDepthHO_ : maxDepthHE_);
0103   maxDepthAll_ = (maxDepthAll_ > maxDepthHF_ ? maxDepthAll_ : maxDepthHF_);
0104 
0105   // Get Phi segmentation from geometry, use the max phi number so that all iphi
0106   // values are included.
0107 
0108   int NphiMax = hcons.getNPhi(0);
0109 
0110   NphiMax = (hcons.getNPhi(1) > NphiMax ? hcons.getNPhi(1) : NphiMax);
0111   NphiMax = (hcons.getNPhi(2) > NphiMax ? hcons.getNPhi(2) : NphiMax);
0112   NphiMax = (hcons.getNPhi(3) > NphiMax ? hcons.getNPhi(3) : NphiMax);
0113 
0114   // Center the iphi bins on the integers
0115   iphi_min_ = 0.5;
0116   iphi_max_ = NphiMax + 0.5;
0117   iphi_bins_ = (int)(iphi_max_ - iphi_min_);
0118 
0119   // Retain classic behavior, all plots have same ieta range.
0120 
0121   int iEtaMax = (hcons.getEtaRange(0).second > hcons.getEtaRange(1).second ? hcons.getEtaRange(0).second
0122                                                                            : hcons.getEtaRange(1).second);
0123   iEtaMax = (iEtaMax > hcons.getEtaRange(2).second ? iEtaMax : hcons.getEtaRange(2).second);
0124   iEtaMax = (iEtaMax > hcons.getEtaRange(3).second ? iEtaMax : hcons.getEtaRange(3).second);
0125 
0126   // Give an empty bin around the subdet ieta range to make it clear that all
0127   // ieta rings have been included
0128   ieta_min_ = -iEtaMax - 1.5;
0129   ieta_max_ = iEtaMax + 1.5;
0130   ieta_bins_ = (int)(ieta_max_ - ieta_min_);
0131 }
0132 
0133 void HcalRecHitsAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0134                                          edm::Run const & /* iRun*/,
0135                                          edm::EventSetup const &)
0136 
0137 {
0138   Char_t histo[200];
0139 
0140   ibooker.setCurrentFolder(topFolderName_);
0141 
0142   // General counters (drawn)
0143 
0144   // Produce both a total per subdetector, and number of rechits per subdetector
0145   // depth
0146   // The bins are 1 unit wide, and the range is determined by the number of
0147   // channels per subdetector
0148 
0149   for (int depth = 0; depth <= maxDepthHB_; depth++) {
0150     if (depth == 0) {
0151       sprintf(histo, "N_HB");
0152     } else {
0153       sprintf(histo, "N_HB_depth%d", depth);
0154     }
0155     int NBins = (int)(nChannels_[1] * 1.1);
0156     Nhb.push_back(ibooker.book1D(histo, histo, NBins, 0., (float)NBins));
0157   }
0158   for (int depth = 0; depth <= maxDepthHE_; depth++) {
0159     if (depth == 0) {
0160       sprintf(histo, "N_HE");
0161     } else {
0162       sprintf(histo, "N_HE_depth%d", depth);
0163     }
0164     int NBins = (int)(nChannels_[2] * 1.1);
0165     Nhe.push_back(ibooker.book1D(histo, histo, NBins, 0., (float)NBins));
0166   }
0167   for (int depth = 0; depth <= maxDepthHO_; depth++) {
0168     if (depth == 0) {
0169       sprintf(histo, "N_HO");
0170     } else {
0171       sprintf(histo, "N_HO_depth%d", depth);
0172     }
0173     int NBins = (int)(nChannels_[3] * 1.1);
0174     Nho.push_back(ibooker.book1D(histo, histo, NBins, 0., (float)NBins));
0175   }
0176   for (int depth = 0; depth <= maxDepthHF_; depth++) {
0177     if (depth == 0) {
0178       sprintf(histo, "N_HF");
0179     } else {
0180       sprintf(histo, "N_HF_depth%d", depth);
0181     }
0182     int NBins = (int)(nChannels_[4] * 1.1);
0183     Nhf.push_back(ibooker.book1D(histo, histo, NBins, 0., (float)NBins));
0184   }
0185 
0186   // ZS
0187   if (subdet_ == 6) {
0188   }
0189 
0190   // ALL others, except ZS
0191   else {
0192     for (int depth = 1; depth <= maxDepthAll_; depth++) {
0193       sprintf(histo, "emap_depth%d", depth);
0194       emap.push_back(ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_));
0195     }
0196     sprintf(histo, "emap_HO");
0197     emap_HO = ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_);
0198 
0199     // The mean energy histos are drawn, but not the RMS or emean seq
0200 
0201     for (int depth = 1; depth <= maxDepthHB_; depth++) {
0202       sprintf(histo, "emean_vs_ieta_HB%d", depth);
0203       emean_vs_ieta_HB.push_back(ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0204 
0205       sprintf(histo, "emean_vs_ieta_M0_HB%d", depth);
0206       emean_vs_ieta_HBM0.push_back(
0207           ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0208 
0209       sprintf(histo, "emean_vs_ieta_M3_HB%d", depth);
0210       emean_vs_ieta_HBM3.push_back(
0211           ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0212     }
0213     for (int depth = 1; depth <= maxDepthHE_; depth++) {
0214       sprintf(histo, "emean_vs_ieta_HE%d", depth);
0215       emean_vs_ieta_HE.push_back(ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0216 
0217       sprintf(histo, "emean_vs_ieta_M0_HE%d", depth);
0218       emean_vs_ieta_HEM0.push_back(
0219           ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0220 
0221       sprintf(histo, "emean_vs_ieta_M3_HE%d", depth);
0222       emean_vs_ieta_HEM3.push_back(
0223           ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0224     }
0225 
0226     if (hep17_) {
0227       for (int depth = 1; depth <= maxDepthHE_; depth++) {
0228         sprintf(histo, "emean_vs_ieta_HEP17_depth%d", depth);
0229         emean_vs_ieta_HEP17.push_back(
0230             ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0231 
0232         sprintf(histo, "emean_vs_ieta_M0_HEP17_depth%d", depth);
0233         emean_vs_ieta_HEP17M0.push_back(
0234             ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0235 
0236         sprintf(histo, "emean_vs_ieta_M3_HEP17_depth%d", depth);
0237         emean_vs_ieta_HEP17M3.push_back(
0238             ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0239       }
0240     }
0241 
0242     for (int depth = 1; depth <= maxDepthHF_; depth++) {
0243       sprintf(histo, "emean_vs_ieta_HF%d", depth);
0244       emean_vs_ieta_HF.push_back(ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " "));
0245     }
0246     sprintf(histo, "emean_vs_ieta_HO");
0247     emean_vs_ieta_HO = ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -10., 2000., " ");
0248 
0249     // The only occupancy histos drawn are occupancy vs. ieta
0250     // but the maps are needed because this is where the latter are filled from
0251 
0252     for (int depth = 1; depth <= maxDepthHB_; depth++) {
0253       sprintf(histo, "occupancy_map_HB%d", depth);
0254       occupancy_map_HB.push_back(
0255           ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_));
0256     }
0257 
0258     for (int depth = 1; depth <= maxDepthHE_; depth++) {
0259       sprintf(histo, "occupancy_map_HE%d", depth);
0260       occupancy_map_HE.push_back(
0261           ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_));
0262     }
0263 
0264     sprintf(histo, "occupancy_map_HO");
0265     occupancy_map_HO = ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_);
0266 
0267     for (int depth = 1; depth <= maxDepthHF_; depth++) {
0268       sprintf(histo, "occupancy_map_HF%d", depth);
0269       occupancy_map_HF.push_back(
0270           ibooker.book2D(histo, histo, ieta_bins_, ieta_min_, ieta_max_, iphi_bins_, iphi_min_, iphi_max_));
0271     }
0272 
0273     // nrechits vs iphi
0274     for (int depth = 1; depth <= maxDepthHB_; depth++) {
0275       sprintf(histo, "occupancy_vs_ieta_HB%d", depth);
0276       occupancy_vs_ieta_HB.push_back(ibooker.book1D(histo, histo, ieta_bins_, ieta_min_, ieta_max_));
0277       sprintf(histo, "nrechits_vs_iphi_HBP_d%d", depth);
0278       nrechits_vs_iphi_HBP.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0279       sprintf(histo, "nrechits_vs_iphi_HBM_d%d", depth);
0280       nrechits_vs_iphi_HBM.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0281     }
0282 
0283     for (int depth = 1; depth <= maxDepthHE_; depth++) {
0284       sprintf(histo, "occupancy_vs_ieta_HE%d", depth);
0285       occupancy_vs_ieta_HE.push_back(ibooker.book1D(histo, histo, ieta_bins_, ieta_min_, ieta_max_));
0286       sprintf(histo, "nrechits_vs_iphi_HEP_d%d", depth);
0287       nrechits_vs_iphi_HEP.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0288       sprintf(histo, "nrechits_vs_iphi_HEM_d%d", depth);
0289       nrechits_vs_iphi_HEM.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0290     }
0291 
0292     sprintf(histo, "occupancy_vs_ieta_HO");
0293     occupancy_vs_ieta_HO = ibooker.book1D(histo, histo, ieta_bins_, ieta_min_, ieta_max_);
0294     sprintf(histo, "nrechits_vs_iphi_HOP");
0295     nrechits_vs_iphi_HOP = ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_);
0296     sprintf(histo, "nrechits_vs_iphi_HOM");
0297     nrechits_vs_iphi_HOM = ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_);
0298 
0299     for (int depth = 1; depth <= maxDepthHF_; depth++) {
0300       sprintf(histo, "occupancy_vs_ieta_HF%d", depth);
0301       occupancy_vs_ieta_HF.push_back(ibooker.book1D(histo, histo, ieta_bins_, ieta_min_, ieta_max_));
0302       sprintf(histo, "nrechits_vs_iphi_HFP_d%d", depth);
0303       nrechits_vs_iphi_HFP.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0304       sprintf(histo, "nrechits_vs_iphi_HFM_d%d", depth);
0305       nrechits_vs_iphi_HFM.push_back(ibooker.book1D(histo, histo, iphi_bins_, iphi_min_, iphi_max_));
0306     }
0307 
0308     // All status word histos except HF67 are drawn
0309     sprintf(histo, "HcalRecHitTask_RecHit_StatusWord_HB");
0310     RecHit_StatusWord_HB = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0311 
0312     sprintf(histo, "HcalRecHitTask_RecHit_StatusWord_HE");
0313     RecHit_StatusWord_HE = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0314 
0315     sprintf(histo, "HcalRecHitTask_RecHit_StatusWord_HF");
0316     RecHit_StatusWord_HF = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0317 
0318     sprintf(histo, "HcalRecHitTask_RecHit_StatusWord_HO");
0319     RecHit_StatusWord_HO = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0320 
0321     // Aux status word histos
0322     sprintf(histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HB");
0323     RecHit_Aux_StatusWord_HB = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0324 
0325     sprintf(histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HE");
0326     RecHit_Aux_StatusWord_HE = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0327 
0328     sprintf(histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HF");
0329     RecHit_Aux_StatusWord_HF = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0330 
0331     sprintf(histo, "HcalRecHitTask_RecHit_Aux_StatusWord_HO");
0332     RecHit_Aux_StatusWord_HO = ibooker.book1DD(histo, histo, 32, -0.5, 31.5);
0333 
0334   }  // end-of (subdet_ =! 6)
0335 
0336   // Status word correlations
0337   sprintf(histo, "HcalRecHitTask_RecHit_StatusWordCorr_HB");
0338   RecHit_StatusWordCorr_HB = ibooker.book2D(histo, histo, 2, -0.5, 1.5, 2, -0.5, 1.5);
0339 
0340   sprintf(histo, "HcalRecHitTask_RecHit_StatusWordCorr_HE");
0341   RecHit_StatusWordCorr_HE = ibooker.book2D(histo, histo, 2, -0.5, 1.5, 2, -0.5, 1.5);
0342 
0343   //======================= Now various cases one by one ===================
0344 
0345   // Histograms drawn for single pion scan
0346   if (subdet_ != 0 && imc != 0) {  // just not for noise
0347     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths");
0348     meEnConeEtaProfile = ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -100., 2000., " ");
0349 
0350     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_E");
0351     meEnConeEtaProfile_E = ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -100., 2000., " ");
0352 
0353     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_EH");
0354     meEnConeEtaProfile_EH = ibooker.bookProfile(histo, histo, ieta_bins_, ieta_min_, ieta_max_, -100., 2000., " ");
0355   }
0356 
0357   // ************** HB **********************************
0358   if (subdet_ == 1 || subdet_ == 5) {
0359     // Only severity level, energy of rechits and overall HB timing histos are
0360     // drawn
0361 
0362     sprintf(histo, "HcalRecHitTask_severityLevel_HB");
0363     sevLvl_HB = ibooker.book1DD(histo, histo, 25, -0.5, 24.5);
0364 
0365     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HB");
0366     meRecHitsEnergyHB = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0367 
0368     sprintf(histo, "HcalRecHitTask_cleaned_energy_of_rechits_HB");
0369     meRecHitsCleanedEnergyHB = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0370 
0371     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M0_HB");
0372     meRecHitsEnergyHBM0 = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0373 
0374     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3_HB");
0375     meRecHitsEnergyHBM3 = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0376 
0377     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M2vM0_HB");
0378     meRecHitsEnergyM2vM0HB = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0379 
0380     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3vM0_HB");
0381     meRecHitsEnergyM3vM0HB = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0382 
0383     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3vM2_HB");
0384     meRecHitsEnergyM3vM2HB = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0385 
0386     sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HB");
0387     meRecHitsM2Chi2HB = ibooker.book1D(histo, histo, 120, -2., 10.);
0388 
0389     sprintf(histo, "HcalRecHitTask_timing_HB");
0390     meTimeHB = ibooker.book1DD(histo, histo, 70, -48., 92.);
0391 
0392     // High, medium and low histograms to reduce RAM usage
0393     sprintf(histo, "HcalRecHitTask_timing_vs_energy_Low_HB");
0394     meTE_Low_HB = ibooker.book2D(histo, histo, 50, -5., 45., 70, -48., 92.);
0395 
0396     sprintf(histo, "HcalRecHitTask_timing_vs_energy_HB");
0397     meTE_HB = ibooker.book2D(histo, histo, 150, -5., 295., 70, -48., 92.);
0398 
0399     sprintf(histo, "HcalRecHitTask_timing_vs_energy_High_HB");
0400     meTE_High_HB = ibooker.book2D(histo, histo, 150, -5., 2995., 70, -48., 92.);
0401 
0402     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HB");
0403     meTEprofileHB_Low = ibooker.bookProfile(histo, histo, 50, -5., 45., -48., 92., " ");
0404 
0405     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HB");
0406     meTEprofileHB = ibooker.bookProfile(histo, histo, 150, -5., 295., -48., 92., " ");
0407 
0408     sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HB");
0409     meLog10Chi2profileHB = ibooker.bookProfile(histo, histo, 150, -5., 295., -2., 9.9, " ");
0410 
0411     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HB");
0412     meTEprofileHB_High = ibooker.bookProfile(histo, histo, 150, -5., 2995., -48., 92., " ");
0413   }
0414 
0415   // ********************** HE ************************************
0416   if (subdet_ == 2 || subdet_ == 5) {
0417     // Only severity level, energy of rechits and overall HB timing histos are
0418     // drawn
0419     sprintf(histo, "HcalRecHitTask_severityLevel_HE");
0420     sevLvl_HE = ibooker.book1DD(histo, histo, 25, -0.5, 24.5);
0421 
0422     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HE");
0423     meRecHitsEnergyHE = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0424 
0425     sprintf(histo, "HcalRecHitTask_cleaned_energy_of_rechits_HE");
0426     meRecHitsCleanedEnergyHE = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0427 
0428     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M0_HE");
0429     meRecHitsEnergyHEM0 = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0430 
0431     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3_HE");
0432     meRecHitsEnergyHEM3 = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0433 
0434     if (hep17_) {
0435       sprintf(histo, "HcalRecHitTask_energy_of_rechits_HEP17");
0436       meRecHitsEnergyHEP17.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0437 
0438       sprintf(histo, "HcalRecHitTask_energy_of_rechits_M0_HEP17");
0439       meRecHitsEnergyHEP17M0.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0440 
0441       sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3_HEP17");
0442       meRecHitsEnergyHEP17M3.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0443       for (int depth = 1; depth <= maxDepthHE_; depth++) {
0444         sprintf(histo, "HcalRecHitTask_energy_of_rechits_HEP17_depth%d", depth);
0445         meRecHitsEnergyHEP17.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0446 
0447         sprintf(histo, "HcalRecHitTask_energy_of_rechits_M0_HEP17_depth%d", depth);
0448         meRecHitsEnergyHEP17M0.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0449 
0450         sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3_HEP17_depth%d", depth);
0451         meRecHitsEnergyHEP17M3.push_back(ibooker.book1D(histo, histo, 2010, -10., 2000.));
0452       }
0453     }
0454 
0455     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M2vM0_HE");
0456     meRecHitsEnergyM2vM0HE = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0457 
0458     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3vM0_HE");
0459     meRecHitsEnergyM3vM0HE = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0460 
0461     sprintf(histo, "HcalRecHitTask_energy_of_rechits_M3vM2_HE");
0462     meRecHitsEnergyM3vM2HE = ibooker.book2D(histo, histo, 42, -10., 200., 42, -10., 200.);
0463 
0464     sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HE");
0465     meRecHitsM2Chi2HE = ibooker.book1D(histo, histo, 120, -2., 10.);
0466 
0467     sprintf(histo, "HcalRecHitTask_timing_HE");
0468     meTimeHE = ibooker.book1DD(histo, histo, 70, -48., 92.);
0469 
0470     sprintf(histo, "HcalRecHitTask_timing_vs_energy_Low_HE");
0471     meTE_Low_HE = ibooker.book2D(histo, histo, 80, -5., 75., 70, -48., 92.);
0472 
0473     sprintf(histo, "HcalRecHitTask_timing_vs_energy_HE");
0474     meTE_HE = ibooker.book2D(histo, histo, 200, -5., 395., 70, -48., 92.);
0475 
0476     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HE");
0477     meTEprofileHE_Low = ibooker.bookProfile(histo, histo, 80, -5., 75., -48., 92., " ");
0478 
0479     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HE");
0480     meTEprofileHE = ibooker.bookProfile(histo, histo, 200, -5., 395., -48., 92., " ");
0481 
0482     sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HE");
0483     meLog10Chi2profileHE = ibooker.bookProfile(histo, histo, 200, -5., 395., -2., 9.9, " ");
0484   }
0485 
0486   // ************** HO ****************************************
0487   if (subdet_ == 3 || subdet_ == 5) {
0488     // Only severity level, energy of rechits and overall HB timing histos are
0489     // drawn
0490 
0491     sprintf(histo, "HcalRecHitTask_severityLevel_HO");
0492     sevLvl_HO = ibooker.book1DD(histo, histo, 25, -0.5, 24.5);
0493 
0494     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HO");
0495     meRecHitsEnergyHO = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0496 
0497     sprintf(histo, "HcalRecHitTask_cleaned_energy_of_rechits_HO");
0498     meRecHitsCleanedEnergyHO = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0499 
0500     sprintf(histo, "HcalRecHitTask_timing_HO");
0501     meTimeHO = ibooker.book1DD(histo, histo, 80, -80., 80.);
0502 
0503     sprintf(histo, "HcalRecHitTask_timing_vs_energy_HO");
0504     meTE_HO = ibooker.book2D(histo, histo, 60, -5., 55., 80, -80., 80.);
0505 
0506     sprintf(histo, "HcalRecHitTask_timing_vs_energy_High_HO");
0507     meTE_High_HO = ibooker.book2D(histo, histo, 100, -5., 995., 80, -80., 80.);
0508 
0509     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HO");
0510     meTEprofileHO = ibooker.bookProfile(histo, histo, 60, -5., 55., -80., 80., " ");
0511 
0512     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HO");
0513     meTEprofileHO_High = ibooker.bookProfile(histo, histo, 100, -5., 995., -80., 80., " ");
0514   }
0515 
0516   // ********************** HF ************************************
0517   if (subdet_ == 4 || subdet_ == 5) {
0518     // Only severity level, energy of rechits and overall HB timing histos are
0519     // drawn
0520 
0521     sprintf(histo, "HcalRecHitTask_severityLevel_HF");
0522     sevLvl_HF = ibooker.book1DD(histo, histo, 25, -0.5, 24.5);
0523 
0524     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HF");
0525     meRecHitsEnergyHF = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0526 
0527     sprintf(histo, "HcalRecHitTask_cleaned_energy_of_rechits_HF");
0528     meRecHitsCleanedEnergyHF = ibooker.book1DD(histo, histo, 2010, -10., 2000.);
0529 
0530     sprintf(histo, "HcalRecHitTask_timing_HF");
0531     meTimeHF = ibooker.book1DD(histo, histo, 70, -48., 92.);
0532 
0533     sprintf(histo, "HcalRecHitTask_timing_vs_energy_Low_HF");
0534     meTE_Low_HF = ibooker.book2D(histo, histo, 100, -5., 195., 70, -48., 92.);
0535 
0536     sprintf(histo, "HcalRecHitTask_timing_vs_energy_HF");
0537     meTE_HF = ibooker.book2D(histo, histo, 200, -5., 995., 70, -48., 92.);
0538 
0539     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HF");
0540     meTEprofileHF_Low = ibooker.bookProfile(histo, histo, 100, -5., 195., -48., 92., " ");
0541 
0542     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HF");
0543     meTEprofileHF = ibooker.bookProfile(histo, histo, 200, -5., 995., -48., 92., " ");
0544   }
0545 }
0546 
0547 void HcalRecHitsAnalyzer::analyze(edm::Event const &ev, edm::EventSetup const &iSetup) {
0548   using namespace edm;
0549 
0550   // cuts for each subdet_ector mimiking  "Scheme B"
0551   //  double cutHB = 0.9, cutHE = 1.4, cutHO = 1.1, cutHFL = 1.2, cutHFS = 1.8;
0552 
0553   // energy in HCAL
0554   double eHcal = 0.;
0555   // Total numbet of RecHits in HCAL, in the cone, above 1 GeV theshold
0556   int nrechits = 0;
0557   int nrechitsThresh = 0;
0558 
0559   // energy in ECAL
0560   double eEcal = 0.;
0561   double eEcalB = 0.;
0562   double eEcalE = 0.;
0563   double eEcalCone = 0.;
0564 
0565   // HCAL energy around MC eta-phi at all depths;
0566   double partR = 0.3;
0567 
0568   // Single particle samples: actual eta-phi position of cluster around
0569   // hottest cell
0570   double etaHot = 99999.;
0571   double phiHot = 99999.;
0572 
0573   //   previously was:  iSetup.get<IdealGeometryRecord>().get (geometry);
0574   geometry = &iSetup.getData(caloGeometryEventToken_);
0575 
0576   // HCAL Topology **************************************************
0577   theHcalTopology = &iSetup.getData(hcalTopologyToken_);
0578 
0579   // HCAL channel status map ****************************************
0580   theHcalChStatus = &iSetup.getData(hcalChannelQualityToken_);
0581 
0582   // Assignment of severity levels **********************************
0583   theHcalSevLvlComputer = &iSetup.getData(hcalSeverityLevelComputerToken_);
0584 
0585   // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
0586   fillRecHitsTmp(subdet_, ev);
0587 
0588   // HB
0589   if (subdet_ == 5 || subdet_ == 1) {
0590     for (unsigned int iv = 0; iv < hcalHBSevLvlVec.size(); iv++) {
0591       sevLvl_HB->Fill(hcalHBSevLvlVec[iv]);
0592     }
0593   }
0594   // HE
0595   if (subdet_ == 5 || subdet_ == 2) {
0596     for (unsigned int iv = 0; iv < hcalHESevLvlVec.size(); iv++) {
0597       sevLvl_HE->Fill(hcalHESevLvlVec[iv]);
0598     }
0599   }
0600   // HO
0601   if (subdet_ == 5 || subdet_ == 3) {
0602     for (unsigned int iv = 0; iv < hcalHOSevLvlVec.size(); iv++) {
0603       sevLvl_HO->Fill(hcalHOSevLvlVec[iv]);
0604     }
0605   }
0606   // HF
0607   if (subdet_ == 5 || subdet_ == 4) {
0608     for (unsigned int iv = 0; iv < hcalHFSevLvlVec.size(); iv++) {
0609       sevLvl_HF->Fill(hcalHFSevLvlVec[iv]);
0610     }
0611   }
0612 
0613   //===========================================================================
0614   // IN ALL other CASES : ieta-iphi maps
0615   //===========================================================================
0616 
0617   // ECAL
0618   if (ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
0619     Handle<EBRecHitCollection> rhitEB;
0620     if (ev.getByToken(tok_EB_, rhitEB)) {
0621       for (const auto &recHit : *(rhitEB.product())) {
0622         double en = recHit.energy();
0623         eEcal += en;
0624         eEcalB += en;
0625       }
0626     }
0627 
0628     Handle<EERecHitCollection> rhitEE;
0629     if (ev.getByToken(tok_EE_, rhitEE)) {
0630       for (const auto &recHit : *(rhitEE.product())) {
0631         double en = recHit.energy();
0632         eEcal += en;
0633         eEcalE += en;
0634       }
0635     }
0636   }  // end of ECAL selection
0637 
0638   // Counting, including ZS items
0639   // Filling HCAL maps  ----------------------------------------------------
0640   //   double maxE = -99999.;
0641 
0642   // element 0: any depth. element 1,2,..: depth 1,2
0643   std::vector<int> nhb_v(maxDepthHB_ + 1, 0);
0644   std::vector<int> nhe_v(maxDepthHE_ + 1, 0);
0645   std::vector<int> nho_v(maxDepthHO_ + 1, 0);
0646   std::vector<int> nhf_v(maxDepthHF_ + 1, 0);
0647 
0648   for (unsigned int i = 0; i < cen.size(); i++) {
0649     int sub = csub[i];
0650     int depth = cdepth[i];
0651     int ieta = cieta[i];
0652     int iphi = ciphi[i];
0653     double en = cen[i];
0654     double enM0 = cenM0[i];
0655     double enM3 = cenM3[i];
0656     //     double eta    = ceta[i];
0657     //     double phi    = cphi[i];
0658     uint32_t stwd = cstwd[i];
0659     uint32_t auxstwd = cauxstwd[i];
0660     //    double z   = cz[i];
0661 
0662     // This will be true if hep17 == "yes" and the rechit is in the hep17 wedge
0663     bool isHEP17 = (sub == 2) && (iphi >= 63) && (iphi <= 66) && (ieta > 0) && (hep17_);
0664 
0665     // Make sure that an invalid depth won't cause an error. We should probably
0666     // report the problem as well.
0667     if (depth < 1)
0668       continue;
0669     if (sub == 1 && depth > maxDepthHB_)
0670       continue;
0671     if (sub == 2 && depth > maxDepthHE_)
0672       continue;
0673     if (sub == 3 && depth > maxDepthHO_)
0674       continue;
0675     if (sub == 4 && depth > maxDepthHF_)
0676       continue;
0677 
0678     if (sub == 1) {
0679       nhb_v[depth]++;
0680       nhb_v[0]++;
0681     }  // element 0: any depth, element 1,2,..: depth 1,2,...
0682     if (sub == 2) {
0683       nhe_v[depth]++;
0684       nhe_v[0]++;
0685     }  //
0686     if (sub == 3) {
0687       nho_v[depth]++;
0688       nho_v[0]++;
0689     }  //
0690     if (sub == 4) {
0691       nhf_v[depth]++;
0692       nhf_v[0]++;
0693     }  //
0694 
0695     if (subdet_ == 6) {  // ZS specific
0696     }
0697 
0698     if (subdet_ != 6) {
0699       int ieta2 = ieta;
0700       int depth2 = depth;
0701       if (sub == 4) {
0702         if (ieta2 < 0)
0703           ieta2--;
0704         else
0705           ieta2++;
0706       }
0707       if (sub == 3)
0708         emap_HO->Fill(double(ieta2), double(iphi), en);  // HO
0709       else
0710         emap[depth2 - 1]->Fill(double(ieta2), double(iphi), en);  // HB+HE+HF
0711 
0712       // to distinguish HE and HF
0713       if (depth == 1 || depth == 2) {
0714         int ieta1 = ieta;
0715         if (sub == 4) {
0716           if (ieta1 < 0)
0717             ieta1--;
0718           else
0719             ieta1++;
0720         }
0721       }
0722 
0723       if (sub == 1) {
0724         emean_vs_ieta_HB[depth - 1]->Fill(double(ieta), en);
0725         emean_vs_ieta_HBM0[depth - 1]->Fill(double(ieta), enM0);
0726         emean_vs_ieta_HBM3[depth - 1]->Fill(double(ieta), enM3);
0727         occupancy_map_HB[depth - 1]->Fill(double(ieta), double(iphi));
0728         if (ieta > 0)
0729           nrechits_vs_iphi_HBP[depth - 1]->Fill(double(iphi));
0730         else
0731           nrechits_vs_iphi_HBM[depth - 1]->Fill(double(iphi));
0732       }
0733       if (sub == 2) {
0734         if (!isHEP17) {
0735           emean_vs_ieta_HE[depth - 1]->Fill(double(ieta), en);
0736           emean_vs_ieta_HEM0[depth - 1]->Fill(double(ieta), enM0);
0737           emean_vs_ieta_HEM3[depth - 1]->Fill(double(ieta), enM3);
0738         } else {
0739           emean_vs_ieta_HEP17[depth - 1]->Fill(double(ieta), en);
0740           emean_vs_ieta_HEP17M0[depth - 1]->Fill(double(ieta), enM0);
0741           emean_vs_ieta_HEP17M3[depth - 1]->Fill(double(ieta), enM3);
0742         }
0743         occupancy_map_HE[depth - 1]->Fill(double(ieta), double(iphi));
0744         if (ieta > 0)
0745           nrechits_vs_iphi_HEP[depth - 1]->Fill(double(iphi));
0746         else
0747           nrechits_vs_iphi_HEM[depth - 1]->Fill(double(iphi));
0748       }
0749       if (sub == 3) {
0750         emean_vs_ieta_HO->Fill(double(ieta), en);
0751         occupancy_map_HO->Fill(double(ieta), double(iphi));
0752         if (ieta > 0)
0753           nrechits_vs_iphi_HOP->Fill(double(iphi));
0754         else
0755           nrechits_vs_iphi_HOM->Fill(double(iphi));
0756       }
0757       if (sub == 4) {
0758         emean_vs_ieta_HF[depth - 1]->Fill(double(ieta), en);
0759         occupancy_map_HF[depth - 1]->Fill(double(ieta), double(iphi));
0760         if (ieta > 0)
0761           nrechits_vs_iphi_HFP[depth - 1]->Fill(double(iphi));
0762         else
0763           nrechits_vs_iphi_HFM[depth - 1]->Fill(double(iphi));
0764       }
0765     }
0766 
0767     // 32-bit status word
0768     uint32_t statadd;
0769     unsigned int isw67 = 0;
0770 
0771     // Statusword correlation
0772     unsigned int sw27 = 27;
0773     unsigned int sw13 = 13;
0774 
0775     uint32_t statadd27 = 0x1 << sw27;
0776     uint32_t statadd13 = 0x1 << sw13;
0777 
0778     float status27 = 0;
0779     float status13 = 0;
0780 
0781     if (stwd & statadd27)
0782       status27 = 1;
0783     if (stwd & statadd13)
0784       status13 = 1;
0785 
0786     if (sub == 1) {
0787       RecHit_StatusWordCorr_HB->Fill(status13, status27);
0788     } else if (sub == 2) {
0789       RecHit_StatusWordCorr_HE->Fill(status13, status27);
0790     }
0791 
0792     for (unsigned int isw = 0; isw < 32; isw++) {
0793       statadd = 0x1 << (isw);
0794       if (stwd & statadd) {
0795         if (sub == 1)
0796           RecHit_StatusWord_HB->Fill(isw);
0797         else if (sub == 2)
0798           RecHit_StatusWord_HE->Fill(isw);
0799         else if (sub == 3)
0800           RecHit_StatusWord_HO->Fill(isw);
0801         else if (sub == 4) {
0802           RecHit_StatusWord_HF->Fill(isw);
0803           if (isw == 6)
0804             isw67 += 1;
0805           if (isw == 7)
0806             isw67 += 2;
0807         }
0808       }
0809     }
0810 
0811     for (unsigned int isw = 0; isw < 32; isw++) {
0812       statadd = 0x1 << (isw);
0813       if (auxstwd & statadd) {
0814         if (sub == 1)
0815           RecHit_Aux_StatusWord_HB->Fill(isw);
0816         else if (sub == 2)
0817           RecHit_Aux_StatusWord_HE->Fill(isw);
0818         else if (sub == 3)
0819           RecHit_Aux_StatusWord_HO->Fill(isw);
0820         else if (sub == 4)
0821           RecHit_Aux_StatusWord_HF->Fill(isw);
0822       }
0823     }
0824   }
0825 
0826   for (int depth = 0; depth <= maxDepthHB_; depth++)
0827     Nhb[depth]->Fill(double(nhb_v[depth]));
0828   for (int depth = 0; depth <= maxDepthHE_; depth++)
0829     Nhe[depth]->Fill(double(nhe_v[depth]));
0830   for (int depth = 0; depth <= maxDepthHO_; depth++)
0831     Nho[depth]->Fill(double(nho_v[depth]));
0832   for (int depth = 0; depth <= maxDepthHF_; depth++)
0833     Nhf[depth]->Fill(double(nhf_v[depth]));
0834 
0835   //===========================================================================
0836   // SUBSYSTEMS,
0837   //===========================================================================
0838 
0839   if ((subdet_ != 6) && (subdet_ != 0)) {
0840     double clusEta = 999.;
0841     double clusPhi = 999.;
0842     double clusEn = 0.;
0843 
0844     double HcalCone = 0.;
0845 
0846     int ietaMax = 9999;
0847     //     double enMax1 = -9999.;
0848     //     double enMax2 = -9999.;
0849     //     double enMax3 = -9999.;
0850     //     double enMax4 = -9999.;
0851     //     double enMax  = -9999.;
0852     //     double etaMax =  9999.;
0853 
0854     //   CYCLE over cells ====================================================
0855 
0856     for (unsigned int i = 0; i < cen.size(); i++) {
0857       int sub = csub[i];
0858       double eta = ceta[i];
0859       double phi = cphi[i];
0860       double ieta = cieta[i];
0861       double iphi = ciphi[i];
0862       double en = cen[i];
0863       double enM0 = cenM0[i];
0864       double enM3 = cenM3[i];
0865       double chi2 = cchi2[i];
0866       double chi2_log10 = 9.99;  // initial value - stay with this value if chi2<0.
0867       if (chi2 > 0.)
0868         chi2_log10 = log10(chi2);
0869       double t = ctime[i];
0870       double depth = cdepth[i];
0871       int sevlev = csevlev[i];
0872 
0873       bool isHEP17 = (sub == 2) && (iphi >= 63) && (iphi <= 66) && (ieta > 0) && (hep17_);
0874 
0875       //       int   ieta = cieta[i];
0876 
0877       double rhot = dR(etaHot, phiHot, eta, phi);
0878       if (rhot < partR && en > 1.) {
0879         clusEta = (clusEta * clusEn + eta * en) / (clusEn + en);
0880         clusPhi = phi12(clusPhi, clusEn, phi, en);
0881         clusEn += en;
0882       }
0883 
0884       nrechits++;
0885       eHcal += en;
0886 
0887       if (en > 1.)
0888         nrechitsThresh++;
0889 
0890       // The energy and overall timing histos are drawn while
0891       // the ones split by depth are not
0892       if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0893         meTimeHB->Fill(t);
0894         meRecHitsEnergyHB->Fill(en);
0895         if (sevlev <= 9)
0896           meRecHitsCleanedEnergyHB->Fill(en);
0897 
0898         meRecHitsEnergyHBM0->Fill(enM0);
0899         meRecHitsEnergyHBM3->Fill(enM3);
0900 
0901         meRecHitsEnergyM2vM0HB->Fill(enM0, en);
0902         meRecHitsEnergyM3vM0HB->Fill(enM0, enM3);
0903         meRecHitsEnergyM3vM2HB->Fill(en, enM3);
0904 
0905         meRecHitsM2Chi2HB->Fill(chi2_log10);
0906         meLog10Chi2profileHB->Fill(en, chi2_log10);
0907 
0908         meTE_Low_HB->Fill(en, t);
0909         meTE_HB->Fill(en, t);
0910         meTE_High_HB->Fill(en, t);
0911         meTEprofileHB_Low->Fill(en, t);
0912         meTEprofileHB->Fill(en, t);
0913         meTEprofileHB_High->Fill(en, t);
0914       }
0915       if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0916         meTimeHE->Fill(t);
0917         if (!isHEP17) {
0918           meRecHitsEnergyHE->Fill(en);
0919           if (sevlev <= 9)
0920             meRecHitsCleanedEnergyHE->Fill(en);
0921 
0922           meRecHitsEnergyHEM0->Fill(enM0);
0923           meRecHitsEnergyHEM3->Fill(enM3);
0924         } else {
0925           meRecHitsEnergyHEP17[0]->Fill(en);
0926           meRecHitsEnergyHEP17M0[0]->Fill(enM0);
0927           meRecHitsEnergyHEP17M3[0]->Fill(enM3);
0928           meRecHitsEnergyHEP17[depth]->Fill(en);
0929           meRecHitsEnergyHEP17M0[depth]->Fill(enM0);
0930           meRecHitsEnergyHEP17M3[depth]->Fill(enM3);
0931         }
0932 
0933         meRecHitsEnergyM2vM0HE->Fill(enM0, en);
0934         meRecHitsEnergyM3vM0HE->Fill(enM0, enM3);
0935         meRecHitsEnergyM3vM2HE->Fill(en, enM3);
0936 
0937         meRecHitsM2Chi2HE->Fill(chi2_log10);
0938         meLog10Chi2profileHE->Fill(en, chi2_log10);
0939 
0940         meTE_Low_HE->Fill(en, t);
0941         meTE_HE->Fill(en, t);
0942         meTEprofileHE_Low->Fill(en, t);
0943         meTEprofileHE->Fill(en, t);
0944       }
0945       if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0946         meTimeHF->Fill(t);
0947         meRecHitsEnergyHF->Fill(en);
0948         if (sevlev <= 9)
0949           meRecHitsCleanedEnergyHF->Fill(en);
0950 
0951         meTE_Low_HF->Fill(en, t);
0952         meTE_HF->Fill(en, t);
0953         meTEprofileHF_Low->Fill(en, t);
0954         meTEprofileHF->Fill(en, t);
0955       }
0956       if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0957         meTimeHO->Fill(t);
0958         meRecHitsEnergyHO->Fill(en);
0959         if (sevlev <= 9)
0960           meRecHitsCleanedEnergyHO->Fill(en);
0961 
0962         meTE_HO->Fill(en, t);
0963         meTE_High_HO->Fill(en, t);
0964         meTEprofileHO->Fill(en, t);
0965         meTEprofileHO_High->Fill(en, t);
0966       }
0967     }
0968 
0969     if (imc != 0) {
0970       // Cone by depth are not drawn, the others are used for pion scan
0971       meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);  //
0972       meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0973       meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0974     }
0975 
0976     // Single particle samples ONLY !  ======================================
0977     // Fill up some histos for "integrated" subsustems.
0978     // These are not drawn
0979   }
0980 
0981   nevtot++;
0982 }
0983 
0984 ///////////////////////////////////////////////////////////////////////////////
0985 void HcalRecHitsAnalyzer::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0986   using namespace edm;
0987 
0988   // initialize data vectors
0989   csub.clear();
0990   cen.clear();
0991   cenM0.clear();
0992   cenM3.clear();
0993   cchi2.clear();
0994   ceta.clear();
0995   cphi.clear();
0996   ctime.clear();
0997   cieta.clear();
0998   ciphi.clear();
0999   cdepth.clear();
1000   cz.clear();
1001   cstwd.clear();
1002   cauxstwd.clear();
1003   csevlev.clear();
1004   hcalHBSevLvlVec.clear();
1005   hcalHESevLvlVec.clear();
1006   hcalHFSevLvlVec.clear();
1007   hcalHOSevLvlVec.clear();
1008 
1009   if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1010     // HBHE
1011     edm::Handle<HBHERecHitCollection> hbhecoll;
1012     if (ev.getByToken(tok_hbhe_, hbhecoll)) {
1013       for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
1014         HcalDetId cell(j->id());
1015         const HcalGeometry *cellGeometry = dynamic_cast<const HcalGeometry *>(geometry->getSubdetectorGeometry(cell));
1016         double eta = cellGeometry->getPosition(cell).eta();
1017         double phi = cellGeometry->getPosition(cell).phi();
1018         double zc = cellGeometry->getPosition(cell).z();
1019         int sub = cell.subdet();
1020         int depth = cell.depth();
1021         int inteta = cell.ieta();
1022         int intphi = cell.iphi();
1023         double en = j->energy();
1024         double enM0 = j->eraw();
1025         double enM3 = j->eaux();
1026         double chi2 = j->chi2();
1027         double t = j->time();
1028         int stwd = j->flags();
1029         int auxstwd = j->aux();
1030 
1031         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
1032         if (cell.subdet() == HcalBarrel) {
1033           hcalHBSevLvlVec.push_back(severityLevel);
1034         } else if (cell.subdet() == HcalEndcap) {
1035           hcalHESevLvlVec.push_back(severityLevel);
1036         }
1037 
1038         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1039           csub.push_back(sub);
1040           cen.push_back(en);
1041           cenM0.push_back(enM0);
1042           cenM3.push_back(enM3);
1043           cchi2.push_back(chi2);
1044           ceta.push_back(eta);
1045           cphi.push_back(phi);
1046           ctime.push_back(t);
1047           cieta.push_back(inteta);
1048           ciphi.push_back(intphi);
1049           cdepth.push_back(depth);
1050           cz.push_back(zc);
1051           cstwd.push_back(stwd);
1052           cauxstwd.push_back(auxstwd);
1053           csevlev.push_back(severityLevel);
1054         }
1055       }
1056     }
1057   }
1058 
1059   if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1060     // HF
1061     edm::Handle<HFRecHitCollection> hfcoll;
1062     if (ev.getByToken(tok_hf_, hfcoll)) {
1063       for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
1064         HcalDetId cell(j->id());
1065         auto cellGeometry = (geometry->getSubdetectorGeometry(cell))->getGeometry(cell);
1066         double eta = cellGeometry->getPosition().eta();
1067         double phi = cellGeometry->getPosition().phi();
1068         double zc = cellGeometry->getPosition().z();
1069         int sub = cell.subdet();
1070         int depth = cell.depth();
1071         int inteta = cell.ieta();
1072         int intphi = cell.iphi();
1073         double en = j->energy();
1074         double enM0 = 0.;
1075         double enM3 = 0.;
1076         double chi2 = 0.;
1077         double t = j->time();
1078         int stwd = j->flags();
1079         int auxstwd = j->aux();
1080 
1081         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
1082         if (cell.subdet() == HcalForward) {
1083           hcalHFSevLvlVec.push_back(severityLevel);
1084         }
1085 
1086         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1087           csub.push_back(sub);
1088           cen.push_back(en);
1089           cenM0.push_back(enM0);
1090           cenM3.push_back(enM3);
1091           cchi2.push_back(chi2);
1092           ceta.push_back(eta);
1093           cphi.push_back(phi);
1094           ctime.push_back(t);
1095           cieta.push_back(inteta);
1096           ciphi.push_back(intphi);
1097           cdepth.push_back(depth);
1098           cz.push_back(zc);
1099           cstwd.push_back(stwd);
1100           cauxstwd.push_back(auxstwd);
1101           csevlev.push_back(severityLevel);
1102         }
1103       }
1104     }
1105   }
1106 
1107   // HO
1108   if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1109     edm::Handle<HORecHitCollection> hocoll;
1110     if (ev.getByToken(tok_ho_, hocoll)) {
1111       for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
1112         HcalDetId cell(j->id());
1113         auto cellGeometry = (geometry->getSubdetectorGeometry(cell))->getGeometry(cell);
1114         double eta = cellGeometry->getPosition().eta();
1115         double phi = cellGeometry->getPosition().phi();
1116         double zc = cellGeometry->getPosition().z();
1117         int sub = cell.subdet();
1118         int depth = cell.depth();
1119         int inteta = cell.ieta();
1120         int intphi = cell.iphi();
1121         double t = j->time();
1122         double en = j->energy();
1123         double enM0 = 0.;
1124         double enM3 = 0.;
1125         double chi2 = 0.;
1126         int stwd = j->flags();
1127         int auxstwd = j->aux();
1128 
1129         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
1130         if (cell.subdet() == HcalOuter) {
1131           hcalHOSevLvlVec.push_back(severityLevel);
1132         }
1133 
1134         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1135           csub.push_back(sub);
1136           cen.push_back(en);
1137           cenM0.push_back(enM0);
1138           cenM3.push_back(enM3);
1139           cchi2.push_back(chi2);
1140           ceta.push_back(eta);
1141           cphi.push_back(phi);
1142           ctime.push_back(t);
1143           cieta.push_back(inteta);
1144           ciphi.push_back(intphi);
1145           cdepth.push_back(depth);
1146           cz.push_back(zc);
1147           cstwd.push_back(stwd);
1148           cauxstwd.push_back(auxstwd);
1149           csevlev.push_back(severityLevel);
1150         }
1151       }
1152     }
1153   }
1154 }
1155 
1156 double HcalRecHitsAnalyzer::dR(double eta1, double phi1, double eta2, double phi2) {
1157   double PI = 3.1415926535898;
1158   double deltaphi = phi1 - phi2;
1159   if (phi2 > phi1) {
1160     deltaphi = phi2 - phi1;
1161   }
1162   if (deltaphi > PI) {
1163     deltaphi = 2. * PI - deltaphi;
1164   }
1165   double deltaeta = eta2 - eta1;
1166   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
1167   return tmp;
1168 }
1169 
1170 double HcalRecHitsAnalyzer::phi12(double phi1, double en1, double phi2, double en2) {
1171   // weighted mean value of phi1 and phi2
1172 
1173   double tmp;
1174   double PI = 3.1415926535898;
1175   double a1 = phi1;
1176   double a2 = phi2;
1177 
1178   if (a1 > 0.5 * PI && a2 < 0.)
1179     a2 += 2 * PI;
1180   if (a2 > 0.5 * PI && a1 < 0.)
1181     a1 += 2 * PI;
1182   tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
1183   if (tmp > PI)
1184     tmp -= 2. * PI;
1185 
1186   return tmp;
1187 }
1188 
1189 double HcalRecHitsAnalyzer::dPhiWsign(double phi1, double phi2) {
1190   // clockwise      phi2 w.r.t phi1 means "+" phi distance
1191   // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
1192 
1193   double PI = 3.1415926535898;
1194   double a1 = phi1;
1195   double a2 = phi2;
1196   double tmp = a2 - a1;
1197   if (a1 * a2 < 0.) {
1198     if (a1 > 0.5 * PI)
1199       tmp += 2. * PI;
1200     if (a2 > 0.5 * PI)
1201       tmp -= 2. * PI;
1202   }
1203   return tmp;
1204 }
1205 
1206 int HcalRecHitsAnalyzer::hcalSevLvl(const CaloRecHit *hit) {
1207   HcalDetId id = hit->detid();
1208   if (theHcalTopology->getMergePositionFlag() && id.subdet() == HcalEndcap) {
1209     id = theHcalTopology->idFront(id);
1210   }
1211 
1212   const uint32_t recHitFlag = hit->flags();
1213   const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1214 
1215   int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1216 
1217   return severityLevel;
1218 }
1219 
1220 DEFINE_FWK_MODULE(HcalRecHitsAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!