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