Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-18 08:23:50

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 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/Utilities/interface/transform.h"
0013 
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016 
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0019 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0020 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0021 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0022 
0023 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0024 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0025 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0026 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0027 
0028 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0029 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0030 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
0031 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0032 
0033 #include <TH1F.h>
0034 #include <TH2F.h>
0035 
0036 #include <memory>
0037 #include <iostream>
0038 #include <fstream>
0039 #include <vector>
0040 #include <map>
0041 #include <string>
0042 
0043 //#define EDM_ML_DEBUG
0044 
0045 class CaloSimHitAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0046 public:
0047   CaloSimHitAnalysis(const edm::ParameterSet& ps);
0048   ~CaloSimHitAnalysis() override {}
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050 
0051 protected:
0052   void beginJob() override {}
0053   void analyze(edm::Event const&, edm::EventSetup const&) override;
0054   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0055   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0056 
0057   void analyzeHits(std::vector<PCaloHit>&, int);
0058   void analyzePassiveHits(std::vector<PassiveHit>& hits);
0059 
0060 private:
0061   const std::string g4Label_;
0062   const std::vector<std::string> hitLab_;
0063   const std::vector<double> timeSliceUnit_;
0064   const double maxEnergy_, maxTime_, tMax_, tScale_, tCut_;
0065   const bool testNumber_, passive_;
0066   const int allSteps_;
0067   const std::vector<std::string> detNames_;
0068   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokGeom_;
0069   const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer> > toks_calo_;
0070   const edm::EDGetTokenT<edm::PassiveHitContainer> tok_passive_;
0071 
0072   const CaloGeometry* caloGeometry_;
0073   const HcalGeometry* hcalGeom_;
0074 
0075   static constexpr int nCalo_ = 6;
0076   TH1F *h_hit_[nCalo_], *h_time_[nCalo_], *h_edep_[nCalo_], *h_edepT_[nCalo_];
0077   TH1F *h_timeT_[nCalo_], *h_edep1_[nCalo_], *h_edepT1_[nCalo_];
0078   TH1F *h_edepEM_[nCalo_], *h_edepHad_[nCalo_], *h_rr_[nCalo_], *h_zz_[nCalo_];
0079   TH1F *h_eta_[nCalo_], *h_phi_[nCalo_], *h_etot_[nCalo_], *h_etotg_[nCalo_];
0080   TH2F *h_rz_, *h_rz1_, *h_etaphi_;
0081   TH1F *h_hitp_, *h_trackp_, *h_edepp_, *h_timep_, *h_stepp_;
0082   std::vector<TH1F*> h_edepTk_, h_timeTk_;
0083   std::vector<TH2F*> h_rzH_;
0084   std::map<int, unsigned int> etaDepth_;
0085 };
0086 
0087 CaloSimHitAnalysis::CaloSimHitAnalysis(const edm::ParameterSet& ps)
0088     : g4Label_(ps.getUntrackedParameter<std::string>("moduleLabel", "g4SimHits")),
0089       hitLab_(ps.getParameter<std::vector<std::string> >("hitCollection")),
0090       timeSliceUnit_(ps.getParameter<std::vector<double> >("timeSliceUnit")),
0091       maxEnergy_(ps.getUntrackedParameter<double>("maxEnergy", 250.0)),
0092       maxTime_(ps.getUntrackedParameter<double>("maxTime", 1000.0)),
0093       tMax_(ps.getUntrackedParameter<double>("timeCut", 100.0)),
0094       tScale_(ps.getUntrackedParameter<double>("timeScale", 1.0)),
0095       tCut_(ps.getUntrackedParameter<double>("timeThreshold", 15.0)),
0096       testNumber_(ps.getUntrackedParameter<bool>("testNumbering", false)),
0097       passive_(ps.getUntrackedParameter<bool>("passiveHits", false)),
0098       allSteps_(ps.getUntrackedParameter<int>("allSteps", 100)),
0099       detNames_(ps.getUntrackedParameter<std::vector<std::string> >("detNames")),
0100       tokGeom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0101       toks_calo_{edm::vector_transform(hitLab_,
0102                                        [this](const std::string& name) {
0103                                          return consumes<edm::PCaloHitContainer>(edm::InputTag{g4Label_, name});
0104                                        })},
0105       tok_passive_(consumes<edm::PassiveHitContainer>(edm::InputTag(g4Label_, "AllPassiveHits"))) {
0106   usesResource(TFileService::kSharedResource);
0107 
0108   edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label_ << "   Hits|timeSliceUnit:";
0109   for (unsigned int i = 0; i < hitLab_.size(); i++)
0110     edm::LogVerbatim("HitStudy") << "[" << i << "] " << hitLab_[i] << " " << timeSliceUnit_[i];
0111   edm::LogVerbatim("HitStudy") << "Passive Hits " << passive_ << " from AllPassiveHits";
0112   edm::LogVerbatim("HitStudy") << "maxEnergy: " << maxEnergy_ << " maxTime: " << maxTime_ << " tMax: " << tMax_;
0113   for (unsigned int k = 0; k < detNames_.size(); ++k)
0114     edm::LogVerbatim("HitStudy") << "Detector[" << k << "] " << detNames_[k];
0115 
0116   edm::Service<TFileService> tfile;
0117   if (!tfile.isAvailable())
0118     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0119                                       << "please add it to config file";
0120   char name[29], title[120];
0121   std::string dets[nCalo_] = {"EB", "EE", "HB", "HE", "HO", "HF"};
0122   for (int i = 0; i < nCalo_; i++) {
0123     sprintf(name, "Hit%d", i);
0124     sprintf(title, "Number of hits in %s", dets[i].c_str());
0125     h_hit_[i] = tfile->make<TH1F>(name, title, 100, 0., 20000.);
0126     h_hit_[i]->GetXaxis()->SetTitle(title);
0127     h_hit_[i]->GetYaxis()->SetTitle("Events");
0128     sprintf(name, "Time%d", i);
0129     sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
0130     h_time_[i] = tfile->make<TH1F>(name, title, 100, 0., 200.);
0131     h_time_[i]->GetXaxis()->SetTitle(title);
0132     h_time_[i]->GetYaxis()->SetTitle("Hits");
0133     sprintf(name, "TimeT%d", i);
0134     sprintf(title, "Time of each hit (ns) in %s", dets[i].c_str());
0135     h_timeT_[i] = tfile->make<TH1F>(name, title, 100, 0., 200.);
0136     h_timeT_[i]->GetXaxis()->SetTitle(title);
0137     h_timeT_[i]->GetYaxis()->SetTitle("Hits");
0138     double ymax = (i > 1) ? 0.01 : 0.1;
0139     double ymx0 = (i > 1) ? 0.0025 : 0.025;
0140     sprintf(name, "Edep%d", i);
0141     sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
0142     h_edep_[i] = tfile->make<TH1F>(name, title, 100, 0., ymax);
0143     h_edep_[i]->GetXaxis()->SetTitle(title);
0144     h_edep_[i]->GetYaxis()->SetTitle("Hits");
0145     sprintf(name, "EdepT%d", i);
0146     sprintf(title, "Energy deposit (GeV) of each hit in %s", dets[i].c_str());
0147     h_edepT_[i] = tfile->make<TH1F>(name, title, 100, 0., ymx0);
0148     h_edepT_[i]->GetXaxis()->SetTitle(title);
0149     h_edepT_[i]->GetYaxis()->SetTitle("Hits");
0150     sprintf(name, "EdepEM%d", i);
0151     sprintf(title, "Energy deposit (GeV) by EM particles in %s", dets[i].c_str());
0152     h_edepEM_[i] = tfile->make<TH1F>(name, title, 100, 0., ymx0);
0153     h_edepEM_[i]->GetXaxis()->SetTitle(title);
0154     h_edepEM_[i]->GetYaxis()->SetTitle("Hits");
0155     sprintf(name, "EdepHad%d", i);
0156     sprintf(title, "Energy deposit (GeV) by hadrons in %s", dets[i].c_str());
0157     h_edepHad_[i] = tfile->make<TH1F>(name, title, 100, 0., ymx0);
0158     h_edepHad_[i]->GetXaxis()->SetTitle(title);
0159     h_edepHad_[i]->GetYaxis()->SetTitle("Hits");
0160     sprintf(name, "Edep15%d", i);
0161     sprintf(title, "Energy deposit (GeV) for T > %4.0f ns in %s", tCut_, dets[i].c_str());
0162     h_edep1_[i] = tfile->make<TH1F>(name, title, 100, 0., ymax);
0163     h_edep1_[i]->GetXaxis()->SetTitle(title);
0164     h_edep1_[i]->GetYaxis()->SetTitle("Hits");
0165     sprintf(name, "EdepT15%d", i);
0166     sprintf(title, "Energy deposit (GeV) of each hit for T > %4.0f in %s", tCut_, dets[i].c_str());
0167     h_edepT1_[i] = tfile->make<TH1F>(name, title, 100, 0., ymx0);
0168     h_edepT1_[i]->GetXaxis()->SetTitle(title);
0169     h_edepT1_[i]->GetYaxis()->SetTitle("Hits");
0170     ymax = (i > 1) ? 1.0 : maxEnergy_;
0171     sprintf(name, "Etot%d", i);
0172     sprintf(title, "Total energy deposit (GeV) in %s", dets[i].c_str());
0173     h_etot_[i] = tfile->make<TH1F>(name, title, 50, 0., ymax);
0174     h_etot_[i]->GetXaxis()->SetTitle(title);
0175     h_etot_[i]->GetYaxis()->SetTitle("Events");
0176     sprintf(name, "EtotG%d", i);
0177     sprintf(title, "Total energy deposit (GeV) in %s (t < 100 ns)", dets[i].c_str());
0178     h_etotg_[i] = tfile->make<TH1F>(name, title, 50, 0., ymax);
0179     h_etotg_[i]->GetXaxis()->SetTitle(title);
0180     h_etotg_[i]->GetYaxis()->SetTitle("Events");
0181     sprintf(name, "rr%d", i);
0182     sprintf(title, "R of hit point (cm) in %s", dets[i].c_str());
0183     h_rr_[i] = tfile->make<TH1F>(name, title, 100, 0., 250.);
0184     h_rr_[i]->GetXaxis()->SetTitle(title);
0185     h_rr_[i]->GetYaxis()->SetTitle("Hits");
0186     sprintf(name, "zz%d", i);
0187     sprintf(title, "z of hit point (cm) in %s", dets[i].c_str());
0188     h_zz_[i] = tfile->make<TH1F>(name, title, 240, -600., 600.);
0189     h_zz_[i]->GetXaxis()->SetTitle(title);
0190     h_zz_[i]->GetYaxis()->SetTitle("Hits");
0191     sprintf(name, "eta%d", i);
0192     sprintf(title, "#eta of hit point in %s", dets[i].c_str());
0193     h_eta_[i] = tfile->make<TH1F>(name, title, 100, -5.0, 5.0);
0194     h_eta_[i]->GetXaxis()->SetTitle(title);
0195     h_eta_[i]->GetYaxis()->SetTitle("Hits");
0196     sprintf(name, "phi%d", i);
0197     sprintf(title, "#phi of hit point in %s", dets[i].c_str());
0198     h_phi_[i] = tfile->make<TH1F>(name, title, 100, -M_PI, M_PI);
0199     h_phi_[i]->GetXaxis()->SetTitle(title);
0200     h_phi_[i]->GetYaxis()->SetTitle("Hits");
0201   }
0202   sprintf(title, "R vs Z of hit point");
0203   h_rz_ = tfile->make<TH2F>("rz", title, 120, 0., 600., 100, 0., 250.);
0204   h_rz_->GetXaxis()->SetTitle("z (cm)");
0205   h_rz_->GetYaxis()->SetTitle("R (cm)");
0206   sprintf(title, "R vs Z of hit point for hits with T > %4.0f ns", tCut_);
0207   h_rz1_ = tfile->make<TH2F>("rz2", title, 120, 0., 600., 100, 0., 250.);
0208   h_rz1_->GetXaxis()->SetTitle("z (cm)");
0209   h_rz1_->GetYaxis()->SetTitle("R (cm)");
0210   sprintf(title, "#phi vs #eta of hit point");
0211   h_etaphi_ = tfile->make<TH2F>("etaphi", title, 100, 0., 5., 100, 0., M_PI);
0212   h_etaphi_->GetXaxis()->SetTitle("#eta");
0213   h_etaphi_->GetYaxis()->SetTitle("#phi");
0214 
0215   if (passive_) {
0216     h_hitp_ = tfile->make<TH1F>("hitp", "All Steps", 100, 0.0, 20000.0);
0217     h_hitp_->GetXaxis()->SetTitle("Hits");
0218     h_hitp_->GetYaxis()->SetTitle("Events");
0219     h_trackp_ = tfile->make<TH1F>("trackp", "All Steps", 100, 0.0, 200000.0);
0220     h_hitp_->GetXaxis()->SetTitle("Tracks");
0221     h_hitp_->GetYaxis()->SetTitle("Events");
0222     h_edepp_ = tfile->make<TH1F>("edepp", "All Steps", 100, 0.0, 50.0);
0223     h_edepp_->GetXaxis()->SetTitle("Energy Deposit (MeV)");
0224     h_edepp_->GetYaxis()->SetTitle("Hits");
0225     h_timep_ = tfile->make<TH1F>("timep", "All Steps", 100, 0.0, 100.0);
0226     h_timep_->GetXaxis()->SetTitle("Hits");
0227     h_timep_->GetYaxis()->SetTitle("Hit Time (ns)");
0228     h_stepp_ = tfile->make<TH1F>("stepp", "All Steps", 1000, 0.0, 100.0);
0229     h_stepp_->GetXaxis()->SetTitle("Hits");
0230     h_stepp_->GetYaxis()->SetTitle("Step length (cm)");
0231     for (unsigned int k = 0; k < detNames_.size(); ++k) {
0232       sprintf(name, "edept%d", k);
0233       sprintf(title, "Energy Deposit (MeV) in %s", detNames_[k].c_str());
0234       h_edepTk_.emplace_back(tfile->make<TH1F>(name, title, 100, 0.0, 1.0));
0235       h_edepTk_.back()->GetYaxis()->SetTitle("Hits");
0236       h_edepTk_.back()->GetXaxis()->SetTitle(title);
0237       sprintf(name, "timet%d", k);
0238       sprintf(title, "Hit Time (ns) in %s", detNames_[k].c_str());
0239       h_timeTk_.emplace_back(tfile->make<TH1F>(name, title, 100, 0.0, 100.0));
0240       h_timeTk_.back()->GetYaxis()->SetTitle("Hits");
0241       h_timeTk_.back()->GetXaxis()->SetTitle(title);
0242     }
0243     if ((allSteps_ / 100) % 10 > 0) {
0244       for (int eta = 1; eta <= 29; ++eta) {
0245         int dmax = (eta < 16) ? 4 : 7;
0246         for (int depth = 1; depth <= dmax; ++depth) {
0247           sprintf(name, "Eta%dDepth%d", eta, depth);
0248           sprintf(title, "R vs Z (#eta = %d, depth = %d)", eta, depth);
0249           etaDepth_[eta * 100 + depth] = h_rzH_.size();
0250           h_rzH_.emplace_back(tfile->make<TH2F>(name, title, 120, 0., 600., 100, 0., 250.));
0251           h_rzH_.back()->GetXaxis()->SetTitle("z (cm)");
0252           h_rzH_.back()->GetYaxis()->SetTitle("R (cm)");
0253         }
0254       }
0255     }
0256   }
0257 }
0258 
0259 void CaloSimHitAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0260   edm::ParameterSetDescription desc;
0261   std::vector<std::string> labels = {"EcalHitsEB1", "EcalHitsEE1", "HcalHits1"};
0262   std::vector<double> times = {1, 1, 1};
0263   desc.addUntracked<std::string>("moduleLabel", "g4SimHits");
0264   desc.add<std::vector<std::string> >("hitCollection", labels);
0265   desc.add<std::vector<double> >("timeSliceUnit", times);
0266   desc.addUntracked<double>("maxEnergy", 250.0);
0267   desc.addUntracked<double>("maxTime", 1000.0);
0268   desc.addUntracked<double>("timeCut", 100.0);
0269   desc.addUntracked<double>("timeScale", 1.0);
0270   desc.addUntracked<double>("timeThreshold", 15.0);
0271   desc.addUntracked<bool>("testNumbering", false);
0272   desc.addUntracked<bool>("passiveHits", false);
0273   std::vector<std::string> names = {"PixelBarrel", "PixelForward", "TIB", "TID", "TOB", "TEC"};
0274   desc.addUntracked<std::vector<std::string> >("detNames", names);
0275   desc.addUntracked<int>("allStep", 100);
0276   descriptions.add("caloSimHitAnalysis", desc);
0277 }
0278 
0279 void CaloSimHitAnalysis::analyze(edm::Event const& e, edm::EventSetup const& set) {
0280   edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis:Run = " << e.id().run() << " Event = " << e.id().event();
0281 
0282   caloGeometry_ = &set.getData(tokGeom_);
0283   hcalGeom_ = static_cast<const HcalGeometry*>(caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0284 
0285   for (unsigned int i = 0; i < toks_calo_.size(); i++) {
0286     const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(toks_calo_[i]);
0287     bool getHits = (hitsCalo.isValid());
0288 #ifdef EDM_ML_DEBUG
0289     edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis: Input flags Hits[" << i << "]: " << getHits;
0290 #endif
0291     if (getHits) {
0292       std::vector<PCaloHit> caloHits;
0293       caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
0294 #ifdef EDM_ML_DEBUG
0295       edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis: Hit buffer [" << i << "] " << caloHits.size();
0296 #endif
0297       analyzeHits(caloHits, i);
0298     }
0299   }
0300 
0301   if (passive_) {
0302     const edm::Handle<edm::PassiveHitContainer>& hitsPassive = e.getHandle(tok_passive_);
0303     bool getHits = (hitsPassive.isValid());
0304 #ifdef EDM_ML_DEBUG
0305     edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis: Passive: " << getHits;
0306 #endif
0307     if (getHits) {
0308       std::vector<PassiveHit> passiveHits;
0309       passiveHits.insert(passiveHits.end(), hitsPassive->begin(), hitsPassive->end());
0310 #ifdef EDM_ML_DEBUG
0311       edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis: Passive Hit buffer  " << passiveHits.size();
0312 #endif
0313       analyzePassiveHits(passiveHits);
0314     }
0315   }
0316 }
0317 
0318 void CaloSimHitAnalysis::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
0319   int nHit = hits.size();
0320   int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nEB = 0, nEE = 0, nBad = 0, iHit = 0;
0321   std::map<CaloHitID, std::pair<double, double> > hitMap;
0322   double etot[nCalo_], etotG[nCalo_];
0323   for (unsigned int k = 0; k < nCalo_; ++k) {
0324     etot[k] = etotG[k] = 0;
0325   }
0326   for (const auto& hit : hits) {
0327     double edep = hit.energy();
0328     double time = tScale_ * hit.time();
0329     uint32_t id = hit.id();
0330     int itra = hit.geantTrackId();
0331     double edepEM = hit.energyEM();
0332     double edepHad = hit.energyHad();
0333     int idx(-1);
0334     if (indx != 2) {
0335       idx = indx;
0336       if (indx == 0)
0337         ++nEB;
0338       else
0339         ++nEE;
0340     } else {
0341       int subdet(0);
0342       if (testNumber_) {
0343         int ieta(0), phi(0), z(0), lay(0), depth(0);
0344         HcalTestNumbering::unpackHcalIndex(id, subdet, z, depth, ieta, phi, lay);
0345         id = HcalDetId(static_cast<HcalSubdetector>(subdet), z * ieta, phi, depth).rawId();
0346       } else {
0347         subdet = HcalDetId(id).subdet();
0348       }
0349       if (subdet == static_cast<int>(HcalBarrel)) {
0350         idx = indx;
0351         nHB++;
0352       } else if (subdet == static_cast<int>(HcalEndcap)) {
0353         idx = indx + 1;
0354         nHE++;
0355       } else if (subdet == static_cast<int>(HcalOuter)) {
0356         idx = indx + 2;
0357         nHO++;
0358       } else if (subdet == static_cast<int>(HcalForward)) {
0359         idx = indx + 3;
0360         nHF++;
0361       }
0362     }
0363 #ifdef EDM_ML_DEBUG
0364     edm::LogVerbatim("HitStudy") << "Hit[" << iHit << ":" << nHit << ":" << idx << "] E " << edep << ":" << edepEM
0365                                  << ":" << edepHad << " T " << time << " itra " << itra << " ID " << std::hex << id
0366                                  << std::dec;
0367 #endif
0368     ++iHit;
0369     if (idx >= 0) {
0370       CaloHitID hid(id, time, itra, 0, timeSliceUnit_[indx]);
0371       auto itr = hitMap.find(hid);
0372       if (itr == hitMap.end())
0373         hitMap[hid] = std::make_pair(time, edep);
0374       else
0375         ((itr->second).second) += edep;
0376       h_edepT_[idx]->Fill(edep);
0377       h_timeT_[idx]->Fill(time);
0378       if (edepEM > 0)
0379         h_edepEM_[idx]->Fill(edepEM);
0380       if (edepHad > 0)
0381         h_edepHad_[idx]->Fill(edepHad);
0382       if (time > tCut_)
0383         h_edepT1_[idx]->Fill(edep);
0384     } else {
0385       ++nBad;
0386     }
0387   }
0388 
0389   //Now make plots
0390   for (auto itr = hitMap.begin(); itr != hitMap.end(); ++itr) {
0391     int idx = -1;
0392     GlobalPoint point;
0393     DetId id((itr->first).unitID());
0394 #ifdef EDM_ML_DEBUG
0395     edm::LogVerbatim("HitStudy") << "Index " << indx << " Geom " << caloGeometry_ << ":" << hcalGeom_ << "  "
0396                                  << std::hex << id.rawId() << std::dec;
0397 #endif
0398     if (indx != 2) {
0399       idx = indx;
0400       point = caloGeometry_->getPosition(id);
0401     } else {
0402       int subdet = id.subdetId();
0403       if (subdet == static_cast<int>(HcalBarrel)) {
0404         idx = indx;
0405       } else if (subdet == static_cast<int>(HcalEndcap)) {
0406         idx = indx + 1;
0407       } else if (subdet == static_cast<int>(HcalOuter)) {
0408         idx = indx + 2;
0409       } else if (subdet == static_cast<int>(HcalForward)) {
0410         idx = indx + 3;
0411       }
0412       point = hcalGeom_->getPosition(id);
0413     }
0414     double edep = (itr->second).second;
0415     double time = (itr->second).first;
0416 #ifdef EDM_ML_DEBUG
0417     edm::LogVerbatim("HitStudy") << "Index " << idx << ":" << nCalo_ << " Point " << point << " E " << edep << " T "
0418                                  << time;
0419 #endif
0420     if (idx >= 0) {
0421       h_time_[idx]->Fill(time);
0422       h_edep_[idx]->Fill(edep);
0423       h_rr_[idx]->Fill(point.perp());
0424       h_zz_[idx]->Fill(point.z());
0425       h_eta_[idx]->Fill(point.eta());
0426       h_phi_[idx]->Fill(point.phi());
0427       h_rz_->Fill(std::abs(point.z()), point.perp());
0428       h_etaphi_->Fill(std::abs(point.eta()), std::abs(point.phi()));
0429       etot[idx] += edep;
0430       if (time < tMax_)
0431         etotG[idx] += edep;
0432       if (time > tCut_) {
0433         h_edep1_[idx]->Fill(edep);
0434         h_rz1_->Fill(std::abs(point.z()), point.perp());
0435       }
0436     }
0437   }
0438 
0439   if (indx < 2) {
0440     h_etot_[indx]->Fill(etot[indx]);
0441     h_etotg_[indx]->Fill(etotG[indx]);
0442     if (indx == 0)
0443       h_hit_[indx]->Fill(double(nEB));
0444     else
0445       h_hit_[indx]->Fill(double(nEE));
0446   } else {
0447     h_hit_[2]->Fill(double(nHB));
0448     h_hit_[3]->Fill(double(nHE));
0449     h_hit_[4]->Fill(double(nHO));
0450     h_hit_[5]->Fill(double(nHF));
0451     for (int idx = 2; idx < nCalo_; idx++) {
0452       h_etot_[idx]->Fill(etot[idx]);
0453       h_etotg_[idx]->Fill(etotG[idx]);
0454     }
0455   }
0456 
0457   edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis::analyzeHits: EB " << nEB << " EE " << nEE << " HB " << nHB
0458                                << " HE " << nHE << " HO " << nHO << " HF " << nHF << " Bad " << nBad << " All " << nHit
0459                                << " Reduced " << hitMap.size();
0460 }
0461 
0462 void CaloSimHitAnalysis::analyzePassiveHits(std::vector<PassiveHit>& hits) {
0463   const std::string active = "Active";
0464   const std::string sensor = "Sensor";
0465   std::map<std::pair<std::string, uint32_t>, int> hitx;
0466   std::map<int, int> tracks;
0467   unsigned int passive1(0), passive2(0);
0468   for (auto& hit : hits) {
0469     std::string name = hit.vname();
0470     std::pair<std::string, uint32_t> volume = std::make_pair(name, (hit.id() % 1000000));
0471     auto itr = hitx.find(volume);
0472     if (itr == hitx.end())
0473       hitx[volume] = 1;
0474     else
0475       ++(itr->second);
0476     auto ktr = tracks.find(hit.trackId());
0477     if (ktr == tracks.end())
0478       tracks[hit.trackId()] = 1;
0479     else
0480       ++(ktr->second);
0481     h_edepp_->Fill(hit.energy());
0482     h_timep_->Fill(hit.time());
0483     h_stepp_->Fill(hit.stepLength());
0484     if ((name.find(active) != std::string::npos) || (name.find(sensor) != std::string::npos)) {
0485       unsigned idet = detNames_.size();
0486       for (unsigned int k = 0; k < detNames_.size(); ++k) {
0487         if (name.find(detNames_[k]) != std::string::npos) {
0488           idet = k;
0489           break;
0490         }
0491       }
0492       if (idet < detNames_.size()) {
0493         h_edepTk_[idet]->Fill(hit.energy());
0494         h_timeTk_[idet]->Fill(hit.time());
0495       }
0496     }
0497 
0498     if ((allSteps_ / 100) % 10 > 0) {
0499       uint32_t id = hit.id();
0500       if (DetId(id).det() == DetId::Hcal) {
0501         HcalDetId hid = HcalDetId(id);
0502         int indx = (100 * hid.ietaAbs() + hid.depth());
0503         auto itr = etaDepth_.find(indx);
0504 #ifdef EDM_ML_DEBUG
0505         edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis::ID: " << hid << " Index " << indx << " Iterator "
0506                                      << (itr != etaDepth_.end());
0507 #endif
0508         ++passive1;
0509         if (itr != etaDepth_.end()) {
0510           uint32_t ipos = itr->second;
0511           double rr = std::sqrt(hit.x() * hit.x() + hit.y() * hit.y());
0512           if (ipos < h_rzH_.size()) {
0513             h_rzH_[ipos]->Fill(hit.z(), rr);
0514             ++passive2;
0515           }
0516         }
0517       }
0518     }
0519   }
0520   h_hitp_->Fill(hitx.size());
0521   h_trackp_->Fill(tracks.size());
0522   edm::LogVerbatim("HitStudy") << "CaloSimHitAnalysis::analyzPassiveHits: Total " << hits.size() << " Cells "
0523                                << hitx.size() << " Tracks " << tracks.size() << " Passive " << passive1 << ":"
0524                                << passive2;
0525 }
0526 
0527 //define this as a plug-in
0528 DEFINE_FWK_MODULE(CaloSimHitAnalysis);