11 #include <beast/beamabort/modules/BeamabortStudyModule.h>
12 #include <beast/beamabort/dataobjects/BeamabortSimHit.h>
13 #include <beast/beamabort/dataobjects/BeamabortHit.h>
14 #include <framework/datastore/StoreArray.h>
15 #include <framework/logging/Logger.h>
16 #include <framework/gearbox/GearDir.h>
27 #include <generators/SAD/dataobjects/SADMetaHit.h>
32 using namespace beamabort;
46 setDescription(
"Study module for Beamaborts (BEAST)");
49 addParam(
"WorkFunction", m_WorkFunction,
"Convert eV to e [e/eV] ", 1.12);
50 addParam(
"FanoFactor", m_FanoFactor,
"e resolution ", 0.1);
53 BeamabortStudyModule::~BeamabortStudyModule()
57 void BeamabortStudyModule::defineHisto()
59 B2INFO(
"BeamabortStudyModule: Initialize");
61 for (
int i = 0; i < 8; i++) {
62 h_dia_dose[i] =
new TH1F(TString::Format(
"dia_dose_%d", i),
"", 10000, 0., 10000.);
63 h_dia_dose[i]->Sumw2();
71 void BeamabortStudyModule::initialize()
78 for (
int i = 0; i < 4; i++) {
79 h_dia_rate[i] =
new TH1F(TString::Format(
"dia_rate_%d", i),
"Count", 8, 0., 8.);
80 h_dia_rate[i]->Sumw2();
82 for (
int i = 0; i < 2; i++) {
83 h_dia_rs_rate[i] =
new TH2F(TString::Format(
"dia_rs_rate_%d", i),
"Count vs ring section", 8, 0., 8., 12, 0., 12.);
84 h_dia_rs_rate[i]->Sumw2();
87 for (
int i = 0; i < 8; i++) {
88 h_dia_doseWeight[i] =
new TH1F(TString::Format(
"dia_doseWeight_%d", i),
"", 10000, 0., 10000.);
89 h_dia_amp[i] =
new TH1F(TString::Format(
"dia_amp_%d", i),
"", 10000, 0., 10000.);
90 h_dia_time[i] =
new TH1F(TString::Format(
"dia_time_%d", i),
"", 1000, 0., 100.);
91 h_dia_vtime[i] =
new TH1F(TString::Format(
"dia_vtime_%d", i),
"", 1000, 0., 100.);
92 h_dia_idose[i] =
new TH1F(TString::Format(
"dia_idose_%d", i),
"", 10000, 0., 10000.);
93 h_dia_idoseWeight[i] =
new TH1F(TString::Format(
"dia_idoseWeight_%d", i),
"", 10000, 0., 10000.);
94 h_dia_rs_idose[i] =
new TH2F(TString::Format(
"dia_rs_idose_%d", i),
"", 10000, 0., 10000., 12, 0., 12.);
95 h_dia_rs_idoseWeight[i] =
new TH2F(TString::Format(
"dia_rs_idoseWeight_%d", i),
"", 10000, 0., 10000., 12, 0., 12.);
96 h_dia_iamp[i] =
new TH1F(TString::Format(
"dia_iamp_%d", i),
"", 10000, 0., 10000.);
97 h_dia_itime[i] =
new TH1F(TString::Format(
"dia_itime_%d", i),
"", 1000, 0., 100.);
98 h_dia_ivtime[i] =
new TH1F(TString::Format(
"dia_ivtime_%d", i),
"", 1000, 0., 100.);
99 h_dia_Amp[i] =
new TH1F(TString::Format(
"dia_Amp_%d", i),
"", 100000, 0., 100000.);
100 h_dia_edep[i] =
new TH1F(TString::Format(
"dia_edep_%d", i),
"", 4000, 0., 4000.);
102 h_dia_doseWeight[i]->Sumw2();
103 h_dia_idose[i]->Sumw2();
104 h_dia_idoseWeight[i]->Sumw2();
105 h_dia_rs_idose[i]->Sumw2();
106 h_dia_rs_idoseWeight[i]->Sumw2();
113 void BeamabortStudyModule::beginRun()
117 void BeamabortStudyModule::event()
122 for (
int i = 0; i < 8; i++) {
137 int ring_section = -1;
138 int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
139 for (
const auto& MetaHit : MetaHits) {
140 rate = MetaHit.getrate();
141 double sad_ssraw = MetaHit.getssraw();
143 if (sad_ssraw >= 0) ssraw = sad_ssraw / 100.;
144 else if (sad_ssraw < 0) ssraw = 3000. + sad_ssraw / 100.;
145 ring_section = section_ordering[(int)((ssraw) / 250.)] - 1;
149 for (
const auto& SimHit : SimHits) {
150 int detNb = SimHit.getCellId();
152 double edep = SimHit.getEnergyDep();
153 double time = SimHit.getFlightTime();
154 double meanEl = edep / m_WorkFunction * 1e9;
155 double sigma = sqrt(m_FanoFactor * meanEl);
156 int NbEle = (int)gRandom->Gaus(meanEl, sigma);
157 double Amp = NbEle / (6.25 * 1e18);
160 h_dia_dose[detNb]->Fill(edep * 1e6);
161 h_dia_doseWeight[detNb]->Fill(edep * 1e6, rate);
162 h_dia_amp[detNb]->Fill(Amp * 1e15);
163 h_dia_time[detNb]->Fill(time);
164 h_dia_vtime[detNb]->Fill(time, Amp * 1e15);
165 h_dia_rate[0]->Fill(detNb);
166 h_dia_rate[1]->Fill(detNb, rate);
170 for (
int i = 0; i < 4; i++) {
171 if (curr[i] > 0 && Edep[i] > 0) {
172 h_dia_Amp[i]->Fill(curr[i] * 1e15);
173 h_dia_edep[i]->Fill(Edep[i] * 1e6);
177 for (
const auto&
Hit : Hits) {
178 int detNb =
Hit.getdetNb();
179 double edep =
Hit.getedep();
180 double current =
Hit.getI();
181 double time =
Hit.gettime();
182 h_dia_idose[detNb]->Fill(edep);
183 h_dia_idoseWeight[detNb]->Fill(edep, rate);
184 h_dia_rs_idose[detNb]->Fill(edep, ring_section);
185 h_dia_rs_idoseWeight[detNb]->Fill(edep, ring_section, rate);
186 h_dia_iamp[detNb]->Fill(current * 1e6);
187 h_dia_itime[detNb]->Fill(time);
188 h_dia_ivtime[detNb]->Fill(time, current * 1e6);
189 h_dia_rate[2]->Fill(detNb);
190 h_dia_rate[3]->Fill(detNb, rate);
192 h_dia_rs_rate[0]->Fill(detNb, ring_section);
193 h_dia_rs_rate[1]->Fill(detNb, ring_section, rate);
199 void BeamabortStudyModule::getXMLData()
201 GearDir content =
GearDir(
"/Detector/DetectorComponent[@name=\"BEAMABORT\"]/Content/");
211 m_WorkFunction = content.getDouble(
"WorkFunction");
212 m_FanoFactor = content.getDouble(
"FanoFactor");
214 B2INFO(
"BeamabortStudy");
217 void BeamabortStudyModule::endRun()
223 void BeamabortStudyModule::terminate()