Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-14 23:41:23 

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