Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:25 

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