File indexing completed on 2024-05-10 02:21:14
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/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.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/HcalDetId/interface/HcalDetId.h"
0019 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0020 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0021 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0022
0023 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0024 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0025
0026 #include <CLHEP/Units/PhysicalConstants.h>
0027 #include <CLHEP/Units/SystemOfUnits.h>
0028
0029 #include <TH1F.h>
0030 #include <TProfile.h>
0031 #include <TProfile2D.h>
0032
0033 #include <fstream>
0034 #include <iostream>
0035 #include <iomanip>
0036 #include <memory>
0037 #include <map>
0038 #include <string>
0039 #include <vector>
0040
0041 class HOSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0042 public:
0043 HOSimHitStudy(const edm::ParameterSet &ps);
0044 ~HOSimHitStudy() override = default;
0045 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0046
0047 protected:
0048 void beginJob() override;
0049 void beginRun(edm::Run const &, edm::EventSetup const &) override {}
0050 void endJob() override {}
0051 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0052 void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0053
0054 void analyzeHits();
0055
0056 private:
0057 const std::string g4Label_;
0058 const std::vector<std::string> hitLab_;
0059 const double maxEnergy_, scaleEB_, scaleHB_, scaleHO_;
0060 const double tcut_;
0061 const bool scheme_, print_;
0062 const edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
0063 const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>> toks_calo_;
0064 std::vector<PCaloHit> ecalHits, hcalHits;
0065 TH1F *hit_[3], *time_[3], *edepTW_[3], *edepTWT_[3];
0066 TH1F *edep_[3], *hitTow_[3], *eneInc_, *etaInc_, *phiInc_;
0067 TH1F *edepT_[3], *eEB_, *eEBHB_, *eEBHBHO_, *eEBHBHOT_;
0068 TH1F *edepZon_[3], *edepZonT_[3], *eEBT_, *eEBHBT_;
0069 TH1F *eHOE17_[15], *eHOE18_[15], *eHOE_[15];
0070 TH1F *eHOE17T_[15], *eHOE18T_[15], *eHOET_[15];
0071 TH1F *eHOEta17_[15], *eHOEta18_[15], *eHOEta_[15];
0072 TH1F *eHOEta17T_[15], *eHOEta18T_[15], *eHOEtaT_[15];
0073 TH1F *nHOE1_[15], *nHOE1T_[15], *nHOEta1_[15], *nHOEta1T_[15];
0074 TProfile *eHO1_, *eHO1T_, *eHO17_, *eHO17T_, *eHO18_, *eHO18T_;
0075 TProfile *nHO1_, *nHO1T_, *nHOE2_[15], *nHOE2T_[15];
0076 TProfile *nHOEta2_[15], *nHOEta2T_[15];
0077 TProfile2D *eHO2_, *eHO2T_, *nHO2_, *nHO2T_;
0078 double eInc, etaInc, phiInc;
0079 };
0080
0081 HOSimHitStudy::HOSimHitStudy(const edm::ParameterSet &ps)
0082 : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0083 hitLab_(ps.getParameter<std::vector<std::string>>("HitCollection")),
0084 maxEnergy_(ps.getUntrackedParameter<double>("MaxEnergy", 200.0)),
0085 scaleEB_(ps.getUntrackedParameter<double>("ScaleEB", 1.0)),
0086 scaleHB_(ps.getUntrackedParameter<double>("ScaleHB", 100.0)),
0087 scaleHO_(ps.getUntrackedParameter<double>("ScaleHO", 2.0)),
0088 tcut_(ps.getUntrackedParameter<double>("TimeCut", 100.0)),
0089 scheme_(ps.getUntrackedParameter<bool>("TestNumbering", false)),
0090 print_(ps.getUntrackedParameter<bool>("PrintExcessEnergy", true)),
0091 tok_evt_(consumes<edm::HepMCProduct>(
0092 edm::InputTag(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")))),
0093 toks_calo_{edm::vector_transform(hitLab_, [this](const std::string &name) {
0094 return consumes<edm::PCaloHitContainer>(edm::InputTag{g4Label_, name});
0095 })} {
0096 usesResource(TFileService::kSharedResource);
0097
0098 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Module Label: " << g4Label_ << " Hits: " << hitLab_[0] << ", "
0099 << hitLab_[1] << " MaxEnergy: " << maxEnergy_ << " Scale factor for EB " << scaleEB_
0100 << ", for HB " << scaleHB_ << " and for HO " << scaleHO_ << " time Cut " << tcut_;
0101 }
0102
0103 void HOSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0104 edm::ParameterSetDescription desc;
0105 std::vector<std::string> labels = {"EcalHitsEB", "HcalHits"};
0106 desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
0107 desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0108 desc.addUntracked<std::vector<std::string>>("HitCollection", labels);
0109 desc.addUntracked<double>("MaxEnergy", 50.0);
0110 desc.addUntracked<double>("ScaleEB", 1.02);
0111 desc.addUntracked<double>("ScaleHB", 104.4);
0112 desc.addUntracked<double>("ScaleHO", 2.33);
0113 desc.addUntracked<double>("TimeCut", 100.0);
0114 desc.addUntracked<bool>("PrintExcessEnergy", true);
0115 desc.addUntracked<bool>("TestNumbering", false);
0116 descriptions.add("HOSimHitStudy", desc);
0117 }
0118
0119 void HOSimHitStudy::beginJob() {
0120 edm::Service<TFileService> tfile;
0121
0122 if (!tfile.isAvailable())
0123 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0124 << "please add it to config file";
0125 std::string dets[3] = {"EB", "HB", "HO"};
0126 char name[60], title[100];
0127 double ymax = maxEnergy_;
0128 sprintf(title, "Incident Energy (GeV)");
0129 eneInc_ = tfile->make<TH1F>("EneInc", title, 1000, 0., ymax);
0130 eneInc_->GetXaxis()->SetTitle(title);
0131 eneInc_->GetYaxis()->SetTitle("Events");
0132 sprintf(title, "Incident #eta");
0133 etaInc_ = tfile->make<TH1F>("EtaInc", title, 200, -5., 5.);
0134 etaInc_->GetXaxis()->SetTitle(title);
0135 etaInc_->GetYaxis()->SetTitle("Events");
0136 sprintf(title, "Incident #phi");
0137 phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
0138 phiInc_->GetXaxis()->SetTitle(title);
0139 phiInc_->GetYaxis()->SetTitle("Events");
0140 int itcut = (int)(tcut_);
0141 for (int i = 0; i < 3; i++) {
0142 sprintf(name, "Hit%d", i);
0143 sprintf(title, "Number of hits in %s", dets[i].c_str());
0144 hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
0145 hit_[i]->GetXaxis()->SetTitle(title);
0146 hit_[i]->GetYaxis()->SetTitle("Events");
0147 sprintf(name, "Time%d", i);
0148 sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
0149 time_[i] = tfile->make<TH1F>(name, title, 1200, 0., 1200.);
0150 time_[i]->GetXaxis()->SetTitle(title);
0151 time_[i]->GetYaxis()->SetTitle("Hits");
0152 if (i > 0)
0153 ymax = 1.0;
0154 else
0155 ymax = 50.0;
0156 sprintf(name, "Edep%d", i);
0157 sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
0158 edep_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0159 edep_[i]->GetXaxis()->SetTitle(title);
0160 edep_[i]->GetYaxis()->SetTitle("Hits");
0161 sprintf(name, "EdepT%d", i);
0162 sprintf(title, "Energy deposit (GeV) in %s for t < %d ns", dets[i].c_str(), itcut);
0163 edepT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0164 edepT_[i]->GetXaxis()->SetTitle(title);
0165 edepT_[i]->GetYaxis()->SetTitle("Hits");
0166 sprintf(name, "HitTow%d", i);
0167 sprintf(title, "Number of towers with hits in %s", dets[i].c_str());
0168 hitTow_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
0169 hitTow_[i]->GetXaxis()->SetTitle(title);
0170 hitTow_[i]->GetYaxis()->SetTitle("Events");
0171 if (i > 0)
0172 ymax = 0.05 * maxEnergy_;
0173 else
0174 ymax = maxEnergy_;
0175 sprintf(name, "EdepTW%d", i);
0176 sprintf(title, "Energy deposit (GeV) in %s Tower", dets[i].c_str());
0177 edepTW_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0178 edepTW_[i]->GetXaxis()->SetTitle(title);
0179 edepTW_[i]->GetYaxis()->SetTitle("Towers");
0180 sprintf(name, "EdepTWT%d", i);
0181 sprintf(title, "Energy deposit (GeV) in %s Tower for t < %d ns", dets[i].c_str(), itcut);
0182 edepTWT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0183 edepTWT_[i]->GetXaxis()->SetTitle(title);
0184 edepTWT_[i]->GetYaxis()->SetTitle("Towers");
0185 sprintf(name, "EdepZone%d", i);
0186 sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
0187 edepZon_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0188 edepZon_[i]->GetXaxis()->SetTitle(title);
0189 edepZon_[i]->GetYaxis()->SetTitle("Events");
0190 sprintf(name, "EdepZoneT%d", i);
0191 sprintf(title, "Energy deposit (GeV) in %s for t < %d ns", dets[i].c_str(), itcut);
0192 edepZonT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
0193 edepZonT_[i]->GetXaxis()->SetTitle(title);
0194 edepZonT_[i]->GetYaxis()->SetTitle("Events");
0195 }
0196 sprintf(title, "Energy Measured in EB (GeV)");
0197 eEB_ = tfile->make<TH1F>("EEB", title, 5000, 0., maxEnergy_);
0198 eEB_->GetXaxis()->SetTitle(title);
0199 eEB_->GetYaxis()->SetTitle("Events");
0200 sprintf(title, "Energy Measured in EB+HB (GeV)");
0201 eEBHB_ = tfile->make<TH1F>("EEBHB", title, 5000, 0., maxEnergy_);
0202 eEBHB_->GetXaxis()->SetTitle(title);
0203 eEBHB_->GetYaxis()->SetTitle("Events");
0204 sprintf(title, "Energy Measured in EB+HB+HO (GeV)");
0205 eEBHBHO_ = tfile->make<TH1F>("EEBHBHO", title, 5000, 0., maxEnergy_);
0206 eEBHBHO_->GetXaxis()->SetTitle(title);
0207 eEBHBHO_->GetYaxis()->SetTitle("Events");
0208 sprintf(title, "Energy Measured in EB (GeV) for t < %d ns", itcut);
0209 eEBT_ = tfile->make<TH1F>("EEBT", title, 5000, 0., maxEnergy_);
0210 eEBT_->GetXaxis()->SetTitle(title);
0211 eEBT_->GetYaxis()->SetTitle("Events");
0212 sprintf(title, "Energy Measured in EB+HB (GeV) for t < %d ns", itcut);
0213 eEBHBT_ = tfile->make<TH1F>("EEBHBT", title, 5000, 0., maxEnergy_);
0214 eEBHBT_->GetXaxis()->SetTitle(title);
0215 eEBHBT_->GetYaxis()->SetTitle("Events");
0216 sprintf(title, "Energy Measured in EB+HB+HO (GeV) for t < %d ns", itcut);
0217 eEBHBHOT_ = tfile->make<TH1F>("EEBHBHOT", title, 5000, 0., maxEnergy_);
0218 eEBHBHOT_->GetXaxis()->SetTitle(title);
0219 eEBHBHOT_->GetYaxis()->SetTitle("Events");
0220 sprintf(title, "SimHit energy in HO");
0221 eHO1_ = tfile->make<TProfile>("EHO1", title, 30, -1.305, 1.305);
0222 eHO1_->GetXaxis()->SetTitle(title);
0223 eHO1_->GetYaxis()->SetTitle("Events");
0224 eHO2_ = tfile->make<TProfile2D>("EHO2", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
0225 eHO2_->GetXaxis()->SetTitle(title);
0226 eHO2_->GetYaxis()->SetTitle("Events");
0227 sprintf(title, "SimHit energy in HO Layer 17");
0228 eHO17_ = tfile->make<TProfile>("EHO17", title, 30, -1.305, 1.305);
0229 eHO17_->GetXaxis()->SetTitle(title);
0230 eHO17_->GetYaxis()->SetTitle("Events");
0231 sprintf(title, "SimHit energy in HO Layer 18");
0232 eHO18_ = tfile->make<TProfile>("EHO18", title, 30, -1.305, 1.305);
0233 eHO18_->GetXaxis()->SetTitle(title);
0234 eHO18_->GetYaxis()->SetTitle("Events");
0235 sprintf(title, "SimHit energy in HO for t < %d ns", itcut);
0236 eHO1T_ = tfile->make<TProfile>("EHO1T", title, 30, -1.305, 1.305);
0237 eHO1T_->GetXaxis()->SetTitle(title);
0238 eHO1T_->GetYaxis()->SetTitle("Events");
0239 eHO2T_ = tfile->make<TProfile2D>("EHO2T", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
0240 eHO2T_->GetXaxis()->SetTitle(title);
0241 eHO2T_->GetYaxis()->SetTitle("Events");
0242 sprintf(title, "SimHit energy in HO Layer 17 for t < %d ns", itcut);
0243 eHO17T_ = tfile->make<TProfile>("EHO17T", title, 30, -1.305, 1.305);
0244 eHO17T_->GetXaxis()->SetTitle(title);
0245 eHO17T_->GetYaxis()->SetTitle("Events");
0246 sprintf(title, "SimHit energy in HO Layer 18 for t < %d ns", itcut);
0247 eHO18T_ = tfile->make<TProfile>("EHO18T", title, 30, -1.305, 1.305);
0248 eHO18T_->GetXaxis()->SetTitle(title);
0249 eHO18T_->GetYaxis()->SetTitle("Events");
0250 sprintf(title, "Number of layers hit in HO");
0251 nHO1_ = tfile->make<TProfile>("NHO1", title, 30, -1.305, 1.305);
0252 nHO1_->GetXaxis()->SetTitle(title);
0253 nHO1_->GetYaxis()->SetTitle("Events");
0254 nHO2_ = tfile->make<TProfile2D>("NHO2", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
0255 nHO2_->GetXaxis()->SetTitle(title);
0256 nHO2_->GetYaxis()->SetTitle("Events");
0257 sprintf(title, "Number of layers hit in HO for t < %d ns", itcut);
0258 nHO1T_ = tfile->make<TProfile>("NHO1T", title, 30, -1.305, 1.305);
0259 nHO1T_->GetXaxis()->SetTitle(title);
0260 nHO1T_->GetYaxis()->SetTitle("Events");
0261 nHO2T_ = tfile->make<TProfile2D>("NHO2T", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
0262 nHO2T_->GetXaxis()->SetTitle(title);
0263 nHO2T_->GetYaxis()->SetTitle("Events");
0264 for (int i = 0; i < 15; i++) {
0265 sprintf(name, "EHOE%d", i + 1);
0266 sprintf(title, "SimHit energy in HO (Beam in #eta=%d bin)", i + 1);
0267 eHOE_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0268 eHOE_[i]->GetXaxis()->SetTitle(title);
0269 eHOE_[i]->GetYaxis()->SetTitle("Events");
0270 sprintf(name, "EHOE17%d", i + 1);
0271 sprintf(title, "SimHit energy in Layer 17 (Beam in #eta=%d bin)", i + 1);
0272 eHOE17_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0273 eHOE17_[i]->GetXaxis()->SetTitle(title);
0274 eHOE17_[i]->GetYaxis()->SetTitle("Events");
0275 sprintf(name, "EHOE18%d", i + 1);
0276 sprintf(title, "SimHit energy in Layer 18 (Beam in #eta=%d bin)", i + 1);
0277 eHOE18_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0278 eHOE18_[i]->GetXaxis()->SetTitle(title);
0279 eHOE18_[i]->GetYaxis()->SetTitle("Events");
0280 sprintf(name, "EHOE%dT", i + 1);
0281 sprintf(title, "SimHit energy in HO (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
0282 eHOET_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0283 eHOET_[i]->GetXaxis()->SetTitle(title);
0284 eHOET_[i]->GetYaxis()->SetTitle("Events");
0285 sprintf(name, "EHOE17%dT", i + 1);
0286 sprintf(title, "SimHit energy in Layer 17 (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
0287 eHOE17T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0288 eHOE17T_[i]->GetXaxis()->SetTitle(title);
0289 eHOE17T_[i]->GetYaxis()->SetTitle("Events");
0290 sprintf(name, "EHOE18%dT", i + 1);
0291 sprintf(title, "SimHit energy in Layer 18 (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
0292 eHOE18T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0293 eHOE18T_[i]->GetXaxis()->SetTitle(title);
0294 eHOE18T_[i]->GetYaxis()->SetTitle("Events");
0295 sprintf(name, "EHOEta%d", i + 1);
0296 sprintf(title, "SimHit energy in HO #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
0297 eHOEta_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0298 eHOEta_[i]->GetXaxis()->SetTitle(title);
0299 eHOEta_[i]->GetYaxis()->SetTitle("Events");
0300 sprintf(name, "EHOEta17%d", i + 1);
0301 sprintf(title, "SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
0302 eHOEta17_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0303 eHOEta17_[i]->GetXaxis()->SetTitle(title);
0304 eHOEta17_[i]->GetYaxis()->SetTitle("Events");
0305 sprintf(name, "EHOEta18%d", i + 1);
0306 sprintf(title, "SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
0307 eHOEta18_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0308 eHOEta18_[i]->GetXaxis()->SetTitle(title);
0309 eHOEta18_[i]->GetYaxis()->SetTitle("Events");
0310 sprintf(name, "EHOEta%dT", i + 1);
0311 sprintf(title, "SimHit energy in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
0312 eHOEtaT_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0313 eHOEtaT_[i]->GetXaxis()->SetTitle(title);
0314 eHOEtaT_[i]->GetYaxis()->SetTitle("Events");
0315 sprintf(name, "EHOEta17%dT", i + 1);
0316 sprintf(title, "SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
0317 eHOEta17T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0318 eHOEta17T_[i]->GetXaxis()->SetTitle(title);
0319 eHOEta17T_[i]->GetYaxis()->SetTitle("Events");
0320 sprintf(name, "EHOEta18%dT", i + 1);
0321 sprintf(title, "SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
0322 eHOEta18T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
0323 eHOEta18T_[i]->GetXaxis()->SetTitle(title);
0324 eHOEta18T_[i]->GetYaxis()->SetTitle("Events");
0325 sprintf(name, "NHOE1%d", i + 1);
0326 sprintf(title, "Number of layers hit in HO (Beam in #eta=%d bin)", i + 1);
0327 nHOE1_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
0328 nHOE1_[i]->GetXaxis()->SetTitle(title);
0329 nHOE1_[i]->GetYaxis()->SetTitle("Events");
0330 sprintf(name, "NHOE2%d", i + 1);
0331 nHOE2_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
0332 nHOE2_[i]->GetXaxis()->SetTitle(title);
0333 nHOE2_[i]->GetYaxis()->SetTitle("Events");
0334 sprintf(name, "NHOE1%dT", i + 1);
0335 sprintf(title, "Number of layers hit in HO (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
0336 nHOE1T_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
0337 nHOE1T_[i]->GetXaxis()->SetTitle(title);
0338 nHOE1T_[i]->GetYaxis()->SetTitle("Events");
0339 sprintf(name, "NHOE2%dT", i + 1);
0340 nHOE2T_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
0341 nHOE2T_[i]->GetXaxis()->SetTitle(title);
0342 nHOE2T_[i]->GetYaxis()->SetTitle("Events");
0343 sprintf(name, "NHOEta1%d", i + 1);
0344 sprintf(title, "Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
0345 nHOEta1_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
0346 nHOEta1_[i]->GetXaxis()->SetTitle(title);
0347 nHOEta1_[i]->GetYaxis()->SetTitle("Events");
0348 sprintf(name, "NHOEta2%d", i + 1);
0349 nHOEta2_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
0350 nHOEta2_[i]->GetXaxis()->SetTitle(title);
0351 nHOEta2_[i]->GetYaxis()->SetTitle("Events");
0352 sprintf(name, "NHOEta1%dT", i + 1);
0353 sprintf(title, "Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
0354 nHOEta1T_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
0355 nHOEta1T_[i]->GetXaxis()->SetTitle(title);
0356 nHOEta1T_[i]->GetYaxis()->SetTitle("Events");
0357 sprintf(name, "NHOEta2%dT", i + 1);
0358 nHOEta2T_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
0359 nHOEta2T_[i]->GetXaxis()->SetTitle(title);
0360 nHOEta2T_[i]->GetYaxis()->SetTitle("Events");
0361 }
0362 }
0363
0364 void HOSimHitStudy::analyze(const edm::Event &e, const edm::EventSetup &) {
0365 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Run = " << e.id().run() << " Event = " << e.id().event();
0366
0367 const edm::Handle<edm::HepMCProduct> &EvtHandle = e.getHandle(tok_evt_);
0368 const HepMC::GenEvent *myGenEvent = EvtHandle->GetEvent();
0369
0370 eInc = etaInc = phiInc = 0;
0371 HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
0372 if (p != myGenEvent->particles_end()) {
0373 eInc = (*p)->momentum().e();
0374 etaInc = (*p)->momentum().eta();
0375 phiInc = (*p)->momentum().phi();
0376 }
0377
0378 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Energy = " << eInc << " Eta = " << etaInc
0379 << " Phi = " << phiInc / CLHEP::deg;
0380
0381 for (int i = 0; i < 2; i++) {
0382 bool getHits = false;
0383 if (i == 0)
0384 ecalHits.clear();
0385 else
0386 hcalHits.clear();
0387 const edm::Handle<edm::PCaloHitContainer> &hitsCalo = e.getHandle(toks_calo_[i]);
0388 if (hitsCalo.isValid())
0389 getHits = true;
0390 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Input flag " << hitLab_[i] << " getHits flag " << getHits;
0391
0392 if (getHits) {
0393 unsigned int isiz;
0394 if (i == 0) {
0395 ecalHits.insert(ecalHits.end(), hitsCalo->begin(), hitsCalo->end());
0396 isiz = ecalHits.size();
0397 } else {
0398 hcalHits.insert(hcalHits.end(), hitsCalo->begin(), hitsCalo->end());
0399 isiz = hcalHits.size();
0400 }
0401 edm::LogVerbatim("HitStudy") << "HOSimHitStudy:: Hit buffer for " << hitLab_[i] << " has " << isiz << " hits";
0402 }
0403 }
0404 analyzeHits();
0405 }
0406
0407 void HOSimHitStudy::analyzeHits() {
0408
0409 int nhit[3];
0410 double etot[3], etotT[3];
0411 std::vector<unsigned int> ebID, hbID, hoID;
0412 std::vector<double> ebEtow, hbEtow, hoEtow;
0413 std::vector<double> ebEtowT, hbEtowT, hoEtowT;
0414 for (int k = 0; k < 3; k++) {
0415 nhit[k] = 0;
0416 etot[k] = 0;
0417 etotT[k] = 0;
0418 }
0419 eneInc_->Fill(eInc);
0420 etaInc_->Fill(etaInc);
0421 phiInc_->Fill(phiInc);
0422
0423 double eHO17 = 0, eHO18 = 0, eHO17T = 0, eHO18T = 0;
0424 double eHOE17[15], eHOE18[15], eHOE17T[15], eHOE18T[15];
0425 for (int k = 0; k < 15; k++) {
0426 eHOE17[k] = eHOE18[k] = eHOE17T[k] = eHOE18T[k] = 0;
0427 }
0428
0429 for (int k = 0; k < 2; k++) {
0430 int nHit;
0431 if (k == 0) {
0432 nHit = ecalHits.size();
0433 } else {
0434 nHit = hcalHits.size();
0435 }
0436 for (int i = 0; i < nHit; i++) {
0437 double edep, time;
0438 int indx = -1;
0439 unsigned int id_;
0440 if (k == 0) {
0441 indx = 0;
0442 edep = ecalHits[i].energy();
0443 time = ecalHits[i].time();
0444 id_ = ecalHits[i].id();
0445 } else {
0446 edep = hcalHits[i].energy();
0447 time = hcalHits[i].time();
0448 id_ = hcalHits[i].id();
0449 int subdet, zside, depth, eta, phi, lay;
0450 if (scheme_) {
0451 HcalTestNumberingScheme::unpackHcalIndex(id_, subdet, zside, depth, eta, phi, lay);
0452 } else {
0453 HcalDetId id = HcalDetId(id_);
0454 subdet = id.subdet();
0455 zside = id.zside();
0456 depth = id.depth();
0457 eta = id.ietaAbs();
0458 phi = id.iphi();
0459 lay = -1;
0460 }
0461 edm::LogVerbatim("HitStudy") << "HOSimHitStudy:: Hit " << k << " Subdet:" << subdet << " zside:" << zside
0462 << " depth:" << depth << " layer:" << lay << " eta:" << eta << " phi:" << phi;
0463 if (subdet == static_cast<int>(HcalBarrel))
0464 indx = 1;
0465 else if (subdet == static_cast<int>(HcalOuter)) {
0466 indx = 2;
0467 if (lay == 18) {
0468 eHO17 += edep;
0469 if (time < tcut_)
0470 eHO17T += edep;
0471 if (eta >= 0 && eta < 15) {
0472 eHOE17[eta] += edep;
0473 if (time < tcut_)
0474 eHOE17T[eta] += edep;
0475 }
0476 } else {
0477 eHO18 += edep;
0478 if (time < tcut_)
0479 eHO18T += edep;
0480 if (eta >= 0 && eta < 15) {
0481 eHOE18[eta] += edep;
0482 if (time < tcut_)
0483 eHOE18T[eta] += edep;
0484 }
0485 }
0486 }
0487 }
0488 if (indx >= 0) {
0489 double edepT = edep;
0490 time_[indx]->Fill(time);
0491 edep_[indx]->Fill(edep);
0492 etot[indx] += edep;
0493 if (time < tcut_) {
0494 etotT[indx] += edep;
0495 edepT_[indx]->Fill(edep);
0496 edepT = 0;
0497 }
0498 nhit[indx]++;
0499 if (indx == 0) {
0500 bool ok = false;
0501 for (unsigned int j = 0; j < ebID.size(); j++) {
0502 if (id_ == ebID[j]) {
0503 ebEtow[j] += edep;
0504 ebEtowT[j] += edepT;
0505 ok = true;
0506 break;
0507 }
0508 }
0509 if (!ok) {
0510 ebID.push_back(id_);
0511 ebEtow.push_back(edep);
0512 ebEtowT.push_back(edepT);
0513 }
0514 } else if (indx == 1) {
0515 bool ok = false;
0516 for (unsigned int j = 0; j < hbID.size(); j++) {
0517 if (id_ == hbID[j]) {
0518 hbEtow[j] += edep;
0519 hbEtowT[j] += edepT;
0520 ok = true;
0521 break;
0522 }
0523 }
0524 if (!ok) {
0525 hbID.push_back(id_);
0526 hbEtow.push_back(edep);
0527 hbEtowT.push_back(edepT);
0528 }
0529 } else {
0530 bool ok = false;
0531 for (unsigned int j = 0; j < hoID.size(); j++) {
0532 if (id_ == hoID[j]) {
0533 hoEtow[j] += edep;
0534 hoEtowT[j] += edepT;
0535 ok = true;
0536 break;
0537 }
0538 }
0539 if (!ok) {
0540 hoID.push_back(id_);
0541 hoEtow.push_back(edep);
0542 hoEtowT.push_back(edepT);
0543 }
0544 }
0545 }
0546 }
0547 }
0548
0549
0550 for (int k = 0; k < 3; k++) {
0551 hit_[k]->Fill(double(nhit[k]));
0552 edepZon_[k]->Fill(etot[k]);
0553 edepZonT_[k]->Fill(etotT[k]);
0554 }
0555 hitTow_[0]->Fill(double(ebEtow.size()));
0556 for (unsigned int i = 0; i < ebEtow.size(); i++) {
0557 edepTW_[0]->Fill(ebEtow[i]);
0558 edepTWT_[0]->Fill(ebEtowT[i]);
0559 }
0560 hitTow_[1]->Fill(double(hbEtow.size()));
0561 for (unsigned int i = 0; i < hbEtow.size(); i++) {
0562 edepTW_[1]->Fill(hbEtow[i]);
0563 edepTWT_[1]->Fill(hbEtowT[i]);
0564 }
0565 hitTow_[2]->Fill(double(hoEtow.size()));
0566 for (unsigned int i = 0; i < hoEtow.size(); i++) {
0567 edepTW_[2]->Fill(hoEtow[i]);
0568 edepTWT_[2]->Fill(hoEtowT[i]);
0569 }
0570 double eEB = scaleEB_ * etot[0];
0571 double eEBHB = eEB + scaleHB_ * etot[1];
0572 double eEBHBHO = eEBHB + scaleHB_ * scaleHO_ * etot[2];
0573 eEB_->Fill(eEB);
0574 eEBHB_->Fill(eEBHB);
0575 eEBHBHO_->Fill(eEBHBHO);
0576 double eEBT = scaleEB_ * etotT[0];
0577 double eEBHBT = eEBT + scaleHB_ * etotT[1];
0578 double eEBHBHOT = eEBHBT + scaleHB_ * scaleHO_ * etotT[2];
0579 eEBT_->Fill(eEBT);
0580 eEBHBT_->Fill(eEBHBT);
0581 eEBHBHOT_->Fill(eEBHBHOT);
0582 eHO1_->Fill(etaInc, eHO17 + eHO18);
0583 eHO2_->Fill(etaInc, phiInc, eHO17 + eHO18);
0584 eHO17_->Fill(etaInc, eHO17);
0585 eHO18_->Fill(etaInc, eHO18);
0586 eHO1T_->Fill(etaInc, eHO17T + eHO18T);
0587 eHO2T_->Fill(etaInc, phiInc, eHO17T + eHO18T);
0588 eHO17T_->Fill(etaInc, eHO17T);
0589 eHO18T_->Fill(etaInc, eHO18T);
0590 int nHO = 0, nHOT = 0;
0591 if (eHO17 > 0)
0592 nHO++;
0593 if (eHO17T > 0)
0594 nHOT++;
0595 if (eHO18 > 0)
0596 nHO++;
0597 if (eHO18T > 0)
0598 nHOT++;
0599 nHO1_->Fill(etaInc, (double)(nHO));
0600 nHO2_->Fill(etaInc, phiInc, (double)(nHO));
0601 nHO1T_->Fill(etaInc, (double)(nHOT));
0602 nHO2T_->Fill(etaInc, phiInc, (double)(nHOT));
0603 int ieta = 15;
0604 for (int k = 0; k < 15; ++k) {
0605 if (std::abs(etaInc) < 0.087 * (k + 1)) {
0606 ieta = k;
0607 break;
0608 }
0609 }
0610 if (ieta >= 0 && ieta < 15) {
0611 eHOE_[ieta]->Fill(eHO17 + eHO18);
0612 eHOE17_[ieta]->Fill(eHO17);
0613 eHOE18_[ieta]->Fill(eHO18);
0614 eHOET_[ieta]->Fill(eHO17T + eHO18T);
0615 eHOE17T_[ieta]->Fill(eHO17T);
0616 eHOE18T_[ieta]->Fill(eHO18T);
0617 eHOEta_[ieta]->Fill(eHOE17[ieta] + eHOE18[ieta]);
0618 eHOEta17_[ieta]->Fill(eHOE17[ieta]);
0619 eHOEta18_[ieta]->Fill(eHOE18[ieta]);
0620 nHOE1_[ieta]->Fill((double)(nHO));
0621 nHOE2_[ieta]->Fill(phiInc, (double)(nHO));
0622 nHOE1T_[ieta]->Fill((double)(nHOT));
0623 nHOE2T_[ieta]->Fill(phiInc, (double)(nHOT));
0624 int nHOE = 0, nHOET = 0;
0625 if (eHOE17[ieta] > 0)
0626 nHOE++;
0627 if (eHOE17T[ieta] > 0)
0628 nHOET++;
0629 if (eHOE18[ieta] > 0)
0630 nHOE++;
0631 if (eHOE18T[ieta] > 0)
0632 nHOET++;
0633 nHOEta1_[ieta]->Fill((double)(nHOE));
0634 nHOEta2_[ieta]->Fill(phiInc, (double)(nHOE));
0635 nHOEta1T_[ieta]->Fill((double)(nHOET));
0636 nHOEta2T_[ieta]->Fill(phiInc, (double)(nHOET));
0637 }
0638
0639 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::analyzeHits: Hits in EB " << nhit[0] << " in " << ebEtow.size()
0640 << " towers with total E " << etot[0] << "|" << etotT[0]
0641 << "\n Hits in HB " << nhit[1] << " in " << hbEtow.size()
0642 << " towers with total E " << etot[1] << "|" << etotT[1]
0643 << "\n Hits in HO " << nhit[2] << " in " << hoEtow.size()
0644 << " towers with total E " << etot[2] << "|" << etotT[2]
0645 << "\n Energy in EB " << eEB << "|" << eEBT << " with HB "
0646 << eEBHB << "|" << eEBHBT << " and with HO " << eEBHBHO << "|" << eEBHBHOT
0647 << "\n E in HO layers " << eHO17 << "|" << eHO17T << " "
0648 << eHO18 << "|" << eHO18T << " number of HO hits " << nHO << "|" << nHOT;
0649
0650 if (eEBHBHO > 0.75 * maxEnergy_ && print_) {
0651 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Event with excessive energy: EB = " << eEB << " EB+HB = " << eEBHB
0652 << " EB+HB+HO = " << eEBHBHO;
0653 const std::string Dets[3] = {"EB", "HB", "HO"};
0654 for (int k = 0; k < 2; k++) {
0655 int nHit;
0656 if (k == 0) {
0657 nHit = ecalHits.size();
0658 } else {
0659 nHit = hcalHits.size();
0660 }
0661 for (int i = 0; i < nHit; i++) {
0662 double edep, time;
0663 int indx = -1;
0664 unsigned int id_;
0665 int ieta, iphi, depth = 0;
0666 if (k == 0) {
0667 indx = 0;
0668 edep = ecalHits[i].energy();
0669 time = ecalHits[i].time();
0670 id_ = ecalHits[i].id();
0671 EBDetId id = EBDetId(id_);
0672 ieta = id.ieta();
0673 iphi = id.iphi();
0674 } else {
0675 indx = -1;
0676 edep = hcalHits[i].energy();
0677 time = hcalHits[i].time();
0678 id_ = hcalHits[i].id();
0679 HcalDetId id = HcalDetId(id_);
0680 int subdet = id.subdet();
0681 if (subdet == static_cast<int>(HcalBarrel))
0682 indx = 1;
0683 else if (subdet == static_cast<int>(HcalOuter))
0684 indx = 2;
0685 ieta = id.ieta();
0686 iphi = id.iphi();
0687 depth = id.depth();
0688 }
0689 if (indx >= 0) {
0690 edm::LogVerbatim("HitStudy") << "HOSimHitStudy::" << Dets[indx] << " " << i << std::hex << id_ << std::dec
0691 << " (" << ieta << "|" << iphi << "|" << depth << ") " << std::setw(8) << edep
0692 << " " << std::setw(8) << time;
0693 }
0694 }
0695 }
0696 }
0697 }
0698
0699
0700 DEFINE_FWK_MODULE(HOSimHitStudy);
