CMS 3D CMS Logo

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