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