Belle II Software  release-06-01-15
MicrotpcStudyModule.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/microtpc/modules/MicrotpcStudyModule.h>
10 #include <beast/microtpc/dataobjects/MicrotpcSimHit.h>
11 #include <beast/microtpc/dataobjects/TPCG4TrackInfo.h>
12 #include <beast/microtpc/dataobjects/MicrotpcHit.h>
13 #include <beast/microtpc/dataobjects/MicrotpcRecoTrack.h>
14 #include <generators/SAD/dataobjects/SADMetaHit.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/logging/Logger.h>
17 #include <framework/gearbox/GearDir.h>
18 #include <framework/gearbox/Const.h>
19 #include <cmath>
20 #include <boost/foreach.hpp>
21 
22 #include <fstream>
23 #include <string>
24 
25 // ROOT
26 #include <TVector3.h>
27 #include <TRandom.h>
28 #include <TH1.h>
29 #include <TH2.h>
30 #include <TH3.h>
31 #include <TMath.h>
32 
33 int eventNum = 0;
34 
35 using namespace std;
36 
37 using namespace Belle2;
38 using namespace microtpc;
39 
40 //-----------------------------------------------------------------
41 // Register the Module
42 //-----------------------------------------------------------------
43 REG_MODULE(MicrotpcStudy)
44 
45 //-----------------------------------------------------------------
46 // Implementation
47 //-----------------------------------------------------------------
48 
49 MicrotpcStudyModule::MicrotpcStudyModule() : HistoModule()
50 {
51  // Set module properties
52  setDescription("Study module for Microtpcs (BEAST)");
53 
54  //Default values are set here. New values can be in MICROTPC.xml.
55  addParam("ChipRowNb", m_ChipRowNb, "Chip number of row", 226);
56  addParam("ChipColumnNb", m_ChipColumnNb, "Chip number of column", 80);
57  addParam("ChipColumnX", m_ChipColumnX, "Chip x dimension in cm / 2", 1.0);
58  addParam("ChipRowY", m_ChipRowY, "Chip y dimension in cm / 2", 0.86);
59  addParam("z_DG", m_z_DG, "Drift gap distance [cm]", 12.0);
60 }
61 
62 MicrotpcStudyModule::~MicrotpcStudyModule()
63 {
64 }
65 
66 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
67 void MicrotpcStudyModule::defineHisto()
68 {
69  for (int i = 0 ; i < 6 ; i++) {
70  h_tpc_rate[i] = new TH1F(TString::Format("h_tpc_rate_%d", i), "detector #", 8, 0., 8.);
71  }
72 
73  h_mctpc_recoil[0] = new TH3F("h_mctpc_recoil_He", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
74  h_mctpc_recoilW[0] = new TH3F("h_mctpc_recoil_w_He", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
75  h_mctpc_recoil[0]->Sumw2();
76  h_mctpc_recoilW[0]->Sumw2();
77 
78  h_mctpc_recoil[1] = new TH3F("h_mctpc_recoil_O", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
79  h_mctpc_recoilW[1] = new TH3F("h_mctpc_recoil_w_O", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
80  h_mctpc_recoil[1]->Sumw2();
81  h_mctpc_recoilW[1]->Sumw2();
82 
83  h_mctpc_recoil[2] = new TH3F("h_mctpc_recoil_C", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
84  h_mctpc_recoilW[2] = new TH3F("h_mctpc_recoil_w_C", "Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
85  h_mctpc_recoil[2]->Sumw2();
86  h_mctpc_recoilW[2]->Sumw2();
87 
88 
89  for (int i = 0 ; i < 12 ; i++) {
90  h_mctpc_kinetic[i] = new TH2F(TString::Format("h_mctpc_kinetic_%d", i), "Neutron kin. energy [GeV]", 8, -0.5, 7.5, 1000, 0., 10.);
91  h_mctpc_kinetic_zoom[i] = new TH2F(TString::Format("h_mctpc_kinetic_zoom_%d", i), "Neutron kin. energy [MeV]", 8, -0.5, 7.5, 1000,
92  0., 10.);
93  h_mctpc_tvp[i] = new TH2F(TString::Format("h_mctpc_tvp_%d", i), "theta v phi", 180, 0., 180., 360, -180., 180.);
94  h_mctpc_tvpW[i] = new TH2F(TString::Format("h_mctpc_tvpW_%d", i), "theta v phi weighted by kin", 180, 0., 180., 360, -180., 180.);
95  h_mctpc_zr[i] = new TH2F(TString::Format("h_mctpc_zr_%d", i), "r v z", 200, -400., 400., 200, 0., 400.);
96  h_mctpc_kinetic[i]->Sumw2();
97  h_mctpc_kinetic_zoom[i]->Sumw2();
98  h_mctpc_tvp[i]->Sumw2();
99  h_mctpc_tvpW[i]->Sumw2();
100  h_mctpc_zr[i]->Sumw2();
101  }
102  for (int i = 0 ; i < 8 ; i++) {
103  for (int j = 0; j < 12; j++) {
104  h_Wtvp1[i][j] = new TH2F(TString::Format("h_Wtvp1_%d_%d", i, j), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
105  h_Wtvp2[i][j] = new TH2F(TString::Format("h_Wtvp2_%d_%d", i, j), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
106  h_Wevtrl1[i][j] = new TH2F(TString::Format("h_Wevtrl1_%d_%d", i, j), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000,
107  200, 0., 6.);
108  h_Wevtrl2[i][j] = new TH2F(TString::Format("h_Wevtrl2_%d_%d", i, j), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000,
109  200, 0., 6.);
110  }
111  }
112 
113  for (int i = 0 ; i < 8 ; i++) {
114  h_z[i] = new TH1F(TString::Format("h_z_%d", i), "Charged density per cm^2", 2000, 0.0, 20.0);
115 
116  h_zr[i] = new TH2F(TString::Format("h_zr_%d", i), "Charged density vs z vs r", 100, 0, 20, 100, 0., 5.);
117 
118  h_xy[i] = new TH2F(TString::Format("h_xy_%d", i), "Charged density vs y vs x", 100, -5., 5., 100, -5., 5.);
119 
120  h_zx[i] = new TH2F(TString::Format("h_zx_%d", i), "Charged density vs x vs r", 100, 0, 20, 100, -5., 5.);
121 
122  h_zy[i] = new TH2F(TString::Format("h_zy_%d", i), "Charged density vs y vs r", 100, 0, 20, 100, -5., 5.);
123 
124  h_evtrl[i] = new TH2F(TString::Format("h_evtrl_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
125  6.);
126  h_evtrlb[i] = new TH2F(TString::Format("h_evtrlb_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
127  6.);
128  h_evtrlc[i] = new TH2F(TString::Format("h_evtrlc_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
129  6.);
130  h_evtrld[i] = new TH2F(TString::Format("h_evtrld_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
131  6.);
132 
133  h_evtrl_p[i] = new TH2F(TString::Format("h_evtrl_p_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
134  6.);
135  h_evtrl_x[i] = new TH2F(TString::Format("h_evtrl_x_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
136  6.);
137  h_evtrl_Hex[i] = new TH2F(TString::Format("h_evtrl_Hex_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
138  0., 6.);
139  h_evtrl_He[i] = new TH2F(TString::Format("h_evtrl_He_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
140  0., 6.);
141  h_evtrl_C[i] = new TH2F(TString::Format("h_evtrl_C_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
142  6.);
143  h_evtrl_O[i] = new TH2F(TString::Format("h_evtrl_O_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
144  6.);
145  h_evtrl_He_pure[i] = new TH2F(TString::Format("h_evtrl_He_pure_%d", i), "Deposited energy [keV] v. track length [cm]", 2000, 0.,
146  2000, 200, 0., 6.);
147 
148 
149  h_tevtrl[i] = new TH2F(TString::Format("h_tevtrl_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
150  6.);
151  h_tevtrl_p[i] = new TH2F(TString::Format("h_tevtrl_p_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
152  0.,
153  6.);
154  h_tevtrl_x[i] = new TH2F(TString::Format("h_tevtrl_x_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
155  0.,
156  6.);
157  h_tevtrl_Hex[i] = new TH2F(TString::Format("h_tevtrl_Hex_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000,
158  200,
159  0., 6.);
160  h_tevtrl_He[i] = new TH2F(TString::Format("h_tevtrl_He_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000,
161  200,
162  0., 6.);
163  h_tevtrl_C[i] = new TH2F(TString::Format("h_tevtrl_C_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
164  0.,
165  6.);
166  h_tevtrl_O[i] = new TH2F(TString::Format("h_tevtrl_O_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
167  0.,
168  6.);
169  h_tevtrl_He_pure[i] = new TH2F(TString::Format("h_tevtrl_He_pure_%d", i), "t: Deposited energy [keV] v. track length [cm]", 2000,
170  0.,
171  2000, 200, 0., 6.);
172 
173  h_tvp[i] = new TH2F(TString::Format("h_tvp_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
174  h_wtvpb[i] = new TH2F(TString::Format("h_wtvpb_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
175  h_wtvpc[i] = new TH2F(TString::Format("h_wtvpc_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
176  h_wtvpd[i] = new TH2F(TString::Format("h_wtvpd_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
177 
178  h_tvpb[i] = new TH2F(TString::Format("h_tvpb_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
179  h_tvpc[i] = new TH2F(TString::Format("h_tvpc_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
180  h_tvpd[i] = new TH2F(TString::Format("h_tvpd_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
181  h_ttvp[i] = new TH2F(TString::Format("h_ttvp_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
182  h_wtvp[i] = new TH2F(TString::Format("h_wtvp_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
183  h_tvp_x[i] = new TH2F(TString::Format("h_tvp_x_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
184  h_ttvp_x[i] = new TH2F(TString::Format("h_ttvp_x_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
185  h_wtvp_x[i] = new TH2F(TString::Format("h_wtvp_x_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
186  h_tvp_p[i] = new TH2F(TString::Format("h_tvp_p_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
187  h_ttvp_p[i] = new TH2F(TString::Format("h_ttvp_p_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
188  h_wtvp_p[i] = new TH2F(TString::Format("h_wtvp_p_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
189  h_tvp_He[i] = new TH2F(TString::Format("h_tvp_He_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
190  h_ttvp_He[i] = new TH2F(TString::Format("h_ttvp_He_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
191  h_wtvp_He[i] = new TH2F(TString::Format("h_wtvp_He_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
192  h_tvp_Hex[i] = new TH2F(TString::Format("h_tvp_Hex_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
193  h_ttvp_Hex[i] = new TH2F(TString::Format("h_ttvp_Hex_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
194  h_wtvp_Hex[i] = new TH2F(TString::Format("h_wtvp_Hex_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180.,
195  180.);
196  h_tvp_He_pure[i] = new TH2F(TString::Format("h_tvp_He_pure_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
197  h_ttvp_He_pure[i] = new TH2F(TString::Format("h_ttvp_He_pure_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180.,
198  180.);
199  h_wtvp_He_pure[i] = new TH2F(TString::Format("h_wtvp_He_pure_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360,
200  -180., 180.);
201  h_twtvp_He_pure[i] = new TH2F(TString::Format("h_twtvp_He_pure_%d", i), "t: Phi [deg] v. theta [deg] - weighted", 180, 0., 180,
202  360,
203  -180., 180.);
204  h_tvp_C[i] = new TH2F(TString::Format("h_tvp_C_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
205  h_ttvp_C[i] = new TH2F(TString::Format("h_ttvp_C_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
206  h_wtvp_C[i] = new TH2F(TString::Format("h_wtvp_C_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
207  h_tvp_O[i] = new TH2F(TString::Format("h_tvp_O_%d", i), "Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
208  h_ttvp_O[i] = new TH2F(TString::Format("h_ttvp_O_%d", i), "t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
209  h_wtvp_O[i] = new TH2F(TString::Format("h_wtvp_O_%d", i), "Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
210 
211  h_tvp[i]->Sumw2();
212  h_tvpb[i]->Sumw2();
213  h_tvpc[i]->Sumw2();
214  h_tvpd[i]->Sumw2();
215 
216  h_wtvpb[i]->Sumw2();
217  h_wtvpc[i]->Sumw2();
218  h_wtvpd[i]->Sumw2();
219 
220  h_ttvp[i]->Sumw2();
221  h_tvp_x[i]->Sumw2();
222  h_ttvp_x[i]->Sumw2();
223  h_tvp_p[i]->Sumw2();
224  h_ttvp_p[i]->Sumw2();
225  h_tvp_He[i]->Sumw2();
226  h_ttvp_He[i]->Sumw2();
227  h_wtvp[i]->Sumw2();
228  h_wtvp_x[i]->Sumw2();
229  h_wtvp_p[i]->Sumw2();
230  h_wtvp_He[i]->Sumw2();
231 
232  }
233 
234 }
235 
236 
237 void MicrotpcStudyModule::initialize()
238 {
239  B2INFO("MicrotpcStudyModule: Initialize");
240 
241  //read microtpc xml file
242  getXMLData();
243 
244  REG_HISTOGRAM
245 
246  //convert sample time into rate in Hz
247  //rateCorrection = m_sampletime / 1e6;
248 }
249 
250 void MicrotpcStudyModule::beginRun()
251 {
252 }
253 
254 void MicrotpcStudyModule::event()
255 {
256  //Here comes the actual event processing
257 
258  StoreArray<MicrotpcSimHit> SimHits;
259  StoreArray<MicrotpcHit> Hits;
260  StoreArray<MicrotpcRecoTrack> Tracks;
261  StoreArray<TPCG4TrackInfo> mcparts;
262  StoreArray<SADMetaHit> sadMetaHits;
263  double rate = 0;
264  int ring_section = 0;
265  int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
266  for (const auto& sadMetaHit : sadMetaHits) {
267  rate = sadMetaHit.getrate();
268  double ss = sadMetaHit.getss() / 100.;
269  if (ss < 0) ss += 3000.;
270  int section = (int)(ss / 250.);
271  if (section >= 0 && section < 12) ring_section = section_ordering[section];
272  }
273 
274  /*
275  StoreArray<MicrotpcDataHit> DataHits;
276  int dentries = DataHits.getEntries();
277  for (int j = 0; j < dentries; j++) {
278  MicrotpcDataHit* aHit = DataHits[j];
279  int detNb = aHit->getdetNb();
280  //int trkID = aHit->gettrkID();
281  int col = aHit->getcolumn();
282  int row = aHit->getrow();
283  int tot = aHit->getTOT();
284  cout << " col " << col << " row " << row << " tot " << tot << " detNb " << detNb << endl;
285  }
286  */
287  //Bool_t EdgeCut[8];
288  double esum[8];
289  //Initialize recoil and hit type counters
290  for (int i = 0; i < 8; i++) {
291  xRec[i] = false;
292  pRec[i] = false;
293  HeRec[i] = false;
294  ORec[i] = false;
295  CRec[i] = false;
296  //ARec[i] = false;
297  pid_old[i] = 0;
298  //EdgeCut[i] = true;
299  esum[i] = 0;
300  }
301 
302  //number of entries in SimHit
303  int nSimHits = SimHits.getEntries();
304 
305  auto phiArray = new vector<double>[8](); //phi
306  auto thetaArray = new vector<double>[8](); //theta
307  auto pidArray = new vector<int>[8](); //PID
308  //auto edgeArray = new vector<int>[8](); // Edge cut
309  auto esumArray = new vector<double>[8](); // esum
310  auto trlArray = new vector<double>[8](); // trl
311 
312  TVector3 EndPoint;
313  //loop over all SimHit entries
314  for (int i = 0; i < nSimHits; i++) {
315  MicrotpcSimHit* aHit = SimHits[i];
316  int detNb = aHit->getdetNb();
317  TVector3 position = aHit->gettkPos();
318  double xpos = position.X() / 100. - TPCCenter[detNb].X();
319  double ypos = position.Y() / 100. - TPCCenter[detNb].Y();
320  double zpos = position.Z() / 100. - TPCCenter[detNb].Z() + m_z_DG / 2.;
321  if (0. < zpos && zpos < m_z_DG) {
322 
323  int PDGid = aHit->gettkPDG();
324  if (PDGid == 1000080160) ORec[detNb] = true;
325  if (PDGid == 1000060120) CRec[detNb] = true;
326  if (PDGid == 1000020040) HeRec[detNb] = true;
327  if (PDGid == Const::proton.getPDGCode()) pRec[detNb] = true;
328  if (fabs(PDGid) == Const::electron.getPDGCode() || PDGid == Const::photon.getPDGCode()) xRec[detNb] = true;
329 
330  double edep = aHit->getEnergyDep();
331  double niel = aHit->getEnergyNiel();
332  double ioni = (edep - niel) * 1e3; //MeV -> keV
333 
334  double r = sqrt(xpos * xpos + ypos * ypos);
335  h_z[detNb]->Fill(zpos, ioni);
336  h_zr[detNb]->Fill(zpos, r, ioni);
337  h_zx[detNb]->Fill(zpos, xpos, ioni);
338  h_xy[detNb]->Fill(xpos, ypos, ioni);
339  h_zy[detNb]->Fill(zpos, ypos, ioni);
340  TVector3 direction = aHit->gettkMomDir();
341  double theta = direction.Theta() * TMath::RadToDeg();
342  double phi = direction.Phi() * TMath::RadToDeg();
343 
344  if ((-m_ChipColumnX < xpos && xpos < m_ChipColumnX) &&
345  (-m_ChipRowY < ypos && ypos < m_ChipRowY) &&
346  (0. < zpos && zpos < m_z_DG)) {
347  //edgeArray].push_back(1);
348  } else {
349  //edgeArray[i].push_back(0);
350  //EdgeCut[detNb] = false;
351  }
352 
353  if (pid_old[detNb] != PDGid) {
354  if (esum[detNb] > 0) {
355  esumArray[detNb].push_back(esum[detNb]);
356  TVector3 BeginPoint;
357  BeginPoint.SetXYZ(xpos, ypos, zpos);
358  double trl0 = BeginPoint * direction.Unit();
359  double trl1 = EndPoint * direction.Unit();
360  trlArray[detNb].push_back(fabs(trl0 - trl1));
361  /*
362  double trl = fabs(trl0 - trl1);
363  double ioniz = esum[detNb];
364  if (PDGid == 1000080160) {
365  h_ttvp_O[detNb]->Fill(theta, phi);
366  h_tevtrl_O[detNb]->Fill(ioniz, trl);
367  }
368  if (PDGid == 1000060120) {
369  h_ttvp_C[detNb]->Fill(theta, phi);
370  h_tevtrl_C[detNb]->Fill(ioniz, trl);
371  }
372  if (PDGid == 1000020040) {
373  h_ttvp_He[detNb]->Fill(theta, phi);
374  h_tevtrl_He[detNb]->Fill(ioniz, trl);
375  }
376  if (PDGid == Const::proton.getPDGCode()) {
377  h_ttvp_p[detNb]->Fill(theta, phi);
378  h_tevtrl_p[detNb]->Fill(ioniz, trl);
379  }
380  if (fabs(PDGid) == Const::electron.getPDGCode() || PDGid == Const::photon.getPDGCode()) {
381  h_ttvp_x[detNb]->Fill(theta, phi);
382  h_tevtrl_x[detNb]->Fill(ioniz, trl);
383  }
384  h_ttvp[detNb]->Fill(theta, phi);
385  h_tevtrl[detNb]->Fill(ioniz, trl);
386 
387  if (EdgeCut[detNb]) {
388  if (PDGid == 1000020040) {
389  h_ttvp_He_pure[detNb]->Fill(theta, phi);
390  h_twtvp_He_pure[detNb]->Fill(theta, phi, ioniz);
391  h_tevtrl_He_pure[detNb]->Fill(ioni, trl);
392  }
393  }
394  */
395 
396  thetaArray[detNb].push_back(theta);
397  phiArray[detNb].push_back(phi);
398  pidArray[detNb].push_back(PDGid);
399 
400  }
401  pid_old[detNb] = PDGid;
402  esum[detNb] = 0;
403 
404  } else {
405  esum[detNb] += ioni;
406  EndPoint.SetXYZ(xpos, ypos, zpos);
407  }
408  }
409  }
410  /*
411  for (int i = 0; i < 8; i++) {
412 
413  for (int j = 0; j < (int)phiArray[i].size(); j++) {
414  //if (EdgeCut[i]) {
415  double phi = phiArray[i][j];
416  double theta = thetaArray[i][j];
417  int PDGid = pidArray[i][j];
418  double ioni = esumArray[i][j];
419  double trl = trlArray[i][j];
420  if (PDGid == 1000080160) {
421  h_ttvp_O[i]->Fill(theta, phi);
422  h_tevtrl_O[i]->Fill(ioni, trl);
423  }
424  if (PDGid == 1000060120) {
425  h_ttvp_C[i]->Fill(theta, phi);
426  h_tevtrl_C[i]->Fill(ioni, trl);
427  }
428  if (PDGid == 1000020040) {
429  h_ttvp_He[i]->Fill(theta, phi);
430  h_tevtrl_He[i]->Fill(ioni, trl);
431  }
432  if (PDGid == Const::proton.getPDGCode()) {
433  h_ttvp_p[i]->Fill(theta, phi);
434  h_tevtrl_p[i]->Fill(ioni, trl);
435  }
436  if (fabs(PDGid) == Const::electron.getPDGCode() || PDGid == Const::photon.getPDGCode()) {
437  h_ttvp_x[i]->Fill(theta, phi);
438  h_tevtrl_x[i]->Fill(ioni, trl);
439  }
440  h_ttvp[i]->Fill(theta, phi);
441  h_tevtrl[i]->Fill(ioni, trl);
442 
443  if (EdgeCut[i]) {
444  //if (PDGid == 1000020040) {
445  h_ttvp_He_pure[i]->Fill(theta, phi);
446  h_twtvp_He_pure[i]->Fill(theta, phi, ioni);
447  h_tevtrl_He_pure[i]->Fill(ioni, trl);
448  //}
449  }
450  //}
451  }
452  }
453  */
454  int trID = 0;
455  for (const auto& mcpart : mcparts) { // start loop over all Tracks
456  const double energy = mcpart.getEnergy();
457  const double mass = mcpart.getMass();
458  double kin = energy - mass;
459  const double PDG = mcpart.getPDG();
460  const TVector3 vtx = mcpart.getProductionVertex();
461  const TVector3 mom = mcpart.getMomentum();
462  double theta = mom.Theta() * TMath::RadToDeg();
463  double phi = mom.Phi() * TMath::RadToDeg();
464  double z = vtx.Z();
465  double r = sqrt(vtx.X() * vtx.X() + vtx.Y() * vtx.Y());
466  int partID[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
467  if (trID == mcpart.getTrackID()) continue;
468  else trID = mcpart.getTrackID();
469  int detNb = -1;
470  int nhit = 0;
471  for (const auto& sHit : SimHits) {
472  if (sHit.gettkID() == trID) {
473  detNb = sHit.getdetNb(); nhit++;
474  kin = sHit.gettkKEnergy() / 1000;
475  }
476  }
477 
478  // the only part that is actually used at the moment // Santelj 28.2.2019
479  if (PDG == Const::neutron.getPDGCode()) {
480  double trlen = abs(2. / TMath::Sin(mom.Theta()));
481  if (trlen > 10) trlen = 10.;
482  int irecoil = 0;
483  for (auto fract : m_maxEnFrac) { // loop over all recoils in beast/microtpc/data/MICROTPC-recoilProb.xml
484  double recoil = gRandom->Uniform(fract) * kin * 1e3; // calculate recoil energy
485  double weight = m_intProb[irecoil]->Eval(kin * 1e3) * trlen; // weight - interaction probability * track lenght
486  if (weight < 0) weight = 0;
487  h_mctpc_recoil[irecoil]->Fill(ring_section, detNb, recoil); // fill recoil energy
488  h_mctpc_recoilW[irecoil]->Fill(ring_section, detNb, recoil, weight); // fill weighted recoil energy
489  // std::cout << ring_section << " " << detNb << " " << recoil << " " << weight << std::endl;
490  irecoil++;
491  }
492  }
493  //------------------------------------------------------
494 
495  if (PDG == Const::electron.getPDGCode()) partID[0] = 1; //positron
496  else if (PDG == -Const::electron.getPDGCode()) partID[1] = 1; //electron
497  else if (PDG == Const::photon.getPDGCode()) partID[2] = 1; //photon
498  else if (PDG == Const::neutron.getPDGCode()) partID[3] = 1; //neutron
499  else if (PDG == Const::proton.getPDGCode()) partID[4] = 1; //proton
500  else if (PDG == 1000080160) partID[5] = 1; // O
501  else if (PDG == 1000060120) partID[6] = 1; // C
502  else if (PDG == 1000020040) partID[7] = 1; // He
503  else partID[8] = 1;
504 
505  if (PDG == Const::neutron.getPDGCode()) {
506 
507  if (r < 10.0) {
508  h_mctpc_kinetic[9]->Fill(detNb, kin);
509  h_mctpc_kinetic_zoom[9]->Fill(detNb, kin * 1e3);
510  h_mctpc_tvp[9]->Fill(theta, phi);
511  h_mctpc_tvpW[9]->Fill(theta, phi, kin);
512  h_mctpc_zr[9]->Fill(z, r);
513  }
514  if (r > 70.0) {
515  h_mctpc_kinetic[10]->Fill(detNb, kin);
516  h_mctpc_kinetic_zoom[10]->Fill(detNb, kin * 1e3);
517  h_mctpc_tvp[10]->Fill(theta, phi);
518  h_mctpc_tvpW[10]->Fill(theta, phi, kin);
519  h_mctpc_zr[10]->Fill(z, r);
520  }
521  }
522 
523  for (int i = 0; i < 9; i++) {
524  if (partID[i] == 1) {
525  h_mctpc_kinetic[i]->Fill(detNb, kin);
526  h_mctpc_kinetic_zoom[i]->Fill(detNb, kin * 1e3);
527  h_mctpc_tvp[i]->Fill(theta, phi);
528  h_mctpc_tvpW[i]->Fill(theta, phi, kin);
529  h_mctpc_zr[i]->Fill(z, r);
530  }
531  }
532  }
533 
534  //number of Tracks
535  //int nTracks = Tracks.getEntries();
536 
537  //loop over all Tracks
538  for (const auto& aTrack : Tracks) { // start loop over all Tracks
539  const int detNb = aTrack.getdetNb();
540  const float phi = aTrack.getphi();
541  const float theta = aTrack.gettheta();
542  const float trl = aTrack.gettrl();
543  const float tesum = aTrack.getesum();
544  const int pixnb = aTrack.getpixnb();
545  //const int time_range = aTrack.gettime_range();
546  int side[16];
547  for (int j = 0; j < 16; j++) {
548  side[j] = 0;
549  side[j] = aTrack.getside()[j];
550  }
551  Bool_t EdgeCuts = false;
552  if (side[0] == 0 && side[1] == 0 && side[2] == 0 && side[3] == 0) EdgeCuts = true;
553  Bool_t Asource = false;
554  if (side[4] == 2 && side[5] == 2) Asource = true;
555  //Bool_t Goodtrk = false;
556  //if (2.015 < trl && trl < 2.03) Goodtrk = true;
557  //Bool_t GoodAngle = false;
558  //if (88.5 < theta && theta < 91.5) GoodAngle = true;
559  int partID[7];
560  partID[0] = 1; //[0] for all events
561  for (int j = 0; j < 6; j++) partID[j + 1] = aTrack.getpartID()[j];
562 
563  if (ORec[detNb] || CRec[detNb] || HeRec[detNb])
564  h_tpc_rate[0]->Fill(detNb);
565  if (pRec[detNb])
566  h_tpc_rate[1]->Fill(detNb);
567  if (xRec[detNb])
568  h_tpc_rate[2]->Fill(detNb);
569 
570  if (EdgeCuts) {
571  if (ORec[detNb] || CRec[detNb] || HeRec[detNb])
572  h_tpc_rate[3]->Fill(detNb);
573  if (pRec[detNb])
574  h_tpc_rate[4]->Fill(detNb);
575  if (xRec[detNb])
576  h_tpc_rate[5]->Fill(detNb);
577  }
578 
579  h_evtrl[detNb]->Fill(tesum, trl);
580  h_tvp[detNb]->Fill(theta, phi);
581  h_wtvp[detNb]->Fill(theta, phi, tesum);
582  h_Wtvp1[detNb][0]->Fill(theta, phi, rate);
583  h_Wevtrl1[detNb][0]->Fill(tesum, trl, rate);
584  h_Wtvp2[detNb][0]->Fill(theta, phi, rate * tesum);
585  //h_Wevtrl1[detNb][0]->Fill(tesum, trl, rate);
586  if (EdgeCuts && pixnb > 10. && tesum > 10.) {
587  h_evtrlb[detNb]->Fill(tesum, trl);
588  h_tvpb[detNb]->Fill(theta, phi);
589  h_wtvpb[detNb]->Fill(theta, phi, tesum);
590  h_Wtvp1[detNb][1]->Fill(theta, phi, rate);
591  h_Wevtrl1[detNb][1]->Fill(tesum, trl, rate);
592  h_Wtvp2[detNb][1]->Fill(theta, phi, rate * tesum);
593  }
594 
595  for (int j = 0; j < 7; j++) {
596  if (j == 3 && !EdgeCuts && (partID[1] == 1 || partID[2] == 1 || partID[4] == 1 || partID[5] == 1 || partID[6] == 1)) partID[j] = 0;
597  if ((j == 4 || j == 5) && !Asource) partID[j] = 0;
598  if (partID[j] == 1) {
599  h_Wtvp1[detNb][2 + j]->Fill(theta, phi, rate);
600  h_Wevtrl1[detNb][2 + j]->Fill(tesum, trl, rate);
601  h_Wtvp2[detNb][2 + j]->Fill(theta, phi, rate * tesum);
602  if (j == 0) {
603  h_evtrlc[detNb]->Fill(tesum, trl);
604  h_tvpc[detNb]->Fill(theta, phi);
605  h_wtvpc[detNb]->Fill(theta, phi, tesum);
606  }
607  if (j == 1) {
608  h_evtrld[detNb]->Fill(tesum, trl);
609  h_tvpd[detNb]->Fill(theta, phi);
610  h_wtvpd[detNb]->Fill(theta, phi, tesum);
611  }
612  if (j == 2) {
613  h_evtrl_x[detNb]->Fill(tesum, trl);
614  h_tvp_x[detNb]->Fill(theta, phi);
615  h_wtvp_x[detNb]->Fill(theta, phi, tesum);
616  }
617  if (j == 3) {
618  h_evtrl_p[detNb]->Fill(tesum, trl);
619  h_tvp_p[detNb]->Fill(theta, phi);
620  h_wtvp_p[detNb]->Fill(theta, phi, tesum);
621  }
622  if (j == 4) {
623  h_evtrl_x[detNb]->Fill(tesum, trl);
624  h_tvp_x[detNb]->Fill(theta, phi);
625  h_wtvp_x[detNb]->Fill(theta, phi, tesum);
626  }
627  if (j == 5) {
628  h_evtrl_He[detNb]->Fill(tesum, trl);
629  h_tvp_He[detNb]->Fill(theta, phi);
630  h_wtvp_He[detNb]->Fill(theta, phi, tesum);
631  }
632  if (j == 6) {
633  h_evtrl_Hex[detNb]->Fill(tesum, trl);
634  h_tvp_Hex[detNb]->Fill(theta, phi);
635  h_wtvp_Hex[detNb]->Fill(theta, phi, tesum);
636  }
637  }
638  }
639  }
640 
641  eventNum++;
642 
643  //delete
644  delete [] phiArray;
645  delete [] thetaArray;
646  delete [] pidArray;
647  //delete [] edgeArray;
648  delete [] esumArray;
649  delete [] trlArray;
650 
651 }
652 //read tube centers, impulse response, and garfield drift data filename from MICROTPC.xml
653 void MicrotpcStudyModule::getXMLData()
654 {
655  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"MICROTPC\"]/Content/");
656 
657  //get the location of the tubes
658  BOOST_FOREACH(const GearDir & activeParams, content.getNodes("Active")) {
659 
660  TPCCenter.push_back(TVector3(activeParams.getLength("TPCpos_x"), activeParams.getLength("TPCpos_y"),
661  activeParams.getLength("TPCpos_z")));
662  nTPC++;
663  }
664 
665  m_ChipColumnNb = content.getInt("ChipColumnNb");
666  m_ChipRowNb = content.getInt("ChipRowNb");
667  m_ChipColumnX = content.getDouble("ChipColumnX");
668  m_ChipRowY = content.getDouble("ChipRowY");
669  m_z_DG = content.getDouble("z_DG");
670 
671  GearDir content2 = GearDir("/Detector/DetectorComponent[@name=\"MICROTPC\"]/Content/RecoilProbability");
672  for (const GearDir& recoil : content2.getNodes("Recoil")) {
673  m_maxEnFrac.push_back(recoil.getDouble("Fraction"));
674  istringstream probstream;
675  double e, prob;
676  probstream.str(recoil.getString("Probability"));
677  TGraph* gr = new TGraph();
678  int i = 0;
679  while (probstream >> e >> prob) {
680  gr->SetPoint(i, e, prob);
681  i++;
682  }
683  m_intProb.push_back(gr);
684  }
685 
686  B2INFO("TpcDigitizer: Aquired tpc locations and gas parameters");
687  B2INFO(" from MICROTPC.xml. There are " << nTPC << " TPCs implemented");
688 
689 }
690 void MicrotpcStudyModule::endRun()
691 {
692 
693  //B2RESULT("MicrotpcStudyModule: # of p recoils: " << npHits);
694  //B2RESULT("MicrotpcStudyModule: # of He recoils: " << nHeHits);
695  //B2RESULT("MicrotpcStudyModule: # of O recoils: " << nOHits);
696  //B2RESULT("MicrotpcStudyModule: # of C recoils: " << nCHits);
697 
698 }
699 
700 void MicrotpcStudyModule::terminate()
701 {
702 }
703 
704 
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::MicrotpcSimHit
ClassMicrotpcSimHit - Geant4 simulated hit for the Microtpc detector.
Definition: MicrotpcSimHit.h:30
Belle2::MicrotpcSimHit::getEnergyDep
float getEnergyDep() const
Return the energy deposition in electrons.
Definition: MicrotpcSimHit.h:65
Belle2::MicrotpcSimHit::gettkPos
TVector3 gettkPos() const
Return the track position.
Definition: MicrotpcSimHit.h:77
Belle2::MicrotpcSimHit::getdetNb
float getdetNb() const
Return the TPC number.
Definition: MicrotpcSimHit.h:73
Belle2::MicrotpcSimHit::getEnergyNiel
float getEnergyNiel() const
Return the non-ionization energy in electrons.
Definition: MicrotpcSimHit.h:67
Belle2::MicrotpcSimHit::gettkMomDir
TVector3 gettkMomDir() const
Return the track momentum direction.
Definition: MicrotpcSimHit.h:81
Belle2::MicrotpcSimHit::gettkPDG
int gettkPDG() const
Return the PDG number of the track.
Definition: MicrotpcSimHit.h:69
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Belle2::gearbox::Interface::getLength
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:259
Belle2::gearbox::Interface::getNodes
std::vector< GearDir > getNodes(const std::string &path="") const
Get vector of GearDirs which point to all the nodes the given path evaluates to.
Definition: Interface.cc:21
Belle2::microtpc::MicrotpcStudyModule
Study module for Microtpcs (BEAST)
Definition: MicrotpcStudyModule.h:36
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