CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HFPMTHitAnalyzer.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/Frameworkfwd.h"
2 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
3 
4 #include "FWCore/Framework/interface/Event.h"
5 #include "FWCore/Framework/interface/EventSetup.h"
6 #include "FWCore/Framework/interface/MakerMacros.h"
7 #include "FWCore/MessageLogger/interface/MessageLogger.h"
8 
9 //SimHits
10 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
11 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
12 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
13 #include "SimDataFormats/Track/interface/SimTrack.h"
14 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
15 
16 //other stuff
17 #include "FWCore/ParameterSet/interface/ParameterSet.h"
18 #include "FWCore/ServiceRegistry/interface/Service.h"
19 #include "FWCore/MessageLogger/interface/MessageLogger.h"
20 #include "CommonTools/UtilAlgos/interface/TFileService.h"
21 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
22 #include "FWCore/Utilities/interface/Exception.h"
23 
24 #include <iostream>
25 #include <fstream>
26 #include <vector>
27 #include <map>
28 #include <string>
29 #include <cmath>
30 
31 #include "TH1F.h"
32 #include "TProfile.h"
33 
34 class HFPMTHitAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
35 public:
36  explicit HFPMTHitAnalyzer(const edm::ParameterSet &);
37  ~HFPMTHitAnalyzer() override {}
38  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
39 
40 private:
41  void beginJob() override;
42  void beginRun(edm::Run const &, edm::EventSetup const &) override {}
43  void analyze(const edm::Event &, const edm::EventSetup &) override;
44  void endJob() override;
45  void endRun(edm::Run const &, edm::EventSetup const &) override {}
46  void analyzeHits(std::vector<PCaloHit> &, const std::vector<SimTrack> &);
47 
48  //user parameters
49  const std::string g4Label_, hcalHits_;
50  const edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
51  const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
52  const edm::EDGetTokenT<edm::SimTrackContainer> tok_track_;
53 
54  int event_no;
55 
56  //root file name and its objects
57  TH1F *h_HFDepHit, *hHF_MC_e, *h_HFEta[3], *h_HFPhi[3];
58 
59  TH1F *hHF_e_1[3], *hHF_em_1[3], *hHF_had_1[3];
60  TH1F *hHF_e_2[3], *hHF_em_2[3], *hHF_had_2[3];
61  TH1F *hHF_e_12[3], *hHF_em_12[3], *hHF_had_12[3];
62 
63  TH1F *hHF1_time[3], *hHF1_time_Ewt[3];
64  TH1F *hHF2_time[3], *hHF2_time_Ewt[3];
65  TH1F *hHF12_time[3], *hHF12_time_Ewt[3];
66 };
67 
68 HFPMTHitAnalyzer::HFPMTHitAnalyzer(const edm::ParameterSet &iConfig)
69  : g4Label_(iConfig.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
70  hcalHits_(iConfig.getUntrackedParameter<std::string>("HitCollection", "HcalHits")),
71  tok_evt_(consumes<edm::HepMCProduct>(
72  edm::InputTag(iConfig.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")))),
73  tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
74  tok_track_(consumes<edm::SimTrackContainer>(edm::InputTag(g4Label_))) {
75  usesResource(TFileService::kSharedResource);
76 }
77 
78 void HFPMTHitAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
79  edm::ParameterSetDescription desc;
80  desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
81  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
82  desc.addUntracked<std::string>("HitCollection", "HcalHits");
83  descriptions.add("HFPMTHitAnalyzer", desc);
84 }
85 
86 void HFPMTHitAnalyzer::beginJob() {
87  event_no = 0;
88  char name[20], title[120], sub[11];
89 
90  edm::Service<TFileService> fs;
91  if (!fs.isAvailable())
92  throw cms::Exception("BadConfig") << "TFileService unavailable: "
93  << "please add it to config file";
94 
95  TFileDirectory HFHitsDir = fs->mkdir("HFPMTHits");
96  h_HFDepHit = HFHitsDir.make<TH1F>("Hit20", "Depths in HF", 50, 0., 50.);
97  h_HFDepHit->GetXaxis()->SetTitle("Depths in HF");
98  for (int i = 0; i < 3; i++) {
99  if (i == 0)
100  sprintf(sub, "(PMT)");
101  else if (i == 1)
102  sprintf(sub, "(Bundle)");
103  else
104  sprintf(sub, "(Jungle)");
105  sprintf(name, "Eta%d", i);
106  sprintf(title, "Eta Index of hits in %s", sub);
107  h_HFEta[i] = HFHitsDir.make<TH1F>(name, title, 100, 0, 100.);
108  h_HFEta[i]->GetXaxis()->SetTitle(title);
109  sprintf(name, "Phi%d", i);
110  sprintf(title, "Phi Index of hits in %s", sub);
111  h_HFPhi[i] = HFHitsDir.make<TH1F>(name, title, 100, 0, 100.);
112  h_HFPhi[i]->GetXaxis()->SetTitle(title);
113  }
114  TFileDirectory HFSourcePart = fs->mkdir("HFMCinfo");
115  hHF_MC_e = HFSourcePart.make<TH1F>("MCEnergy", "Energy of Generated Particle", 1000, 0., 500.);
116  hHF_MC_e->GetXaxis()->SetTitle("Energy of Generated Particle");
117 
118  //energy Histograms
119  TFileDirectory HFPCaloHitEnergyDir = fs->mkdir("HFPCaloHitEnergy");
120  for (int i = 0; i < 3; i++) {
121  if (i == 0)
122  sprintf(sub, "(Absorber)");
123  else if (i == 1)
124  sprintf(sub, "(PMT)");
125  else
126  sprintf(sub, "(All)");
127  sprintf(name, "Energy1%d", i);
128  sprintf(title, "Energy in depth 1 %s", sub);
129  hHF_e_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
130  hHF_e_1[i]->GetXaxis()->SetTitle(title);
131  sprintf(name, "Energy2%d", i);
132  sprintf(title, "Energy in depth 2 %s", sub);
133  hHF_e_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
134  hHF_e_2[i]->GetXaxis()->SetTitle(title);
135  sprintf(name, "Energy12%d", i);
136  sprintf(title, "Energy in depths 1,2 %s", sub);
137  hHF_e_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
138  hHF_e_12[i]->GetXaxis()->SetTitle(title);
139  sprintf(name, "Em1%d", i);
140  sprintf(title, "EM energy in depth 1 %s", sub);
141  hHF_em_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
142  hHF_em_1[i]->GetXaxis()->SetTitle(title);
143  sprintf(name, "Em2%d", i);
144  sprintf(title, "EM energy in depth 2 %s", sub);
145  hHF_em_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
146  hHF_em_2[i]->GetXaxis()->SetTitle(title);
147  sprintf(name, "Em12%d", i);
148  sprintf(title, "EM energy in depths 1,2 %s", sub);
149  hHF_em_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 500.);
150  hHF_em_12[i]->GetXaxis()->SetTitle(title);
151  sprintf(name, "Had1%d", i);
152  sprintf(title, "Had energy in depth 1 %s", sub);
153  hHF_had_1[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
154  hHF_had_1[i]->GetXaxis()->SetTitle(title);
155  sprintf(name, "Had2%d", i);
156  sprintf(title, "Had energy in depth 2 %s", sub);
157  hHF_had_2[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
158  hHF_had_2[i]->GetXaxis()->SetTitle(title);
159  sprintf(name, "Had12%d", i);
160  sprintf(title, "Had energy in depths 1,2 %s", sub);
161  hHF_had_12[i] = HFPCaloHitEnergyDir.make<TH1F>(name, title, 1000, 0., 0.1);
162  hHF_had_12[i]->GetXaxis()->SetTitle(title);
163  }
164 
165  //Timing Histograms
166  TFileDirectory HFPCaloHitTimeDir = fs->mkdir("HFPCaloHitTime");
167  for (int i = 0; i < 3; i++) {
168  if (i == 0)
169  sprintf(sub, "(Absorber)");
170  else if (i == 1)
171  sprintf(sub, "(PMT)");
172  else
173  sprintf(sub, "(All)");
174  sprintf(name, "Time1Ewt%d", i);
175  sprintf(title, "Time (energy weighted) in depth 1 %s", sub);
176  hHF1_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
177  hHF1_time_Ewt[i]->GetXaxis()->SetTitle(title);
178  sprintf(name, "Time2Ewt%d", i);
179  sprintf(title, "Time (energy weighted) in depth 2 %s", sub);
180  hHF2_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
181  hHF2_time_Ewt[i]->GetXaxis()->SetTitle(title);
182  sprintf(name, "Time12Ewt%d", i);
183  sprintf(title, "Time (energy weighted) in depths 1,2 %s", sub);
184  hHF12_time_Ewt[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
185  hHF12_time_Ewt[i]->GetXaxis()->SetTitle(title);
186  sprintf(name, "Time1%d", i);
187  sprintf(title, "Time in depth 1 %s", sub);
188  hHF1_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
189  hHF1_time[i]->GetXaxis()->SetTitle(title);
190  sprintf(name, "Time2%d", i);
191  sprintf(title, "Time in depth 2 %s", sub);
192  hHF2_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
193  hHF2_time[i]->GetXaxis()->SetTitle(title);
194  sprintf(name, "Time12%d", i);
195  sprintf(title, "Time in depths 1,2 %s", sub);
196  hHF12_time[i] = HFPCaloHitTimeDir.make<TH1F>(name, title, 400, 0., 400.);
197  hHF12_time[i]->GetXaxis()->SetTitle(title);
198  }
199 }
200 
201 void HFPMTHitAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
202  ++event_no;
203  if (event_no % 500 == 0)
204  edm::LogVerbatim("HcalSim") << "Event # " << event_no << " processed.";
205 
206  const edm::Handle<edm::PCaloHitContainer> &hitsHcal = iEvent.getHandle(tok_calo_);
207  const edm::Handle<edm::SimTrackContainer> &Tracks = iEvent.getHandle(tok_track_);
208  const edm::Handle<edm::HepMCProduct> EvtHandle = iEvent.getHandle(tok_evt_);
209  const HepMC::GenEvent *myGenEvent = EvtHandle->GetEvent();
210 
211  float orig_energy = 0;
212  for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
213  ++p) {
214  orig_energy = (*p)->momentum().e();
215  hHF_MC_e->Fill(orig_energy);
216  }
217 
218  std::vector<PCaloHit> caloHits;
219  caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
220  analyzeHits(caloHits, *Tracks);
221 }
222 
223 void HFPMTHitAnalyzer::analyzeHits(std::vector<PCaloHit> &hits, const std::vector<SimTrack> &tracks1) {
224  int nHit = hits.size();
225  float energy1[3], energy2[3], energy12[3];
226  float em1[3], had1[3], em2[3], had2[3], em12[3], had12[3];
227  for (int i = 0; i < 3; i++) {
228  energy1[i] = energy2[i] = energy12[i] = 0;
229  em1[i] = em2[i] = em12[i] = 0;
230  had1[i] = had2[i] = had12[i] = 0;
231  }
232  edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser::Entry " << event_no << " with " << nHit << " hits";
233  for (int i = 0; i < nHit; i++) {
234  double energy = hits[i].energy();
235  double em = hits[i].energyEM();
236  double had = hits[i].energyHad();
237  double time = hits[i].time();
238  uint32_t id_ = hits[i].id();
239  uint16_t pmtHit = hits[i].depth();
240  uint16_t depthX = pmtHit;
241  edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser::Hit " << i << " ID " << std::hex << id_ << std::dec << " PMT "
242  << pmtHit << " E (e|h|t) " << em << " " << had << " " << energy << " Time " << time;
243 
244  HcalDetId id = HcalDetId(id_);
245  int det = id.det();
246  int subdet = id.subdet();
247  int depth = id.depth();
248  if (pmtHit != 0)
249  pmtHit = 1;
250 
251  if (det == 4) {
252  if (subdet == static_cast<int>(HcalForward)) {
253  h_HFDepHit->Fill(double(depth + 10 * depthX));
254  if (depthX > 0) {
255  int ieta = id.ietaAbs();
256  int iphi = id.iphi();
257  if (depth != 1) {
258  ieta += 50;
259  iphi += 50;
260  }
261  if (depthX <= 3) {
262  h_HFEta[depthX - 1]->Fill(double(ieta));
263  h_HFPhi[depthX - 1]->Fill(double(iphi));
264  }
265  }
266  if (depth == 1) {
267  energy1[pmtHit] += energy;
268  energy12[pmtHit] += energy;
269  em1[pmtHit] += em;
270  em12[pmtHit] += em;
271  had1[pmtHit] += had;
272  had12[pmtHit] += had;
273  energy1[2] += energy;
274  energy12[2] += energy;
275  em1[2] += em;
276  em12[2] += em;
277  had1[2] += had;
278  had12[2] += had;
279  hHF1_time[pmtHit]->Fill(time);
280  hHF1_time[2]->Fill(time);
281  hHF1_time_Ewt[pmtHit]->Fill(time, energy);
282  hHF1_time_Ewt[2]->Fill(time, energy);
283  }
284  if (depth == 2) {
285  energy2[pmtHit] += energy;
286  energy12[pmtHit] += energy;
287  em2[pmtHit] += em;
288  em12[pmtHit] += em;
289  had2[pmtHit] += had;
290  had12[pmtHit] += had;
291  energy2[2] += energy;
292  energy12[2] += energy;
293  em2[2] += em;
294  em12[2] += em;
295  had2[2] += had;
296  had12[2] += had;
297  hHF2_time[pmtHit]->Fill(time);
298  hHF2_time[2]->Fill(time);
299  hHF2_time_Ewt[pmtHit]->Fill(time, energy);
300  hHF2_time_Ewt[2]->Fill(time, energy);
301  }
302  }
303  }
304  }
305  for (int i = 0; i < 3; i++) {
306  hHF_e_1[i]->Fill(energy1[i]);
307  hHF_e_2[i]->Fill(energy2[i]);
308  hHF_e_12[i]->Fill(energy12[i]);
309  hHF_em_1[i]->Fill(em1[i]);
310  hHF_em_2[i]->Fill(em2[i]);
311  hHF_em_12[i]->Fill(em12[i]);
312  hHF_had_1[i]->Fill(had1[i]);
313  hHF_had_2[i]->Fill(had2[i]);
314  hHF_had_12[i]->Fill(had12[i]);
315  edm::LogVerbatim("HFShower") << "HFPMTHitAnalyser:: Type " << i << " Energy 1|2| " << energy1[i] << " "
316  << energy2[i] << " " << energy12[i] << " EM Energy 1|2| " << em1[i] << " " << em2[i]
317  << " " << em12[i] << " Had Energy 1|2| " << had1[i] << " " << had2[i] << " "
318  << had12[i];
319  }
320 }
321 
322 void HFPMTHitAnalyzer::endJob() {}
323 
324 DEFINE_FWK_MODULE(HFPMTHitAnalyzer);
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
edm::PCaloHitContainer
std::vector< PCaloHit > PCaloHitContainer
Definition: PCaloHitContainer.h:8
mps_fire.i
i
Definition: mps_fire.py:428
edm::ParameterSetDescription::addUntracked
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:100
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HFPMTHitAnalyzer::hHF12_time_Ewt
TH1F * hHF12_time_Ewt[3]
Definition: HFPMTHitAnalyzer.cc:65
HFPMTHitAnalyzer::~HFPMTHitAnalyzer
~HFPMTHitAnalyzer() override
Definition: HFPMTHitAnalyzer.cc:37
HLT_FULL_cff.Tracks
tuple Tracks
Definition: HLT_FULL_cff.py:65533
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HFPMTHitAnalyzer::beginJob
void beginJob() override
Definition: HFPMTHitAnalyzer.cc:86
HFPMTHitAnalyzer::tok_calo_
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_calo_
Definition: HFPMTHitAnalyzer.cc:51
edm::Event::getHandle
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:563
HFPMTHitAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: HFPMTHitAnalyzer.cc:201
iEvent
int iEvent
Definition: GenABIO.cc:224
HFPMTHitAnalyzer::hcalHits_
const std::string hcalHits_
Definition: HFPMTHitAnalyzer.cc:49
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
TFileDirectory::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileDirectory.h:53
HFPMTHitAnalyzer::HFPMTHitAnalyzer
HFPMTHitAnalyzer(const edm::ParameterSet &)
Definition: HFPMTHitAnalyzer.cc:68
HLT_FULL_cff.depth
tuple depth
Definition: HLT_FULL_cff.py:12091
HFPMTHitAnalyzer::tok_evt_
const edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
Definition: HFPMTHitAnalyzer.cc:50
HFPMTHitAnalyzer::tok_track_
const edm::EDGetTokenT< edm::SimTrackContainer > tok_track_
Definition: HFPMTHitAnalyzer.cc:52
TFileService::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
HFPMTHitAnalyzer::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: HFPMTHitAnalyzer.cc:42
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:71427
HFPMTHitAnalyzer::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: HFPMTHitAnalyzer.cc:45
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
edm::SimTrackContainer
std::vector< SimTrack > SimTrackContainer
Definition: SimTrackContainer.h:12
HFPMTHitAnalyzer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HFPMTHitAnalyzer.cc:78
edm::Run
Definition: Run.h:45
HFPMTHitAnalyzer::analyzeHits
void analyzeHits(std::vector< PCaloHit > &, const std::vector< SimTrack > &)
Definition: HFPMTHitAnalyzer.cc:223
HFPMTHitAnalyzer::endJob
void endJob() override
Definition: HFPMTHitAnalyzer.cc:322
HFPMTHitAnalyzer::g4Label_
const std::string g4Label_
Definition: HFPMTHitAnalyzer.cc:49
Generated for CMSSW Reference Manual by doxygen 1.8.5