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
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
0629 DEFINE_FWK_MODULE(CaloSimHitStudy);
