Belle II Software  release-05-01-25
CsIStudyModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: beaulieu *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/core/HistoModule.h>
12 #include <framework/datastore/StoreArray.h>
13 #include <framework/datastore/StoreObjPtr.h>
14 #include <framework/dataobjects/EventMetaData.h>
15 #include <beast/csi/modules/CsiModule.h>
16 #include <beast/csi/dataobjects/CsiSimHit.h>
17 #include <beast/csi/dataobjects/CsiHit.h>
18 #include <beast/csi/dataobjects/CsiDigiHit.h>
19 #include <beast/csi/modules/CsIStudyModule.h>
20 
21 // ROOT
22 #include <TH1.h>
23 #include <TFile.h>
24 
25 using namespace std;
26 using namespace Belle2;
27 using namespace csi;
28 
29 //-----------------------------------------------------------------
30 // Register the Module
31 //-----------------------------------------------------------------
32 REG_MODULE(CsIStudy)
33 
34 //-----------------------------------------------------------------
35 // Implementation
36 //-----------------------------------------------------------------
37 
39  fWF(NULL),
40  h_CrystalEdep(NULL),
41  h_CrystalSpectrum(NULL),
42  h_CrystalRadDose(NULL),
43  h_CrystalRadDoseSH(NULL),
44  h_NhitCrystal(NULL),
45  h_LightYieldCrystal(NULL),
46  h_Waveform(NULL),
47  h_Gate(NULL),
48  h_Charge(NULL),
49  h_TrueEdep(NULL),
50  h_Height(NULL),
51  h_Weight(NULL)
52 
53 {
54  // Set module properties
55  setDescription("Analyze simulations of CsI readings in BEAST. Requires HistoManager module.");
56 
57  // Parameter definitions
58  addParam("nWFSamples", m_nWFSamples, "Number of samples in the saved waveforms", 8000);
59  addParam("paramTemplate", m_paramTemplate, "Template of an input parameter. Noop for now.", 0.0);
60 
61  string rootfilepath = "output/waveforms.root";
62  addParam("waveformFilename", m_waveformFilename, "Path to the root file to save waveforms", rootfilepath);
63 
64 }
65 
66 CsIStudyModule::~CsIStudyModule()
67 {
68 }
69 
70 void CsIStudyModule::defineHisto()
71 {
72  h_CrystalEdep = new TH1F("Crystal_Edep", "Energy distribution in each crystal;CellID", 16, -0.5, 15.5);
73  h_CrystalSpectrum = new TH1F("Crystal_Edist", "Photon energy distribution in all crystals;", 100, 0, 0.1);
74  h_NhitCrystal = new TH1F("Crystal_HitRate", "Number of hit per crystal;CellID; hit/s", 16, -0.5, 15.5);
75  h_CrystalRadDose = new TH1F("Crystal_RadDose", "Crystal Radiation Dose;CellID;Gy/yr", 16, -0.5, 15.5);
76  h_CrystalRadDoseSH = new TH1F("Crystal_RadDose_SH", "Crystal Radiation Dose from SimHits;CellID;Gy/yr", 16, -0.5, 15.5);
77  h_LightYieldCrystal = new TH1F("Crystal_g_yield", "Light yield each crystal;CellID;gamma/sample", 16, -0.5, 15.5);
78  h_Waveform = new TH1S("Waveform", "Recorded waveform;Time index;ADC bits", m_nWFSamples, 0, m_nWFSamples - 1);
79  h_Gate = new TH1C("Gate", "Recorded gate;Time index;ADC bits", m_nWFSamples, 0, m_nWFSamples - 1);
80  h_Charge = new TH1F("Charge", "log_{10} of integrated Charge", 200, 0, 10);
81  h_TrueEdep = new TH1F("TrueEdep", "log_{10} of true deposited energy/GeV", 200, -6, 0);
82  h_Height = new TH1F("Height", "log10 of pulse height", 200, 0, 6);
83  h_Weight = new TH1F("Weight", "Distribution of the event weights form the generator", 1000, 0, 10);
84 }
85 
86 void CsIStudyModule::initialize()
87 {
88  REG_HISTOGRAM // required to register histograms to HistoManager
89 
90  m_hits.isOptional();
91  m_simhits.isRequired();
92  m_digihits.isOptional();
93 
94  fWF = new TFile(m_waveformFilename.c_str(), "recreate");
95 }
96 
97 void CsIStudyModule::beginRun()
98 {
99 }
100 
101 void CsIStudyModule::endRun()
102 {
103  // h_LightYieldCrystal->Divide( h_CrystalEdep );
104  // h_LightYieldCrystal->Scale( 1e-6 );
105 }
106 
107 
108 void CsIStudyModule::event()
109 {
110 
111  StoreObjPtr<EventMetaData> eventMetaDataPtr;
112  double eventWeight = eventMetaDataPtr->getGeneratedWeight();
113 
114  h_Weight->Fill(eventWeight);
115 
116  //Loop over CsiHits
117  if (m_hits.getEntries() > 0) {
118  int hitNum = m_hits.getEntries();
120  double E_tmp[16] = {0};
122  double Mass = 5;
123  double edeptodose = GeVtoJ / Mass * usInYr / Sampletime;
125  for (int i = 0; i < hitNum; i++) {
127  CsiHit* aCsIHit = m_hits[i];
128  int m_cellID = aCsIHit->getCellId();
129  double edep = aCsIHit->getEnergyDep();
130  //double hitTime = aCsIHit->getTimeAve(); /**< Time of the hit*/
131 
132  E_tmp[m_cellID] += edep;
133 
134  // Fill histograms
135  h_CrystalSpectrum->Fill(edep, eventWeight);
136  h_CrystalEdep->Fill(m_cellID, edep * eventWeight);
137  h_CrystalRadDose->Fill(m_cellID, edep * edeptodose * eventWeight);
138 
139  //Number of hits per second
140  h_NhitCrystal->Fill(m_cellID, eventWeight * 1.0e9 / Sampletime);
141  }
142  }
143 
144  //Loop over CsiSimHits
145  if (m_simhits.getEntries() > 0) {
146  int hitNum = m_simhits.getEntries();
148  double Mass = 5;
149  double edeptodose = GeVtoJ / Mass * usInYr / Sampletime;
151  for (int i = 0; i < hitNum; i++) {
153  CsiSimHit* aCsIHit = m_simhits[i];
154  int m_cellID = aCsIHit->getCellId();
155  double edep = aCsIHit->getEnergyDep();
157  // Fill histograms
158  h_CrystalRadDoseSH->Fill(m_cellID, eventWeight * edep * edeptodose);
159 
160  // To get the Number of photons per GeV (divide by total edep before plotting)
161  if (22 == aCsIHit->getPDGCode()) {
162  h_LightYieldCrystal->Fill(m_cellID, eventWeight);
163  }
164  }
165  }
166 
167  //Loop over Digihits
168  if (m_digihits.getEntries() > 0) {
169  int hitNum = m_digihits.getEntries();
171  for (int i = 0; i < hitNum; i++) {
173  CsiDigiHit* aCsIDigiHit = m_digihits[i];
174 
175  int cellID = aCsIDigiHit->getCellId();
176  // cast to double before taking the log
177  double charge = (double) aCsIDigiHit->getCharge();
178  double height = (double) aCsIDigiHit->getMaxVal();
179  double baseline = (double) aCsIDigiHit->getBaseline();
180  double trueEdep = aCsIDigiHit->getTrueEdep();
181  vector<uint16_t>* waveform = aCsIDigiHit->getWaveform();
182  vector<uint8_t>* status = aCsIDigiHit->getStatusBits();
183 
184  h_TrueEdep->Fill(log10(trueEdep), eventWeight);
185  h_Charge->Fill(log10(charge - baseline), eventWeight);
186  h_Height->Fill(log10(height), eventWeight);
187 
188  // Write the wavforms. All have the same number in the TFile, but
189  // we can access them by their cycle number:
190  // root[] TFile *f1 = new TFile("filename.root")
191  // root[] TH1F *Waveform4; f1->GetObject("h_Waveform;4",Waveform4);
192 
193  if (waveform->size()) {
194  h_Gate->Reset();
195  h_Waveform->Reset();
196 
197  char histoTitle[80];
198  sprintf(histoTitle, "Waveform Hit No. %i, Cell No %i", i, cellID);
199  h_Waveform->SetTitle(histoTitle);
200 
201 
202  for (uint iBin = 0; iBin < waveform->size(); iBin++) {
203  bool* statusBits = readDPPStatusBits(status->at(iBin));
204  h_Gate->Fill(iBin + 1, (int) statusBits[1]);
205  h_Waveform->Fill(iBin + 1, waveform->at(iBin));
206  }
207 
208  fWF->cd();
209  h_Gate->Write();
210  h_Waveform->Write();
211  }
212  }
213  }
214 }
215 
216 
217 void CsIStudyModule::terminate()
218 {
219  fWF->Close();
220 }
221 
222 
223 bool* CsIStudyModule::readDPPStatusBits(char data)
224 {
225 
226  bool* bit = new bool[8];
227 
228  // fill data
229  for (int i = 0; i < 8; i++) {
230  bit[i] = ((data >> i) & 0x01);
231  }
232 
233  return bit;
234 }
Belle2::CsiSimHit::getEnergyDep
double getEnergyDep() const
Get Deposit energy.
Definition: CsiSimHit.h:124
Belle2::CsIStudyModule
Analyze simulations of CsI readings in BEAST.
Definition: CsIStudyModule.h:51
Belle2::CsiDigiHit::getCellId
uint8_t getCellId() const
Get Cell ID.
Definition: CsiDigiHit.h:73
Belle2::CsiDigiHit::getBaseline
uint16_t getBaseline() const
Get Baseline.
Definition: CsiDigiHit.h:99
Belle2::CsiDigiHit
Class to store Csi digitized hits (output of CsIDigitizer) relation to CsiHit filled in beast/csi/mod...
Definition: CsiDigiHit.h:38
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::CsiSimHit
ClassCsiSimHit - Geant4 simulated hits in CsI crystals in BEAST.
Definition: CsiSimHit.h:41
Belle2::CsiDigiHit::getCharge
uint32_t getCharge() const
Get Integrated Charge.
Definition: CsiDigiHit.h:79
Belle2::CsiHit::getCellId
int getCellId() const
Get Cell ID.
Definition: CsiHit.h:85
Belle2::CsiDigiHit::getMaxVal
uint16_t getMaxVal() const
Get Maximal Value.
Definition: CsiDigiHit.h:107
Belle2::CsiHit::getEnergyDep
double getEnergyDep() const
Get deposit energy.
Definition: CsiHit.h:90
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CsiDigiHit::getWaveform
std::vector< uint16_t > * getWaveform()
Get Waveform array.
Definition: CsiDigiHit.h:117
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::CsiSimHit::getCellId
int getCellId() const
Get Cell ID.
Definition: CsiSimHit.h:104
Belle2::CsiDigiHit::getTrueEdep
double getTrueEdep() const
Get true deposited energy.
Definition: CsiDigiHit.h:89
Belle2::CsiSimHit::getPDGCode
int getPDGCode() const
Get Particle PDG (can be one of secondaries)
Definition: CsiSimHit.h:114
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::CsiDigiHit::getStatusBits
std::vector< uint8_t > * getStatusBits()
Get Status bits array.
Definition: CsiDigiHit.h:126
Belle2::CsiHit
Class to store simulated hits which equate to average of CsiSimHit on crystals Input for digitization...
Definition: CsiHit.h:38