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