Belle II Software  release-05-01-25
He3tubeStudyModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Sam de Jong *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <beast/he3tube/modules/He3tubeStudyModule.h>
12 #include <beast/he3tube/dataobjects/He3tubeSimHit.h>
13 #include <beast/he3tube/dataobjects/He3tubeHit.h>
14 #include <beast/he3tube/dataobjects/HE3G4TrackInfo.h>
15 #include <generators/SAD/dataobjects/SADMetaHit.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/logging/Logger.h>
18 #include <cmath>
19 
20 #include <iostream>
21 #include <fstream>
22 #include <sstream>
23 #include <string>
24 
25 // ROOT
26 #include <TVector3.h>
27 #include <TH1.h>
28 #include <TH2.h>
29 #include <TMath.h>
30 
31 int eventNum = 0;
32 
33 using namespace std;
34 
35 using namespace Belle2;
36 using namespace he3tube;
37 
38 //-----------------------------------------------------------------
39 // Register the Module
40 //-----------------------------------------------------------------
41 REG_MODULE(He3tubeStudy)
42 
43 //-----------------------------------------------------------------
44 // Implementation
45 //-----------------------------------------------------------------
46 
47 He3tubeStudyModule::He3tubeStudyModule() : HistoModule()
48 {
49  // Set module properties
50  setDescription("Study module for He3tubes (BEAST)");
51 
52  // Parameter definitions
53  addParam("sampleTime", m_sampletime, "Length of sample, in us", 10000);
54 
55 }
56 
57 He3tubeStudyModule::~He3tubeStudyModule()
58 {
59 }
60 
61 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
62 void He3tubeStudyModule::defineHisto()
63 {
64  for (int i = 0 ; i < 9 ; i++) {
65  h_mche3_kinetic[i] = new TH1F(TString::Format("h_mche3_kinetic_%d", i), "MC kin. energy [GeV]", 1000, 0., 10.);
66  h_mche3_kinetic_zoom[i] = new TH1F(TString::Format("h_mche3_kinetic_zoom_%d", i), "MC kin. energy [MeV]", 1000, 0., 10.);
67  h_mche3_tvp[i] = new TH2F(TString::Format("h_mche3_tvp_%d", i), "theta v phi", 180, 0., 180., 360, -180., 180.);
68  h_mche3_tvpW[i] = new TH2F(TString::Format("h_mche3_tvpW_%d", i), "theta v phi weighted by kin", 180, 0., 180., 360, -180., 180.);
69  h_mche3_zr[i] = new TH2F(TString::Format("h_mche3_zr_%d", i), "r v z", 200, -400., 400., 200, 0., 400.);
70  h_mche3_kinetic[i]->Sumw2();
71  h_mche3_kinetic_zoom[i]->Sumw2();
72  h_mche3_tvp[i]->Sumw2();
73  h_mche3_tvpW[i]->Sumw2();
74  h_mche3_zr[i]->Sumw2();
75  }
76 
77  h_NeutronHits = new TH1F("NeutronHits", "Neutron Hits;Tube ", 4, -0.5, 3.5);
78  h_NeutronHitsWeighted = new TH1F("NeutronHitsWeighted", "Neutron Hits;Tube ", 4, -0.5, 3.5);
79  h_DefNeutronHits = new TH1F("DefNeutronHits", "Definite Neutron Hits;Tube ", 4, -0.5, 3.5);
80  h_DefNeutronHitsWeighted = new TH1F("DefNeutronHitsWeighted", "Definite Neutron Hits;Tube ", 4, -0.5, 3.5);
81  h_NeutronRate = new TH1F("NeutronRate", "Neutron Hits per second;Tube; Rate (Hz)", 4, -0.5, 3.5);
82  h_DefNeutronRate = new TH1F("DefNeutronRate", "Neutron Hits per second;Tube; Rate (Hz)", 4, -0.5, 3.5);
83 
84  h_NeutronHitsVrs = new TH2F("NeutronHitsVrs", "Neutron Hits;Tube ", 4, -0.5, 3.5, 12, 0., 12.);
85  h_NeutronHitsWeightedVrs = new TH2F("NeutronHitsWeightedVrs", "Neutron Hits;Tube ", 4, -0.5, 3.5, 12, 0., 12.);
86  h_DefNeutronHitsVrs = new TH2F("DefNeutronHitsVrs", "Definite Neutron Hits;Tube ", 4, -0.5, 3.5, 12, 0., 12.);
87  h_DefNeutronHitsWeightedVrs = new TH2F("DefNeutronHitsWeightedVrs", "Definite Neutron Hits;Tube ", 4, -0.5, 3.5, 12, 0., 12.);
88  h_NeutronRateVrs = new TH2F("NeutronRateVrs", "Neutron Hits per second;Tube; Rate (Hz)", 4, -0.5, 3.5, 12, 0., 12.);
89  h_DefNeutronRateVrs = new TH2F("DefNeutronRateVrs", "Neutron Hits per second;Tube; Rate (Hz)", 4, -0.5, 3.5, 12, 0., 12.);
90 
91  h_Edep1H3H = new TH1F("Edep1H3H" , "Energy deposited by Proton and Tritium; MeV", 100, 0.7, 0.8);
92  h_Edep1H3H_detNB = new TH1F("Edep1H3H_tube", "Energy deposited by Proton and Tritium in each tube;Tube Num; MeV", 4, -0.5, 3.5);
93  h_Edep1H = new TH1F("Edep1H" , "Energy deposited by Protons;MeV", 100, 0, 0.7);
94  h_Edep3H = new TH1F("Edep3H" , "Energy deposited by Tritiums;MeV", 100, 0, 0.4);
95  h_TotEdep = new TH1F("TotEdep" , "Total energy deposited;MeV", 100, 0, 1.5);
96  h_DetN_Edep = new TH1F("DetN_Edep" , "Energy deposited vs detector number;Tube Num;MeV", 4, -0.5, 3.5);
97 
98  h_PulseHeights_NotNeutron = new TH1F("PulseHeights_NotNeutron", "Pulse height of waveforms from non-neutron events", 100, 0, 18000);
99  h_PulseHeights_Neutron = new TH1F("PulseHeights_Neutron" , "Pulse height of waveforms from neutron events", 100, 0, 18000);
100  h_PulseHeights_DefNeutron = new TH1F("PulseHeights_DefNeutron" , "Pulse height of waveforms from definite neutron events", 100,
101  0, 18000);
102  h_PulseHeights_All = new TH1F("PulseHeights_All" , "Pulse height of waveforms from all events", 100, 0, 18000);
103 
104  h_DefNeutronHits->Sumw2();
105  h_DefNeutronRate->Sumw2();
106  h_DefNeutronHitsWeighted->Sumw2();
107  h_NeutronHits->Sumw2();
108  h_NeutronRate->Sumw2();
109  h_NeutronHitsWeighted->Sumw2();
110 
111  h_DefNeutronHitsVrs->Sumw2();
112  h_DefNeutronRateVrs->Sumw2();
113  h_DefNeutronHitsWeightedVrs->Sumw2();
114  h_NeutronHitsVrs->Sumw2();
115  h_NeutronRateVrs->Sumw2();
116  h_NeutronHitsWeightedVrs->Sumw2();
117 
118 }
119 
120 
121 void He3tubeStudyModule::initialize()
122 {
123  B2INFO("He3tubeStudyModule: Initialize");
124 
125  REG_HISTOGRAM
126 
127  //convert sample time into rate in Hz
128  rateCorrection = m_sampletime / 1e6;
129 }
130 
131 void He3tubeStudyModule::beginRun()
132 {
133 }
134 
135 void He3tubeStudyModule::event()
136 {
137  //Here comes the actual event processing
138 
139  StoreArray<He3tubeSimHit> simHits;
140  StoreArray<He3tubeHit> Hits;
141  StoreArray<HE3G4TrackInfo> mcparts;
142  StoreArray<SADMetaHit> MetaHits;
143 
144  //Look at the meta data to extract IR rate and scattering ring section
145  double rate = 0;
146  int ring_section = -1;
147  for (const auto& MetaHit : MetaHits) {
148  rate = MetaHit.getrate();
149  ring_section = MetaHit.getring_section() - 1;
150  }
151 
152  for (const auto& mcpart : mcparts) { // start loop over all Tracks
153  const double energy = mcpart.getEnergy();
154  const double mass = mcpart.getMass();
155  double kin = energy - mass;
156  const double PDG = mcpart.getPDG();
157  const TVector3 vtx = mcpart.getProductionVertex();
158  const TVector3 mom = mcpart.getMomentum();
159  double theta = mom.Theta() * TMath::RadToDeg();
160  double phi = mom.Phi() * TMath::RadToDeg();
161  double z = vtx.Z();
162  double r = sqrt(vtx.X() * vtx.X() + vtx.Y() * vtx.Y()); int partID[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
163 
164  if (PDG == 11) partID[0] = 1; //positron
165  else if (PDG == -11) partID[1] = 1; //electron
166  else if (PDG == 22) partID[2] = 1; //photon
167  else if (PDG == 2112) partID[3] = 1; //neutron
168  else if (PDG == 2212) partID[4] = 1; //proton
169  else if (PDG == 1000080160) partID[5] = 1; // O
170  else if (PDG == 1000060120) partID[6] = 1; // C
171  else if (PDG == 1000020040) partID[7] = 1; // He
172  else partID[8] = 1;
173  for (int i = 0; i < 9; i++) {
174  if (partID[i] == 1) {
175  h_mche3_kinetic[i]->Fill(kin);
176  h_mche3_kinetic_zoom[i]->Fill(kin * 1e3);
177  h_mche3_tvp[i]->Fill(theta, phi);
178  h_mche3_tvpW[i]->Fill(theta, phi, kin);
179  h_mche3_zr[i]->Fill(z, r);
180  }
181  }
182  }
183 
184  //initalize various counters
185  double edepSum = 0;
186  double edepSum_1H = 0;
187  double edepSum_3H = 0;
188  double totedepSum = 0;
189  double totedepDet[4] = {0};
190 
191  bool ContainsP[4] = {false};
192  bool Contains3H[4] = {false};
193 
194  for (int i = 0; i < simHits.getEntries(); i++) {
195  He3tubeSimHit* aHit = simHits[i];
196  double edep = aHit->getEnergyDep();
197  int detNB = aHit->getdetNb();
198  int PID = aHit->gettkPDG();
199 
200  if (detNB > 3) {
201  B2WARNING("Hit registered in undefined tube, ignoring");
202  continue;
203 
204  }
205 
206 
207  h_DetN_Edep->AddBinContent(detNB + 1, edep);
208  totedepSum = totedepSum + edep;
209  totedepDet[detNB] = totedepDet[detNB] + edep;
210 
211  //energy deposited by protons
212  if (PID == 2212) {
213  ContainsP[detNB] = true;
214  //nPhits++;
215  h_Edep1H3H_detNB->AddBinContent(detNB + 1, edep);
216  edepSum_1H = edepSum_1H + edep;
217  edepSum = edepSum + edep;
218 
219  }
220  //energy deposited by tritiums
221  if (PID == 1000010030) {
222  Contains3H[detNB] = true;
223  //n3Hhits++;
224  h_Edep1H3H_detNB->AddBinContent(detNB + 1, edep);
225  edepSum_3H = edepSum_3H + edep;
226  edepSum = edepSum + edep;
227  }
228 
229 
230  }
231 
232  if (totedepSum != 0) h_TotEdep->Fill(totedepSum);
233  if (edepSum != 0) h_Edep1H3H->Fill(edepSum);
234  if (edepSum_1H != 0) h_Edep1H->Fill(edepSum_1H);
235  if (edepSum_3H != 0) h_Edep3H->Fill(edepSum_3H);
236 
237 
238  int neutronStatus = 0;
239  int tubeNum = -1;
240 
241  //pulse heights of digitized waveforms
242  for (int i = 0; i < Hits.getEntries(); i++) {
243  He3tubeHit* aHit = Hits[i];
244  if (aHit->definiteNeutron()) { //if this is true, this hit was definitely caused by a neutron.
245  nDefiniteNeutron++;
246  h_DefNeutronHits->Fill(aHit->getdetNb());
247  h_DefNeutronHitsWeighted->Fill(aHit->getdetNb(), rate);
248  h_DefNeutronRate->Fill(aHit->getdetNb(), 1 / rateCorrection);
249  h_DefNeutronHitsVrs->Fill(aHit->getdetNb(), ring_section);
250  h_DefNeutronHitsWeightedVrs->Fill(aHit->getdetNb(), ring_section, rate);
251  h_DefNeutronRateVrs->Fill(aHit->getdetNb(), ring_section, 1 / rateCorrection);
252  h_PulseHeights_DefNeutron->Fill(aHit->getPeakV());
253  neutronStatus = 1;
254  tubeNum = aHit->getdetNb();
255  }
256 
257  if (ContainsP[aHit->getdetNb()] && Contains3H[aHit->getdetNb()]) {
258  if (neutronStatus == 0) neutronStatus = 2;
259  tubeNum = aHit->getdetNb();
260  h_PulseHeights_Neutron->Fill(aHit->getPeakV());
261  h_NeutronHits->Fill(aHit->getdetNb());
262  h_NeutronRate->Fill(aHit->getdetNb(), 1 / rateCorrection);
263  h_NeutronHitsWeighted->Fill(aHit->getdetNb(), rate);
264  h_NeutronHitsVrs->Fill(aHit->getdetNb(), ring_section);
265  h_NeutronRateVrs->Fill(aHit->getdetNb(), ring_section, 1 / rateCorrection);
266  h_NeutronHitsWeightedVrs->Fill(aHit->getdetNb(), ring_section, rate);
267  nNeutronHits++;
268  } else {
269  h_PulseHeights_NotNeutron->Fill(aHit->getPeakV());
270  }
271  h_PulseHeights_All->Fill(aHit->getPeakV());
272 
273  }
274 
275  if (neutronStatus == 1) B2DEBUG(80, "He3tubeStudyModule: Definite Neutron in tube #" << tubeNum);
276  else if (neutronStatus == 2) B2DEBUG(80, "He3tubeStudyModule: Possible Neutron in tube #" << tubeNum);
277 
278  for (int i = 0; i < 4; i++) {
279  if (ContainsP[i]) nPhits++;
280  if (Contains3H[i]) n3Hhits++;
281  }
282 
283  eventNum++;
284 
285 
286 
287 }
288 
289 void He3tubeStudyModule::endRun()
290 {
291 
292  B2RESULT("He3tubeStudyModule: # of neutrons: " << nNeutronHits);
293  B2RESULT(" # of definite neutrons " << nDefiniteNeutron);
294  B2RESULT(" # of 3H hits: " << n3Hhits);
295  B2RESULT(" # of H hits: " << nPhits);
296 
297 
298 
299 }
300 
301 void He3tubeStudyModule::terminate()
302 {
303 }
304 
305 
Belle2::He3tubeHit::getPeakV
double getPeakV() const
Return peak.
Definition: He3tubeHit.h:62
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::He3tubeHit::getdetNb
int getdetNb() const
Return the tube number.
Definition: He3tubeHit.h:60
Belle2::He3tubeSimHit::getEnergyDep
float getEnergyDep() const
Return the energy deposition in electrons.
Definition: He3tubeSimHit.h:65
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::He3tubeSimHit::gettkPDG
int gettkPDG() const
Return the PDG number of the track.
Definition: He3tubeSimHit.h:69
Belle2::he3tube::He3tubeStudyModule
Study module for He3tubes (BEAST)
Definition: He3tubeStudyModule.h:42
Belle2::He3tubeSimHit
ClassHe3tubeSimHit - Geant4 simulated hit for the He3tube detector.
Definition: He3tubeSimHit.h:40
Belle2::He3tubeSimHit::getdetNb
float getdetNb() const
Return the He3tube number.
Definition: He3tubeSimHit.h:73
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::He3tubeHit
ClassHe3Hit - digitization simulated hit for the He3tube detector.
Definition: He3tubeHit.h:36
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::He3tubeHit::definiteNeutron
bool definiteNeutron() const
true if this is definitely a neutron event
Definition: He3tubeHit.h:66