11 #include <beast/dosi/modules/DosiStudyModule.h>
12 #include <beast/dosi/dataobjects/DosiSimHit.h>
13 #include <beast/dosi/dataobjects/DosiHit.h>
15 #include <mdst/dataobjects/MCParticle.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/datastore/RelationIndex.h>
18 #include <framework/logging/Logger.h>
30 #include <generators/SAD/dataobjects/SADMetaHit.h>
48 setDescription(
"Study module for Dosis (BEAST)");
51 DosiStudyModule::~DosiStudyModule()
56 void DosiStudyModule::defineHisto()
58 for (
int i = 0; i < 5; i++) {
76 h_dosi_edep0[i] =
new TH1F(TString::Format(
"dosi_edep0_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
77 h_dosi_edep1[i] =
new TH1F(TString::Format(
"dosi_edep1_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
78 h_dosi_edep2[i] =
new TH1F(TString::Format(
"dosi_edep2_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
79 h_dosi_edep3[i] =
new TH1F(TString::Format(
"dosi_edep3_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
80 h_dosi_edep4[i] =
new TH1F(TString::Format(
"dosi_edep4_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
81 h_dosi_edep5[i] =
new TH1F(TString::Format(
"dosi_edep5_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
82 h_dosi_edep6[i] =
new TH1F(TString::Format(
"dosi_edep6_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
83 h_dosi_edep7[i] =
new TH1F(TString::Format(
"dosi_edep7_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
84 h_dosi_edep8[i] =
new TH1F(TString::Format(
"dosi_edep8_%d", i),
"Energy deposited [MeV]", 50000, 0., 400.);
86 h_dosi_edep0[i]->Sumw2();
87 h_dosi_edep1[i]->Sumw2();
88 h_dosi_edep2[i]->Sumw2();
89 h_dosi_edep3[i]->Sumw2();
90 h_dosi_edep4[i]->Sumw2();
91 h_dosi_edep5[i]->Sumw2();
92 h_dosi_edep6[i]->Sumw2();
93 h_dosi_edep7[i]->Sumw2();
94 h_dosi_edep8[i]->Sumw2();
96 h_dosi_rs_edep0[i] =
new TH2F(TString::Format(
"dosi_rs_edep0_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
97 h_dosi_rs_edep1[i] =
new TH2F(TString::Format(
"dosi_rs_edep1_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
98 h_dosi_rs_edep2[i] =
new TH2F(TString::Format(
"dosi_rs_edep2_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
99 h_dosi_rs_edep3[i] =
new TH2F(TString::Format(
"dosi_rs_edep3_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
100 h_dosi_rs_edep4[i] =
new TH2F(TString::Format(
"dosi_rs_edep4_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
101 h_dosi_rs_edep5[i] =
new TH2F(TString::Format(
"dosi_rs_edep5_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
102 h_dosi_rs_edep6[i] =
new TH2F(TString::Format(
"dosi_rs_edep6_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
103 h_dosi_rs_edep7[i] =
new TH2F(TString::Format(
"dosi_rs_edep7_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
104 h_dosi_rs_edep8[i] =
new TH2F(TString::Format(
"dosi_rs_edep8_%d", i),
"Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
106 h_dosi_rs_edep0[i]->Sumw2();
107 h_dosi_rs_edep1[i]->Sumw2();
108 h_dosi_rs_edep2[i]->Sumw2();
109 h_dosi_rs_edep3[i]->Sumw2();
110 h_dosi_rs_edep4[i]->Sumw2();
111 h_dosi_rs_edep5[i]->Sumw2();
112 h_dosi_rs_edep6[i]->Sumw2();
113 h_dosi_rs_edep7[i]->Sumw2();
114 h_dosi_rs_edep8[i]->Sumw2();
120 void DosiStudyModule::initialize()
122 B2INFO(
"DosiStudyModule: Initialize");
134 void DosiStudyModule::beginRun()
138 void DosiStudyModule::event()
149 int ring_section = -1;
150 int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
151 for (
const auto& MetaHit : MetaHits) {
153 double sad_ssraw = MetaHit.getssraw();
155 if (sad_ssraw >= 0) ssraw = sad_ssraw / 100.;
156 else if (sad_ssraw < 0) ssraw = 3000. + sad_ssraw / 100.;
157 ring_section = section_ordering[(int)((ssraw) / 250.)] - 1;
164 for (
const auto& SimHit : SimHits) {
165 const int detNB = SimHit.getCellId();
166 const double Edep = SimHit.getEnergyDep() * 1e3;
167 h_dosi_edep0[detNB]->Fill(Edep);
168 h_dosi_rs_edep0[detNB]->Fill(Edep, ring_section);
171 const float Mass = particle->getMass();
172 const float Energy = particle->getEnergy();
173 const float Kinetic = (Energy - Mass) * 1e3;
174 const int pdg = particle->getPDG();
175 h_dosi_edep1[detNB]->Fill(Edep);
176 if (0.005 <= Kinetic && Kinetic <= 10.0)
177 h_dosi_edep2[detNB]->Fill(Edep);
178 if (0.005 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == 11))
179 h_dosi_edep3[detNB]->Fill(Edep);
180 if (0.1 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == 11))
181 h_dosi_edep4[detNB]->Fill(Edep);
182 if (0.001 <= Kinetic && Kinetic <= 0.050)
183 h_dosi_edep5[detNB]->Fill(Edep);
185 h_dosi_edep6[detNB]->Fill(Edep);
186 if (pdg == 2112 && (0.001 <= Kinetic && Kinetic <= 10.0))
187 h_dosi_edep7[detNB]->Fill(Edep);
190 h_dosi_rs_edep1[detNB]->Fill(Edep, ring_section);
191 if (0.005 <= Kinetic && Kinetic <= 10.0)
192 h_dosi_rs_edep2[detNB]->Fill(Edep, ring_section);
193 if (0.005 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == 11))
194 h_dosi_rs_edep3[detNB]->Fill(Edep, ring_section);
195 if (0.1 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == 11))
196 h_dosi_rs_edep4[detNB]->Fill(Edep, ring_section);
197 if (0.001 <= Kinetic && Kinetic <= 0.050)
198 h_dosi_rs_edep5[detNB]->Fill(Edep, ring_section);
200 h_dosi_rs_edep6[detNB]->Fill(Edep, ring_section);
201 if (pdg == 2112 && (0.001 <= Kinetic && Kinetic <= 10.0))
202 h_dosi_rs_edep7[detNB]->Fill(Edep, ring_section);
217 void DosiStudyModule::endRun()
223 void DosiStudyModule::terminate()