Belle II Software  release-05-02-19
QcsmonitorStudyModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Igal Jaegle *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <beast/qcsmonitor/modules/QcsmonitorStudyModule.h>
12 #include <beast/qcsmonitor/dataobjects/QcsmonitorSimHit.h>
13 #include <beast/qcsmonitor/dataobjects/QcsmonitorHit.h>
14 #include <generators/SAD/dataobjects/SADMetaHit.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/gearbox/GearDir.h>
17 #include <framework/logging/Logger.h>
18 #include <cmath>
19 
20 #include <fstream>
21 #include <string>
22 
23 // ROOT
24 #include <TH1.h>
25 #include <TH2.h>
26 
27 int eventNum = 0;
28 
29 using namespace std;
30 
31 using namespace Belle2;
32 using namespace qcsmonitor;
33 
34 //-----------------------------------------------------------------
35 // Register the Module
36 //-----------------------------------------------------------------
37 REG_MODULE(QcsmonitorStudy)
38 
39 //-----------------------------------------------------------------
40 // Implementation
41 //-----------------------------------------------------------------
42 
43 QcsmonitorStudyModule::QcsmonitorStudyModule() : HistoModule()
44 {
45  // Set module properties
46  setDescription("Study module for Qcsmonitors (BEAST)");
47 
48  addParam("Ethres", m_Ethres, "Energy threshold in MeV", 0.0);
49 }
50 
51 QcsmonitorStudyModule::~QcsmonitorStudyModule()
52 {
53 }
54 
55 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
56 void QcsmonitorStudyModule::defineHisto()
57 {
58  for (int i = 0; i < 40; i++) {
59  h_qcss_Evtof1[i] = new TH2F(TString::Format("qcss_Evtof1_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 500, 0., 1000.,
60  100, 0., 10.);
61  h_qcss_Evtof2[i] = new TH2F(TString::Format("qcss_Evtof2_%d", i), "Energy deposited [MeV] vs TOF [ns] - only photons", 500, 0.,
62  100., 1000, 0., 10.);
63  h_qcss_Evtof3[i] = new TH2F(TString::Format("qcss_Evtof3_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 500, 0.,
64  100., 1000, 0., 10.);
65  h_qcss_Evtof4[i] = new TH2F(TString::Format("qcss_Evtof4_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 500, 0.,
66  100., 1000, 0., 10.);
67  h_qcss_edep[i] = new TH1F(TString::Format("qcss_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 10.);
68  h_Wqcss_edep[i] = new TH1F(TString::Format("Wqcss_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 10.);
69  }
70 
71  h_qcss_hitrate0 = new TH1F("qcss_hitrate0", "Hit distributions", 100, 0., 100.);
72  h_qcss_hitrate1 = new TH1F("qcss_hitrate1", "Hit distributions", 100, 0., 100.);
73  h_qcss_hitrate2 = new TH1F("qcss_hitrate2", "Hit distributions", 100, 0., 100.);
74  h_qcss_hitrate1W = new TH1F("qcss_hitrate1W", "Hit distributions", 100, 0., 100.);
75  h_qcss_hitrate2W = new TH1F("qcss_hitrate2W", "Hit distributions", 100, 0., 100.);
76 
77  h_qcss_hitrate0->Sumw2();
78  h_qcss_hitrate1->Sumw2();
79  h_qcss_hitrate1W->Sumw2();
80  h_qcss_hitrate2->Sumw2();
81  h_qcss_hitrate2W->Sumw2();
82 
83  h_qcss_rs_hitrate1 = new TH2F("qcss_rs_hitrate1", "Hit distributions vs rs", 100, 0., 100., 12, 0., 12.);
84  h_qcss_rs_hitrate2 = new TH2F("qcss_rs_hitrate2", "Hit distributions vs rs", 100, 0., 100., 12, 0., 12.);
85  h_qcss_rs_hitrate1W = new TH2F("qcss_rs_hitrate1W", "Hit distributions vs rs", 100, 0., 100., 12, 0., 12.);
86  h_qcss_rs_hitrate2W = new TH2F("qcss_rs_hitrate2W", "Hit distributions vs rs", 100, 0., 100., 12, 0., 12.);
87 
88  h_qcss_rs_hitrate1->Sumw2();
89  h_qcss_rs_hitrate1W->Sumw2();
90  h_qcss_rs_hitrate2->Sumw2();
91  h_qcss_rs_hitrate2W->Sumw2();
92 
93  for (int i = 0; i < 40; i++) {
94  h_qcss_rate1[i] = new TH1F(TString::Format("qcss_rate1_%d", i), "PE distributions", 500, 0., 500.);
95  h_qcss_rate2[i] = new TH1F(TString::Format("qcss_rate2_%d", i), "PE distributions", 500, 0., 500.);
96  h_qcss_rate1W[i] = new TH1F(TString::Format("qcss_rate1W_%d", i), "PE distributions", 500, 0., 500.);
97  h_qcss_rate2W[i] = new TH1F(TString::Format("qcss_rate2W_%d", i), "PE distributions", 500, 0., 500.);
98  h_qcss_pe1[i] = new TH2F(TString::Format("qcss_pe1_%d", i), "PE distributions", 500, 0., 500., 100, 0., 1000.);
99  h_qcss_pe2[i] = new TH2F(TString::Format("qcss_pe2_%d", i), "PE distributions", 500, 0., 500., 100, 0., 1000.);
100  h_qcss_pe1W[i] = new TH2F(TString::Format("qcss_pe1W_%d", i), "PE distributions", 500, 0., 500., 100, 0., 1000.);
101  h_qcss_pe2W[i] = new TH2F(TString::Format("qcss_pe2W_%d", i), "PE distributions", 500, 0., 500., 100, 0., 1000.);
102 
103  h_qcss_rs_rate1[i] = new TH2F(TString::Format("qcss_rs_rate1_%d", i), "PE distributions", 500, 0., 500., 12, 0., 12.);
104  h_qcss_rs_rate2[i] = new TH2F(TString::Format("qcss_rs_rate2_%d", i), "PE distributions", 500, 0., 500., 12, 0., 12.);
105  h_qcss_rs_rate1W[i] = new TH2F(TString::Format("qcss_rs_rate1W_%d", i), "PE distributions", 500, 0., 500., 12, 0., 12.);
106  h_qcss_rs_rate2W[i] = new TH2F(TString::Format("qcss_rs_rate2W_%d", i), "PE distributions", 500, 0., 500., 12, 0., 12.);
107 
108  h_qcss_rate1[i]->Sumw2();
109  h_qcss_rate2[i]->Sumw2();
110  h_qcss_rate1W[i]->Sumw2();
111  h_qcss_rate2W[i]->Sumw2();
112  h_qcss_rs_rate1[i]->Sumw2();
113  h_qcss_rs_rate2[i]->Sumw2();
114  h_qcss_rs_rate1W[i]->Sumw2();
115  h_qcss_rs_rate2W[i]->Sumw2();
116  h_qcss_pe1[i]->Sumw2();
117  h_qcss_pe2[i]->Sumw2();
118  h_qcss_pe1W[i]->Sumw2();
119  h_qcss_pe2W[i]->Sumw2();
120  }
121 }
122 
123 
124 void QcsmonitorStudyModule::initialize()
125 {
126  B2INFO("QcsmonitorStudyModule: Initialize");
127 
128  REG_HISTOGRAM
129 
130  //get QCSMONITORS paramters ie energy threshold
131  getXMLData();
132 }
133 
134 void QcsmonitorStudyModule::beginRun()
135 {
136 }
137 
138 void QcsmonitorStudyModule::event()
139 {
140  //Here comes the actual event processing
141  StoreArray<QcsmonitorSimHit> SimHits;
142  StoreArray<QcsmonitorHit> Hits;
143  StoreArray<SADMetaHit> MetaHits;
144 
145  double rate = 0;
146  int ring_section = -1;
147  for (const auto& MetaHit : MetaHits) {
148  rate = MetaHit.getrate();
149  ring_section = MetaHit.getring_section() - 1;
150  }
151 
152  //number of entries in SimHits
153  int nSimHits = SimHits.getEntries();
154 
155  //loop over all SimHit entries
156  for (int i = 0; i < nSimHits; i++) {
157  QcsmonitorSimHit* aHit = SimHits[i];
158  int detNB = aHit->getCellId();
159  if (detNB < 40) {
160  //int trkID = aHit->getTrackId();
161  int pdg = aHit->getPDGCode();
162  double Edep = aHit->getEnergyDep() * 1e3; //GeV -> MeV
163  double tof = aHit->getFlightTime(); //ns
164 
165  h_qcss_hitrate0->Fill(detNB);
166  h_qcss_Evtof1[detNB]->Fill(tof, Edep);
167  if (pdg == 22) h_qcss_Evtof2[detNB]->Fill(tof, Edep);
168  else if (fabs(pdg) == 11) h_qcss_Evtof3[detNB]->Fill(tof, Edep);
169  else h_qcss_Evtof4[detNB]->Fill(tof, Edep);
170  if (Edep > m_Ethres) {
171  h_qcss_edep[detNB]->Fill(Edep);
172  h_Wqcss_edep[detNB]->Fill(Edep, rate);
173  }
174  }
175  }
176 
177  for (const auto& Hit : Hits) {
178  const int detNb = Hit.getdetNb();
179  if (detNb < 40) {
180  const int timebin = Hit.gettime();
181  const float edep = Hit.getedep();
182  const float pe = Hit.getPE();
183  h_qcss_hitrate1->Fill(detNb);
184  h_qcss_hitrate1W->Fill(detNb, rate);
185  h_qcss_rate1[detNb]->Fill(pe);
186  h_qcss_rate1W[detNb]->Fill(pe, rate);
187  h_qcss_rs_rate1[detNb]->Fill(pe, ring_section);
188  h_qcss_rs_rate1W[detNb]->Fill(pe, ring_section, rate);
189  h_qcss_rs_hitrate1->Fill(detNb, ring_section);
190  h_qcss_rs_hitrate1W->Fill(detNb, ring_section, rate);
191  h_qcss_pe1[detNb]->Fill(timebin, pe);
192  h_qcss_pe1W[detNb]->Fill(timebin, pe, rate);
193  if (edep > m_Ethres) {
194  h_qcss_hitrate2->Fill(detNb);
195  h_qcss_hitrate2W->Fill(detNb, rate);
196  h_qcss_rate2[detNb]->Fill(pe);
197  h_qcss_rate2W[detNb]->Fill(pe, rate);
198  h_qcss_rs_rate2[detNb]->Fill(pe, ring_section);
199  h_qcss_rs_rate2W[detNb]->Fill(pe, ring_section, rate);
200  h_qcss_rs_hitrate2->Fill(detNb, ring_section);
201  h_qcss_rs_hitrate2W->Fill(detNb, ring_section, rate);
202  h_qcss_pe2[detNb]->Fill(timebin, pe);
203  h_qcss_pe2W[detNb]->Fill(timebin, pe, rate);
204  }
205  }
206  }
207 
208 }
209 
210 //read energy threshold from QCSMONITOR.xml
211 void QcsmonitorStudyModule::getXMLData()
212 {
213  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"QCSMONITOR\"]/Content/");
214  m_Ethres = content.getDouble("Ethres");
215 
216  B2INFO("QcsmonitorStudy");
217 }
218 
219 
220 void QcsmonitorStudyModule::endRun()
221 {
222 
223 
224 
225 }
226 
227 void QcsmonitorStudyModule::terminate()
228 {
229 }
230 
231 
Belle2::QcsmonitorSimHit::getCellId
int getCellId() const
Get Cell ID.
Definition: QcsmonitorSimHit.h:107
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::QcsmonitorSimHit
ClassQcsmonitorSimHit - Geant4 simulated hit for the Qcsmonitor crystal in beast.
Definition: QcsmonitorSimHit.h:41
Belle2::QcsmonitorSimHit::getFlightTime
double getFlightTime() const
Get Flight time from IP.
Definition: QcsmonitorSimHit.h:122
Belle2::qcsmonitor::QcsmonitorStudyModule
Study module for Qcsmonitor (BEAST)
Definition: QcsmonitorStudyModule.h:39
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::QcsmonitorSimHit::getPDGCode
int getPDGCode() const
Get Particle PDG (can be one of secondaries)
Definition: QcsmonitorSimHit.h:117
Belle2::QcsmonitorSimHit::getEnergyDep
double getEnergyDep() const
Get Deposit energy.
Definition: QcsmonitorSimHit.h:127
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::HistoModule
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29
Belle2::Hit
Structure to hold some of the calpulse data.
Definition: TOPTimeBaseCalibratorModule.h:40