Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-13 23:39: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   // Total numbet of RecHits in HCAL, in the cone, above 1 GeV theshold
0554   int nrechits = 0;
0555   int nrechitsThresh = 0;
0556 
0557   // energy in ECAL
0558   double eEcalCone = 0.;
0559 
0560   // HCAL energy around MC eta-phi at all depths;
0561   double partR = 0.3;
0562 
0563   // Single particle samples: actual eta-phi position of cluster around
0564   // hottest cell
0565   double etaHot = 99999.;
0566   double phiHot = 99999.;
0567 
0568   //   previously was:  iSetup.get<IdealGeometryRecord>().get (geometry);
0569   geometry = &iSetup.getData(caloGeometryEventToken_);
0570 
0571   // HCAL Topology **************************************************
0572   theHcalTopology = &iSetup.getData(hcalTopologyToken_);
0573 
0574   // HCAL channel status map ****************************************
0575   theHcalChStatus = &iSetup.getData(hcalChannelQualityToken_);
0576 
0577   // Assignment of severity levels **********************************
0578   theHcalSevLvlComputer = &iSetup.getData(hcalSeverityLevelComputerToken_);
0579 
0580   // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
0581   fillRecHitsTmp(subdet_, ev);
0582 
0583   // HB
0584   if (subdet_ == 5 || subdet_ == 1) {
0585     for (unsigned int iv = 0; iv < hcalHBSevLvlVec.size(); iv++) {
0586       sevLvl_HB->Fill(hcalHBSevLvlVec[iv]);
0587     }
0588   }
0589   // HE
0590   if (subdet_ == 5 || subdet_ == 2) {
0591     for (unsigned int iv = 0; iv < hcalHESevLvlVec.size(); iv++) {
0592       sevLvl_HE->Fill(hcalHESevLvlVec[iv]);
0593     }
0594   }
0595   // HO
0596   if (subdet_ == 5 || subdet_ == 3) {
0597     for (unsigned int iv = 0; iv < hcalHOSevLvlVec.size(); iv++) {
0598       sevLvl_HO->Fill(hcalHOSevLvlVec[iv]);
0599     }
0600   }
0601   // HF
0602   if (subdet_ == 5 || subdet_ == 4) {
0603     for (unsigned int iv = 0; iv < hcalHFSevLvlVec.size(); iv++) {
0604       sevLvl_HF->Fill(hcalHFSevLvlVec[iv]);
0605     }
0606   }
0607 
0608   // Counting, including ZS items
0609   // Filling HCAL maps  ----------------------------------------------------
0610   //   double maxE = -99999.;
0611 
0612   // element 0: any depth. element 1,2,..: depth 1,2
0613   std::vector<int> nhb_v(maxDepthHB_ + 1, 0);
0614   std::vector<int> nhe_v(maxDepthHE_ + 1, 0);
0615   std::vector<int> nho_v(maxDepthHO_ + 1, 0);
0616   std::vector<int> nhf_v(maxDepthHF_ + 1, 0);
0617 
0618   for (unsigned int i = 0; i < cen.size(); i++) {
0619     int sub = csub[i];
0620     int depth = cdepth[i];
0621     int ieta = cieta[i];
0622     int iphi = ciphi[i];
0623     double en = cen[i];
0624     double enM0 = cenM0[i];
0625     double enM3 = cenM3[i];
0626     //     double eta    = ceta[i];
0627     //     double phi    = cphi[i];
0628     uint32_t stwd = cstwd[i];
0629     uint32_t auxstwd = cauxstwd[i];
0630     //    double z   = cz[i];
0631 
0632     // This will be true if hep17 == "yes" and the rechit is in the hep17 wedge
0633     bool isHEP17 = (sub == 2) && (iphi >= 63) && (iphi <= 66) && (ieta > 0) && (hep17_);
0634 
0635     // Make sure that an invalid depth won't cause an error. We should probably
0636     // report the problem as well.
0637     if (depth < 1)
0638       continue;
0639     if (sub == 1 && depth > maxDepthHB_)
0640       continue;
0641     if (sub == 2 && depth > maxDepthHE_)
0642       continue;
0643     if (sub == 3 && depth > maxDepthHO_)
0644       continue;
0645     if (sub == 4 && depth > maxDepthHF_)
0646       continue;
0647 
0648     if (sub == 1) {
0649       nhb_v[depth]++;
0650       nhb_v[0]++;
0651     }  // element 0: any depth, element 1,2,..: depth 1,2,...
0652     if (sub == 2) {
0653       nhe_v[depth]++;
0654       nhe_v[0]++;
0655     }  //
0656     if (sub == 3) {
0657       nho_v[depth]++;
0658       nho_v[0]++;
0659     }  //
0660     if (sub == 4) {
0661       nhf_v[depth]++;
0662       nhf_v[0]++;
0663     }  //
0664 
0665     if (subdet_ == 6) {  // ZS specific
0666     }
0667 
0668     if (subdet_ != 6) {
0669       int ieta2 = ieta;
0670       int depth2 = depth;
0671       if (sub == 4) {
0672         if (ieta2 < 0)
0673           ieta2--;
0674         else
0675           ieta2++;
0676       }
0677       if (sub == 3)
0678         emap_HO->Fill(double(ieta2), double(iphi), en);  // HO
0679       else
0680         emap[depth2 - 1]->Fill(double(ieta2), double(iphi), en);  // HB+HE+HF
0681 
0682       // to distinguish HE and HF
0683       if (depth == 1 || depth == 2) {
0684         int ieta1 = ieta;
0685         if (sub == 4) {
0686           if (ieta1 < 0)
0687             ieta1--;
0688           else
0689             ieta1++;
0690         }
0691       }
0692 
0693       if (sub == 1) {
0694         emean_vs_ieta_HB[depth - 1]->Fill(double(ieta), en);
0695         emean_vs_ieta_HBM0[depth - 1]->Fill(double(ieta), enM0);
0696         emean_vs_ieta_HBM3[depth - 1]->Fill(double(ieta), enM3);
0697         occupancy_map_HB[depth - 1]->Fill(double(ieta), double(iphi));
0698         if (ieta > 0)
0699           nrechits_vs_iphi_HBP[depth - 1]->Fill(double(iphi));
0700         else
0701           nrechits_vs_iphi_HBM[depth - 1]->Fill(double(iphi));
0702       }
0703       if (sub == 2) {
0704         if (!isHEP17) {
0705           emean_vs_ieta_HE[depth - 1]->Fill(double(ieta), en);
0706           emean_vs_ieta_HEM0[depth - 1]->Fill(double(ieta), enM0);
0707           emean_vs_ieta_HEM3[depth - 1]->Fill(double(ieta), enM3);
0708         } else {
0709           emean_vs_ieta_HEP17[depth - 1]->Fill(double(ieta), en);
0710           emean_vs_ieta_HEP17M0[depth - 1]->Fill(double(ieta), enM0);
0711           emean_vs_ieta_HEP17M3[depth - 1]->Fill(double(ieta), enM3);
0712         }
0713         occupancy_map_HE[depth - 1]->Fill(double(ieta), double(iphi));
0714         if (ieta > 0)
0715           nrechits_vs_iphi_HEP[depth - 1]->Fill(double(iphi));
0716         else
0717           nrechits_vs_iphi_HEM[depth - 1]->Fill(double(iphi));
0718       }
0719       if (sub == 3) {
0720         emean_vs_ieta_HO->Fill(double(ieta), en);
0721         occupancy_map_HO->Fill(double(ieta), double(iphi));
0722         if (ieta > 0)
0723           nrechits_vs_iphi_HOP->Fill(double(iphi));
0724         else
0725           nrechits_vs_iphi_HOM->Fill(double(iphi));
0726       }
0727       if (sub == 4) {
0728         emean_vs_ieta_HF[depth - 1]->Fill(double(ieta), en);
0729         occupancy_map_HF[depth - 1]->Fill(double(ieta), double(iphi));
0730         if (ieta > 0)
0731           nrechits_vs_iphi_HFP[depth - 1]->Fill(double(iphi));
0732         else
0733           nrechits_vs_iphi_HFM[depth - 1]->Fill(double(iphi));
0734       }
0735     }
0736 
0737     // 32-bit status word
0738     uint32_t statadd;
0739 
0740     // Statusword correlation
0741     unsigned int sw27 = 27;
0742     unsigned int sw13 = 13;
0743 
0744     uint32_t statadd27 = 0x1 << sw27;
0745     uint32_t statadd13 = 0x1 << sw13;
0746 
0747     float status27 = 0;
0748     float status13 = 0;
0749 
0750     if (stwd & statadd27)
0751       status27 = 1;
0752     if (stwd & statadd13)
0753       status13 = 1;
0754 
0755     if (sub == 1) {
0756       RecHit_StatusWordCorr_HB->Fill(status13, status27);
0757     } else if (sub == 2) {
0758       RecHit_StatusWordCorr_HE->Fill(status13, status27);
0759     }
0760 
0761     for (unsigned int isw = 0; isw < 32; isw++) {
0762       statadd = 0x1 << (isw);
0763       if (stwd & statadd) {
0764         if (sub == 1)
0765           RecHit_StatusWord_HB->Fill(isw);
0766         else if (sub == 2)
0767           RecHit_StatusWord_HE->Fill(isw);
0768         else if (sub == 3)
0769           RecHit_StatusWord_HO->Fill(isw);
0770         else if (sub == 4) {
0771           RecHit_StatusWord_HF->Fill(isw);
0772         }
0773       }
0774     }
0775 
0776     for (unsigned int isw = 0; isw < 32; isw++) {
0777       statadd = 0x1 << (isw);
0778       if (auxstwd & statadd) {
0779         if (sub == 1)
0780           RecHit_Aux_StatusWord_HB->Fill(isw);
0781         else if (sub == 2)
0782           RecHit_Aux_StatusWord_HE->Fill(isw);
0783         else if (sub == 3)
0784           RecHit_Aux_StatusWord_HO->Fill(isw);
0785         else if (sub == 4)
0786           RecHit_Aux_StatusWord_HF->Fill(isw);
0787       }
0788     }
0789   }
0790 
0791   for (int depth = 0; depth <= maxDepthHB_; depth++)
0792     Nhb[depth]->Fill(double(nhb_v[depth]));
0793   for (int depth = 0; depth <= maxDepthHE_; depth++)
0794     Nhe[depth]->Fill(double(nhe_v[depth]));
0795   for (int depth = 0; depth <= maxDepthHO_; depth++)
0796     Nho[depth]->Fill(double(nho_v[depth]));
0797   for (int depth = 0; depth <= maxDepthHF_; depth++)
0798     Nhf[depth]->Fill(double(nhf_v[depth]));
0799 
0800   //===========================================================================
0801   // SUBSYSTEMS,
0802   //===========================================================================
0803 
0804   if ((subdet_ != 6) && (subdet_ != 0)) {
0805     double clusEta = 999.;
0806     double clusPhi = 999.;
0807     double clusEn = 0.;
0808 
0809     double HcalCone = 0.;
0810 
0811     int ietaMax = 9999;
0812     //     double enMax1 = -9999.;
0813     //     double enMax2 = -9999.;
0814     //     double enMax3 = -9999.;
0815     //     double enMax4 = -9999.;
0816     //     double enMax  = -9999.;
0817     //     double etaMax =  9999.;
0818 
0819     //   CYCLE over cells ====================================================
0820 
0821     for (unsigned int i = 0; i < cen.size(); i++) {
0822       int sub = csub[i];
0823       double eta = ceta[i];
0824       double phi = cphi[i];
0825       double ieta = cieta[i];
0826       double iphi = ciphi[i];
0827       double en = cen[i];
0828       double enM0 = cenM0[i];
0829       double enM3 = cenM3[i];
0830       double chi2 = cchi2[i];
0831       double chi2_log10 = 9.99;  // initial value - stay with this value if chi2<0.
0832       if (chi2 > 0.)
0833         chi2_log10 = log10(chi2);
0834       double t = ctime[i];
0835       double depth = cdepth[i];
0836       int sevlev = csevlev[i];
0837 
0838       bool isHEP17 = (sub == 2) && (iphi >= 63) && (iphi <= 66) && (ieta > 0) && (hep17_);
0839 
0840       //       int   ieta = cieta[i];
0841 
0842       double rhot = dR(etaHot, phiHot, eta, phi);
0843       if (rhot < partR && en > 1.) {
0844         clusEta = (clusEta * clusEn + eta * en) / (clusEn + en);
0845         clusPhi = phi12(clusPhi, clusEn, phi, en);
0846         clusEn += en;
0847       }
0848 
0849       nrechits++;
0850 
0851       if (en > 1.)
0852         nrechitsThresh++;
0853 
0854       // The energy and overall timing histos are drawn while
0855       // the ones split by depth are not
0856       if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0857         meTimeHB->Fill(t);
0858         meRecHitsEnergyHB->Fill(en);
0859         if (sevlev <= 9)
0860           meRecHitsCleanedEnergyHB->Fill(en);
0861 
0862         meRecHitsEnergyHBM0->Fill(enM0);
0863         meRecHitsEnergyHBM3->Fill(enM3);
0864 
0865         meRecHitsEnergyM2vM0HB->Fill(enM0, en);
0866         meRecHitsEnergyM3vM0HB->Fill(enM0, enM3);
0867         meRecHitsEnergyM3vM2HB->Fill(en, enM3);
0868 
0869         meRecHitsM2Chi2HB->Fill(chi2_log10);
0870         meLog10Chi2profileHB->Fill(en, chi2_log10);
0871 
0872         meTE_Low_HB->Fill(en, t);
0873         meTE_HB->Fill(en, t);
0874         meTE_High_HB->Fill(en, t);
0875         meTEprofileHB_Low->Fill(en, t);
0876         meTEprofileHB->Fill(en, t);
0877         meTEprofileHB_High->Fill(en, t);
0878       }
0879       if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0880         meTimeHE->Fill(t);
0881         if (!isHEP17) {
0882           meRecHitsEnergyHE->Fill(en);
0883           if (sevlev <= 9)
0884             meRecHitsCleanedEnergyHE->Fill(en);
0885 
0886           meRecHitsEnergyHEM0->Fill(enM0);
0887           meRecHitsEnergyHEM3->Fill(enM3);
0888         } else {
0889           meRecHitsEnergyHEP17[0]->Fill(en);
0890           meRecHitsEnergyHEP17M0[0]->Fill(enM0);
0891           meRecHitsEnergyHEP17M3[0]->Fill(enM3);
0892           meRecHitsEnergyHEP17[depth]->Fill(en);
0893           meRecHitsEnergyHEP17M0[depth]->Fill(enM0);
0894           meRecHitsEnergyHEP17M3[depth]->Fill(enM3);
0895         }
0896 
0897         meRecHitsEnergyM2vM0HE->Fill(enM0, en);
0898         meRecHitsEnergyM3vM0HE->Fill(enM0, enM3);
0899         meRecHitsEnergyM3vM2HE->Fill(en, enM3);
0900 
0901         meRecHitsM2Chi2HE->Fill(chi2_log10);
0902         meLog10Chi2profileHE->Fill(en, chi2_log10);
0903 
0904         meTE_Low_HE->Fill(en, t);
0905         meTE_HE->Fill(en, t);
0906         meTEprofileHE_Low->Fill(en, t);
0907         meTEprofileHE->Fill(en, t);
0908       }
0909       if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0910         meTimeHF->Fill(t);
0911         meRecHitsEnergyHF->Fill(en);
0912         if (sevlev <= 9)
0913           meRecHitsCleanedEnergyHF->Fill(en);
0914 
0915         meTE_Low_HF->Fill(en, t);
0916         meTE_HF->Fill(en, t);
0917         meTEprofileHF_Low->Fill(en, t);
0918         meTEprofileHF->Fill(en, t);
0919       }
0920       if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0921         meTimeHO->Fill(t);
0922         meRecHitsEnergyHO->Fill(en);
0923         if (sevlev <= 9)
0924           meRecHitsCleanedEnergyHO->Fill(en);
0925 
0926         meTE_HO->Fill(en, t);
0927         meTE_High_HO->Fill(en, t);
0928         meTEprofileHO->Fill(en, t);
0929         meTEprofileHO_High->Fill(en, t);
0930       }
0931     }
0932 
0933     if (imc != 0) {
0934       // Cone by depth are not drawn, the others are used for pion scan
0935       meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);  //
0936       meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0937       meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0938     }
0939 
0940     // Single particle samples ONLY !  ======================================
0941     // Fill up some histos for "integrated" subsustems.
0942     // These are not drawn
0943   }
0944 
0945   nevtot++;
0946 }
0947 
0948 ///////////////////////////////////////////////////////////////////////////////
0949 void HcalRecHitsAnalyzer::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0950   using namespace edm;
0951 
0952   // initialize data vectors
0953   csub.clear();
0954   cen.clear();
0955   cenM0.clear();
0956   cenM3.clear();
0957   cchi2.clear();
0958   ceta.clear();
0959   cphi.clear();
0960   ctime.clear();
0961   cieta.clear();
0962   ciphi.clear();
0963   cdepth.clear();
0964   cz.clear();
0965   cstwd.clear();
0966   cauxstwd.clear();
0967   csevlev.clear();
0968   hcalHBSevLvlVec.clear();
0969   hcalHESevLvlVec.clear();
0970   hcalHFSevLvlVec.clear();
0971   hcalHOSevLvlVec.clear();
0972 
0973   if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0974     // HBHE
0975     edm::Handle<HBHERecHitCollection> hbhecoll;
0976     if (ev.getByToken(tok_hbhe_, hbhecoll)) {
0977       for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
0978         HcalDetId cell(j->id());
0979         const HcalGeometry *cellGeometry = dynamic_cast<const HcalGeometry *>(geometry->getSubdetectorGeometry(cell));
0980         double eta = cellGeometry->getPosition(cell).eta();
0981         double phi = cellGeometry->getPosition(cell).phi();
0982         double zc = cellGeometry->getPosition(cell).z();
0983         int sub = cell.subdet();
0984         int depth = cell.depth();
0985         int inteta = cell.ieta();
0986         int intphi = cell.iphi();
0987         double en = j->energy();
0988         double enM0 = j->eraw();
0989         double enM3 = j->eaux();
0990         double chi2 = j->chi2();
0991         double t = j->time();
0992         int stwd = j->flags();
0993         int auxstwd = j->aux();
0994 
0995         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
0996         if (cell.subdet() == HcalBarrel) {
0997           hcalHBSevLvlVec.push_back(severityLevel);
0998         } else if (cell.subdet() == HcalEndcap) {
0999           hcalHESevLvlVec.push_back(severityLevel);
1000         }
1001 
1002         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1003           csub.push_back(sub);
1004           cen.push_back(en);
1005           cenM0.push_back(enM0);
1006           cenM3.push_back(enM3);
1007           cchi2.push_back(chi2);
1008           ceta.push_back(eta);
1009           cphi.push_back(phi);
1010           ctime.push_back(t);
1011           cieta.push_back(inteta);
1012           ciphi.push_back(intphi);
1013           cdepth.push_back(depth);
1014           cz.push_back(zc);
1015           cstwd.push_back(stwd);
1016           cauxstwd.push_back(auxstwd);
1017           csevlev.push_back(severityLevel);
1018         }
1019       }
1020     }
1021   }
1022 
1023   if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1024     // HF
1025     edm::Handle<HFRecHitCollection> hfcoll;
1026     if (ev.getByToken(tok_hf_, hfcoll)) {
1027       for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
1028         HcalDetId cell(j->id());
1029         auto cellGeometry = (geometry->getSubdetectorGeometry(cell))->getGeometry(cell);
1030         double eta = cellGeometry->getPosition().eta();
1031         double phi = cellGeometry->getPosition().phi();
1032         double zc = cellGeometry->getPosition().z();
1033         int sub = cell.subdet();
1034         int depth = cell.depth();
1035         int inteta = cell.ieta();
1036         int intphi = cell.iphi();
1037         double en = j->energy();
1038         double enM0 = 0.;
1039         double enM3 = 0.;
1040         double chi2 = 0.;
1041         double t = j->time();
1042         int stwd = j->flags();
1043         int auxstwd = j->aux();
1044 
1045         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
1046         if (cell.subdet() == HcalForward) {
1047           hcalHFSevLvlVec.push_back(severityLevel);
1048         }
1049 
1050         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1051           csub.push_back(sub);
1052           cen.push_back(en);
1053           cenM0.push_back(enM0);
1054           cenM3.push_back(enM3);
1055           cchi2.push_back(chi2);
1056           ceta.push_back(eta);
1057           cphi.push_back(phi);
1058           ctime.push_back(t);
1059           cieta.push_back(inteta);
1060           ciphi.push_back(intphi);
1061           cdepth.push_back(depth);
1062           cz.push_back(zc);
1063           cstwd.push_back(stwd);
1064           cauxstwd.push_back(auxstwd);
1065           csevlev.push_back(severityLevel);
1066         }
1067       }
1068     }
1069   }
1070 
1071   // HO
1072   if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
1073     edm::Handle<HORecHitCollection> hocoll;
1074     if (ev.getByToken(tok_ho_, hocoll)) {
1075       for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
1076         HcalDetId cell(j->id());
1077         auto cellGeometry = (geometry->getSubdetectorGeometry(cell))->getGeometry(cell);
1078         double eta = cellGeometry->getPosition().eta();
1079         double phi = cellGeometry->getPosition().phi();
1080         double zc = cellGeometry->getPosition().z();
1081         int sub = cell.subdet();
1082         int depth = cell.depth();
1083         int inteta = cell.ieta();
1084         int intphi = cell.iphi();
1085         double t = j->time();
1086         double en = j->energy();
1087         double enM0 = 0.;
1088         double enM3 = 0.;
1089         double chi2 = 0.;
1090         int stwd = j->flags();
1091         int auxstwd = j->aux();
1092 
1093         int severityLevel = hcalSevLvl((CaloRecHit *)&*j);
1094         if (cell.subdet() == HcalOuter) {
1095           hcalHOSevLvlVec.push_back(severityLevel);
1096         }
1097 
1098         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
1099           csub.push_back(sub);
1100           cen.push_back(en);
1101           cenM0.push_back(enM0);
1102           cenM3.push_back(enM3);
1103           cchi2.push_back(chi2);
1104           ceta.push_back(eta);
1105           cphi.push_back(phi);
1106           ctime.push_back(t);
1107           cieta.push_back(inteta);
1108           ciphi.push_back(intphi);
1109           cdepth.push_back(depth);
1110           cz.push_back(zc);
1111           cstwd.push_back(stwd);
1112           cauxstwd.push_back(auxstwd);
1113           csevlev.push_back(severityLevel);
1114         }
1115       }
1116     }
1117   }
1118 }
1119 
1120 double HcalRecHitsAnalyzer::dR(double eta1, double phi1, double eta2, double phi2) {
1121   double PI = 3.1415926535898;
1122   double deltaphi = phi1 - phi2;
1123   if (phi2 > phi1) {
1124     deltaphi = phi2 - phi1;
1125   }
1126   if (deltaphi > PI) {
1127     deltaphi = 2. * PI - deltaphi;
1128   }
1129   double deltaeta = eta2 - eta1;
1130   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
1131   return tmp;
1132 }
1133 
1134 double HcalRecHitsAnalyzer::phi12(double phi1, double en1, double phi2, double en2) {
1135   // weighted mean value of phi1 and phi2
1136 
1137   double tmp;
1138   double PI = 3.1415926535898;
1139   double a1 = phi1;
1140   double a2 = phi2;
1141 
1142   if (a1 > 0.5 * PI && a2 < 0.)
1143     a2 += 2 * PI;
1144   if (a2 > 0.5 * PI && a1 < 0.)
1145     a1 += 2 * PI;
1146   tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
1147   if (tmp > PI)
1148     tmp -= 2. * PI;
1149 
1150   return tmp;
1151 }
1152 
1153 double HcalRecHitsAnalyzer::dPhiWsign(double phi1, double phi2) {
1154   // clockwise      phi2 w.r.t phi1 means "+" phi distance
1155   // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
1156 
1157   double PI = 3.1415926535898;
1158   double a1 = phi1;
1159   double a2 = phi2;
1160   double tmp = a2 - a1;
1161   if (a1 * a2 < 0.) {
1162     if (a1 > 0.5 * PI)
1163       tmp += 2. * PI;
1164     if (a2 > 0.5 * PI)
1165       tmp -= 2. * PI;
1166   }
1167   return tmp;
1168 }
1169 
1170 int HcalRecHitsAnalyzer::hcalSevLvl(const CaloRecHit *hit) {
1171   HcalDetId id = hit->detid();
1172   if (theHcalTopology->getMergePositionFlag() && id.subdet() == HcalEndcap) {
1173     id = theHcalTopology->idFront(id);
1174   }
1175 
1176   const uint32_t recHitFlag = hit->flags();
1177   const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1178 
1179   int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1180 
1181   return severityLevel;
1182 }
1183 
1184 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!