Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-17 03:23:34

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "FWCore/Utilities/interface/transform.h"
0014 
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0017 
0018 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0019 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0020 
0021 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0022 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0023 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0024 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0025 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0026 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0027 
0028 #include <TH1F.h>
0029 
0030 #include <memory>
0031 #include <iostream>
0032 #include <fstream>
0033 #include <vector>
0034 #include <map>
0035 #include <string>
0036 
0037 //#define EDM_ML_DEBUG
0038 
0039 class CaloSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0040 public:
0041   CaloSimHitStudy(const edm::ParameterSet& ps);
0042   ~CaloSimHitStudy() override {}
0043   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044 
0045 protected:
0046   void beginJob() override {}
0047   void analyze(edm::Event const&, edm::EventSetup const&) override;
0048   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0049   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0050 
0051   void analyzeHits(std::vector<PCaloHit>&, int);
0052   void analyzeHits(const edm::Handle<edm::PSimHitContainer>&, int);
0053   void analyzeHits(std::vector<PCaloHit>&, std::vector<PCaloHit>&, std::vector<PCaloHit>&);
0054 
0055 private:
0056   const std::string g4Label_;
0057   const std::vector<std::string> hitLab_;
0058   const double maxEnergy_, tmax_, eMIP_;
0059   const bool storeRL_, testNumber_;
0060   edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
0061   std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>> toks_calo_;
0062   std::vector<edm::EDGetTokenT<edm::PSimHitContainer>> toks_track_;
0063   std::vector<edm::EDGetTokenT<edm::PSimHitContainer>> toks_tkHigh_;
0064   std::vector<edm::EDGetTokenT<edm::PSimHitContainer>> toks_tkLow_;
0065 
0066   const std::vector<std::string> muonLab_ = {"MuonRPCHits", "MuonCSCHits", "MuonDTHits", "MuonGEMHits"};
0067   const std::vector<std::string> tkHighLab_ = {"TrackerHitsPixelBarrelHighTof",
0068                                                "TrackerHitsPixelEndcapHighTof",
0069                                                "TrackerHitsTECHighTof",
0070                                                "TrackerHitsTIBHighTof",
0071                                                "TrackerHitsTIDHighTof",
0072                                                "TrackerHitsTOBHighTof"};
0073   const std::vector<std::string> tkLowLab_ = {"TrackerHitsPixelBarrelLowTof",
0074                                               "TrackerHitsPixelEndcapLowTof",
0075                                               "TrackerHitsTECLowTof",
0076                                               "TrackerHitsTIBLowTof",
0077                                               "TrackerHitsTIDLowTof",
0078                                               "TrackerHitsTOBLowTof"};
0079 
0080   TH1F *hit_[9], *time_[9], *edepEM_[9], *edepHad_[9], *edep_[9];
0081   TH1F *etot_[9], *etotg_[9], *timeAll_[9], *hitMu, *hitHigh;
0082   TH1F *hitLow, *eneInc_, *etaInc_, *phiInc_, *ptInc_;
0083   TH1F *hitTk_[16], *edepTk_[16], *tofTk_[16], *edepC_[9], *edepT_[9];
0084 };
0085 
0086 CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps)
0087     : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel")),
0088       hitLab_(ps.getUntrackedParameter<std::vector<std::string>>("CaloCollection")),
0089       maxEnergy_(ps.getUntrackedParameter<double>("MaxEnergy", 200.0)),
0090       tmax_(ps.getUntrackedParameter<double>("TimeCut", 100.0)),
0091       eMIP_(ps.getUntrackedParameter<double>("MIPCut", 0.70)),
0092       storeRL_(ps.getUntrackedParameter<bool>("StoreRL", false)),
0093       testNumber_(ps.getUntrackedParameter<bool>("TestNumbering", true)) {
0094   usesResource(TFileService::kSharedResource);
0095 
0096   tok_evt_ =
0097       consumes<edm::HepMCProduct>(edm::InputTag(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")));
0098   for (unsigned i = 0; i != hitLab_.size(); i++)
0099     toks_calo_.emplace_back(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hitLab_[i])));
0100   for (unsigned i = 0; i != muonLab_.size(); i++)
0101     toks_track_.emplace_back(consumes<edm::PSimHitContainer>(edm::InputTag(g4Label_, muonLab_[i])));
0102   for (unsigned i = 0; i != tkHighLab_.size(); i++)
0103     toks_tkHigh_.emplace_back(consumes<edm::PSimHitContainer>(edm::InputTag(g4Label_, tkHighLab_[i])));
0104   for (unsigned i = 0; i != tkLowLab_.size(); i++)
0105     toks_tkLow_.emplace_back(consumes<edm::PSimHitContainer>(edm::InputTag(g4Label_, tkLowLab_[i])));
0106 
0107   edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label_ << "   Hits: " << hitLab_[0] << ", " << hitLab_[1]
0108                                << ", " << hitLab_[2] << ", " << hitLab_[3] << "   MaxEnergy: " << maxEnergy_
0109                                << "  Tmax: " << tmax_ << "  MIP Cut: " << eMIP_;
0110 
0111   edm::Service<TFileService> tfile;
0112 
0113   if (!tfile.isAvailable())
0114     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0115                                       << "please add it to config file";
0116   char name[20], title[120];
0117   sprintf(title, "Incident PT (GeV)");
0118   ptInc_ = tfile->make<TH1F>("PtInc", title, 1000, 0., maxEnergy_);
0119   ptInc_->GetXaxis()->SetTitle(title);
0120   ptInc_->GetYaxis()->SetTitle("Events");
0121   sprintf(title, "Incident Energy (GeV)");
0122   eneInc_ = tfile->make<TH1F>("EneInc", title, 1000, 0., maxEnergy_);
0123   eneInc_->GetXaxis()->SetTitle(title);
0124   eneInc_->GetYaxis()->SetTitle("Events");
0125   sprintf(title, "Incident #eta");
0126   etaInc_ = tfile->make<TH1F>("EtaInc", title, 200, -5., 5.);
0127   etaInc_->GetXaxis()->SetTitle(title);
0128   etaInc_->GetYaxis()->SetTitle("Events");
0129   sprintf(title, "Incident #phi");
0130   phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
0131   phiInc_->GetXaxis()->SetTitle(title);
0132   phiInc_->GetYaxis()->SetTitle("Events");
0133 #ifdef EDM_ML_DEBUG
0134   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for incident particle";
0135 #endif
0136   std::string dets[9] = {"EB", "EB(APD)", "EB(ATJ)", "EE", "ES", "HB", "HE", "HO", "HF"};
0137   double nhcMax[9] = {40000., 2000., 2000., 40000., 10000., 10000., 10000., 2000., 10000.};
0138   for (int i = 0; i < 9; i++) {
0139     sprintf(name, "Hit%d", i);
0140     sprintf(title, "Number of hits in %s", dets[i].c_str());
0141     hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., nhcMax[i]);
0142     hit_[i]->GetXaxis()->SetTitle(title);
0143     hit_[i]->GetYaxis()->SetTitle("Events");
0144     sprintf(name, "Time%d", i);
0145     sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
0146     time_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
0147     time_[i]->GetXaxis()->SetTitle(title);
0148     time_[i]->GetYaxis()->SetTitle("Hits");
0149     sprintf(name, "TimeAll%d", i);
0150     sprintf(title, "Hit time (ns) in %s with no check on Track ID", dets[i].c_str());
0151     timeAll_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
0152     timeAll_[i]->GetXaxis()->SetTitle(title);
0153     timeAll_[i]->GetYaxis()->SetTitle("Hits");
0154     double ymax = maxEnergy_;
0155     if (i == 1 || i == 2 || i == 4)
0156       ymax = 1.0;
0157     else if (i > 4 && i < 8)
0158       ymax = 10.0;
0159     sprintf(name, "Edep%d", i);
0160     sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
0161     edep_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0162     edep_[i]->GetXaxis()->SetTitle(title);
0163     edep_[i]->GetYaxis()->SetTitle("Hits");
0164     sprintf(name, "EdepEM%d", i);
0165     sprintf(title, "Energy deposit (GeV) by EM particles in %s", dets[i].c_str());
0166     edepEM_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0167     edepEM_[i]->GetXaxis()->SetTitle(title);
0168     edepEM_[i]->GetYaxis()->SetTitle("Hits");
0169     sprintf(name, "EdepHad%d", i);
0170     sprintf(title, "Energy deposit (GeV) by hadrons in %s", dets[i].c_str());
0171     edepHad_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0172     edepHad_[i]->GetXaxis()->SetTitle(title);
0173     edepHad_[i]->GetYaxis()->SetTitle("Hits");
0174     sprintf(name, "Etot%d", i);
0175     sprintf(title, "Total energy deposit (GeV) in %s", dets[i].c_str());
0176     etot_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0177     etot_[i]->GetXaxis()->SetTitle(title);
0178     etot_[i]->GetYaxis()->SetTitle("Events");
0179     sprintf(name, "EtotG%d", i);
0180     sprintf(title, "Total energy deposit (GeV) in %s (t < 100 ns)", dets[i].c_str());
0181     etotg_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0182     etotg_[i]->GetXaxis()->SetTitle(title);
0183     etotg_[i]->GetYaxis()->SetTitle("Events");
0184   }
0185 #ifdef EDM_ML_DEBUG
0186   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for first level of Calorimeter";
0187 #endif
0188   std::string detx[9] = {"EB/EE (MIP)",
0189                          "HB/HE (MIP)",
0190                          "HB/HE/HO (MIP)",
0191                          "EB/EE (no MIP)",
0192                          "HB/HE (no MIP)",
0193                          "HB/HE/HO (no MIP)",
0194                          "EB/EE (All)",
0195                          "HB/HE (All)",
0196                          "HB/HE/HO (All)"};
0197   for (int i = 0; i < 9; i++) {
0198     double ymax = 1.0;
0199     if (i == 0 || i == 3 || i == 6)
0200       ymax = maxEnergy_;
0201     sprintf(name, "EdepCal%d", i);
0202     sprintf(title, "Energy deposit in %s", detx[i].c_str());
0203     edepC_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0204     edepC_[i]->GetXaxis()->SetTitle(title);
0205     edepC_[i]->GetYaxis()->SetTitle("Events");
0206     sprintf(name, "EdepCalT%d", i);
0207     sprintf(title, "Energy deposit (t < %f ns) in %s", tmax_, detx[i].c_str());
0208     edepT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0209     edepT_[i]->GetXaxis()->SetTitle(title);
0210     edepT_[i]->GetYaxis()->SetTitle("Events");
0211   }
0212 #ifdef EDM_ML_DEBUG
0213   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for second level of Calorimeter";
0214 #endif
0215   hitLow = tfile->make<TH1F>("HitLow", "Number of hits in Track (Low)", 1000, 0, 20000.);
0216   hitLow->GetXaxis()->SetTitle("Number of hits in Track (Low)");
0217   hitLow->GetYaxis()->SetTitle("Events");
0218   hitHigh = tfile->make<TH1F>("HitHigh", "Number of hits in Track (High)", 1000, 0, 5000.);
0219   hitHigh->GetXaxis()->SetTitle("Number of hits in Track (High)");
0220   hitHigh->GetYaxis()->SetTitle("Events");
0221   hitMu = tfile->make<TH1F>("HitMu", "Number of hits in Track (Muon)", 1000, 0, 2000.);
0222   hitMu->GetXaxis()->SetTitle("Number of hits in Muon");
0223   hitMu->GetYaxis()->SetTitle("Events");
0224 #ifdef EDM_ML_DEBUG
0225   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for general tracking hits";
0226 #endif
0227   std::string dett[16] = {"Pixel Barrel (High)",
0228                           "Pixel Endcap (High)",
0229                           "TEC (High)",
0230                           "TIB (High)",
0231                           "TID (High)",
0232                           "TOB (High)",
0233                           "Pixel Barrel (Low)",
0234                           "Pixel Endcap (Low)",
0235                           "TEC (Low)",
0236                           "TIB (Low)",
0237                           "TID (Low)",
0238                           "TOB (Low)",
0239                           "RPC",
0240                           "CSC",
0241                           "DT",
0242                           "GEM"};
0243   double nhtMax[16] = {
0244       500., 500., 1000., 1000., 500., 1000., 5000., 2000., 10000., 5000., 2000., 5000., 500., 1000., 1000., 500.};
0245   for (int i = 0; i < 16; i++) {
0246     sprintf(name, "HitTk%d", i);
0247     sprintf(title, "Number of hits in %s", dett[i].c_str());
0248     hitTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., nhtMax[i]);
0249     hitTk_[i]->GetXaxis()->SetTitle(title);
0250     hitTk_[i]->GetYaxis()->SetTitle("Events");
0251     sprintf(name, "TimeTk%d", i);
0252     sprintf(title, "Time of the hit (ns) in %s", dett[i].c_str());
0253     tofTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., 200.);
0254     tofTk_[i]->GetXaxis()->SetTitle(title);
0255     tofTk_[i]->GetYaxis()->SetTitle("Hits");
0256     sprintf(name, "EdepTk%d", i);
0257     sprintf(title, "Energy deposit (GeV) in %s", dett[i].c_str());
0258     double ymax = (i < 12) ? 0.25 : 0.005;
0259     edepTk_[i] = tfile->make<TH1F>(name, title, 250, 0., ymax);
0260     edepTk_[i]->GetXaxis()->SetTitle(title);
0261     edepTk_[i]->GetYaxis()->SetTitle("Hits");
0262   }
0263 #ifdef EDM_ML_DEBUG
0264   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for SimHit objects";
0265 #endif
0266 }
0267 
0268 void CaloSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   edm::ParameterSetDescription desc;
0270   std::vector<std::string> calonames = {"EcalHitsEB", "EcalHitsEE", "EcalHitsES", "HcalHits"};
0271   desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
0272   desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0273   desc.addUntracked<std::vector<std::string>>("CaloCollection", calonames);
0274   desc.addUntracked<double>("MaxEnergy", 200.0);
0275   desc.addUntracked<double>("TimeCut", 100.0);
0276   desc.addUntracked<double>("MIPCut", 0.70);
0277   desc.addUntracked<bool>("StoreRL", false);
0278   desc.addUntracked<bool>("TestNumbering", true);
0279   descriptions.add("CaloSimHitStudy", desc);
0280 }
0281 
0282 void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) {
0283   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy:Run = " << e.id().run() << " Event = " << e.id().event();
0284 
0285   const edm::Handle<edm::HepMCProduct> EvtHandle = e.getHandle(tok_evt_);
0286   const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent();
0287 
0288   double eInc = 0, etaInc = 0, phiInc = 0;
0289   HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
0290   if (p != myGenEvent->particles_end()) {
0291     eInc = (*p)->momentum().e();
0292     etaInc = (*p)->momentum().eta();
0293     phiInc = (*p)->momentum().phi();
0294   }
0295   double ptInc = eInc / std::cosh(etaInc);
0296   ptInc_->Fill(ptInc);
0297   eneInc_->Fill(eInc);
0298   etaInc_->Fill(etaInc);
0299   phiInc_->Fill(phiInc);
0300 
0301   std::vector<PCaloHit> ebHits, eeHits, hcHits;
0302   for (unsigned int i = 0; i < toks_calo_.size(); i++) {
0303     bool getHits = false;
0304     const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(toks_calo_[i]);
0305     if (hitsCalo.isValid())
0306       getHits = true;
0307 #ifdef EDM_ML_DEBUG
0308     edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Input flags Hits " << getHits;
0309 #endif
0310     if (getHits) {
0311       std::vector<PCaloHit> caloHits;
0312       caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
0313       if (i == 0)
0314         ebHits.insert(ebHits.end(), hitsCalo->begin(), hitsCalo->end());
0315       else if (i == 1)
0316         eeHits.insert(eeHits.end(), hitsCalo->begin(), hitsCalo->end());
0317       else if (i == 3)
0318         hcHits.insert(hcHits.end(), hitsCalo->begin(), hitsCalo->end());
0319 #ifdef EDM_ML_DEBUG
0320       edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Hit buffer " << caloHits.size();
0321 #endif
0322       analyzeHits(caloHits, i);
0323     }
0324   }
0325   analyzeHits(ebHits, eeHits, hcHits);
0326 
0327   std::vector<PSimHit> muonHits;
0328   for (unsigned int i = 0; i < toks_track_.size(); i++) {
0329     const edm::Handle<edm::PSimHitContainer>& hitsTrack = e.getHandle(toks_track_[i]);
0330     if (hitsTrack.isValid()) {
0331       muonHits.insert(muonHits.end(), hitsTrack->begin(), hitsTrack->end());
0332       analyzeHits(hitsTrack, i + 12);
0333     }
0334   }
0335   unsigned int nhmu = muonHits.size();
0336   hitMu->Fill(double(nhmu));
0337   std::vector<PSimHit> tkHighHits;
0338   for (unsigned int i = 0; i < toks_tkHigh_.size(); i++) {
0339     const edm::Handle<edm::PSimHitContainer>& hitsTrack = e.getHandle(toks_tkHigh_[i]);
0340     if (hitsTrack.isValid()) {
0341       tkHighHits.insert(tkHighHits.end(), hitsTrack->begin(), hitsTrack->end());
0342       analyzeHits(hitsTrack, i);
0343     }
0344   }
0345   unsigned int nhtkh = tkHighHits.size();
0346   hitHigh->Fill(double(nhtkh));
0347   std::vector<PSimHit> tkLowHits;
0348   for (unsigned int i = 0; i < toks_tkLow_.size(); i++) {
0349     const edm::Handle<edm::PSimHitContainer>& hitsTrack = e.getHandle(toks_tkLow_[i]);
0350     if (hitsTrack.isValid()) {
0351       tkLowHits.insert(tkLowHits.end(), hitsTrack->begin(), hitsTrack->end());
0352       analyzeHits(hitsTrack, i + 6);
0353     }
0354   }
0355   unsigned int nhtkl = tkLowHits.size();
0356   hitLow->Fill(double(nhtkl));
0357 }
0358 
0359 void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
0360   int nHit = hits.size();
0361   int nHB(0), nHE(0), nHO(0), nHF(0), nEB(0), nEBAPD(0), nEBATJ(0), nEE(0), nES(0);
0362 #ifdef EDM_ML_DEBUG
0363   int nBad(0);
0364 #endif
0365   std::map<unsigned int, double> hitMap;
0366   std::vector<double> etot(9, 0), etotG(9, 0);
0367   for (int i = 0; i < nHit; i++) {
0368     double edep = hits[i].energy();
0369     double time = hits[i].time();
0370     unsigned int id_ = hits[i].id();
0371     double edepEM = hits[i].energyEM();
0372     double edepHad = hits[i].energyHad();
0373     if (indx == 0) {
0374       int dep = (hits[i].depth()) & PCaloHit::kEcalDepthIdMask;
0375       if (dep == 1)
0376         id_ |= 0x20000;
0377       else if (dep == 2)
0378         id_ |= 0x40000;
0379     }
0380     std::map<unsigned int, double>::const_iterator it = hitMap.find(id_);
0381     if (it == hitMap.end()) {
0382       hitMap.insert(std::pair<unsigned int, double>(id_, time));
0383     }
0384     int idx = -1;
0385     if (indx != 3) {
0386       if (indx == 0) {
0387         if (storeRL_)
0388           idx = 0;
0389         else
0390           idx = ((hits[i].depth()) & PCaloHit::kEcalDepthIdMask);
0391       } else
0392         idx = indx + 2;
0393       time_[idx]->Fill(time);
0394       edep_[idx]->Fill(edep);
0395       edepEM_[idx]->Fill(edepEM);
0396       edepHad_[idx]->Fill(edepHad);
0397       if (idx == 0)
0398         nEB++;
0399       else if (idx == 1)
0400         nEBAPD++;
0401       else if (idx == 2)
0402         nEBATJ++;
0403       else if (idx == 3)
0404         nEE++;
0405       else if (idx == 4)
0406         nES++;
0407 #ifdef EDM_ML_DEBUG
0408       else
0409         nBad++;
0410 #endif
0411       if (indx >= 0 && indx < 3) {
0412         etot[idx] += edep;
0413         if (time < 100)
0414           etotG[idx] += edep;
0415       }
0416     } else {
0417       int subdet(0);
0418       if (testNumber_) {
0419         int ieta(0), phi(0), z(0), lay(0), depth(0);
0420         HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay);
0421       } else {
0422         HcalDetId id = HcalDetId(id_);
0423         subdet = id.subdet();
0424       }
0425       if (subdet == static_cast<int>(HcalBarrel)) {
0426         idx = indx + 2;
0427         nHB++;
0428       } else if (subdet == static_cast<int>(HcalEndcap)) {
0429         idx = indx + 3;
0430         nHE++;
0431       } else if (subdet == static_cast<int>(HcalOuter)) {
0432         idx = indx + 4;
0433         nHO++;
0434       } else if (subdet == static_cast<int>(HcalForward)) {
0435         idx = indx + 5;
0436         nHF++;
0437       }
0438       if (idx > 0) {
0439         time_[idx]->Fill(time);
0440         edep_[idx]->Fill(edep);
0441         edepEM_[idx]->Fill(edepEM);
0442         edepHad_[idx]->Fill(edepHad);
0443         etot[idx] += edep;
0444         if (time < 100)
0445           etotG[idx] += edep;
0446       } else {
0447 #ifdef EDM_ML_DEBUG
0448         nBad++;
0449 #endif
0450       }
0451     }
0452   }
0453   if (indx < 3) {
0454     etot_[indx + 2]->Fill(etot[indx + 2]);
0455     etotg_[indx + 2]->Fill(etotG[indx + 2]);
0456     if (indx == 0) {
0457       etot_[indx]->Fill(etot[indx]);
0458       etotg_[indx]->Fill(etotG[indx]);
0459       etot_[indx + 1]->Fill(etot[indx + 1]);
0460       etotg_[indx + 1]->Fill(etotG[indx + 1]);
0461       hit_[indx]->Fill(double(nEB));
0462       hit_[indx + 1]->Fill(double(nEBAPD));
0463       hit_[indx + 2]->Fill(double(nEBATJ));
0464     } else {
0465       hit_[indx + 2]->Fill(double(nHit));
0466     }
0467   } else if (indx == 3) {
0468     hit_[5]->Fill(double(nHB));
0469     hit_[6]->Fill(double(nHE));
0470     hit_[7]->Fill(double(nHO));
0471     hit_[8]->Fill(double(nHF));
0472     for (int idx = 5; idx < 9; idx++) {
0473       etot_[idx]->Fill(etot[idx]);
0474       etotg_[idx]->Fill(etotG[idx]);
0475     }
0476   }
0477 
0478 #ifdef EDM_ML_DEBUG
0479   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: EB " << nEB << ", " << nEBAPD << ", " << nEBATJ
0480                                << " EE " << nEE << " ES " << nES << " HB " << nHB << " HE " << nHE << " HO " << nHO
0481                                << " HF " << nHF << " Bad " << nBad << " All " << nHit << " Reduced " << hitMap.size();
0482 #endif
0483   std::map<unsigned int, double>::const_iterator it = hitMap.begin();
0484   for (; it != hitMap.end(); it++) {
0485     double time = it->second;
0486     unsigned int id_ = (it->first);
0487     int idx = -1;
0488     if (indx < 3) {
0489       if (indx == 0) {
0490         if ((id_ & 0x20000) != 0)
0491           idx = indx + 1;
0492         else if ((id_ & 0x40000) != 0)
0493           idx = indx + 1;
0494         else
0495           idx = indx;
0496       } else {
0497         idx = indx + 2;
0498       }
0499       if (idx >= 0 && idx < 5)
0500         timeAll_[idx]->Fill(time);
0501     } else if (indx == 3) {
0502       int idx(-1), subdet(0);
0503       if (testNumber_) {
0504         int ieta(0), phi(0), z(0), lay(0), depth(0);
0505         HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay);
0506       } else {
0507         HcalDetId id = HcalDetId(id_);
0508         subdet = id.subdet();
0509       }
0510       if (subdet == static_cast<int>(HcalBarrel)) {
0511         idx = indx + 2;
0512       } else if (subdet == static_cast<int>(HcalEndcap)) {
0513         idx = indx + 3;
0514       } else if (subdet == static_cast<int>(HcalOuter)) {
0515         idx = indx + 4;
0516       } else if (subdet == static_cast<int>(HcalForward)) {
0517         idx = indx + 5;
0518       }
0519       if (idx > 0) {
0520         timeAll_[idx]->Fill(time);
0521       }
0522     }
0523   }
0524 }
0525 
0526 void CaloSimHitStudy::analyzeHits(const edm::Handle<edm::PSimHitContainer>& hits, int indx) {
0527   int nHit = 0;
0528   edm::PSimHitContainer::const_iterator ihit;
0529   std::string label(" ");
0530   if (indx >= 0 && indx < 6)
0531     label = tkHighLab_[indx];
0532   else if (indx >= 6 && indx < 12)
0533     label = tkLowLab_[indx - 6];
0534   else if (indx >= 12 && indx < 16)
0535     label = muonLab_[indx - 12];
0536   for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
0537     edepTk_[indx]->Fill(ihit->energyLoss());
0538     tofTk_[indx]->Fill(ihit->timeOfFlight());
0539     nHit++;
0540   }
0541   hitTk_[indx]->Fill(float(nHit));
0542 #ifdef EDM_ML_DEBUG
0543   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: for " << label << " Index " << indx << " # of Hits "
0544                                << nHit;
0545 #endif
0546 }
0547 
0548 void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& ebHits,
0549                                   std::vector<PCaloHit>& eeHits,
0550                                   std::vector<PCaloHit>& hcHits) {
0551   double edepEB = 0, edepEBT = 0;
0552   for (unsigned int i = 0; i < ebHits.size(); i++) {
0553     double edep = ebHits[i].energy();
0554     double time = ebHits[i].time();
0555     if (((ebHits[i].depth()) & PCaloHit::kEcalDepthIdMask) == 0) {
0556       edepEB += edep;
0557       if (time < tmax_)
0558         edepEBT += edep;
0559     }
0560   }
0561   double edepEE = 0, edepEET = 0;
0562   for (unsigned int i = 0; i < eeHits.size(); i++) {
0563     double edep = eeHits[i].energy();
0564     double time = eeHits[i].time();
0565     edepEE += edep;
0566     if (time < tmax_)
0567       edepEET += edep;
0568   }
0569   double edepH = 0, edepHT = 0, edepHO = 0, edepHOT = 0;
0570   for (unsigned int i = 0; i < hcHits.size(); i++) {
0571     double edep = hcHits[i].energy();
0572     double time = hcHits[i].time();
0573     int subdet(0);
0574     if (testNumber_) {
0575       int ieta(0), phi(0), z(0), lay(0), depth(0);
0576       HcalTestNumbering::unpackHcalIndex(hcHits[i].id(), subdet, z, depth, ieta, phi, lay);
0577     } else {
0578       HcalDetId id = HcalDetId(hcHits[i].id());
0579       subdet = id.subdet();
0580     }
0581     if (subdet == static_cast<int>(HcalBarrel) || subdet == static_cast<int>(HcalEndcap)) {
0582       edepH += edep;
0583       if (time < tmax_)
0584         edepHT += edep;
0585     } else if (subdet == static_cast<int>(HcalOuter)) {
0586       edepHO += edep;
0587       if (time < tmax_)
0588         edepHOT += edep;
0589     }
0590   }
0591   double edepE = edepEB + edepEE;
0592   double edepET = edepEBT + edepEET;
0593   double edepHC = edepH + edepHO;
0594   double edepHCT = edepHT + edepHOT;
0595 #ifdef EDM_ML_DEBUG
0596   edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::energy in EB " << edepEB << " (" << edepEBT << ") from "
0597                                << ebHits.size() << " hits; "
0598                                << " energy in EE " << edepEE << " (" << edepEET << ") from " << eeHits.size()
0599                                << " hits; energy in HC " << edepH << ", " << edepHO << " (" << edepHT << ", " << edepHOT
0600                                << ") from " << hcHits.size() << " hits";
0601 #endif
0602   edepC_[6]->Fill(edepE);
0603   edepT_[6]->Fill(edepET);
0604   edepC_[7]->Fill(edepH);
0605   edepT_[7]->Fill(edepHT);
0606   edepC_[8]->Fill(edepHC);
0607   edepT_[8]->Fill(edepHCT);
0608   if (edepE < eMIP_) {
0609     edepC_[0]->Fill(edepE);
0610     edepC_[1]->Fill(edepH);
0611     edepC_[2]->Fill(edepHC);
0612   } else {
0613     edepC_[3]->Fill(edepE);
0614     edepC_[4]->Fill(edepH);
0615     edepC_[5]->Fill(edepHC);
0616   }
0617   if (edepET < eMIP_) {
0618     edepT_[0]->Fill(edepET);
0619     edepT_[1]->Fill(edepHT);
0620     edepT_[2]->Fill(edepHCT);
0621   } else {
0622     edepT_[3]->Fill(edepET);
0623     edepT_[4]->Fill(edepHT);
0624     edepT_[5]->Fill(edepHCT);
0625   }
0626 }
0627 
0628 //define this as a plug-in
0629 DEFINE_FWK_MODULE(CaloSimHitStudy);