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