Belle II Software  release-06-00-14
CsiStudy_v2Module.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <beast/csi/modules/CsiStudy_v2Module.h>
10 #include <framework/logging/Logger.h>
11 #include <framework/gearbox/GearDir.h>
12 #include <framework/gearbox/Const.h>
13 #include <cmath>
14 
15 #include <fstream>
16 #include <string>
17 
18 // ROOT
19 #include <TH1.h>
20 #include <TH2.h>
21 
22 using namespace std;
23 
24 using namespace Belle2;
25 using namespace csi;
26 
27 //-----------------------------------------------------------------
28 // Register the Module
29 //-----------------------------------------------------------------
30 REG_MODULE(CsiStudy_v2)
31 
32 //-----------------------------------------------------------------
33 // Implementation
34 //-----------------------------------------------------------------
35 
36 CsiStudy_v2Module::CsiStudy_v2Module() : HistoModule()
37 {
38  // Set module properties
39  setDescription("Study_v2 module for Csis (BEAST)");
40 
41  addParam("Ethres", m_Ethres, "Energy threshold in MeV", 0.0);
42 }
43 
44 CsiStudy_v2Module::~CsiStudy_v2Module()
45 {
46 }
47 
48 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
49 void CsiStudy_v2Module::defineHisto()
50 {
51  for (int i = 0; i < 153; i++) {
52  h_csi_drate[i] = new TH1F(TString::Format("csi_rate_%d", i), "count", 18, 0., 18.);
53  h_csi_drate[i]->Sumw2();
54 
55  h_csi_rs_drate[i] = new TH2F(TString::Format("csi_rs_drate_%d", i), "count v. ring section", 18, 0., 18., 12, 0., 12.);
56  h_csi_rs_drate[i]->Sumw2();
57 
58  for (int j = 0; j < 18; j++) {
59 
60  h_csi_dedep[j][i] = new TH1F(TString::Format("csi_dedep_%d_%d", j, i), "Energy deposited [MeV]", 5000, 0., 400.);
61  h_csi_rs_dedep[j][i] = new TH2F(TString::Format("csi_rs_dedep_%d_%d", j, i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
62 
63  h_csi_denergy[j][i] = new TH1F(TString::Format("csi_denergy_%d_%d", j, i), "Energy deposited [MeV]", 200, 0.1, 3.);
64  h_csi_rs_denergy[j][i] = new TH2F(TString::Format("csi_rs_denergy_%d_%d", j, i), "Energy deposited [MeV] per section", 200, 0.1, 3.,
65  12, 0.,
66  12.);
67 
68  h_csi_dedep[j][i]->Sumw2();
69  h_csi_rs_dedep[j][i]->Sumw2();
70  h_csi_denergy[j][i]->Sumw2();
71  h_csi_rs_denergy[j][i]->Sumw2();
72  }
73  }
74  /*
75  for (int i = 0; i < 2; i++) {
76  h_csi_rate[i] = new TH1F(TString::Format("csi_rate_%d", i), "count", 18, 0., 18.);
77  h_csi_rate[i]->Sumw2();
78 
79  h_csi_rs_rate[i] = new TH2F(TString::Format("csi_rs_rate_%d", i), "count v. ring section", 18, 0., 18., 12, 0., 12.);
80  h_csi_rs_rate[i]->Sumw2();
81  }
82  for (int i = 0; i < 18; i++) {
83  h_csi_energyVrs1[i] = new TH2F(TString::Format("csi_energyVrs1_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
84  12.);
85  h_csi_energyVrs2[i] = new TH2F(TString::Format("csi_energyVrs2_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
86  12.);
87  h_csi_energyVrs3[i] = new TH2F(TString::Format("csi_energyVrs3_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
88  12.);
89 
90  h_csi_energyVrs1W[i] = new TH2F(TString::Format("csi_energyVrs1W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
91  0.,
92  12.);
93  h_csi_energyVrs2W[i] = new TH2F(TString::Format("csi_energyVrs2W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
94  0.,
95  12.);
96  h_csi_energyVrs3W[i] = new TH2F(TString::Format("csi_energyVrs3W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
97  0.,
98  12.);
99 
100  h_csi_energy1[i] = new TH1F(TString::Format("csi_energy1_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
101  h_csi_energy2[i] = new TH1F(TString::Format("csi_energy2_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
102  h_csi_energy3[i] = new TH1F(TString::Format("csi_energy3_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
103 
104  h_csi_energy1W[i] = new TH1F(TString::Format("csi_energy1W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
105  h_csi_energy2W[i] = new TH1F(TString::Format("csi_energy2W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
106  h_csi_energy3W[i] = new TH1F(TString::Format("csi_energy3W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
107 
108  h_csi_Evtof[i] = new TH2F(TString::Format("csi_Evtof_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
109  1000, 0., 400.);
110  h_csi_Evtof1[i] = new TH2F(TString::Format("csi_Evtof1_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
111  1000, 0., 400.);
112  h_csi_Evtof2[i] = new TH2F(TString::Format("csi_Evtof2_%d", i), "Energy deposited [MeV] vs TOF [ns] - only photons", 5000, 0.,
113  1000., 1000, 0., 400.);
114  h_csi_Evtof3[i] = new TH2F(TString::Format("csi_Evtof3_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 5000, 0.,
115  1000., 1000, 0., 400.);
116 
117  h_csi_edep[i] = new TH1F(TString::Format("csi_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
118  h_csi_edep1[i] = new TH1F(TString::Format("csi_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
119  h_csi_edep2[i] = new TH1F(TString::Format("csi_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
120  h_csi_edep1Weight[i] = new TH1F(TString::Format("csi_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
121  h_csi_edep2Weight[i] = new TH1F(TString::Format("csi_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
122 
123  h_csi_rs_edep1[i] = new TH2F(TString::Format("csi_rs_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
124  h_csi_rs_edep2[i] = new TH2F(TString::Format("csi_rs_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
125  h_csi_rs_edep1Weight[i] = new TH2F(TString::Format("csi_rs_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
126  12.);
127  h_csi_rs_edep2Weight[i] = new TH2F(TString::Format("csi_rs_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
128  12.);
129 
130  h_csi_energy1[i]->Sumw2();
131  h_csi_energy2[i]->Sumw2();
132  h_csi_energy3[i]->Sumw2();
133  h_csi_energy1W[i]->Sumw2();
134  h_csi_energy2W[i]->Sumw2();
135  h_csi_energy3W[i]->Sumw2();
136  h_csi_energyVrs1[i]->Sumw2();
137  h_csi_energyVrs2[i]->Sumw2();
138  h_csi_energyVrs3[i]->Sumw2();
139  h_csi_energyVrs1W[i]->Sumw2();
140  h_csi_energyVrs2W[i]->Sumw2();
141  h_csi_energyVrs3W[i]->Sumw2();
142  h_csi_edep[i]->Sumw2();
143  h_csi_edep1[i]->Sumw2();
144  h_csi_edep2[i]->Sumw2();
145  h_csi_edep1Weight[i]->Sumw2();
146  h_csi_edep2Weight[i]->Sumw2();
147  h_csi_rs_edep1[i]->Sumw2();
148  h_csi_rs_edep2[i]->Sumw2();
149  h_csi_rs_edep1Weight[i]->Sumw2();
150  h_csi_rs_edep2Weight[i]->Sumw2();
151  }
152  */
153 }
154 
155 
156 void CsiStudy_v2Module::initialize()
157 {
158  B2INFO("CsiStudy_v2Module: Initialize");
159 
160  //load
161  MetaHits.isRequired();
162  SimHits.isOptional();
163  Hits.isOptional();
164 
165  //read csi xml file
166  getXMLData();
167 
168  REG_HISTOGRAM
169 
170 }
171 
172 void CsiStudy_v2Module::beginRun()
173 {
174 }
175 
176 void CsiStudy_v2Module::event()
177 {
178  //Skip events with no Hits
179  if (SimHits.getEntries() == 0) {
180  return;
181  }
182 
183  //Look at the meta data to extract IR rate and scattering ring section
184  //double rate = 0;
185  int ring_section = -1;
186  int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
187  for (const auto& MetaHit : MetaHits) {
188  //rate = MetaHit.getrate();
189  double sad_ssraw = MetaHit.getssraw();
190  double ssraw = 0;
191  if (sad_ssraw >= 0) ssraw = sad_ssraw / 100.;
192  else ssraw = 3000. + sad_ssraw / 100.;
193  //else if (sad_ssraw < 0) ssraw = 3000. + sad_ssraw / 100.;
194  ring_section = section_ordering[(int)((ssraw) / 250.)] - 1;
195  //ring_section = MetaHit.getring_section() - 1;
196  }
197  /*
198  //Loop over SimHit
199  for (const auto& SimHit : SimHits) {
200  int detNB = SimHit.getCellId();
201  int pdg = SimHit.getPDGCode();
202  double Edep = SimHit.getEnergyDep() * 1e3; //GeV -> MeV
203  double tof = SimHit.getFlightTime(); //ns
204 
205  h_csi_edep[detNB]->Fill(Edep);
206  h_csi_Evtof[detNB]->Fill(tof, Edep);
207  if (pdg == Const::photon.getPDGCode()) h_csi_Evtof1[detNB]->Fill(tof, Edep);
208  else if (fabs(pdg) == Const::electron.getPDGCode()) h_csi_Evtof2[detNB]->Fill(tof, Edep);
209  h_csi_energy1[detNB]->Fill(log10(Edep));
210  h_csi_energy1W[detNB]->Fill(log10(Edep), rate);
211  h_csi_energyVrs1[detNB]->Fill(log10(Edep), ring_section);
212  h_csi_energyVrs1W[detNB]->Fill(log10(Edep), ring_section, rate);
213  }
214  */
215  //Loop over DigiHit
216  for (const auto& Hit : Hits) {
217  int detNB = Hit.getCellId();
218  double Edep = Hit.getEnergyDep() * 1e3; //GeV -> MeV
219  double RecEdep = Edep;//Hit.getEnergyRecDep() * 1e3; //GeV -> MeV
220  //double tof = Hit.getFlightTime(); //ns
221  /*
222  h_csi_rate[0]->Fill(detNB);
223  h_csi_rate[1]->Fill(detNB, rate);
224  h_csi_edep1[detNB]->Fill(Edep);
225  h_csi_edep2[detNB]->Fill(RecEdep);
226  h_csi_edep1Weight[detNB]->Fill(Edep, rate);
227  h_csi_edep2Weight[detNB]->Fill(RecEdep, rate);
228  h_csi_Evtof3[detNB]->Fill(tof, RecEdep);
229  h_csi_rs_rate[0]->Fill(detNB, ring_section);
230  h_csi_rs_rate[1]->Fill(detNB, ring_section, rate);
231  h_csi_rs_edep1[detNB]->Fill(Edep, ring_section);
232  h_csi_rs_edep2[detNB]->Fill(RecEdep, ring_section);
233  h_csi_rs_edep1Weight[detNB]->Fill(Edep, ring_section, rate);
234  h_csi_rs_edep2Weight[detNB]->Fill(RecEdep, ring_section, rate);
235  h_csi_energy2[detNB]->Fill(log10(Edep));
236  h_csi_energy3[detNB]->Fill(log10(RecEdep));
237  h_csi_energy2W[detNB]->Fill(log10(Edep), rate);
238  h_csi_energy3W[detNB]->Fill(log10(RecEdep), rate);
239  h_csi_energyVrs2[detNB]->Fill(log10(Edep), ring_section);
240  h_csi_energyVrs3[detNB]->Fill(log10(RecEdep), ring_section);
241  h_csi_energyVrs2W[detNB]->Fill(log10(Edep), ring_section, rate);
242  h_csi_energyVrs3W[detNB]->Fill(log10(RecEdep), ring_section, rate);
243  */
244  for (int i = 0; i < 153; i ++) {
245  /*cout << " thres1 " << m_Thres_hitRate[detNB][i]
246  << " thres2 " << m_Thres_sumE[detNB][i]
247  << " edep " << RecEdep << " true Edep " << Edep << endl;*/
248  if (RecEdep >= m_Thres_hitRate[detNB][i] && m_Thres_hitRate[detNB][i] > 0) {
249  h_csi_drate[i]->Fill(detNB);
250  h_csi_rs_drate[i]->Fill(detNB, ring_section);
251  }
252  if (RecEdep >= m_Thres_sumE[detNB][i] && m_Thres_sumE[detNB][i] > 0) {
253  h_csi_dedep[detNB][i]->Fill(RecEdep);
254  h_csi_rs_dedep[detNB][i]->Fill(RecEdep, ring_section);
255 
256  h_csi_denergy[detNB][i]->Fill(log10(RecEdep));
257  h_csi_rs_denergy[detNB][i]->Fill(log10(RecEdep), ring_section);
258  }
259  }
260  }
261 }
262 //read tube centers, impulse response, and garfield drift data filename from CSI.xml
263 void CsiStudy_v2Module::getXMLData()
264 {
265  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"CSI\"]/Content/");
266  m_Ethres = content.getDouble("Ethres");
267 
268 
269  for (int i = 0; i < 18; i ++) {
270  int iThres_hitRate = 0;
271  for (double Thres : content.getArray(TString::Format("thres_hitRate_det%d", i).Data(), {0})) {
272  m_Thres_hitRate[i][iThres_hitRate] = Thres;
273  iThres_hitRate++;
274  }
275  int iThres_sumE = 0;
276  for (double Thres : content.getArray(TString::Format("thres_sumE_det%d", i).Data(), {0})) {
277  m_Thres_sumE[i][iThres_sumE] = Thres;
278  iThres_sumE++;
279  }
280  }
281 
282  B2INFO("CsiStudy_v2");
283 
284 }
285 void CsiStudy_v2Module::endRun()
286 {
287 
288 
289 }
290 
291 void CsiStudy_v2Module::terminate()
292 {
293 }
294 
295 
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:31
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::csi::CsiStudy_v2Module
Study module for Csis (BEAST)
Definition: CsiStudy_v2Module.h:33
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
Belle2::Hit
Structure to hold some of the calpulse data.
Definition: TOPTimeBaseCalibratorModule.h:30