Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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 //define this as a plug-in
0673 DEFINE_FWK_MODULE(CaloSimHitStudy);