Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-20 23:14:39

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/GlobalPhysicalConstants.h>
0027 #include <CLHEP/Units/GlobalSystemOfUnits.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   //initialize
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   // Loop over containers
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   // Now for towers and total energy deposits
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 //define this as a plug-in
0700 DEFINE_FWK_MODULE(HOSimHitStudy);