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