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