File indexing completed on 2023-10-25 10:04:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06Histo.h"
0015 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06BeamSD.h"
0016
0017
0018 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0019 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0020 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022
0023 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0024 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0025
0026 #include "FWCore/PluginManager/interface/ModuleDef.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/Run.h"
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036
0037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0038
0039 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
0040 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0041
0042 #include "Randomize.hh"
0043 #include "globals.hh"
0044
0045 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0046 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0047
0048
0049 #include <memory>
0050 #include <string>
0051 #include <vector>
0052
0053
0054
0055 class HcalTB06Analysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0056 public:
0057 explicit HcalTB06Analysis(const edm::ParameterSet& p);
0058 ~HcalTB06Analysis() override = default;
0059
0060 void beginJob() override;
0061 void endJob() override;
0062 void analyze(const edm::Event& e, const edm::EventSetup& c) override;
0063
0064 HcalTB06Analysis(const HcalTB06Analysis&) = delete;
0065 const HcalTB06Analysis& operator=(const HcalTB06Analysis&) = delete;
0066
0067 private:
0068 const bool m_ECAL;
0069 const double m_eta;
0070 const double m_phi;
0071 const double m_ener;
0072 const std::vector<int> m_PDG;
0073
0074 edm::EDGetTokenT<edm::PCaloHitContainer> m_EcalToken;
0075 const edm::EDGetTokenT<edm::PCaloHitContainer> m_HcalToken;
0076 const edm::EDGetTokenT<edm::PCaloHitContainer> m_BeamToken;
0077
0078 const edm::ParameterSet m_ptb;
0079 const double m_timeLimit;
0080 const double m_widthEcal;
0081 const double m_widthHcal;
0082 const double m_factEcal;
0083 const double m_factHcal;
0084 const double m_eMIP;
0085
0086 int count;
0087 int m_idxetaEcal;
0088 int m_idxphiEcal;
0089 int m_idxetaHcal;
0090 int m_idxphiHcal;
0091
0092 std::unique_ptr<HcalTB06Histo> m_histo;
0093 };
0094
0095 HcalTB06Analysis::HcalTB06Analysis(const edm::ParameterSet& p)
0096 : m_ECAL(p.getParameter<bool>("ECAL")),
0097 m_eta(p.getParameter<double>("MinEta")),
0098 m_phi(p.getParameter<double>("MinPhi")),
0099 m_ener(p.getParameter<double>("MinE")),
0100 m_PDG(p.getParameter<std::vector<int> >("PartID")),
0101 m_HcalToken(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"))),
0102 m_BeamToken(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalTB06BeamHits"))),
0103 m_ptb(p.getParameter<edm::ParameterSet>("TestBeamAnalysis")),
0104 m_timeLimit(m_ptb.getParameter<double>("TimeLimit")),
0105 m_widthEcal(m_ptb.getParameter<double>("EcalWidth")),
0106 m_widthHcal(m_ptb.getParameter<double>("HcalWidth")),
0107 m_factEcal(m_ptb.getParameter<double>("EcalFactor")),
0108 m_factHcal(m_ptb.getParameter<double>("HcalFactor")),
0109 m_eMIP(m_ptb.getParameter<double>("MIP")),
0110 count(0) {
0111 usesResource("TFileService");
0112 if (m_ECAL)
0113 m_EcalToken = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
0114 double minEta = p.getParameter<double>("MinEta");
0115 double maxEta = p.getParameter<double>("MaxEta");
0116 double minPhi = p.getParameter<double>("MinPhi");
0117 double maxPhi = p.getParameter<double>("MaxPhi");
0118 double beamEta = (maxEta + minEta) * 0.5;
0119 double beamPhi = (maxPhi + minPhi) * 0.5;
0120 if (beamPhi < 0) {
0121 beamPhi += twopi;
0122 }
0123
0124 m_idxetaEcal = 13;
0125 m_idxphiEcal = 13;
0126
0127 m_idxetaHcal = static_cast<int>(beamEta / 0.087) + 1;
0128 m_idxphiHcal = static_cast<int>(beamPhi / 0.087) + 6;
0129 if (m_idxphiHcal > 72) {
0130 m_idxphiHcal -= 73;
0131 }
0132
0133 edm::LogInfo("HcalTB06Analysis") << "Beam parameters: E(GeV)= " << m_ener << " pdgID= " << m_PDG[0]
0134 << "\n eta= " << m_eta << " idx_etaEcal= " << m_idxetaEcal
0135 << " idx_etaHcal= " << m_idxetaHcal << " phi= " << m_phi
0136 << " idx_phiEcal= " << m_idxphiEcal << " idx_phiHcal= " << m_idxphiHcal
0137 << "\n EcalFactor= " << m_factEcal << " EcalWidth= " << m_widthEcal << " GeV"
0138 << "\n HcalFactor= " << m_factHcal << " HcalWidth= " << m_widthHcal << " GeV"
0139 << " MIP= " << m_eMIP << " GeV"
0140 << "\n TimeLimit= " << m_timeLimit << " ns"
0141 << "\n";
0142 m_histo = std::make_unique<HcalTB06Histo>(m_ptb);
0143 }
0144
0145 void HcalTB06Analysis::beginJob() { edm::LogInfo("HcalTB06Analysis") << " =====> Begin of Run"; }
0146
0147 void HcalTB06Analysis::endJob() {
0148 edm::LogInfo("HcalTB06Analysis") << " =====> End of Run; Total number of events: " << count;
0149 }
0150
0151 void HcalTB06Analysis::analyze(const edm::Event& evt, const edm::EventSetup&) {
0152 ++count;
0153
0154
0155 m_histo->fillPrimary(m_ener, m_eta, m_phi);
0156
0157 std::vector<double> eCalo(6, 0), eTrig(7, 0);
0158
0159 const std::vector<PCaloHit>* EcalHits = nullptr;
0160 if (m_ECAL) {
0161 const edm::Handle<edm::PCaloHitContainer>& Ecal = evt.getHandle(m_EcalToken);
0162 EcalHits = Ecal.product();
0163 }
0164 const edm::Handle<edm::PCaloHitContainer>& Hcal = evt.getHandle(m_HcalToken);
0165 const std::vector<PCaloHit>* HcalHits = Hcal.product();
0166 const edm::Handle<edm::PCaloHitContainer>& Beam = evt.getHandle(m_BeamToken);
0167 const std::vector<PCaloHit>* BeamHits = Beam.product();
0168
0169
0170 double eecals = 0.;
0171 double ehcals = 0.;
0172
0173 unsigned int ne = 0;
0174 unsigned int nh = 0;
0175 if (m_ECAL) {
0176 ne = EcalHits->size();
0177 for (unsigned int i = 0; i < ne; ++i) {
0178 EBDetId ecalid((*EcalHits)[i].id());
0179 #ifdef EDM_ML_DEBUG
0180 edm::LogVerbatim("HcalTBSim") << "EB " << i << " " << ecalid.ieta() << ":" << m_idxetaEcal << " "
0181 << ecalid.iphi() << ":" << m_idxphiEcal << " " << (*EcalHits)[i].time() << ":"
0182 << m_timeLimit << " " << (*EcalHits)[i].energy();
0183 #endif
0184
0185 if (std::abs(m_idxetaEcal - ecalid.ieta()) <= 3 && std::abs(m_idxphiEcal - ecalid.iphi()) <= 3 &&
0186 (*EcalHits)[i].time() < m_timeLimit) {
0187 eCalo[0] += (*EcalHits)[i].energy();
0188 }
0189 }
0190 if (m_widthEcal > 0.0) {
0191 eCalo[1] = G4RandGauss::shoot(0.0, m_widthEcal);
0192 }
0193 eecals = m_factEcal * (eCalo[0] + eCalo[1]);
0194 }
0195 if (HcalHits) {
0196 nh = HcalHits->size();
0197 for (unsigned int i = 0; i < nh; ++i) {
0198 HcalDetId hcalid((*HcalHits)[i].id());
0199 #ifdef EDM_ML_DEBUG
0200 edm::LogVerbatim("HcalTBSim") << "HC " << i << " " << hcalid.subdet() << " " << hcalid.ieta() << ":"
0201 << m_idxetaHcal << " " << hcalid.iphi() << ":" << m_idxphiHcal << " "
0202 << (*HcalHits)[i].time() << ":" << m_timeLimit << " " << (*HcalHits)[i].energy();
0203 #endif
0204
0205 if (std::abs(m_idxetaHcal - hcalid.ieta()) <= 1 && std::abs(m_idxphiHcal - hcalid.iphi()) <= 1 &&
0206 (*HcalHits)[i].time() < m_timeLimit) {
0207 if (hcalid.subdet() != HcalOuter) {
0208 eCalo[2] += (*HcalHits)[i].energy();
0209 } else {
0210 eCalo[4] += (*HcalHits)[i].energy();
0211 }
0212 }
0213 }
0214 if (m_widthHcal > 0.0) {
0215 eCalo[3] = G4RandGauss::shoot(0.0, m_widthHcal);
0216 eCalo[5] = G4RandGauss::shoot(0.0, m_widthHcal);
0217 }
0218 ehcals = m_factHcal * eCalo[2] + eCalo[3];
0219 }
0220 double etots = eecals + ehcals;
0221
0222 edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Etot(MeV)= " << etots << " E(Ecal)= " << eecals
0223 << " E(Hcal)= " << ehcals << " Nhits(ECAL)= " << ne << " Nhits(HCAL)= " << nh;
0224 m_histo->fillEdep(etots, eecals, ehcals);
0225
0226 if (BeamHits) {
0227 for (unsigned int i = 0; i < BeamHits->size(); ++i) {
0228 unsigned int id = ((*BeamHits)[i].id());
0229 int det, lay, ix, iy;
0230 HcalTestBeamNumbering::unpackIndex(id, det, lay, ix, iy);
0231 if ((det == 1) && ((*BeamHits)[i].time() < m_timeLimit)) {
0232 if (lay > 0 && lay <= 4) {
0233 eTrig[lay - 1] += (*BeamHits)[i].energy();
0234 } else if (lay == 7 || lay == 8) {
0235 eTrig[lay - 2] += (*BeamHits)[i].energy();
0236 } else if (lay >= 11 && lay <= 14) {
0237 eTrig[4] += (*BeamHits)[i].energy();
0238 }
0239 }
0240 }
0241 }
0242
0243 edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Trigger Info: " << eTrig[0] << ":" << eTrig[1] << ":" << eTrig[2]
0244 << ":" << eTrig[3] << ":" << eTrig[4] << ":" << eTrig[5] << ":" << eTrig[6];
0245
0246 m_histo->fillTree(eCalo, eTrig);
0247 }
0248
0249 #include "FWCore/PluginManager/interface/ModuleDef.h"
0250
0251 DEFINE_FWK_MODULE(HcalTB06Analysis);
