Belle II Software  release-05-02-19
CsiStudy_v2Module.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/csi/modules/CsiStudy_v2Module.h>
12 #include <framework/logging/Logger.h>
13 #include <framework/gearbox/GearDir.h>
14 #include <cmath>
15 
16 #include <fstream>
17 #include <string>
18 
19 // ROOT
20 #include <TH1.h>
21 #include <TH2.h>
22 
23 using namespace std;
24 
25 using namespace Belle2;
26 using namespace csi;
27 
28 //-----------------------------------------------------------------
29 // Register the Module
30 //-----------------------------------------------------------------
31 REG_MODULE(CsiStudy_v2)
32 
33 //-----------------------------------------------------------------
34 // Implementation
35 //-----------------------------------------------------------------
36 
38 {
39  // Set module properties
40  setDescription("Study_v2 module for Csis (BEAST)");
41 
42  addParam("Ethres", m_Ethres, "Energy threshold in MeV", 0.0);
43 }
44 
45 CsiStudy_v2Module::~CsiStudy_v2Module()
46 {
47 }
48 
49 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
50 void CsiStudy_v2Module::defineHisto()
51 {
52  for (int i = 0; i < 153; i++) {
53  h_csi_drate[i] = new TH1F(TString::Format("csi_rate_%d", i), "count", 18, 0., 18.);
54  h_csi_drate[i]->Sumw2();
55 
56  h_csi_rs_drate[i] = new TH2F(TString::Format("csi_rs_drate_%d", i), "count v. ring section", 18, 0., 18., 12, 0., 12.);
57  h_csi_rs_drate[i]->Sumw2();
58 
59  for (int j = 0; j < 18; j++) {
60 
61  h_csi_dedep[j][i] = new TH1F(TString::Format("csi_dedep_%d_%d", j, i), "Energy deposited [MeV]", 5000, 0., 400.);
62  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.);
63 
64  h_csi_denergy[j][i] = new TH1F(TString::Format("csi_denergy_%d_%d", j, i), "Energy deposited [MeV]", 200, 0.1, 3.);
65  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.,
66  12, 0.,
67  12.);
68 
69  h_csi_dedep[j][i]->Sumw2();
70  h_csi_rs_dedep[j][i]->Sumw2();
71  h_csi_denergy[j][i]->Sumw2();
72  h_csi_rs_denergy[j][i]->Sumw2();
73  }
74  }
75  /*
76  for (int i = 0; i < 2; i++) {
77  h_csi_rate[i] = new TH1F(TString::Format("csi_rate_%d", i), "count", 18, 0., 18.);
78  h_csi_rate[i]->Sumw2();
79 
80  h_csi_rs_rate[i] = new TH2F(TString::Format("csi_rs_rate_%d", i), "count v. ring section", 18, 0., 18., 12, 0., 12.);
81  h_csi_rs_rate[i]->Sumw2();
82  }
83  for (int i = 0; i < 18; i++) {
84  h_csi_energyVrs1[i] = new TH2F(TString::Format("csi_energyVrs1_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
85  12.);
86  h_csi_energyVrs2[i] = new TH2F(TString::Format("csi_energyVrs2_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
87  12.);
88  h_csi_energyVrs3[i] = new TH2F(TString::Format("csi_energyVrs3_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12, 0.,
89  12.);
90 
91  h_csi_energyVrs1W[i] = new TH2F(TString::Format("csi_energyVrs1W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
92  0.,
93  12.);
94  h_csi_energyVrs2W[i] = new TH2F(TString::Format("csi_energyVrs2W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
95  0.,
96  12.);
97  h_csi_energyVrs3W[i] = new TH2F(TString::Format("csi_energyVrs3W_%d", i), "Energy deposited [MeV] per section", 200, 0.1, 3., 12,
98  0.,
99  12.);
100 
101  h_csi_energy1[i] = new TH1F(TString::Format("csi_energy1_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
102  h_csi_energy2[i] = new TH1F(TString::Format("csi_energy2_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
103  h_csi_energy3[i] = new TH1F(TString::Format("csi_energy3_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
104 
105  h_csi_energy1W[i] = new TH1F(TString::Format("csi_energy1W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
106  h_csi_energy2W[i] = new TH1F(TString::Format("csi_energy2W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
107  h_csi_energy3W[i] = new TH1F(TString::Format("csi_energy3W_%d", i), "Energy deposited [MeV]", 200, 0.1, 3.);
108 
109  h_csi_Evtof[i] = new TH2F(TString::Format("csi_Evtof_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
110  1000, 0., 400.);
111  h_csi_Evtof1[i] = new TH2F(TString::Format("csi_Evtof1_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
112  1000, 0., 400.);
113  h_csi_Evtof2[i] = new TH2F(TString::Format("csi_Evtof2_%d", i), "Energy deposited [MeV] vs TOF [ns] - only photons", 5000, 0.,
114  1000., 1000, 0., 400.);
115  h_csi_Evtof3[i] = new TH2F(TString::Format("csi_Evtof3_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 5000, 0.,
116  1000., 1000, 0., 400.);
117 
118  h_csi_edep[i] = new TH1F(TString::Format("csi_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
119  h_csi_edep1[i] = new TH1F(TString::Format("csi_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
120  h_csi_edep2[i] = new TH1F(TString::Format("csi_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
121  h_csi_edep1Weight[i] = new TH1F(TString::Format("csi_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
122  h_csi_edep2Weight[i] = new TH1F(TString::Format("csi_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
123 
124  h_csi_rs_edep1[i] = new TH2F(TString::Format("csi_rs_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
125  h_csi_rs_edep2[i] = new TH2F(TString::Format("csi_rs_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
126  h_csi_rs_edep1Weight[i] = new TH2F(TString::Format("csi_rs_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
127  12.);
128  h_csi_rs_edep2Weight[i] = new TH2F(TString::Format("csi_rs_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
129  12.);
130 
131  h_csi_energy1[i]->Sumw2();
132  h_csi_energy2[i]->Sumw2();
133  h_csi_energy3[i]->Sumw2();
134  h_csi_energy1W[i]->Sumw2();
135  h_csi_energy2W[i]->Sumw2();
136  h_csi_energy3W[i]->Sumw2();
137  h_csi_energyVrs1[i]->Sumw2();
138  h_csi_energyVrs2[i]->Sumw2();
139  h_csi_energyVrs3[i]->Sumw2();
140  h_csi_energyVrs1W[i]->Sumw2();
141  h_csi_energyVrs2W[i]->Sumw2();
142  h_csi_energyVrs3W[i]->Sumw2();
143  h_csi_edep[i]->Sumw2();
144  h_csi_edep1[i]->Sumw2();
145  h_csi_edep2[i]->Sumw2();
146  h_csi_edep1Weight[i]->Sumw2();
147  h_csi_edep2Weight[i]->Sumw2();
148  h_csi_rs_edep1[i]->Sumw2();
149  h_csi_rs_edep2[i]->Sumw2();
150  h_csi_rs_edep1Weight[i]->Sumw2();
151  h_csi_rs_edep2Weight[i]->Sumw2();
152  }
153  */
154 }
155 
156 
157 void CsiStudy_v2Module::initialize()
158 {
159  B2INFO("CsiStudy_v2Module: Initialize");
160 
161  //load
162  MetaHits.isRequired();
163  SimHits.isOptional();
164  Hits.isOptional();
165 
166  //read csi xml file
167  getXMLData();
168 
169  REG_HISTOGRAM
170 
171 }
172 
173 void CsiStudy_v2Module::beginRun()
174 {
175 }
176 
177 void CsiStudy_v2Module::event()
178 {
179  //Skip events with no Hits
180  if (SimHits.getEntries() == 0) {
181  return;
182  }
183 
184  //Look at the meta data to extract IR rate and scattering ring section
185  //double rate = 0;
186  int ring_section = -1;
187  int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
188  for (const auto& MetaHit : MetaHits) {
189  //rate = MetaHit.getrate();
190  double sad_ssraw = MetaHit.getssraw();
191  double ssraw = 0;
192  if (sad_ssraw >= 0) ssraw = 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 == 22) h_csi_Evtof1[detNB]->Fill(tof, Edep);
208  else if (fabs(pdg) == 11) 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::csi::CsiStudy_v2Module
Study module for Csis (BEAST)
Definition: CsiStudy_v2Module.h:43
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
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::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