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