Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-29 02:41:39 

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