Belle II Software  release-06-00-14
DosiStudyModule.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/dosi/modules/DosiStudyModule.h>
10 #include <beast/dosi/dataobjects/DosiSimHit.h>
11 #include <beast/dosi/dataobjects/DosiHit.h>
12 
13 #include <mdst/dataobjects/MCParticle.h>
14 #include <framework/datastore/StoreArray.h>
15 #include <framework/datastore/RelationIndex.h>
16 #include <framework/logging/Logger.h>
17 #include <framework/gearbox/Const.h>
18 
19 //c++
20 #include <cmath>
21 #include <string>
22 #include <fstream>
23 
24 // ROOT
25 #include <TString.h>
26 #include <TH1.h>
27 #include <TH2.h>
28 
29 #include <generators/SAD/dataobjects/SADMetaHit.h>
30 
31 using namespace std;
32 using namespace Belle2;
33 using namespace dosi;
34 
35 //-----------------------------------------------------------------
36 // Register the Module
37 //-----------------------------------------------------------------
38 REG_MODULE(DosiStudy)
39 
40 //-----------------------------------------------------------------
41 // Implementation
42 //-----------------------------------------------------------------
43 
44 DosiStudyModule::DosiStudyModule() : HistoModule()
45 {
46  // Set module properties
47  setDescription("Study module for Dosis (BEAST)");
48 }
49 
50 DosiStudyModule::~DosiStudyModule()
51 {
52 }
53 
54 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
55 void DosiStudyModule::defineHisto()
56 {
57  for (int i = 0; i < 5; i++) {
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75  h_dosi_edep0[i] = new TH1F(TString::Format("dosi_edep0_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
76  h_dosi_edep1[i] = new TH1F(TString::Format("dosi_edep1_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
77  h_dosi_edep2[i] = new TH1F(TString::Format("dosi_edep2_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
78  h_dosi_edep3[i] = new TH1F(TString::Format("dosi_edep3_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
79  h_dosi_edep4[i] = new TH1F(TString::Format("dosi_edep4_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
80  h_dosi_edep5[i] = new TH1F(TString::Format("dosi_edep5_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
81  h_dosi_edep6[i] = new TH1F(TString::Format("dosi_edep6_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
82  h_dosi_edep7[i] = new TH1F(TString::Format("dosi_edep7_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
83  h_dosi_edep8[i] = new TH1F(TString::Format("dosi_edep8_%d", i), "Energy deposited [MeV]", 50000, 0., 400.);
84 
85  h_dosi_edep0[i]->Sumw2();
86  h_dosi_edep1[i]->Sumw2();
87  h_dosi_edep2[i]->Sumw2();
88  h_dosi_edep3[i]->Sumw2();
89  h_dosi_edep4[i]->Sumw2();
90  h_dosi_edep5[i]->Sumw2();
91  h_dosi_edep6[i]->Sumw2();
92  h_dosi_edep7[i]->Sumw2();
93  h_dosi_edep8[i]->Sumw2();
94 
95  h_dosi_rs_edep0[i] = new TH2F(TString::Format("dosi_rs_edep0_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
96  h_dosi_rs_edep1[i] = new TH2F(TString::Format("dosi_rs_edep1_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
97  h_dosi_rs_edep2[i] = new TH2F(TString::Format("dosi_rs_edep2_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
98  h_dosi_rs_edep3[i] = new TH2F(TString::Format("dosi_rs_edep3_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
99  h_dosi_rs_edep4[i] = new TH2F(TString::Format("dosi_rs_edep4_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
100  h_dosi_rs_edep5[i] = new TH2F(TString::Format("dosi_rs_edep5_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
101  h_dosi_rs_edep6[i] = new TH2F(TString::Format("dosi_rs_edep6_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
102  h_dosi_rs_edep7[i] = new TH2F(TString::Format("dosi_rs_edep7_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
103  h_dosi_rs_edep8[i] = new TH2F(TString::Format("dosi_rs_edep8_%d", i), "Energy deposited [MeV]", 50000, 0., 400., 12, 0., 12.);
104 
105  h_dosi_rs_edep0[i]->Sumw2();
106  h_dosi_rs_edep1[i]->Sumw2();
107  h_dosi_rs_edep2[i]->Sumw2();
108  h_dosi_rs_edep3[i]->Sumw2();
109  h_dosi_rs_edep4[i]->Sumw2();
110  h_dosi_rs_edep5[i]->Sumw2();
111  h_dosi_rs_edep6[i]->Sumw2();
112  h_dosi_rs_edep7[i]->Sumw2();
113  h_dosi_rs_edep8[i]->Sumw2();
114  }
115 
116 }
117 
118 
119 void DosiStudyModule::initialize()
120 {
121  B2INFO("DosiStudyModule: Initialize");
122 
123  //read dosi xml file
124  //getXMLData();
125  //StoreArray<MCParticle> mcParticles;
126  //StoreArray<DosiSimHit> SimHits;
127  //RelationArray relMCSimHit(mcParticles, SimHits);
128 
129  REG_HISTOGRAM
130 
131 }
132 
133 void DosiStudyModule::beginRun()
134 {
135 }
136 
137 void DosiStudyModule::event()
138 {
139  //Here comes the actual event processing
140  StoreArray<MCParticle> mcParticles;
141  StoreArray<DosiSimHit> SimHits;
142  StoreArray<DosiHit> Hits;
143  StoreArray<SADMetaHit> MetaHits;
144 
145  //Look at the meta data to extract IR rate and scattering ring section
146 
147  //double rate = 0;
148  int ring_section = -1;
149  int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
150  for (const auto& MetaHit : MetaHits) {
151  //rate = MetaHit.getrate();
152  double sad_ssraw = MetaHit.getssraw();
153  double ssraw = 0;
154  if (sad_ssraw >= 0) ssraw = sad_ssraw / 100.;
155  else ssraw = 3000. + sad_ssraw / 100.;
156  //else if (sad_ssraw < 0) ssraw = 3000. + sad_ssraw / 100.;
157  ring_section = section_ordering[(int)((ssraw) / 250.)] - 1;
158  //ring_section = MetaHit.getring_section() - 1;
159  }
160 
161  RelationIndex<MCParticle, DosiSimHit> relMCSimHit;
162  typedef RelationIndex<MCParticle, DosiSimHit>::Element RelationElement;
163 
164  for (const auto& SimHit : SimHits) {
165  const int detNB = SimHit.getCellId();
166  const double Edep = SimHit.getEnergyDep() * 1e3; //GeV -> MeV
167  h_dosi_edep0[detNB]->Fill(Edep);
168  h_dosi_rs_edep0[detNB]->Fill(Edep, ring_section);
169  for (const RelationElement& rel : relMCSimHit.getElementsTo(SimHit)) {
170  const MCParticle* particle = rel.from;
171  const float Mass = particle->getMass();
172  const float Energy = particle->getEnergy();
173  const float Kinetic = (Energy - Mass) * 1e3; //GeV to MeV
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) == Const::electron.getPDGCode()))
179  h_dosi_edep3[detNB]->Fill(Edep);
180  if (0.1 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == Const::electron.getPDGCode()))
181  h_dosi_edep4[detNB]->Fill(Edep);
182  if (0.001 <= Kinetic && Kinetic <= 0.050)
183  h_dosi_edep5[detNB]->Fill(Edep);
184  if (pdg == Const::neutron.getPDGCode())
185  h_dosi_edep6[detNB]->Fill(Edep);
186  if (pdg == Const::neutron.getPDGCode() && (0.001 <= Kinetic && Kinetic <= 10.0))
187  h_dosi_edep7[detNB]->Fill(Edep);
188 
189 
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) == Const::electron.getPDGCode()))
194  h_dosi_rs_edep3[detNB]->Fill(Edep, ring_section);
195  if (0.1 <= Kinetic && Kinetic <= 10.0 && (fabs(pdg) == Const::electron.getPDGCode()))
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);
199  if (pdg == Const::neutron.getPDGCode())
200  h_dosi_rs_edep6[detNB]->Fill(Edep, ring_section);
201  if (pdg == Const::neutron.getPDGCode() && (0.001 <= Kinetic && Kinetic <= 10.0))
202  h_dosi_rs_edep7[detNB]->Fill(Edep, ring_section);
203  }
204  }
205 }
206 /*
207 //read tube centers, impulse response, and garfield drift data filename from DOSI.xml
208 void DosiStudyModule::getXMLData()
209 {
210  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"DOSI\"]/Content/");
211  m_Ethres = content.getDouble("Ethres");
212 
213  B2INFO("DosiStudy");
214 
215 }
216 */
217 void DosiStudyModule::endRun()
218 {
219 
220 
221 }
222 
223 void DosiStudyModule::terminate()
224 {
225 }
226 
227 
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::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Belle2::RelationElement
Class to store a single element of a relation.
Definition: RelationElement.h:25
Belle2::RelationIndex
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
Definition: RelationIndex.h:76
Belle2::RelationIndex::getElementsTo
range_to getElementsTo(const TO *to) const
Return a range of all elements pointing to the given object.
Definition: RelationIndex.h:171
Belle2::dosi::DosiStudyModule
Study module for Dosis (BEAST)
Definition: DosiStudyModule.h:32
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17