Belle II Software  release-05-02-19
BeamabortStudyModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Igal Jaegle *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
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>
17 #include <cmath>
18 
19 #include <fstream>
20 #include <string>
21 
22 // ROOT
23 #include <TRandom.h>
24 #include <TH1.h>
25 #include <TH2.h>
26 
27 #include <generators/SAD/dataobjects/SADMetaHit.h>
28 
29 using namespace std;
30 
31 using namespace Belle2;
32 using namespace beamabort;
33 
34 //-----------------------------------------------------------------
35 // Register the Module
36 //-----------------------------------------------------------------
37 REG_MODULE(BeamabortStudy)
38 
39 //-----------------------------------------------------------------
40 // Implementation
41 //-----------------------------------------------------------------
42 
44 {
45  // Set module properties
46  setDescription("Study module for Beamaborts (BEAST)");
47 
48  //Default values are set here. New values can be in PINDIODE.xml.
49  addParam("WorkFunction", m_WorkFunction, "Convert eV to e [e/eV] ", 1.12);
50  addParam("FanoFactor", m_FanoFactor, "e resolution ", 0.1);
51 }
52 
53 BeamabortStudyModule::~BeamabortStudyModule()
54 {
55 }
56 
57 void BeamabortStudyModule::defineHisto()
58 {
59  B2INFO("BeamabortStudyModule: Initialize");
60 
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();
64  }
65 
66  //read beamabort xml file
67 
68 }
69 
70 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
71 void BeamabortStudyModule::initialize()
72 {
73  getXMLData();
74 
75 
76 
77  //Default values are set here. New values can be in BEAMABORT.xml.
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();
81  }
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();
85  }
86 
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.);
101 
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();
107  }
108 
109  REG_HISTOGRAM
110 }
111 
112 
113 void BeamabortStudyModule::beginRun()
114 {
115 }
116 
117 void BeamabortStudyModule::event()
118 {
119  //Here comes the actual event processing
120  double curr[8];
121  double Edep[8];
122  for (int i = 0; i < 8; i++) {
123  curr[i] = 0;
124  Edep[i] = 0;
125  }
128  StoreArray<SADMetaHit> MetaHits;
129 
130  //Skip events with no Hits
131  if (SimHits.getEntries() == 0) {
132  return;
133  }
134 
135  //Look at the meta data to extract IR rate and scattering ring section
136  double rate = 0;
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();
142  double ssraw = 0;
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;
146  //ring_section = MetaHit.getring_section() - 1;
147  }
148 
149  for (const auto& SimHit : SimHits) {
150  int detNb = SimHit.getCellId();
151  if (detNb < 8) {
152  double edep = SimHit.getEnergyDep();
153  double time = SimHit.getFlightTime();
154  double meanEl = edep / m_WorkFunction * 1e9; //GeV to eV
155  double sigma = sqrt(m_FanoFactor * meanEl); //sigma in electron
156  int NbEle = (int)gRandom->Gaus(meanEl, sigma); //electron number
157  double Amp = NbEle / (6.25 * 1e18); // A x s
158  Edep[detNb] += edep;
159  curr[detNb] += Amp;
160  h_dia_dose[detNb]->Fill(edep * 1e6); //GeV to keV
161  h_dia_doseWeight[detNb]->Fill(edep * 1e6, rate); //GeV to keV
162  h_dia_amp[detNb]->Fill(Amp * 1e15); //A x s -> mA x s
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);
167  }
168  }
169 
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);
174  }
175  }
176 
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); //keV
183  h_dia_idoseWeight[detNb]->Fill(edep, rate); //keV
184  h_dia_rs_idose[detNb]->Fill(edep, ring_section); //keV
185  h_dia_rs_idoseWeight[detNb]->Fill(edep, ring_section, rate); //keV
186  h_dia_iamp[detNb]->Fill(current * 1e6); //I to uI
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);
191 
192  h_dia_rs_rate[0]->Fill(detNb, ring_section);
193  h_dia_rs_rate[1]->Fill(detNb, ring_section, rate);
194  }
195 
196  //}
197 }
198 //read tube centers, impulse response, and garfield drift data filename from BEAMABORT.xml
199 void BeamabortStudyModule::getXMLData()
200 {
201  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"BEAMABORT\"]/Content/");
202  /*
203  //get the location of the tubes
204  BOOST_FOREACH(const GearDir & activeParams, content.getNodes("Active")) {
205 
206  PINCenter.push_back(TVector3(activeParams.getLength("z_beamabort"), activeParams.getLength("r_beamabort"),
207  activeParams.getLength("Phi")));
208  nPIN++;
209  }
210  */
211  m_WorkFunction = content.getDouble("WorkFunction");
212  m_FanoFactor = content.getDouble("FanoFactor");
213 
214  B2INFO("BeamabortStudy");
215 
216 }
217 void BeamabortStudyModule::endRun()
218 {
219 
220 
221 }
222 
223 void BeamabortStudyModule::terminate()
224 {
225 }
226 
227 
Belle2::beamabort::BeamabortStudyModule
Study module for Beamaborts (BEAST)
Definition: BeamabortStudyModule.h:39
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
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::Hit
Structure to hold some of the calpulse data.
Definition: TOPTimeBaseCalibratorModule.h:40