Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:26 

0001 // -*- C++ -*-
0002 //
0003 // Package:    HcalDigisValidation
0004 // Class:      HcalDigisValidation
0005 //
0006 /**\class HcalDigisValidation HcalDigisValidation.cc Validation/HcalDigis/src/HcalDigisValidation.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012  */
0013 //
0014 // Original Author:  Ali Fahim,22 R-013,+41227672649,
0015 //         Created:  Wed Mar 23 11:42:34 CET 2011
0016 //
0017 //
0018 
0019 #include "Validation/HcalDigis/interface/HcalDigisValidation.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0022 
0023 HcalDigisValidation::HcalDigisValidation(const edm::ParameterSet& iConfig) {
0024   subdet_ = iConfig.getUntrackedParameter<std::string>("subdetector", "all");
0025   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "");
0026   //    inputLabel_ = iConfig.getParameter<std::string > ("digiLabel");
0027   inputTag_ = iConfig.getParameter<edm::InputTag>("digiTag");
0028   QIE10inputTag_ = iConfig.getParameter<edm::InputTag>("QIE10digiTag");
0029   QIE11inputTag_ = iConfig.getParameter<edm::InputTag>("QIE11digiTag");
0030   emulTPsTag_ = iConfig.getParameter<edm::InputTag>("emulTPs");
0031   dataTPsTag_ = iConfig.getParameter<edm::InputTag>("dataTPs");
0032   mc_ = iConfig.getUntrackedParameter<std::string>("mc", "no");
0033   mode_ = iConfig.getUntrackedParameter<std::string>("mode", "multi");
0034   dirName_ = iConfig.getUntrackedParameter<std::string>("dirName", "HcalDigisV/HcalDigiTask");
0035   testNumber_ = iConfig.getParameter<bool>("TestNumber");
0036   hep17_ = iConfig.getParameter<bool>("hep17");
0037   HEPhase1_ = iConfig.getParameter<bool>("HEPhase1");
0038   HBPhase1_ = iConfig.getParameter<bool>("HBPhase1");
0039   Plot_TP_ver_ = iConfig.getParameter<bool>("Plot_TP_ver");
0040 
0041   // register for data access
0042   if (iConfig.exists("simHits")) {
0043     tok_mc_ = consumes<edm::PCaloHitContainer>(iConfig.getUntrackedParameter<edm::InputTag>("simHits"));
0044   }
0045   tok_hbhe_ = consumes<HBHEDigiCollection>(inputTag_);
0046   tok_ho_ = consumes<HODigiCollection>(inputTag_);
0047   tok_hf_ = consumes<HFDigiCollection>(inputTag_);
0048   tok_emulTPs_ = consumes<HcalTrigPrimDigiCollection>(emulTPsTag_);
0049   if (dataTPsTag_ == edm::InputTag(""))
0050     skipDataTPs = true;
0051   else {
0052     skipDataTPs = false;
0053     tok_dataTPs_ = consumes<HcalTrigPrimDigiCollection>(dataTPsTag_);
0054   }
0055 
0056   tok_qie10_hf_ = consumes<QIE10DigiCollection>(QIE10inputTag_);
0057   tok_qie11_hbhe_ = consumes<QIE11DigiCollection>(QIE11inputTag_);
0058 
0059   tok_HRNDC_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>();
0060   tok_Geom_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
0061   tok_Decoder_ = esConsumes<CaloTPGTranscoder, CaloTPGRecord>();
0062   tok_TPGeom_ = esConsumes<HcalTrigTowerGeometry, CaloGeometryRecord>();
0063   tok_Topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0064   tok_Cond_ = esConsumes<HcalDbService, HcalDbRecord>();
0065 
0066   nevent1 = 0;
0067   nevent2 = 0;
0068   nevent3 = 0;
0069   nevent4 = 0;
0070   nevtot = 0;
0071 
0072   msm_ = new std::map<std::string, MonitorElement*>();
0073 
0074   if (!outputFile_.empty())
0075     edm::LogVerbatim("OutputInfo") << " Hcal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
0076   else
0077     edm::LogVerbatim("OutputInfo") << " Hcal Digi Task histograms will NOT be saved";
0078 }
0079 
0080 HcalDigisValidation::~HcalDigisValidation() { delete msm_; }
0081 
0082 void HcalDigisValidation::dqmBeginRun(const edm::Run& run, const edm::EventSetup& es) {
0083   hcons_ = &es.getData(tok_HRNDC_);
0084 
0085   maxDepth_[1] = hcons_->getMaxDepth(0);  // HB
0086   maxDepth_[2] = hcons_->getMaxDepth(1);  // HE
0087   maxDepth_[3] = hcons_->getMaxDepth(3);  // HO
0088   maxDepth_[4] = hcons_->getMaxDepth(2);  // HF
0089   maxDepth_[0] = (maxDepth_[1] > maxDepth_[2] ? maxDepth_[1] : maxDepth_[2]);
0090   maxDepth_[0] = (maxDepth_[0] > maxDepth_[3] ? maxDepth_[0] : maxDepth_[3]);
0091   maxDepth_[0] = (maxDepth_[0] > maxDepth_[4] ? maxDepth_[0] : maxDepth_[4]);  // any of HB/HE/HO/HF
0092 
0093   const CaloGeometry* geo = &es.getData(tok_Geom_);
0094   const HcalGeometry* gHB = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0095   const HcalGeometry* gHE = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0096   const HcalGeometry* gHO = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalOuter));
0097   const HcalGeometry* gHF = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalForward));
0098 
0099   nChannels_[1] = gHB->getHxSize(1);
0100   nChannels_[2] = gHE->getHxSize(2);
0101   nChannels_[3] = gHO->getHxSize(3);
0102   nChannels_[4] = gHF->getHxSize(4);
0103 
0104   nChannels_[0] = nChannels_[1] + nChannels_[2] + nChannels_[3] + nChannels_[4];
0105 }
0106 
0107 void HcalDigisValidation::bookHistograms(DQMStore::IBooker& ib, edm::Run const& run, edm::EventSetup const& es) {
0108   ib.setCurrentFolder(dirName_);
0109 
0110   // book
0111   book1D(ib, "nevtot", 1, 0, 1);
0112   int bnoise = 0;
0113   int bmc = 0;
0114   if (subdet_ == "noise")
0115     bnoise = 1;
0116   if (mc_ == "yes")
0117     bmc = 1;
0118   if (subdet_ == "noise" || subdet_ == "all") {
0119     booking(ib, "HB", bnoise, bmc);
0120     booking(ib, "HO", bnoise, bmc);
0121     booking(ib, "HF", bnoise, bmc);
0122     booking(ib, "HE", bnoise, bmc);
0123   } else {
0124     booking(ib, subdet_, 0, bmc);
0125   }
0126 
0127   if (skipDataTPs)
0128     return;
0129 
0130   HistLim tp_hl_et(260, -10, 250);
0131   HistLim tp_hl_ntp(640, -20, 3180);
0132   HistLim tp_hl_ntp_sub(404, -20, 2000);
0133   HistLim tp_hl_ieta(85, -42.5, 42.5);
0134   HistLim tp_hl_iphi(74, -0.5, 73.5);
0135 
0136   book1D(ib, "HcalDigiTask_tp_et", tp_hl_et);
0137   book1D(ib, "HcalDigiTask_tp_et_HB", tp_hl_et);
0138   book1D(ib, "HcalDigiTask_tp_et_HE", tp_hl_et);
0139   book1D(ib, "HcalDigiTask_tp_et_HF", tp_hl_et);
0140   book1D(ib, "HcalDigiTask_tp_ntp", tp_hl_ntp);
0141   book1D(ib, "HcalDigiTask_tp_ntp_HB", tp_hl_ntp_sub);
0142   book1D(ib, "HcalDigiTask_tp_ntp_HE", tp_hl_ntp_sub);
0143   book1D(ib, "HcalDigiTask_tp_ntp_HF", tp_hl_ntp_sub);
0144   book1D(ib, "HcalDigiTask_tp_ntp_ieta", tp_hl_ieta);
0145   book1D(ib, "HcalDigiTask_tp_ntp_iphi", tp_hl_iphi);
0146   book1D(ib, "HcalDigiTask_tp_ntp_10_ieta", tp_hl_ieta);
0147   book2D(ib, "HcalDigiTask_tp_et_ieta", tp_hl_ieta, tp_hl_et);
0148   book2D(ib, "HcalDigiTask_tp_ieta_iphi", tp_hl_ieta, tp_hl_iphi);
0149   bookPf(ib, "HcalDigiTask_tp_ave_et_ieta", tp_hl_ieta, tp_hl_et, " ");
0150   if (Plot_TP_ver_) {
0151     book1D(ib, "HcalDigiTask_tp_et_v0", tp_hl_et);
0152     book1D(ib, "HcalDigiTask_tp_et_v1", tp_hl_et);
0153     book1D(ib, "HcalDigiTask_tp_et_HF_v0", tp_hl_et);
0154     book1D(ib, "HcalDigiTask_tp_et_HF_v1", tp_hl_et);
0155     book1D(ib, "HcalDigiTask_tp_ntp_v0", tp_hl_ntp);
0156     book1D(ib, "HcalDigiTask_tp_ntp_v1", tp_hl_ntp);
0157     book1D(ib, "HcalDigiTask_tp_ntp_HF_v0", tp_hl_ntp_sub);
0158     book1D(ib, "HcalDigiTask_tp_ntp_HF_v1", tp_hl_ntp_sub);
0159     book1D(ib, "HcalDigiTask_tp_ntp_ieta_v0", tp_hl_ieta);
0160     book1D(ib, "HcalDigiTask_tp_ntp_ieta_v1", tp_hl_ieta);
0161     book1D(ib, "HcalDigiTask_tp_ntp_iphi_v0", tp_hl_iphi);
0162     book1D(ib, "HcalDigiTask_tp_ntp_iphi_v1", tp_hl_iphi);
0163     book1D(ib, "HcalDigiTask_tp_ntp_10_ieta_v0", tp_hl_ieta);
0164     book1D(ib, "HcalDigiTask_tp_ntp_10_ieta_v1", tp_hl_ieta);
0165     book2D(ib, "HcalDigiTask_tp_et_ieta_v0", tp_hl_ieta, tp_hl_et);
0166     book2D(ib, "HcalDigiTask_tp_et_ieta_v1", tp_hl_ieta, tp_hl_et);
0167     book2D(ib, "HcalDigiTask_tp_ieta_iphi_v0", tp_hl_ieta, tp_hl_iphi);
0168     book2D(ib, "HcalDigiTask_tp_ieta_iphi_v1", tp_hl_ieta, tp_hl_iphi);
0169     bookPf(ib, "HcalDigiTask_tp_ave_et_ieta_v0", tp_hl_ieta, tp_hl_et, " ");
0170     bookPf(ib, "HcalDigiTask_tp_ave_et_ieta_v1", tp_hl_ieta, tp_hl_et, " ");
0171   }
0172 }
0173 
0174 void HcalDigisValidation::booking(DQMStore::IBooker& ib, const std::string bsubdet, int bnoise, int bmc) {
0175   // Adjust/Optimize binning (JR Dittmann, 16-JUL-2015)
0176 
0177   HistLim Ndigis(2600, 0., 2600.);
0178   HistLim sime(200, 0., 1.0);
0179 
0180   HistLim digiAmp(360, -100., 7100.);
0181   HistLim digiAmpWide(2410, -3000., 720000.);  //300 fC binning
0182   HistLim ratio(2000, -100., 3900.);
0183   HistLim sumAmp(100, -500., 1500.);
0184 
0185   HistLim nbin(11, -0.5, 10.5);
0186 
0187   HistLim pedestal(80, -1.0, 15.);
0188   HistLim pedestalfC(400, -10., 30.);
0189 
0190   HistLim frac(80, -0.20, 1.40);
0191 
0192   HistLim pedLim(80, 0., 8.);
0193   HistLim pedWidthLim(100, 0., 2.);
0194 
0195   HistLim gainLim(120, 0., 0.6);
0196   HistLim gainWidthLim(160, 0., 0.32);
0197 
0198   HistLim ietaLim(85, -42.5, 42.5);
0199   HistLim iphiLim(74, -0.5, 73.5);
0200 
0201   HistLim depthLim(15, -0.5, 14.5);
0202 
0203   //...TDC
0204   HistLim tdcLim(250, 0., 250.);
0205   HistLim adcLim(256, 0., 256.);
0206 
0207   if (bsubdet == "HB") {
0208     Ndigis = HistLim(((int)(nChannels_[1] / 100) + 1) * 100, 0., (float)((int)(nChannels_[1] / 100) + 1) * 100);
0209   } else if (bsubdet == "HE") {
0210     sime = HistLim(200, 0., 1.0);
0211     Ndigis = HistLim(((int)(nChannels_[2] / 100) + 1) * 100, 0., (float)((int)(nChannels_[2] / 100) + 1) * 100);
0212   } else if (bsubdet == "HF") {
0213     sime = HistLim(100, 0., 100.);
0214     pedLim = HistLim(100, 0., 20.);
0215     pedWidthLim = HistLim(100, 0., 5.);
0216     frac = HistLim(400, -4.00, 4.00);
0217     Ndigis = HistLim(((int)(nChannels_[4] / 100) + 1) * 100, 0., (float)((int)(nChannels_[4] / 100) + 1) * 100);
0218   } else if (bsubdet == "HO") {
0219     sime = HistLim(200, 0., 1.0);
0220     gainLim = HistLim(160, 0., 1.6);
0221     Ndigis = HistLim(((int)(nChannels_[3] / 100) + 1) * 100, 0., (float)((int)(nChannels_[3] / 100) + 1) * 100);
0222   }
0223 
0224   int isubdet = 0;
0225   if (bsubdet == "HB")
0226     isubdet = 1;
0227   else if (bsubdet == "HE")
0228     isubdet = 2;
0229   else if (bsubdet == "HO")
0230     isubdet = 3;
0231   else if (bsubdet == "HF")
0232     isubdet = 4;
0233   else
0234     edm::LogWarning("HcalDigisValidation") << "HcalDigisValidation Warning: not HB/HE/HF/HO " << bsubdet << std::endl;
0235 
0236   Char_t histo[100];
0237   const char* sub = bsubdet.c_str();
0238   if (bnoise == 0) {
0239     // number of digis in each subdetector
0240     sprintf(histo, "HcalDigiTask_Ndigis_%s", sub);
0241     book1D(ib, histo, Ndigis);
0242 
0243     // maps of occupancies
0244     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0245       sprintf(histo, "HcalDigiTask_ieta_iphi_occupancy_map_depth%d_%s", depth, sub);
0246       book2D(ib, histo, ietaLim, iphiLim);
0247     }
0248 
0249     //Depths
0250     sprintf(histo, "HcalDigiTask_depths_%s", sub);
0251     book1D(ib, histo, depthLim);
0252 
0253     // occupancies vs ieta
0254     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0255       sprintf(histo, "HcalDigiTask_occupancy_vs_ieta_depth%d_%s", depth, sub);
0256       book1D(ib, histo, ietaLim);
0257     }
0258 
0259     // just 1D of all cells' amplitudes
0260     sprintf(histo, "HcalDigiTask_sum_all_amplitudes_%s", sub);
0261     if ((HBPhase1_ && bsubdet == "HB") || (HEPhase1_ && bsubdet == "HE"))
0262       book1D(ib, histo, digiAmpWide);
0263     else
0264       book1D(ib, histo, digiAmp);
0265 
0266     sprintf(histo, "HcalDigiTask_number_of_amplitudes_above_10fC_%s", sub);
0267     book1D(ib, histo, Ndigis);
0268 
0269     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0270       sprintf(histo, "HcalDigiTask_ADC0_adc_depth%d_%s", depth, sub);
0271       book1D(ib, histo, pedestal);
0272     }
0273 
0274     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0275       sprintf(histo, "HcalDigiTask_ADC0_fC_depth%d_%s", depth, sub);
0276       book1D(ib, histo, pedestalfC);
0277     }
0278 
0279     sprintf(histo, "HcalDigiTask_signal_amplitude_%s", sub);
0280     if ((HBPhase1_ && bsubdet == "HB") || (HEPhase1_ && bsubdet == "HE"))
0281       book1D(ib, histo, digiAmpWide);
0282     else
0283       book1D(ib, histo, digiAmp);
0284 
0285     if (hep17_ && bsubdet == "HE") {
0286       sprintf(histo, "HcalDigiTask_signal_amplitude_HEP17");
0287       book1D(ib, histo, digiAmpWide);
0288     }
0289     //
0290     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0291       sprintf(histo, "HcalDigiTask_signal_amplitude_depth%d_%s", depth, sub);
0292       if ((HBPhase1_ && bsubdet == "HB") || (HEPhase1_ && bsubdet == "HE"))
0293         book1D(ib, histo, digiAmpWide);
0294       else
0295         book1D(ib, histo, digiAmp);
0296       if (hep17_ && bsubdet == "HE") {
0297         sprintf(histo, "HcalDigiTask_signal_amplitude_depth%d_HEP17", depth);
0298         book1D(ib, histo, digiAmpWide);
0299       }
0300     }
0301 
0302     sprintf(histo, "HcalDigiTask_signal_amplitude_vs_bin_all_depths_%s", sub);
0303     if ((HBPhase1_ && bsubdet == "HB") || (HEPhase1_ && bsubdet == "HE"))
0304       book2D(ib, histo, nbin, digiAmpWide);
0305     else
0306       book2D(ib, histo, nbin, digiAmp);
0307     if (hep17_ && bsubdet == "HE") {
0308       sprintf(histo, "HcalDigiTask_signal_amplitude_vs_bin_all_depths_HEP17");
0309       book2D(ib, histo, nbin, digiAmpWide);
0310     }
0311 
0312     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0313       sprintf(histo, "HcalDigiTask_all_amplitudes_vs_bin_1D_depth%d_%s", depth, sub);
0314       book1D(ib, histo, nbin);
0315       if (hep17_ && bsubdet == "HE") {
0316         sprintf(histo, "HcalDigiTask_all_amplitudes_vs_bin_1D_depth%d_HEP17", depth);
0317         book1D(ib, histo, nbin);
0318       }
0319     }
0320 
0321     sprintf(histo, "HcalDigiTask_SOI_frac_%s", sub);
0322     book1D(ib, histo, frac);
0323     sprintf(histo, "HcalDigiTask_postSOI_frac_%s", sub);
0324     book1D(ib, histo, frac);
0325 
0326     if (bmc == 1) {
0327       sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_%s", sub);
0328       book2D(ib, histo, sime, digiAmp);
0329       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0330         sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_depth%d_%s", depth, sub);
0331         book2D(ib, histo, sime, digiAmp);
0332       }
0333 
0334       sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_profile_%s", sub);
0335       bookPf(ib, histo, sime, digiAmp);
0336       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0337         sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_profile_depth%d_%s", depth, sub);
0338         bookPf(ib, histo, sime, digiAmp);
0339       }
0340 
0341       sprintf(histo, "HcalDigiTask_ratio_amplitude_vs_simhits_%s", sub);
0342       book1D(ib, histo, ratio);
0343       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0344         sprintf(histo, "HcalDigiTask_ratio_amplitude_vs_simhits_depth%d_%s", depth, sub);
0345         book1D(ib, histo, ratio);
0346       }
0347 
0348       //...TDC
0349       if (bsubdet == "HB" || bsubdet == "HE") {
0350         sprintf(histo, "HcalDigiTask_TDCtime_%s", sub);
0351         book1D(ib, histo, tdcLim);
0352 
0353         sprintf(histo, "HcalDigiTask_TDCtime_vs_ADC_%s", sub);
0354         book2D(ib, histo, adcLim, tdcLim);
0355       }
0356 
0357     }  //mc only
0358 
0359   } else {  // noise only
0360 
0361     // EVENT "1" distributions of all cells properties
0362 
0363     //KH
0364     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0365       sprintf(histo, "HcalDigiTask_gain_capId0_Depth%d_%s", depth, sub);
0366       book1D(ib, histo, gainLim);
0367       sprintf(histo, "HcalDigiTask_gain_capId1_Depth%d_%s", depth, sub);
0368       book1D(ib, histo, gainLim);
0369       sprintf(histo, "HcalDigiTask_gain_capId2_Depth%d_%s", depth, sub);
0370       book1D(ib, histo, gainLim);
0371       sprintf(histo, "HcalDigiTask_gain_capId3_Depth%d_%s", depth, sub);
0372       book1D(ib, histo, gainLim);
0373 
0374       sprintf(histo, "HcalDigiTask_gainWidth_capId0_Depth%d_%s", depth, sub);
0375       book1D(ib, histo, gainWidthLim);
0376       sprintf(histo, "HcalDigiTask_gainWidth_capId1_Depth%d_%s", depth, sub);
0377       book1D(ib, histo, gainWidthLim);
0378       sprintf(histo, "HcalDigiTask_gainWidth_capId2_Depth%d_%s", depth, sub);
0379       book1D(ib, histo, gainWidthLim);
0380       sprintf(histo, "HcalDigiTask_gainWidth_capId3_Depth%d_%s", depth, sub);
0381       book1D(ib, histo, gainWidthLim);
0382 
0383       sprintf(histo, "HcalDigiTask_pedestal_capId0_Depth%d_%s", depth, sub);
0384       book1D(ib, histo, pedLim);
0385       sprintf(histo, "HcalDigiTask_pedestal_capId1_Depth%d_%s", depth, sub);
0386       book1D(ib, histo, pedLim);
0387       sprintf(histo, "HcalDigiTask_pedestal_capId2_Depth%d_%s", depth, sub);
0388       book1D(ib, histo, pedLim);
0389       sprintf(histo, "HcalDigiTask_pedestal_capId3_Depth%d_%s", depth, sub);
0390       book1D(ib, histo, pedLim);
0391 
0392       sprintf(histo, "HcalDigiTask_pedestal_width_capId0_Depth%d_%s", depth, sub);
0393       book1D(ib, histo, pedWidthLim);
0394       sprintf(histo, "HcalDigiTask_pedestal_width_capId1_Depth%d_%s", depth, sub);
0395       book1D(ib, histo, pedWidthLim);
0396       sprintf(histo, "HcalDigiTask_pedestal_width_capId2_Depth%d_%s", depth, sub);
0397       book1D(ib, histo, pedWidthLim);
0398       sprintf(histo, "HcalDigiTask_pedestal_width_capId3_Depth%d_%s", depth, sub);
0399       book1D(ib, histo, pedWidthLim);
0400     }
0401 
0402     //KH
0403     for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0404       sprintf(histo, "HcalDigiTask_gainMap_Depth%d_%s", depth, sub);
0405       book2D(ib, histo, ietaLim, iphiLim);
0406       sprintf(histo, "HcalDigiTask_pwidthMap_Depth%d_%s", depth, sub);
0407       book2D(ib, histo, ietaLim, iphiLim);
0408     }
0409 
0410   }  //end of noise-only
0411 }  //book
0412 
0413 void HcalDigisValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0414   conditions_ = &iSetup.getData(tok_Cond_);
0415 
0416   //TP Code
0417   const auto& decoder = &iSetup.getData(tok_Decoder_);
0418   const auto& tp_geometry = &iSetup.getData(tok_TPGeom_);
0419   htopo_ = &iSetup.getData(tok_Topo_);
0420 
0421   //Get all handles
0422   edm::Handle<HcalTrigPrimDigiCollection> emulTPs;
0423   iEvent.getByToken(tok_emulTPs_, emulTPs);
0424 
0425   edm::Handle<HcalTrigPrimDigiCollection> dataTPs;
0426   if (!skipDataTPs)
0427     iEvent.getByToken(tok_dataTPs_, dataTPs);
0428   //iEvent.getByLabel("hcalDigis", dataTPs);
0429 
0430   //~TP Code
0431 
0432   if (subdet_ != "all") {
0433     noise_ = 0;
0434     if (subdet_ == "HB") {
0435       reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0436       reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0437     }
0438     if (subdet_ == "HE") {
0439       reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0440       reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0441     }
0442     if (subdet_ == "HO")
0443       reco<HODataFrame>(iEvent, iSetup, tok_ho_);
0444     if (subdet_ == "HF") {
0445       reco<HFDataFrame>(iEvent, iSetup, tok_hf_);
0446       reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
0447     }
0448 
0449     if (subdet_ == "noise") {
0450       noise_ = 1;
0451       subdet_ = "HB";
0452       reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0453       reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0454       subdet_ = "HE";
0455       reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0456       reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0457       subdet_ = "HO";
0458       reco<HODataFrame>(iEvent, iSetup, tok_ho_);
0459       subdet_ = "HF";
0460       reco<HFDataFrame>(iEvent, iSetup, tok_hf_);
0461       reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
0462       subdet_ = "noise";
0463     }
0464   }  // all subdetectors
0465   else {
0466     noise_ = 0;
0467 
0468     subdet_ = "HB";
0469     reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0470     reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0471     subdet_ = "HE";
0472     reco<HBHEDataFrame>(iEvent, iSetup, tok_hbhe_);
0473     reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
0474     subdet_ = "HO";
0475     reco<HODataFrame>(iEvent, iSetup, tok_ho_);
0476     subdet_ = "HF";
0477     reco<HFDataFrame>(iEvent, iSetup, tok_hf_);
0478     reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
0479     subdet_ = "all";
0480   }
0481 
0482   fill1D("nevtot", 0);
0483   nevtot++;
0484 
0485   //TP Code
0486   //Counters
0487   int c = 0, chb = 0, che = 0, chf = 0, cv0 = 0, cv1 = 0, chfv0 = 0, chfv1 = 0;
0488 
0489   if (skipDataTPs)
0490     return;
0491 
0492   for (HcalTrigPrimDigiCollection::const_iterator itr = dataTPs->begin(); itr != dataTPs->end(); ++itr) {
0493     int ieta = itr->id().ieta();
0494     int iphi = itr->id().iphi();
0495 
0496     HcalSubdetector subdet = (HcalSubdetector)0;
0497 
0498     if (abs(ieta) <= 16)
0499       subdet = HcalSubdetector::HcalBarrel;
0500     else if (abs(ieta) < tp_geometry->firstHFTower(itr->id().version()))
0501       subdet = HcalSubdetector::HcalEndcap;
0502     else if (abs(ieta) <= 42)
0503       subdet = HcalSubdetector::HcalForward;
0504 
0505     //Right now, the only case where version matters is in HF
0506     //If the subdetector is not HF, set version to -1
0507     int tpVersion = (subdet == HcalSubdetector::HcalForward ? itr->id().version() : -1);
0508 
0509     float en = decoder->hcaletValue(itr->id(), itr->t0());
0510 
0511     if (en < 0.00001)
0512       continue;
0513 
0514     //Plot the variables
0515     //Retain classic behavior (include all tps)
0516     //Additional plots that only include HF 3x2 or HF 1x1
0517 
0518     //Classics
0519     fill1D("HcalDigiTask_tp_et", en);
0520     fill2D("HcalDigiTask_tp_et_ieta", ieta, en);
0521     fill2D("HcalDigiTask_tp_ieta_iphi", ieta, iphi);
0522     fillPf("HcalDigiTask_tp_ave_et_ieta", ieta, en);
0523     fill1D("HcalDigiTask_tp_ntp_ieta", ieta);
0524     fill1D("HcalDigiTask_tp_ntp_iphi", iphi);
0525     if (en > 10.)
0526       fill1D("HcalDigiTask_tp_ntp_10_ieta", ieta);
0527 
0528     //3x2 Trig Primitives (tpVersion == 0)
0529     if ((subdet != HcalSubdetector::HcalForward || tpVersion == 0) && Plot_TP_ver_) {
0530       fill1D("HcalDigiTask_tp_et_v0", en);
0531       fill2D("HcalDigiTask_tp_et_ieta_v0", ieta, en);
0532       fill2D("HcalDigiTask_tp_ieta_iphi_v0", ieta, iphi);
0533       fillPf("HcalDigiTask_tp_ave_et_ieta_v0", ieta, en);
0534       fill1D("HcalDigiTask_tp_ntp_ieta_v0", ieta);
0535       fill1D("HcalDigiTask_tp_ntp_iphi_v0", iphi);
0536       if (en > 10.)
0537         fill1D("HcalDigiTask_tp_ntp_10_ieta_v0", ieta);
0538     }
0539 
0540     //1x1 Trig Primitives (tpVersion == 1)
0541     if ((subdet != HcalSubdetector::HcalForward || tpVersion == 1) && Plot_TP_ver_) {
0542       fill1D("HcalDigiTask_tp_et_v1", en);
0543       fill2D("HcalDigiTask_tp_et_ieta_v1", ieta, en);
0544       fill2D("HcalDigiTask_tp_ieta_iphi_v1", ieta, iphi);
0545       fillPf("HcalDigiTask_tp_ave_et_ieta_v1", ieta, en);
0546       fill1D("HcalDigiTask_tp_ntp_ieta_v1", ieta);
0547       fill1D("HcalDigiTask_tp_ntp_iphi_v1", iphi);
0548       if (en > 10.)
0549         fill1D("HcalDigiTask_tp_ntp_10_ieta_v1", ieta);
0550     }
0551 
0552     ++c;
0553     if (subdet == HcalSubdetector::HcalBarrel) {
0554       fill1D("HcalDigiTask_tp_et_HB", en);
0555       ++chb;
0556       if (Plot_TP_ver_) {
0557         ++cv0;
0558         ++cv1;
0559       }
0560     }
0561     if (subdet == HcalSubdetector::HcalEndcap) {
0562       fill1D("HcalDigiTask_tp_et_HE", en);
0563       ++che;
0564       if (Plot_TP_ver_) {
0565         ++cv0;
0566         ++cv1;
0567       }
0568     }
0569     if (subdet == HcalSubdetector::HcalForward) {
0570       fill1D("HcalDigiTask_tp_et_HF", en);
0571       ++chf;
0572 
0573       if (tpVersion == 0 && Plot_TP_ver_) {
0574         fill1D("HcalDigiTask_tp_et_HF_v0", en);
0575         ++chfv0;
0576         ++cv0;
0577       }
0578 
0579       if (tpVersion == 1 && Plot_TP_ver_) {
0580         fill1D("HcalDigiTask_tp_et_HF_v1", en);
0581         ++chfv1;
0582         ++cv1;
0583       }
0584     }
0585 
0586   }  //end data TP collection
0587 
0588   fill1D("HcalDigiTask_tp_ntp", c);
0589   fill1D("HcalDigiTask_tp_ntp_HB", chb);
0590   fill1D("HcalDigiTask_tp_ntp_HE", che);
0591   fill1D("HcalDigiTask_tp_ntp_HF", chf);
0592   if (Plot_TP_ver_) {
0593     fill1D("HcalDigiTask_tp_ntp_v0", cv0);
0594     fill1D("HcalDigiTask_tp_ntp_v1", cv1);
0595     fill1D("HcalDigiTask_tp_ntp_HF_v0", chfv0);
0596     fill1D("HcalDigiTask_tp_ntp_HF_v1", chfv1);
0597   }
0598 
0599   //~TP Code
0600 }
0601 
0602 template <class Digi>
0603 void HcalDigisValidation::reco(const edm::Event& iEvent,
0604                                const edm::EventSetup& iSetup,
0605                                const edm::EDGetTokenT<edm::SortedCollection<Digi> >& tok) {
0606   // HistLim =============================================================
0607 
0608   std::string strtmp;
0609 
0610   // ======================================================================
0611   typename edm::Handle<edm::SortedCollection<Digi> > digiCollection;
0612   typename edm::SortedCollection<Digi>::const_iterator digiItr;
0613 
0614   // ADC2fC
0615   CaloSamples tool;
0616   iEvent.getByToken(tok, digiCollection);
0617   if (!digiCollection.isValid())
0618     return;
0619   int isubdet = 0;
0620   if (subdet_ == "HB")
0621     isubdet = 1;
0622   if (subdet_ == "HE")
0623     isubdet = 2;
0624   if (subdet_ == "HO")
0625     isubdet = 3;
0626   if (subdet_ == "HF")
0627     isubdet = 4;
0628 
0629   if (isubdet == 1)
0630     nevent1++;
0631   if (isubdet == 2)
0632     nevent2++;
0633   if (isubdet == 3)
0634     nevent3++;
0635   if (isubdet == 4)
0636     nevent4++;
0637 
0638   int indigis = 0;
0639   //  amplitude for signal cell at diff. depths
0640   std::vector<double> v_ampl_c(maxDepth_[isubdet] + 1, 0);
0641 
0642   // is set to 1 if "seed" SimHit is found
0643   int seedSimHit = 0;
0644 
0645   int ieta_Sim = 9999;
0646   int iphi_Sim = 9999;
0647   double emax_Sim = -9999.;
0648 
0649   // SimHits MC only
0650   if (mc_ == "yes") {
0651     edm::Handle<edm::PCaloHitContainer> hcalHits;
0652     iEvent.getByToken(tok_mc_, hcalHits);
0653     const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0654 
0655     if (isubdet != 0 && noise_ == 0) {  // signal only SimHits
0656 
0657       for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0658            ++simhits) {
0659         unsigned int id = simhits->id();
0660         int sub, ieta, iphi;
0661         HcalDetId hid;
0662         if (testNumber_)
0663           hid = HcalHitRelabeller::relabel(id, hcons_);
0664         else
0665           hid = HcalDetId(id);
0666         sub = hid.subdet();
0667         ieta = hid.ieta();
0668         iphi = hid.iphi();
0669 
0670         double en = simhits->energy();
0671 
0672         if (en > emax_Sim && sub == isubdet) {
0673           emax_Sim = en;
0674           ieta_Sim = ieta;
0675           iphi_Sim = iphi;
0676           // to limit "seed" SimHit energy in case of "multi" event
0677           if (mode_ == "multi" && ((sub == 4 && en < 100. && en > 1.) || ((sub != 4) && en < 1. && en > 0.02))) {
0678             seedSimHit = 1;
0679             break;
0680           }
0681         }
0682 
0683       }  // end of SimHits cycle
0684 
0685       // found highest-energy SimHit for single-particle
0686       if (mode_ != "multi" && emax_Sim > 0.)
0687         seedSimHit = 1;
0688     }  // end of SimHits
0689   }    // end of mc_ == "yes"
0690 
0691   // CYCLE OVER CELLS ========================================================
0692   int Ndig = 0;
0693 
0694   for (digiItr = digiCollection->begin(); digiItr != digiCollection->end(); digiItr++) {
0695     HcalDetId cell(digiItr->id());
0696     int depth = cell.depth();
0697     int iphi = cell.iphi();
0698     int ieta = cell.ieta();
0699     int sub = cell.subdet();
0700 
0701     if (depth > maxDepth_[isubdet] && sub == isubdet) {
0702       edm::LogWarning("HcalDetId") << "HcalDetID presents conflicting information. Depth: " << depth
0703                                    << ", iphi: " << iphi << ", ieta: " << ieta
0704                                    << ". Max depth from geometry is: " << maxDepth_[isubdet]
0705                                    << ". TestNumber = " << testNumber_;
0706       continue;
0707     }
0708 
0709     //  amplitude for signal cell at diff. depths
0710     std::vector<double> v_ampl(maxDepth_[isubdet] + 1, 0);
0711 
0712     // Gains, pedestals (once !) and only for "noise" case
0713     if (((nevent1 == 1 && isubdet == 1) || (nevent2 == 1 && isubdet == 2) || (nevent3 == 1 && isubdet == 3) ||
0714          (nevent4 == 1 && isubdet == 4)) &&
0715         noise_ == 1 && sub == isubdet) {
0716       HcalGenericDetId hcalGenDetId(digiItr->id());
0717       const HcalPedestal* pedestal = conditions_->getPedestal(hcalGenDetId);
0718       const HcalGain* gain = conditions_->getGain(hcalGenDetId);
0719       const HcalGainWidth* gainWidth = conditions_->getGainWidth(hcalGenDetId);
0720       const HcalPedestalWidth* pedWidth = conditions_->getPedestalWidth(hcalGenDetId);
0721 
0722       for (int i = 0; i < 4; i++) {
0723         fill1D("HcalDigiTask_gain_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gain->getValue(i));
0724         fill1D("HcalDigiTask_gainWidth_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gainWidth->getValue(i));
0725         fill1D("HcalDigiTask_pedestal_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedestal->getValue(i));
0726         fill1D("HcalDigiTask_pedestal_width_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_,
0727                pedWidth->getWidth(i));
0728       }
0729 
0730       fill2D("HcalDigiTask_gainMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), gain->getValue(0));
0731       fill2D("HcalDigiTask_pwidthMap_Depth" + str(depth) + "_" + subdet_,
0732              double(ieta),
0733              double(iphi),
0734              pedWidth->getWidth(0));
0735 
0736     }  // end of event #1
0737 
0738     if (sub == isubdet)
0739       Ndig++;  // subdet number of digi
0740 
0741     // No-noise case, only single  subdet selected  ===========================
0742 
0743     if (sub == isubdet && noise_ == 0) {
0744       HcalCalibrations calibrations = conditions_->getHcalCalibrations(cell);
0745 
0746       const HcalQIECoder* channelCoder = conditions_->getHcalCoder(cell);
0747       const HcalQIEShape* shape = conditions_->getHcalShape(channelCoder);
0748       HcalCoderDb coder(*channelCoder, *shape);
0749       coder.adc2fC(*digiItr, tool);
0750 
0751       // for dynamic digi time sample analysis
0752       int soi = tool.presamples();
0753       int lastbin = tool.size() - 1;
0754 
0755       double noiseADC = (*digiItr)[0].adc();
0756       double noisefC = tool[0];
0757       // noise evaluations from "pre-samples"
0758       fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
0759       fill1D("HcalDigiTask_ADC0_fC_depth" + str(depth) + "_" + subdet_, noisefC);
0760 
0761       // OCCUPANCY maps fill
0762       fill2D("HcalDigiTask_ieta_iphi_occupancy_map_depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi));
0763 
0764       fill1D("HcalDigiTask_depths_" + subdet_, double(depth));
0765 
0766       // Cycle on time slices
0767       // - for each Digi
0768       // - for one Digi with max SimHits E in subdet
0769 
0770       int closen = 0;  // =1 if 1) seedSimHit = 1 and 2) the cell is the same
0771       if (ieta == ieta_Sim && iphi == iphi_Sim)
0772         closen = seedSimHit;
0773 
0774       for (int ii = 0; ii < tool.size(); ii++) {
0775         int capid = (*digiItr)[ii].capid();
0776         // single ts amplitude
0777         double val = (tool[ii] - calibrations.pedestal(capid));
0778 
0779         if (val > 100.) {
0780           fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
0781           strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
0782           fill1D(strtmp, double(ii), val);
0783         }
0784 
0785         if (closen == 1) {
0786           strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
0787           fill2D(strtmp, double(ii), val);
0788         }
0789 
0790         // all detectors
0791         if (ii >= soi && ii <= lastbin) {
0792           v_ampl[0] += val;
0793           v_ampl[depth] += val;
0794 
0795           if (closen == 1) {
0796             v_ampl_c[0] += val;
0797             v_ampl_c[depth] += val;
0798           }
0799         }
0800       }
0801       // end of time bucket sample
0802 
0803       // maps of sum of amplitudes (sum lin.digis(4,5,6,7) - ped) all depths
0804       // just 1D of all cells' amplitudes
0805       strtmp = "HcalDigiTask_sum_all_amplitudes_" + subdet_;
0806       fill1D(strtmp, v_ampl[0]);
0807 
0808       std::vector<int> v_ampl_sub(v_ampl.begin() + 1, v_ampl.end());  // remove element 0, which is the sum of any depth
0809       double ampl_max = *std::max_element(v_ampl_sub.begin(), v_ampl_sub.end());
0810       if (ampl_max > 10.)
0811         indigis++;
0812       //KH if (ampl1 > 10. || ampl2 > 10. || ampl3 > 10. || ampl4 > 10.) indigis++;
0813 
0814       // fraction 5,6 bins if ampl. is big.
0815       if (v_ampl[depth] > 30.) {
0816         double fbinSOI = tool[soi] - calibrations.pedestal((*digiItr)[soi].capid());
0817         double fbinPS = 0;
0818 
0819         for (int j = soi + 1; j <= lastbin; j++)
0820           fbinPS += tool[j] - calibrations.pedestal((*digiItr)[j].capid());
0821 
0822         fbinSOI /= v_ampl[depth];
0823         fbinPS /= v_ampl[depth];
0824         strtmp = "HcalDigiTask_SOI_frac_" + subdet_;
0825         fill1D(strtmp, fbinSOI);
0826         strtmp = "HcalDigiTask_postSOI_frac_" + subdet_;
0827         fill1D(strtmp, fbinPS);
0828       }
0829 
0830       strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
0831       fill1D(strtmp, v_ampl[0]);
0832       strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
0833       fill1D(strtmp, v_ampl[depth]);
0834     }
0835   }  // End of CYCLE OVER CELLS =============================================
0836 
0837   if (isubdet != 0 && noise_ == 0) {  // signal only, once per event
0838     strtmp = "HcalDigiTask_number_of_amplitudes_above_10fC_" + subdet_;
0839     fill1D(strtmp, indigis);
0840 
0841     // SimHits once again !!!
0842     double eps = 1.e-3;
0843     std::vector<double> v_ehits(maxDepth_[isubdet] + 1, 0);
0844 
0845     if (mc_ == "yes") {
0846       edm::Handle<edm::PCaloHitContainer> hcalHits;
0847       iEvent.getByToken(tok_mc_, hcalHits);
0848       const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0849       for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0850            ++simhits) {
0851         unsigned int id = simhits->id();
0852         int sub, depth, ieta, iphi;
0853         HcalDetId hid;
0854         if (testNumber_)
0855           hid = HcalHitRelabeller::relabel(id, hcons_);
0856         else
0857           hid = HcalDetId(id);
0858         sub = hid.subdet();
0859         depth = hid.depth();
0860         ieta = hid.ieta();
0861         iphi = hid.iphi();
0862 
0863         if (depth > maxDepth_[isubdet] && sub == isubdet) {
0864           edm::LogWarning("HcalDetId") << "HcalDetID(SimHit) presents conflicting information. Depth: " << depth
0865                                        << ", iphi: " << iphi << ", ieta: " << ieta
0866                                        << ". Max depth from geometry is: " << maxDepth_[isubdet]
0867                                        << ". TestNumber = " << testNumber_;
0868           continue;
0869         }
0870 
0871         // take cell already found to be max energy in a particular subdet
0872         if (sub == isubdet && ieta == ieta_Sim && iphi == iphi_Sim) {
0873           double en = simhits->energy();
0874 
0875           v_ehits[0] += en;
0876           v_ehits[depth] += en;
0877         }
0878       }  // simhit loop
0879 
0880       strtmp = "HcalDigiTask_amplitude_vs_simhits_" + subdet_;
0881       if (v_ehits[0] > eps)
0882         fill2D(strtmp, v_ehits[0], v_ampl_c[0]);
0883       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0884         strtmp = "HcalDigiTask_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
0885         if (v_ehits[depth] > eps)
0886           fill2D(strtmp, v_ehits[depth], v_ampl_c[depth]);
0887       }
0888 
0889       strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_" + subdet_;
0890       if (v_ehits[0] > eps)
0891         fillPf(strtmp, v_ehits[0], v_ampl_c[0]);
0892       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0893         strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
0894         if (v_ehits[depth] > eps)
0895           fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
0896       }
0897 
0898       strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_" + subdet_;
0899       if (v_ehits[0] > eps)
0900         fill1D(strtmp, v_ampl_c[0] / v_ehits[0]);
0901       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
0902         strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
0903         if (v_ehits[depth] > eps)
0904           fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
0905         strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
0906         if (v_ehits[depth] > eps)
0907           fill1D(strtmp, v_ampl_c[depth] / v_ehits[depth]);
0908       }
0909 
0910     }  // end of if(mc_ == "yes")
0911 
0912     strtmp = "HcalDigiTask_Ndigis_" + subdet_;
0913     fill1D(strtmp, double(Ndig));
0914 
0915   }  //  end of if( subdet != 0 && noise_ == 0) { // signal only
0916 }
0917 template <class dataFrameType>
0918 void HcalDigisValidation::reco(const edm::Event& iEvent,
0919                                const edm::EventSetup& iSetup,
0920                                const edm::EDGetTokenT<HcalDataFrameContainer<dataFrameType> >& tok) {
0921   // HistLim =============================================================
0922 
0923   std::string strtmp;
0924 
0925   // ======================================================================
0926   typename edm::Handle<HcalDataFrameContainer<dataFrameType> > digiHandle;
0927   //typename HcalDataFrameContainer<dataFrameType>::const_iterator digiItr;
0928 
0929   // ADC2fC
0930   CaloSamples tool;
0931   iEvent.getByToken(tok, digiHandle);
0932   if (!digiHandle.isValid())
0933     return;
0934   const HcalDataFrameContainer<dataFrameType>* digiCollection = digiHandle.product();
0935   int isubdet = 0;
0936   if (subdet_ == "HB")
0937     isubdet = 1;
0938   if (subdet_ == "HE")
0939     isubdet = 2;
0940   if (subdet_ == "HO")
0941     isubdet = 3;
0942   if (subdet_ == "HF")
0943     isubdet = 4;
0944 
0945   if (isubdet == 1)
0946     nevent1++;
0947   if (isubdet == 2)
0948     nevent2++;
0949   if (isubdet == 3)
0950     nevent3++;
0951   if (isubdet == 4)
0952     nevent4++;
0953 
0954   int indigis = 0;
0955   //  amplitude for signal cell at diff. depths
0956   std::vector<double> v_ampl_c(maxDepth_[isubdet] + 1, 0);
0957 
0958   // is set to 1 if "seed" SimHit is found
0959   int seedSimHit = 0;
0960 
0961   int ieta_Sim = 9999;
0962   int iphi_Sim = 9999;
0963   double emax_Sim = -9999.;
0964 
0965   // SimHits MC only
0966   if (mc_ == "yes") {
0967     edm::Handle<edm::PCaloHitContainer> hcalHits;
0968     iEvent.getByToken(tok_mc_, hcalHits);
0969     const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0970 
0971     if (isubdet != 0 && noise_ == 0) {  // signal only SimHits
0972 
0973       for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0974            ++simhits) {
0975         unsigned int id = simhits->id();
0976         int sub, ieta, iphi;
0977         HcalDetId hid;
0978         if (testNumber_)
0979           hid = HcalHitRelabeller::relabel(id, hcons_);
0980         else
0981           hid = HcalDetId(id);
0982         sub = hid.subdet();
0983         ieta = hid.ieta();
0984         iphi = hid.iphi();
0985 
0986         double en = simhits->energy();
0987 
0988         if (en > emax_Sim && sub == isubdet) {
0989           emax_Sim = en;
0990           ieta_Sim = ieta;
0991           iphi_Sim = iphi;
0992           // to limit "seed" SimHit energy in case of "multi" event
0993           if (mode_ == "multi" && ((sub == 4 && en < 100. && en > 1.) || ((sub != 4) && en < 1. && en > 0.02))) {
0994             seedSimHit = 1;
0995             break;
0996           }
0997         }
0998 
0999       }  // end of SimHits cycle
1000 
1001       // found highest-energy SimHit for single-particle
1002       if (mode_ != "multi" && emax_Sim > 0.)
1003         seedSimHit = 1;
1004     }  // end of SimHits
1005   }    // end of mc_ == "yes"
1006 
1007   // CYCLE OVER CELLS ========================================================
1008   int Ndig = 0;
1009 
1010   for (typename HcalDataFrameContainer<dataFrameType>::const_iterator digiItr = digiCollection->begin();
1011        digiItr != digiCollection->end();
1012        digiItr++) {
1013     dataFrameType dataFrame = *digiItr;
1014 
1015     HcalDetId cell(digiItr->id());
1016     int depth = cell.depth();
1017     int iphi = cell.iphi();
1018     int ieta = cell.ieta();
1019     int sub = cell.subdet();
1020 
1021     //Is this in HEP17
1022     bool isHEP17 = (iphi >= 63) && (iphi <= 66) && (ieta > 0) && (sub == 2);
1023 
1024     if (depth > maxDepth_[isubdet] && sub == isubdet) {
1025       edm::LogWarning("HcalDetId") << "HcalDetID presents conflicting information. Depth: " << depth
1026                                    << ", iphi: " << iphi << ", ieta: " << ieta
1027                                    << ". Max depth from geometry is: " << maxDepth_[isubdet]
1028                                    << ". TestNumber = " << testNumber_;
1029       continue;
1030     }
1031 
1032     //  amplitude for signal cell at diff. depths
1033     std::vector<double> v_ampl(maxDepth_[isubdet] + 1, 0);
1034 
1035     // Gains, pedestals (once !) and only for "noise" case
1036     if (((nevent1 == 1 && isubdet == 1) || (nevent2 == 1 && isubdet == 2) || (nevent3 == 1 && isubdet == 3) ||
1037          (nevent4 == 1 && isubdet == 4)) &&
1038         noise_ == 1 && sub == isubdet) {
1039       HcalGenericDetId hcalGenDetId(digiItr->id());
1040       const HcalPedestal* pedestal = conditions_->getPedestal(hcalGenDetId);
1041       const HcalGain* gain = conditions_->getGain(hcalGenDetId);
1042       const HcalGainWidth* gainWidth = conditions_->getGainWidth(hcalGenDetId);
1043       const HcalPedestalWidth* pedWidth = conditions_->getPedestalWidth(hcalGenDetId);
1044 
1045       for (int i = 0; i < 4; i++) {
1046         fill1D("HcalDigiTask_gain_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gain->getValue(i));
1047         fill1D("HcalDigiTask_gainWidth_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gainWidth->getValue(i));
1048         fill1D("HcalDigiTask_pedestal_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedestal->getValue(i));
1049         fill1D("HcalDigiTask_pedestal_width_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_,
1050                pedWidth->getWidth(i));
1051       }
1052 
1053       fill2D("HcalDigiTask_gainMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), gain->getValue(0));
1054       fill2D("HcalDigiTask_pwidthMap_Depth" + str(depth) + "_" + subdet_,
1055              double(ieta),
1056              double(iphi),
1057              pedWidth->getWidth(0));
1058 
1059     }  // end of event #1
1060     //edm::LogVerbatim("OutputInfo") << "==== End of event noise block in cell cycle";
1061 
1062     if (sub == isubdet)
1063       Ndig++;  // subdet number of digi
1064 
1065     // No-noise case, only single  subdet selected  ===========================
1066 
1067     if (sub == isubdet && noise_ == 0) {
1068       HcalCalibrations calibrations = conditions_->getHcalCalibrations(cell);
1069 
1070       const HcalQIECoder* channelCoder = conditions_->getHcalCoder(cell);
1071       const HcalQIEShape* shape = conditions_->getHcalShape(channelCoder);
1072       HcalCoderDb coder(*channelCoder, *shape);
1073       coder.adc2fC(dataFrame, tool);
1074 
1075       // for dynamic digi time sample analysis
1076       int soi = tool.presamples();
1077       int lastbin = tool.size() - 1;
1078 
1079       double noiseADC = (dataFrame)[0].adc();
1080       double noisefC = tool[0];
1081       // noise evaluations from "pre-samples"
1082       fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
1083       fill1D("HcalDigiTask_ADC0_fC_depth" + str(depth) + "_" + subdet_, noisefC);
1084 
1085       // OCCUPANCY maps fill
1086       fill2D("HcalDigiTask_ieta_iphi_occupancy_map_depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi));
1087 
1088       fill1D("HcalDigiTask_depths_" + subdet_, double(depth));
1089 
1090       // Cycle on time slices
1091       // - for each Digi
1092       // - for one Digi with max SimHits E in subdet
1093 
1094       int closen = 0;  // =1 if 1) seedSimHit = 1 and 2) the cell is the same
1095       if (ieta == ieta_Sim && iphi == iphi_Sim)
1096         closen = seedSimHit;
1097 
1098       for (int ii = 0; ii < tool.size(); ii++) {
1099         int capid = (dataFrame)[ii].capid();
1100         // single ts amplitude
1101         double val = (tool[ii] - calibrations.pedestal(capid));
1102 
1103         if (val > 100.) {
1104           fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
1105           if (hep17_) {
1106             if (!isHEP17) {
1107               strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
1108               fill1D(strtmp, double(ii), val);
1109             } else {
1110               strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_HEP17";
1111               fill1D(strtmp, double(ii), val);
1112             }
1113           } else {
1114             strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
1115             fill1D(strtmp, double(ii), val);
1116           }
1117         }
1118 
1119         if (closen == 1) {
1120           if (hep17_) {
1121             if (!isHEP17) {
1122               strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
1123               fill2D(strtmp, double(ii), val);
1124             } else {
1125               strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_HEP17";
1126               fill2D(strtmp, double(ii), val);
1127             }
1128           } else {
1129             strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
1130             fill2D(strtmp, double(ii), val);
1131           }
1132         }
1133 
1134         // all detectors
1135         if (ii >= soi && ii <= lastbin) {
1136           v_ampl[0] += val;
1137           v_ampl[depth] += val;
1138 
1139           if (closen == 1) {
1140             v_ampl_c[0] += val;
1141             v_ampl_c[depth] += val;
1142           }
1143         }
1144 
1145         //...TDC
1146 
1147         if ((HBPhase1_ && sub == 1) || (HEPhase1_ && sub == 2)) {
1148           double digiADC = (dataFrame)[ii].adc();
1149           const QIE11DataFrame dataFrameHBHE = static_cast<const QIE11DataFrame>(*digiItr);
1150           double digiTDC = (dataFrameHBHE)[ii].tdc();
1151           if (digiTDC < 50) {
1152             double time = ii * 25. + (digiTDC * 0.5);
1153             strtmp = "HcalDigiTask_TDCtime_" + subdet_;
1154             fill1D(strtmp, time);
1155 
1156             strtmp = "HcalDigiTask_TDCtime_vs_ADC_" + subdet_;
1157             fill2D(strtmp, digiADC, time);
1158           }
1159         }
1160       }
1161       // end of time bucket sample
1162 
1163       // just 1D of all cells' amplitudes
1164       strtmp = "HcalDigiTask_sum_all_amplitudes_" + subdet_;
1165       fill1D(strtmp, v_ampl[0]);
1166 
1167       std::vector<int> v_ampl_sub(v_ampl.begin() + 1, v_ampl.end());  // remove element 0, which is the sum of any depth
1168       double ampl_max = *std::max_element(v_ampl_sub.begin(), v_ampl_sub.end());
1169       if (ampl_max > 10.)
1170         indigis++;
1171       //KH if (ampl1 > 10. || ampl2 > 10. || ampl3 > 10. || ampl4 > 10.) indigis++;
1172 
1173       // fraction 5,6 bins if ampl. is big.
1174       //histogram names have not been changed, but it should be understood that bin_5 is soi, and bin_6_7 is latter TS'
1175       if ((v_ampl[depth] > 30. && (subdet_ != "HE" || subdet_ != "HB")) ||
1176           (v_ampl[depth] > 300.)) {  //300 fC cut for QIE-11 HB & HE
1177         double fbinSOI = tool[soi] - calibrations.pedestal((dataFrame)[soi].capid());
1178         double fbinPS = 0;
1179 
1180         for (int j = soi + 1; j <= lastbin; j++)
1181           fbinPS += tool[j] - calibrations.pedestal((dataFrame)[j].capid());
1182 
1183         fbinSOI /= v_ampl[depth];
1184         fbinPS /= v_ampl[depth];
1185         strtmp = "HcalDigiTask_SOI_frac_" + subdet_;
1186         fill1D(strtmp, fbinSOI);
1187         strtmp = "HcalDigiTask_postSOI_frac_" + subdet_;
1188         fill1D(strtmp, fbinPS);
1189       }
1190 
1191       if (hep17_) {
1192         if (!isHEP17) {
1193           strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
1194           fill1D(strtmp, v_ampl[0]);
1195           strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
1196           fill1D(strtmp, v_ampl[depth]);
1197         } else {
1198           strtmp = "HcalDigiTask_signal_amplitude_HEP17";
1199           fill1D(strtmp, v_ampl[0]);
1200           strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_HEP17";
1201           fill1D(strtmp, v_ampl[depth]);
1202         }
1203       } else {
1204         strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
1205         fill1D(strtmp, v_ampl[0]);
1206         strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
1207         fill1D(strtmp, v_ampl[depth]);
1208       }
1209     }
1210   }  // End of CYCLE OVER CELLS =============================================
1211 
1212   if (isubdet != 0 && noise_ == 0) {  // signal only, once per event
1213     strtmp = "HcalDigiTask_number_of_amplitudes_above_10fC_" + subdet_;
1214     fill1D(strtmp, indigis);
1215 
1216     // SimHits once again !!!
1217     double eps = 1.e-3;
1218     std::vector<double> v_ehits(maxDepth_[isubdet] + 1, 0);
1219 
1220     if (mc_ == "yes") {
1221       edm::Handle<edm::PCaloHitContainer> hcalHits;
1222       iEvent.getByToken(tok_mc_, hcalHits);
1223       const edm::PCaloHitContainer* simhitResult = hcalHits.product();
1224       for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
1225            ++simhits) {
1226         unsigned int id = simhits->id();
1227         int sub, depth, ieta, iphi;
1228         HcalDetId hid;
1229         if (testNumber_)
1230           hid = HcalHitRelabeller::relabel(id, hcons_);
1231         else
1232           hid = HcalDetId(id);
1233         sub = hid.subdet();
1234         depth = hid.depth();
1235         ieta = hid.ieta();
1236         iphi = hid.iphi();
1237 
1238         if (depth > maxDepth_[isubdet] && sub == isubdet) {
1239           edm::LogWarning("HcalDetId") << "HcalDetID(SimHit) presents conflicting information. Depth: " << depth
1240                                        << ", iphi: " << iphi << ", ieta: " << ieta
1241                                        << ". Max depth from geometry is: " << maxDepth_[isubdet]
1242                                        << ". TestNumber = " << testNumber_;
1243           continue;
1244         }
1245 
1246         // take cell already found to be max energy in a particular subdet
1247         if (sub == isubdet && ieta == ieta_Sim && iphi == iphi_Sim) {
1248           double en = simhits->energy();
1249 
1250           v_ehits[0] += en;
1251           v_ehits[depth] += en;
1252         }
1253       }  // simhit loop
1254 
1255       strtmp = "HcalDigiTask_amplitude_vs_simhits_" + subdet_;
1256       if (v_ehits[0] > eps)
1257         fill2D(strtmp, v_ehits[0], v_ampl_c[0]);
1258       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1259         strtmp = "HcalDigiTask_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
1260         if (v_ehits[depth] > eps)
1261           fill2D(strtmp, v_ehits[depth], v_ampl_c[depth]);
1262       }
1263 
1264       strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_" + subdet_;
1265       if (v_ehits[0] > eps)
1266         fillPf(strtmp, v_ehits[0], v_ampl_c[0]);
1267       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1268         strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
1269         if (v_ehits[depth] > eps)
1270           fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
1271       }
1272 
1273       strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_" + subdet_;
1274       if (v_ehits[0] > eps)
1275         fill1D(strtmp, v_ampl_c[0] / v_ehits[0]);
1276       for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1277         strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
1278         if (v_ehits[depth] > eps)
1279           fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
1280         strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
1281         if (v_ehits[depth] > eps)
1282           fill1D(strtmp, v_ampl_c[depth] / v_ehits[depth]);
1283       }
1284 
1285     }  // end of if(mc_ == "yes")
1286 
1287     strtmp = "HcalDigiTask_Ndigis_" + subdet_;
1288     fill1D(strtmp, double(Ndig));
1289 
1290   }  //  end of if( subdet != 0 && noise_ == 0) { // signal only
1291 
1292 }  //HcalDataFrameContainer
1293 
1294 void HcalDigisValidation::book1D(DQMStore::IBooker& ib, std::string name, int n, double min, double max) {
1295   if (!msm_->count(name))
1296     (*msm_)[name] = ib.book1D(name.c_str(), name.c_str(), n, min, max);
1297 }
1298 
1299 void HcalDigisValidation::book1D(DQMStore::IBooker& ib, std::string name, const HistLim& limX) {
1300   if (!msm_->count(name))
1301     (*msm_)[name] = ib.book1D(name.c_str(), name.c_str(), limX.n, limX.min, limX.max);
1302 }
1303 
1304 void HcalDigisValidation::fill1D(std::string name, double X, double weight) {
1305   msm_->find(name)->second->Fill(X, weight);
1306 }
1307 
1308 void HcalDigisValidation::book2D(DQMStore::IBooker& ib, std::string name, const HistLim& limX, const HistLim& limY) {
1309   if (!msm_->count(name))
1310     (*msm_)[name] = ib.book2D(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max);
1311 }
1312 
1313 void HcalDigisValidation::fill2D(std::string name, double X, double Y, double weight) {
1314   msm_->find(name)->second->Fill(X, Y, weight);
1315 }
1316 
1317 void HcalDigisValidation::bookPf(DQMStore::IBooker& ib, std::string name, const HistLim& limX, const HistLim& limY) {
1318   if (!msm_->count(name))
1319     (*msm_)[name] = ib.bookProfile(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max);
1320 }
1321 
1322 void HcalDigisValidation::bookPf(
1323     DQMStore::IBooker& ib, std::string name, const HistLim& limX, const HistLim& limY, const char* option) {
1324   if (!msm_->count(name))
1325     (*msm_)[name] =
1326         ib.bookProfile(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max, option);
1327 }
1328 
1329 void HcalDigisValidation::fillPf(std::string name, double X, double Y) { msm_->find(name)->second->Fill(X, Y); }
1330 
1331 HcalDigisValidation::MonitorElement* HcalDigisValidation::monitor(std::string name) {
1332   if (!msm_->count(name))
1333     return nullptr;
1334   else
1335     return msm_->find(name)->second;
1336 }
1337 
1338 std::string HcalDigisValidation::str(int x) {
1339   std::stringstream out;
1340   out << x;
1341   return out.str();
1342 }
1343 
1344 //define this as a plug-in
1345 DEFINE_FWK_MODULE(HcalDigisValidation);
[ 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!