Belle II Software  release-05-02-19
TPCStudyModule.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/TPCStudyModule.h>
12 #include <beast/microtpc/dataobjects/MicrotpcSimHit.h>
13 #include <framework/datastore/StoreArray.h>
14 #include <framework/logging/Logger.h>
15 #include <framework/gearbox/GearDir.h>
16 #include <boost/foreach.hpp>
17 
18 #include <iostream>
19 #include <string>
20 
21 // ROOT
22 #include <TVector3.h>
23 #include <TH1.h>
24 #include <TH2.h>
25 
26 int ctr = 0;
27 int co_ctr[4] = {0, 0, 0, 0};
28 int coe_ctr[4] = {0, 0, 0, 0};
29 int h1_ctr[4] = {0, 0, 0, 0};
30 int n_ctr[4] = {0, 0, 0, 0};
31 int he4_ctr[4] = {0, 0, 0, 0};
32 int o16_ctr[4] = {0, 0, 0, 0};
33 int c12_ctr[4] = {0, 0, 0, 0};
34 int ctr_ele[4] = {0, 0, 0, 0};
35 int ctr_pos[4] = {0, 0, 0, 0};
36 int ctr_pro[4] = {0, 0, 0, 0};
37 int ctr_neu[4] = {0, 0, 0, 0};
38 int ctr_good_neu[4] = {0, 0, 0, 0};
39 int ctr_bad_neu[4] = {0, 0, 0, 0};
40 int ctr_bak[4] = {0, 0, 0, 0};
41 
42 using namespace std;
43 
44 using namespace Belle2;
45 using namespace microtpc;
46 
47 //-----------------------------------------------------------------
48 // Register the Module
49 //-----------------------------------------------------------------
50 REG_MODULE(TPCStudy)
51 
52 //-----------------------------------------------------------------
53 // Implementation
54 //-----------------------------------------------------------------
55 
56 TPCStudyModule::TPCStudyModule() : HistoModule()
57 {
58  // Set module properties
59  setDescription("Study module for Microtpcs (BEAST)");
60 
61  //Default values are set here. New values can be in MICROTPC.xml.
62  addParam("ChipRowNb", m_ChipRowNb, "Chip number of row", 226);
63  addParam("ChipColumnNb", m_ChipColumnNb, "Chip number of column", 80);
64  addParam("ChipColumnX", m_ChipColumnX, "Chip x dimension in cm / 2", 1.0);
65  addParam("ChipRowY", m_ChipRowY, "Chip y dimension in cm / 2", 0.86);
66  addParam("z_DG", m_z_DG, "Drift gap distance [cm]", 12.0);
67 }
68 
69 TPCStudyModule::~TPCStudyModule()
70 {
71 }
72 
73 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
74 void TPCStudyModule::defineHisto()
75 {
76  for (int i = 0; i < 11; i ++) {
77  h_tpc_kin[i] = new TH1F(TString::Format("tpc_kin_%d", i), "", 1000, 0., 100.);
78  h_tpc_xy[i] = new TH2F(TString::Format("tpc_xy_%d", i), "", 300, -3., 3., 300, -3., 3.);
79  }
80 
81 }
82 
83 void TPCStudyModule::initialize()
84 {
85  B2INFO("TPCStudyModule: Initialize");
86 
87  //read microtpc xml file
88  getXMLData();
89 
90  REG_HISTOGRAM
91 
92 }
93 
94 void TPCStudyModule::beginRun()
95 {
96 }
97 
98 void TPCStudyModule::event()
99 {
100  //Here comes the actual event processing
101 
102  StoreArray<MicrotpcSimHit> SimHits;
103  //StoreArray<MicrotpcHit> Hits;
104  //StoreArray<MicrotpcRecoTrack> Tracks;
105  //StoreArray<TPCG4TrackInfo> mcparts;
106  //StoreArray<SADMetaHit> sadMetaHits;
107 
108  int old_trkID[4] = { -1, -1, -1, -1};
109  int old_trkID_h1[4] = { -1, -1, -1, -1};
110  int old_trkID_he4[4] = { -1, -1, -1, -1};
111  int old_trkID_c12[4] = { -1, -1, -1, -1};
112  int old_trkID_o16[4] = { -1, -1, -1, -1};
113  bool aneu[4] = {false, false, false, false};
114  bool apro[4] = {false, false, false, false};
115  bool atrk[4] = {false, false, false, false};
116  bool aele[4] = {false, false, false, false};
117  bool apos[4] = {false, false, false, false};
118 
119  vector <double> RecoilE;
120  vector <int> Pdg[4];
121  int NbofPart[4] = {0, 0, 0, 0};
122 
123  cout << "Look at evt " << ctr << endl;
124 
125  for (const auto& MicrotpcSimHit : SimHits) {
126  int detNb = MicrotpcSimHit.getdetNb();
127  int pdg = MicrotpcSimHit.gettkPDG();
128  int trkID = MicrotpcSimHit.gettkID();
129  TVector3 position = MicrotpcSimHit.gettkPos();
130  TVector3 direction = MicrotpcSimHit.gettkMomDir();
131  double xpos = position.X() / 100. - TPCCenter[detNb].X();
132  double ypos = position.Y() / 100. - TPCCenter[detNb].Y();
133  double zpos = position.Z() / 100. - TPCCenter[detNb].Z() + m_z_DG / 2.;
134  //if (ctr == 95 && detNb == 0 && (MicrotpcSimHit.getEnergyDep() * 1e6) > 250.0) {
135  if (ctr == 95 && detNb == 0 && MicrotpcSimHit.getEnergyDep() > 0) {
136  cout << " pdg " << pdg << " Energy deposited " << MicrotpcSimHit.getEnergyDep() << " zpos " << zpos << endl;
137  h_tpc_xy[7]->Fill(xpos, ypos);
138  if (pdg == 1000020040) h_tpc_xy[8]->Fill(xpos, ypos);
139  if (pdg == 2212) h_tpc_xy[9]->Fill(xpos, ypos);
140  if (pdg == 11) h_tpc_xy[10]->Fill(xpos, ypos);
141  }
142  if (old_trkID[detNb] != trkID && MicrotpcSimHit.gettkKEnergy() > 0 && MicrotpcSimHit.getEnergyDep() > 0) {
143  old_trkID[detNb] = trkID;
144  Pdg[detNb].push_back(pdg);
145  NbofPart[detNb] ++;
146  }
147 
148  if (pdg == 1000020040) {
149  //cout << "He4 detNb " << detNb << " trID " << trkID << endl;
150  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
151  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
152  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
153  if (old_trkID_he4[detNb] != trkID && MicrotpcSimHit.gettkKEnergy() > 0 && MicrotpcSimHit.getEnergyDep() > 0) {
154  old_trkID_he4[detNb] = trkID;
155  //cout << "Output alpha track direction and vertex and momentum"<<endl;
156  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
157  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
158  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
159  atrk[detNb] = true;
160  RecoilE.push_back(MicrotpcSimHit.gettkKEnergy());
161  h_tpc_kin[0]->Fill(MicrotpcSimHit.gettkKEnergy());
162  he4_ctr[detNb] ++;
163  }
164  h_tpc_xy[0]->Fill(xpos, ypos);
165  }
166 
167  if (pdg == 1000060120) {
168  //cout << "C 12 detNb " << detNb << " trID " << trkID << endl;
169  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
170  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
171  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
172  if (old_trkID_c12[detNb] != trkID && MicrotpcSimHit.gettkKEnergy() > 0 && MicrotpcSimHit.getEnergyDep() > 0) {
173  old_trkID_c12[detNb] = trkID;
174  //cout << "Output alpha track direction and vertex and momentum"<<endl;
175  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
176  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
177  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
178  RecoilE.push_back(MicrotpcSimHit.gettkKEnergy());
179  atrk[detNb] = true;
180  h_tpc_kin[1]->Fill(MicrotpcSimHit.gettkKEnergy());
181  c12_ctr[detNb] ++;
182  }
183  h_tpc_xy[1]->Fill(xpos, ypos);
184  }
185 
186  if (pdg == 1000080160) {
187  //cout << "O 16 detNb " << detNb << " trID " << trkID << endl;
188  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
189  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
190  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
191  if (old_trkID_o16[detNb] != trkID && MicrotpcSimHit.gettkKEnergy() > 0 && MicrotpcSimHit.getEnergyDep() > 0) {
192  old_trkID_o16[detNb] = trkID;
193  //cout << "Output alpha track direction and vertex and momentum"<<endl;
194  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
195  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
196  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
197  RecoilE.push_back(MicrotpcSimHit.gettkKEnergy());
198  atrk[detNb] = true;
199  h_tpc_kin[2]->Fill(MicrotpcSimHit.gettkKEnergy());
200  o16_ctr[detNb] ++;
201  }
202  h_tpc_xy[2]->Fill(xpos, ypos);
203  }
204 
205  if (pdg == 2212) {
206  //atrk[detNb] = true;
207  //cout << "He4 detNb " << detNb << " trID " << trkID << endl;
208  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
209  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
210  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
211  if (old_trkID_h1[detNb] != trkID && MicrotpcSimHit.gettkKEnergy() > 0 && MicrotpcSimHit.getEnergyDep() > 0) {
212  old_trkID_h1[detNb] = trkID;
213  //cout << "Output alpha track direction and vertex and momentum"<<endl;
214  //cout << "Direction x " << direction.X() << " y " << direction.Y() << " z " << direction.Z() << endl;
215  //cout << "Vertex x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
216  //cout << "Kinetic energy " << MicrotpcSimHit.gettkKEnergy() << endl;
217  apro[detNb] = true;
218  h_tpc_kin[3]->Fill(MicrotpcSimHit.gettkKEnergy());
219  h1_ctr[detNb] ++;
220  }
221  h_tpc_xy[3]->Fill(xpos, ypos);
222  }
223 
224  if (pdg == 2112) {
225  aneu[detNb] = true;
226  h_tpc_kin[4]->Fill(MicrotpcSimHit.gettkKEnergy());
227  h_tpc_xy[4]->Fill(xpos, ypos);
228  n_ctr[detNb] ++;
229  }
230 
231  if (pdg == 11) {
232  aele[detNb] = true;
233  h_tpc_kin[5]->Fill(MicrotpcSimHit.gettkKEnergy());
234  h_tpc_xy[5]->Fill(xpos, ypos);
235  }
236 
237  if (pdg == -11) {
238  apos[detNb] = true;
239  h_tpc_kin[6]->Fill(MicrotpcSimHit.gettkKEnergy());
240  h_tpc_xy[6]->Fill(xpos, ypos);
241  }
242  }
243  for (int i = 0; i < 4; i ++) {
244  if (apro[i] && atrk[i] && aneu[i]) co_ctr[i]++;
245  if (apro[i] && atrk[i] && aneu[i] && aele[i]) coe_ctr[i]++;
246  if (apro[i] && atrk[i]) {
247  for (int j = 0; j < (int) RecoilE.size(); j ++) {
248  h_tpc_kin[7]->Fill(RecoilE[j]);
249  }
250  ctr_pro[i] ++;
251  }
252  if (apro[i] && aneu[i]) ctr_bak[i] ++;
253  if (atrk[i] && aneu[i]) {
254  for (int j = 0; j < (int) RecoilE.size(); j ++) {
255  h_tpc_kin[8]->Fill(RecoilE[j]);
256  }
257  ctr_neu[i] ++;
258  if (NbofPart[i] == 1) ctr_good_neu[i] ++;
259  if (NbofPart[i] > 1) ctr_bad_neu[i] ++;
260  }
261  if (atrk[i] && aele[i]) {
262  for (int j = 0; j < (int) RecoilE.size(); j ++) {
263  h_tpc_kin[9]->Fill(RecoilE[j]);
264  }
265  ctr_ele[i] ++;
266  }
267  if (atrk[i] && apos[i]) {
268  for (int j = 0; j < (int) RecoilE.size(); j ++) {
269  h_tpc_kin[10]->Fill(RecoilE[j]);
270  }
271  ctr_pos[i] ++;
272  }
273  }
274 
275  for (int i = 0; i < 4; i ++) {
276  if (NbofPart[i] > 0) {
277  cout << "det # " << i << " number of particle " << NbofPart[i] << endl;
278  for (int j = 0; j < (int) Pdg[i].size(); j ++) {
279  cout << " pdg " << Pdg[i][j] << endl;
280  }
281  }
282  }
283 
284  ctr ++;
285 }
286 //read tube centers, impulse response, and garfield drift data filename from MICROTPC.xml
287 void TPCStudyModule::getXMLData()
288 {
289  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"MICROTPC\"]/Content/");
290 
291  //get the location of the tubes
292  BOOST_FOREACH(const GearDir & activeParams, content.getNodes("Active")) {
293 
294  TPCCenter.push_back(TVector3(activeParams.getLength("TPCpos_x"), activeParams.getLength("TPCpos_y"),
295  activeParams.getLength("TPCpos_z")));
296  nTPC++;
297  }
298 
299  m_ChipColumnNb = content.getInt("ChipColumnNb");
300  m_ChipRowNb = content.getInt("ChipRowNb");
301  m_ChipColumnX = content.getDouble("ChipColumnX");
302  m_ChipRowY = content.getDouble("ChipRowY");
303  m_z_DG = content.getDouble("z_DG");
304 
305  B2INFO("TpcDigitizer: Aquired tpc locations and gas parameters");
306  B2INFO(" from MICROTPC.xml. There are " << nTPC << " TPCs implemented");
307 
308 }
309 void TPCStudyModule::endRun()
310 {
311 
312  //B2RESULT("TPCStudyModule: # of p recoils: " << npHits);
313  //B2RESULT("TPCStudyModule: # of He recoils: " << nHeHits);
314  //B2RESULT("TPCStudyModule: # of O recoils: " << nOHits);
315  //B2RESULT("TPCStudyModule: # of C recoils: " << nCHits);
316  cout << " Total nb of evts " << ctr << endl;
317  for (int i = 0; i < 4; i ++) {
318  cout << "n " << n_ctr[i] << " n-recoil-p " << co_ctr[i] << " n-recoil " << ctr_neu[i] << " p-n " << ctr_bak[i] << " p-He4 " <<
319  ctr_pro[i] << " H1 " << h1_ctr[i] << " He4 " << he4_ctr[i] << " C12 " << c12_ctr[i] << " O16 " << o16_ctr[i] << " good n-recoil " <<
320  ctr_good_neu[i] << " bad n-recoil " << ctr_bad_neu[i];
321  cout << " w/ e- " << ctr_ele[i] << " w/ e+ " << ctr_pos[i] << " n-p-recoil-e- " << coe_ctr[i] << endl;
322  }
323 }
324 
325 void TPCStudyModule::terminate()
326 {
327 }
328 
329 
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
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::MicrotpcSimHit::gettkID
int gettkID() const
Return track ID.
Definition: MicrotpcSimHit.h:64
Belle2::microtpc::TPCStudyModule
Study module for TPCs (BEAST)
Definition: TPCStudyModule.h:41
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::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::MicrotpcSimHit::gettkKEnergy
float gettkKEnergy() const
Return the kinetic energy of the track.
Definition: MicrotpcSimHit.h:72
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::HistoModule
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29