CMS 3D CMS Logo
HcalTB06Analysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HcalTestBeam
4 // Class : HcalTB06Analysis
5 //
6 // Implementation:
7 // Main analysis class for Hcal Test Beam 2006 Analysis
8 //
9 // Original Author:
10 // Created: 19 November 2015
11 //
12 
13 // user include files
14 #include "HcalTB06Histo.h"
15 #include "HcalTB06BeamSD.h"
16 
17 // to retreive hits
18 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
19 #include "DataFormats/EcalDetId/interface/EBDetId.h"
20 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
21 #include "DataFormats/Common/interface/Handle.h"
22 
23 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
24 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
25 
26 #include "FWCore/Framework/interface/ESHandle.h"
27 #include "FWCore/PluginManager/interface/ModuleDef.h"
28 #include "FWCore/MessageLogger/interface/MessageLogger.h"
29 #include "FWCore/Framework/interface/Frameworkfwd.h"
30 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
31 #include "FWCore/Framework/interface/Event.h"
32 #include "FWCore/Framework/interface/MakerMacros.h"
33 #include "FWCore/Framework/interface/EventSetup.h"
34 #include "FWCore/Framework/interface/Run.h"
35 #include "FWCore/ParameterSet/interface/ParameterSet.h"
36 
37 #include "CommonTools/UtilAlgos/interface/TFileService.h"
38 
39 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
40 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
41 
42 #include "CLHEP/Units/GlobalSystemOfUnits.h"
43 #include "CLHEP/Units/GlobalPhysicalConstants.h"
44 #include "globals.hh"
45 #include "Randomize.hh"
46 
47 // system include files
48 #include <iostream>
49 #include <iomanip>
50 #include <memory>
51 #include <vector>
52 #include <string>
53 
54 //#define EDM_ML_DEBUG
55 
56 class HcalTB06Analysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
57 public:
58  explicit HcalTB06Analysis(const edm::ParameterSet& p);
59  ~HcalTB06Analysis() override;
60 
61  void beginJob() override;
62  void endJob() override;
63  void analyze(const edm::Event& e, const edm::EventSetup& c) override;
64 
65  HcalTB06Analysis(const HcalTB06Analysis&) = delete;
66  const HcalTB06Analysis& operator=(const HcalTB06Analysis&) = delete;
67 
68 private:
69  edm::EDGetTokenT<edm::PCaloHitContainer> m_EcalToken;
70  edm::EDGetTokenT<edm::PCaloHitContainer> m_HcalToken;
71  edm::EDGetTokenT<edm::PCaloHitContainer> m_BeamToken;
72  bool m_ECAL;
73 
74  int count;
75  int m_idxetaEcal;
76  int m_idxphiEcal;
77  int m_idxetaHcal;
78  int m_idxphiHcal;
79 
80  double m_eta;
81  double m_phi;
82  double m_ener;
83  double m_timeLimit;
84  double m_widthEcal;
85  double m_widthHcal;
86  double m_factEcal;
87  double m_factHcal;
88  std::vector<int> m_PDG;
89 
90  HcalTB06Histo* m_histo;
91 };
92 
93 HcalTB06Analysis::HcalTB06Analysis(const edm::ParameterSet& p) : count(0) {
94  usesResource("TFileService");
95 
96  m_ECAL = p.getParameter<bool>("ECAL");
97  if (m_ECAL) {
98  m_EcalToken = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
99  }
100  m_HcalToken = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
101  m_BeamToken = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalTB06BeamHits"));
102  m_eta = p.getParameter<double>("MinEta");
103  m_phi = p.getParameter<double>("MinPhi");
104  m_ener = p.getParameter<double>("MinE");
105  m_PDG = p.getParameter<std::vector<int> >("PartID");
106 
107  double minEta = p.getParameter<double>("MinEta");
108  double maxEta = p.getParameter<double>("MaxEta");
109  double minPhi = p.getParameter<double>("MinPhi");
110  double maxPhi = p.getParameter<double>("MaxPhi");
111  double beamEta = (maxEta + minEta) * 0.5;
112  double beamPhi = (maxPhi + minPhi) * 0.5;
113  if (beamPhi < 0) {
114  beamPhi += twopi;
115  }
116 
117  m_idxetaEcal = 13;
118  m_idxphiEcal = 13;
119 
120  m_idxetaHcal = (int)(beamEta / 0.087) + 1;
121  m_idxphiHcal = (int)(beamPhi / 0.087) + 6;
122  if (m_idxphiHcal > 72) {
123  m_idxphiHcal -= 73;
124  }
125 
126  edm::ParameterSet ptb = p.getParameter<edm::ParameterSet>("TestBeamAnalysis");
127  m_timeLimit = ptb.getParameter<double>("TimeLimit");
128  m_widthEcal = ptb.getParameter<double>("EcalWidth");
129  m_widthHcal = ptb.getParameter<double>("HcalWidth");
130  m_factEcal = ptb.getParameter<double>("EcalFactor");
131  m_factHcal = ptb.getParameter<double>("HcalFactor");
132  double eMIP = ptb.getParameter<double>("MIP");
133 
134  edm::LogInfo("HcalTB06Analysis") << "Beam parameters: E(GeV)= " << m_ener << " pdgID= " << m_PDG[0]
135  << "\n eta= " << m_eta << " idx_etaEcal= " << m_idxetaEcal
136  << " idx_etaHcal= " << m_idxetaHcal << " phi= " << m_phi
137  << " idx_phiEcal= " << m_idxphiEcal << " idx_phiHcal= " << m_idxphiHcal
138  << "\n EcalFactor= " << m_factEcal << " EcalWidth= " << m_widthEcal << " GeV"
139  << "\n HcalFactor= " << m_factHcal << " HcalWidth= " << m_widthHcal << " GeV"
140  << " MIP= " << eMIP << " GeV"
141  << "\n TimeLimit= " << m_timeLimit << " ns"
142  << "\n";
143  m_histo = new HcalTB06Histo(ptb);
144 }
145 
146 HcalTB06Analysis::~HcalTB06Analysis() { delete m_histo; }
147 
148 void HcalTB06Analysis::beginJob() { edm::LogInfo("HcalTB06Analysis") << " =====> Begin of Run"; }
149 
150 void HcalTB06Analysis::endJob() {
151  edm::LogInfo("HcalTB06Analysis") << " =====> End of Run; Total number of events: " << count;
152 }
153 
154 void HcalTB06Analysis::analyze(const edm::Event& evt, const edm::EventSetup&) {
155  ++count;
156 
157  //Beam Information
158  m_histo->fillPrimary(m_ener, m_eta, m_phi);
159 
160  edm::Handle<edm::PCaloHitContainer> Ecal;
161  edm::Handle<edm::PCaloHitContainer> Hcal;
162  edm::Handle<edm::PCaloHitContainer> Beam;
163  std::vector<double> eCalo(6, 0), eTrig(7, 0);
164 
165  const std::vector<PCaloHit>* EcalHits = nullptr;
166  if (m_ECAL) {
167  evt.getByToken(m_EcalToken, Ecal);
168  EcalHits = Ecal.product();
169  }
170  evt.getByToken(m_HcalToken, Hcal);
171  const std::vector<PCaloHit>* HcalHits = Hcal.product();
172  evt.getByToken(m_BeamToken, Beam);
173  const std::vector<PCaloHit>* BeamHits = Beam.product();
174 
175  // Total Energy
176  double eecals = 0.;
177  double ehcals = 0.;
178 
179  unsigned int ne = 0;
180  unsigned int nh = 0;
181  if (m_ECAL) {
182  ne = EcalHits->size();
183  for (unsigned int i = 0; i < ne; ++i) {
184  EBDetId ecalid((*EcalHits)[i].id());
185 #ifdef EDM_ML_DEBUG
186  std::cout << "EB " << i << " " << ecalid.ieta() << ":" << m_idxetaEcal << " " << ecalid.iphi() << ":"
187  << m_idxphiEcal << " " << (*EcalHits)[i].time() << ":" << m_timeLimit << " "
188  << (*EcalHits)[i].energy() << std::endl;
189 #endif
190  // 7x7 crystal selection
191  if (std::abs(m_idxetaEcal - ecalid.ieta()) <= 3 && std::abs(m_idxphiEcal - ecalid.iphi()) <= 3 &&
192  (*EcalHits)[i].time() < m_timeLimit) {
193  eCalo[0] += (*EcalHits)[i].energy();
194  }
195  }
196  if (m_widthEcal > 0.0) {
197  eCalo[1] = G4RandGauss::shoot(0.0, m_widthEcal);
198  }
199  eecals = m_factEcal * (eCalo[0] + eCalo[1]);
200  }
201  if (HcalHits) {
202  nh = HcalHits->size();
203  for (unsigned int i = 0; i < nh; ++i) {
204  HcalDetId hcalid((*HcalHits)[i].id());
205 #ifdef EDM_ML_DEBUG
206  std::cout << "HC " << i << " " << hcalid.subdet() << " " << hcalid.ieta() << ":" << m_idxetaHcal << " "
207  << hcalid.iphi() << ":" << m_idxphiHcal << " " << (*HcalHits)[i].time() << ":" << m_timeLimit << " "
208  << (*HcalHits)[i].energy() << std::endl;
209 #endif
210  // 3x3 towers selection
211  if (std::abs(m_idxetaHcal - hcalid.ieta()) <= 1 && std::abs(m_idxphiHcal - hcalid.iphi()) <= 1 &&
212  (*HcalHits)[i].time() < m_timeLimit) {
213  if (hcalid.subdet() != HcalOuter) {
214  eCalo[2] += (*HcalHits)[i].energy();
215  } else {
216  eCalo[4] += (*HcalHits)[i].energy();
217  }
218  }
219  }
220  if (m_widthHcal > 0.0) {
221  eCalo[3] = G4RandGauss::shoot(0.0, m_widthHcal);
222  eCalo[5] = G4RandGauss::shoot(0.0, m_widthHcal);
223  }
224  ehcals = m_factHcal * eCalo[2] + eCalo[3];
225  }
226  double etots = eecals + ehcals;
227 
228  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Etot(MeV)= " << etots << " E(Ecal)= " << eecals
229  << " E(Hcal)= " << ehcals << " Nhits(ECAL)= " << ne << " Nhits(HCAL)= " << nh;
230  m_histo->fillEdep(etots, eecals, ehcals);
231 
232  if (BeamHits) {
233  for (unsigned int i = 0; i < BeamHits->size(); ++i) {
234  unsigned int id = ((*BeamHits)[i].id());
235  int det, lay, ix, iy;
236  HcalTestBeamNumbering::unpackIndex(id, det, lay, ix, iy);
237  if ((det == 1) && ((*BeamHits)[i].time() < m_timeLimit)) {
238  if (lay > 0 && lay <= 4) {
239  eTrig[lay - 1] += (*BeamHits)[i].energy();
240  } else if (lay == 7 || lay == 8) {
241  eTrig[lay - 2] += (*BeamHits)[i].energy();
242  } else if (lay >= 11 && lay <= 14) {
243  eTrig[4] += (*BeamHits)[i].energy();
244  }
245  }
246  }
247  }
248 
249  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Trigger Info: " << eTrig[0] << ":" << eTrig[1] << ":" << eTrig[2]
250  << ":" << eTrig[3] << ":" << eTrig[4] << ":" << eTrig[5] << ":" << eTrig[6];
251 
252  m_histo->fillTree(eCalo, eTrig);
253 }
254 
255 #include "FWCore/PluginManager/interface/ModuleDef.h"
256 
257 DEFINE_FWK_MODULE(HcalTB06Analysis);
HcalTB06Analysis::m_idxphiHcal
int m_idxphiHcal
Definition: HcalTB06Analysis.cc:78
HcalTB06Analysis::m_idxetaEcal
int m_idxetaEcal
Definition: HcalTB06Analysis.cc:75
Handle.h
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
EDAnalyzer.h
mps_fire.i
i
Definition: mps_fire.py:428
MessageLogger.h
HcalTB06Analysis::~HcalTB06Analysis
~HcalTB06Analysis() override
Definition: HcalTB06Analysis.cc:146
edm::Handle::product
T const * product() const
Definition: Handle.h:70
ESHandle.h
HcalTB06Analysis::m_ener
double m_ener
Definition: HcalTB06Analysis.cc:82
HcalDetId::iphi
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
CaloG4HitCollection.h
HcalTB06Analysis::m_idxphiEcal
int m_idxphiEcal
Definition: HcalTB06Analysis.cc:76
edm::EDGetTokenT< edm::PCaloHitContainer >
EBDetId
Definition: EBDetId.h:17
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HcalTB06Histo
Definition: HcalTB06Histo.h:31
gather_cfg.cout
cout
Definition: gather_cfg.py:144
EBDetId.h
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
edm::one::EDAnalyzer
Definition: EDAnalyzer.h:30
HcalTB06Analysis::count
int count
Definition: HcalTB06Analysis.cc:74
edm::Handle
Definition: AssociativeIterator.h:50
HcalTB06Histo::fillTree
void fillTree(std::vector< double > &ecalo, std::vector< double > &etrig)
Definition: HcalTB06Histo.cc:108
HcalTB06Analysis::m_EcalToken
edm::EDGetTokenT< edm::PCaloHitContainer > m_EcalToken
Definition: HcalTB06Analysis.cc:69
HcalTB06Analysis::m_idxetaHcal
int m_idxetaHcal
Definition: HcalTB06Analysis.cc:77
HcalTB06Analysis::m_factHcal
double m_factHcal
Definition: HcalTB06Analysis.cc:87
HLT_FULL_cff.maxPhi
maxPhi
Definition: HLT_FULL_cff.py:52976
HcalTB06Analysis::endJob
void endJob() override
Definition: HcalTB06Analysis.cc:150
HcalTB06Analysis::m_widthEcal
double m_widthEcal
Definition: HcalTB06Analysis.cc:84
MakerMacros.h
HcalTB06Histo.h
HcalTB06Analysis
Definition: HcalTB06Analysis.cc:56
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HcalTB06Analysis::HcalTB06Analysis
HcalTB06Analysis(const edm::ParameterSet &p)
Definition: HcalTB06Analysis.cc:93
HcalTestBeamNumbering.h
HcalTB06BeamSD.h
cms::cuda::nh
uint32_t nh
Definition: HistoContainer.h:23
maxEta
double maxEta
Definition: PFJetBenchmarkAnalyzer.cc:76
Run.h
HcalTB06Analysis::operator=
const HcalTB06Analysis & operator=(const HcalTB06Analysis &)=delete
submitPVResolutionJobs.count
count
Definition: submitPVResolutionJobs.py:352
HcalOuter
Definition: HcalAssistant.h:35
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:531
HcalTB06Analysis::m_PDG
std::vector< int > m_PDG
Definition: HcalTB06Analysis.cc:88
TFileService.h
HcalTB06Histo::fillPrimary
void fillPrimary(double energy, double eta, double phi)
Definition: HcalTB06Histo.cc:81
HcalTB06Analysis::m_ECAL
bool m_ECAL
Definition: HcalTB06Analysis.cc:72
HcalDetId::ieta
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
edm::ParameterSet
Definition: ParameterSet.h:47
HcalTB06Analysis::analyze
void analyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: HcalTB06Analysis.cc:154
Event.h
HcalDetId.h
HcalTestBeamNumbering::unpackIndex
static void unpackIndex(const uint32_t &idx, int &det, int &lay, int &x, int &y)
Definition: HcalTestBeamNumbering.cc:31
HcalTB06Analysis::m_widthHcal
double m_widthHcal
Definition: HcalTB06Analysis.cc:85
HcalDetId::subdet
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
HcalDetId
Definition: HcalDetId.h:12
ModuleDef.h
createfilelist.int
int
Definition: createfilelist.py:10
edm::EventSetup
Definition: EventSetup.h:57
HcalSubdetector.h
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
HcalTB06Analysis::m_eta
double m_eta
Definition: HcalTB06Analysis.cc:80
HLT_FULL_cff.minPhi
minPhi
Definition: HLT_FULL_cff.py:52964
HcalTB06Analysis::m_factEcal
double m_factEcal
Definition: HcalTB06Analysis.cc:86
HcalTB06Analysis::m_phi
double m_phi
Definition: HcalTB06Analysis.cc:81
HcalTB06Analysis::beginJob
void beginJob() override
Definition: HcalTB06Analysis.cc:148
Frameworkfwd.h
HcalTB06Analysis::m_BeamToken
edm::EDGetTokenT< edm::PCaloHitContainer > m_BeamToken
Definition: HcalTB06Analysis.cc:71
EventSetup.h
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
PCaloHitContainer.h
HcalTB06Analysis::m_HcalToken
edm::EDGetTokenT< edm::PCaloHitContainer > m_HcalToken
Definition: HcalTB06Analysis.cc:70
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
ntuplemaker.time
time
Definition: ntuplemaker.py:310
edm::Event
Definition: Event.h:73
EgHLTOffEleSelection_cfi.minEta
minEta
Definition: EgHLTOffEleSelection_cfi.py:11
HcalTB06Analysis::m_timeLimit
double m_timeLimit
Definition: HcalTB06Analysis.cc:83
CaloG4Hit.h
edm::InputTag
Definition: InputTag.h:15
HcalTB06Histo::fillEdep
void fillEdep(double etots, double eecals, double ehcals)
Definition: HcalTB06Histo.cc:92
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HcalTB06Analysis::m_histo
HcalTB06Histo * m_histo
Definition: HcalTB06Analysis.cc:90
Generated for CMSSW Reference Manual by doxygen 1.8.16