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