Belle II Software  release-08-01-10
BgoStudyModule.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/bgo/modules/BgoStudyModule.h>
10 #include <beast/bgo/dataobjects/BgoSimHit.h>
11 #include <beast/bgo/dataobjects/BgoHit.h>
12 #include <framework/datastore/StoreArray.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
15 #include <cmath>
16 
17 #include <fstream>
18 #include <string>
19 
20 // ROOT
21 #include <TH1.h>
22 #include <TH2.h>
23 
24 #include <generators/SAD/dataobjects/SADMetaHit.h>
25 
26 using namespace std;
27 
28 using namespace Belle2;
29 using namespace bgo;
30 
31 //-----------------------------------------------------------------
32 // Register the Module
33 //-----------------------------------------------------------------
34 REG_MODULE(BgoStudy);
35 
36 //-----------------------------------------------------------------
37 // Implementation
38 //-----------------------------------------------------------------
39 
40 BgoStudyModule::BgoStudyModule() : HistoModule()
41 {
42  // Set module properties
43  setDescription("Study module for Bgos (BEAST)");
44 
45  addParam("Ethres", m_Ethres, "Energy threshold in MeV");
46 }
47 
49 {
50 }
51 
52 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
54 {
55  for (int i = 0; i < 2; i++) {
56  h_bgo_rate[i] = new TH1F(TString::Format("bgo_rate_%d", i), "count", 18, 0., 18.);
57  h_bgo_rate[i]->Sumw2();
58 
59  h_bgo_rs_rate[i] = new TH2F(TString::Format("bgo_rs_rate_%d", i), "count v. ring section", 18, 0., 18., 12, 0., 12.);
60  h_bgo_rs_rate[i]->Sumw2();
61  }
62  for (int i = 0; i < 18; i++) {
63  h_bgo_Evtof[i] = new TH2F(TString::Format("bgo_Evtof_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
64  1000, 0., 400.);
65  h_bgo_Evtof1[i] = new TH2F(TString::Format("bgo_Evtof1_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
66  1000, 0., 400.);
67  h_bgo_Evtof2[i] = new TH2F(TString::Format("bgo_Evtof2_%d", i), "Energy deposited [MeV] vs TOF [ns] - only photons", 5000, 0.,
68  1000., 1000, 0., 400.);
69  h_bgo_Evtof3[i] = new TH2F(TString::Format("bgo_Evtof3_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 5000, 0.,
70  1000., 1000, 0., 400.);
71 
72  h_bgo_edep[i] = new TH1F(TString::Format("bgo_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
73  h_bgo_edep1[i] = new TH1F(TString::Format("bgo_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
74  h_bgo_edep2[i] = new TH1F(TString::Format("bgo_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
75  h_bgo_edep1Weight[i] = new TH1F(TString::Format("bgo_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
76  h_bgo_edep2Weight[i] = new TH1F(TString::Format("bgo_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400.);
77 
78  h_bgo_rs_edep1[i] = new TH2F(TString::Format("bgo_rs_edep1_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
79  h_bgo_rs_edep2[i] = new TH2F(TString::Format("bgo_rs_edep2_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0., 12.);
80  h_bgo_rs_edep1Weight[i] = new TH2F(TString::Format("bgo_rs_edep1Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
81  12.);
82  h_bgo_rs_edep2Weight[i] = new TH2F(TString::Format("bgo_rs_edep2Weight_%d", i), "Energy deposited [MeV]", 5000, 0., 400., 12, 0.,
83  12.);
84 
85  h_bgo_edep[i]->Sumw2();
86  h_bgo_edep1[i]->Sumw2();
87  h_bgo_edep2[i]->Sumw2();
88  h_bgo_edep1Weight[i]->Sumw2();
89  h_bgo_edep2Weight[i]->Sumw2();
90  h_bgo_rs_edep1[i]->Sumw2();
91  h_bgo_rs_edep2[i]->Sumw2();
92  h_bgo_rs_edep1Weight[i]->Sumw2();
93  h_bgo_rs_edep2Weight[i]->Sumw2();
94  }
95 
96 }
97 
98 
100 {
101  B2INFO("BgoStudyModule: Initialize");
102 
103  //read bgo xml file
104  //getXMLData();
105 
106  REG_HISTOGRAM
107 
108 }
109 
111 {
112 }
113 
115 {
116  //Here comes the actual event processing
117  StoreArray<BgoSimHit> SimHits;
118  StoreArray<BgoHit> Hits;
119  StoreArray<SADMetaHit> MetaHits;
120 
121  //Skip events with no Hits
122  if (SimHits.getEntries() == 0) {
123  return;
124  }
125 
126  //Look at the meta data to extract IR rate and scattering ring section
127  double rate = 0;
128  int ring_section = -1;
129  const int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
130  for (const auto& MetaHit : MetaHits) {
131  rate = MetaHit.getrate();
132  double sad_ssraw = MetaHit.getssraw();
133  double ssraw = 0;
134  if (sad_ssraw >= 0) ssraw = sad_ssraw / 100.;
135  else ssraw = 3000. + sad_ssraw / 100.;
136  //else if (sad_ssraw < 0) ssraw = 3000. + sad_ssraw / 100.;
137  ring_section = section_ordering[(int)((ssraw) / 250.)] - 1;
138  //ring_section = MetaHit.getring_section() - 1;
139  }
140  //Loop over SimHit
141  for (const auto& SimHit : SimHits) {
142  int detNB = SimHit.getCellId();
143  int pdg = SimHit.getPDGCode();
144  double Edep = SimHit.getEnergyDep() * 1e3; //GeV -> MeV
145  double tof = SimHit.getFlightTime(); //ns
146 
147  h_bgo_edep[detNB]->Fill(Edep);
148  h_bgo_Evtof[detNB]->Fill(tof, Edep);
149  if (pdg == Const::photon.getPDGCode()) h_bgo_Evtof1[detNB]->Fill(tof, Edep);
150  else if (fabs(pdg) == Const::electron.getPDGCode()) h_bgo_Evtof2[detNB]->Fill(tof, Edep);
151  double RecEdep = Edep;
152  h_bgo_rate[0]->Fill(detNB);
153  h_bgo_rate[1]->Fill(detNB, rate);
154  h_bgo_edep1[detNB]->Fill(Edep);
155  h_bgo_edep2[detNB]->Fill(RecEdep);
156  h_bgo_edep1Weight[detNB]->Fill(Edep, rate);
157  h_bgo_edep2Weight[detNB]->Fill(RecEdep, rate);
158  h_bgo_Evtof3[detNB]->Fill(tof, RecEdep);
159  h_bgo_rs_rate[0]->Fill(detNB, ring_section);
160  h_bgo_rs_rate[1]->Fill(detNB, ring_section, rate);
161  h_bgo_rs_edep1[detNB]->Fill(Edep, ring_section);
162  h_bgo_rs_edep2[detNB]->Fill(RecEdep, ring_section);
163  h_bgo_rs_edep1Weight[detNB]->Fill(Edep, ring_section, rate);
164  h_bgo_rs_edep2Weight[detNB]->Fill(RecEdep, ring_section, rate);
165  }
166  /*
167  //Loop over DigiHit
168  for (const auto& Hit : Hits) {
169  int detNB = Hit.getCellId();
170  double Edep = Hit.getEnergyDep() * 1e3; //GeV -> MeV
171  double RecEdep = Hit.getEnergyRecDep() * 1e3; //GeV -> MeV
172  double tof = Hit.getFlightTime(); //ns
173  h_bgo_rate[0]->Fill(detNB);
174  h_bgo_rate[1]->Fill(detNB, rate);
175  h_bgo_edep1[detNB]->Fill(Edep);
176  h_bgo_edep2[detNB]->Fill(RecEdep);
177  h_bgo_edep1Weight[detNB]->Fill(Edep, rate);
178  h_bgo_edep2Weight[detNB]->Fill(RecEdep, rate);
179  h_bgo_Evtof3[detNB]->Fill(tof, RecEdep);
180  h_bgo_rs_rate[0]->Fill(detNB, ring_section);
181  h_bgo_rs_rate[1]->Fill(detNB, ring_section, rate);
182  h_bgo_rs_edep1[detNB]->Fill(Edep, ring_section);
183  h_bgo_rs_edep2[detNB]->Fill(RecEdep, ring_section);
184  h_bgo_rs_edep1Weight[detNB]->Fill(Edep, ring_section, rate);
185  h_bgo_rs_edep2Weight[detNB]->Fill(RecEdep, ring_section, rate);
186  }
187  */
188 }
189 /*
190 //read tube centers, impulse response, and garfield drift data filename from BGO.xml
191 void BgoStudyModule::getXMLData()
192 {
193  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"BGO\"]/Content/");
194  m_Ethres = content.getDouble("Ethres");
195 
196  B2INFO("BgoStudy");
197 
198 }
199 */
201 {
202 
203 
204 }
205 
207 {
208 }
209 
210 
int getPDGCode() const
PDG code.
Definition: Const.h:464
static const ParticleType photon
photon particle
Definition: Const.h:664
static const ChargedStable electron
electron particle
Definition: Const.h:650
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29
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
std::vector< Double_t > m_Ethres
reads data from BGO.xml: tube location, drift data filename, sigma of impulse response function
virtual ~BgoStudyModule()
Destructor.
TH1F * h_bgo_edep[18]
Energy deposited.
TH1F * h_bgo_edep1[18]
Energy deposited.
virtual void initialize() override
Initialize the Module.
TH2F * h_bgo_rs_edep2[18]
Energy deposited.
virtual void event() override
Event processor.
TH2F * h_bgo_rs_edep1Weight[18]
Energy deposited.
virtual void endRun() override
End-of-run action.
TH2F * h_bgo_Evtof[18]
Energy deposited vs TOF.
virtual void terminate() override
Termination action.
TH1F * h_bgo_edep2[18]
Energy deposited.
TH2F * h_bgo_rs_edep1[18]
Energy deposited.
virtual void beginRun() override
Called when entering a new run.
TH1F * h_bgo_edep1Weight[18]
Energy deposited.
TH2F * h_bgo_Evtof1[18]
Energy deposited vs TOF.
TH1F * h_bgo_edep2Weight[18]
Energy deposited.
TH2F * h_bgo_Evtof2[18]
Energy deposited vs TOF.
TH2F * h_bgo_Evtof3[18]
Energy deposited vs TOF.
TH2F * h_bgo_rs_edep2Weight[18]
Energy deposited.
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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.