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>
41 h_CrystalSpectrum(NULL),
42 h_CrystalRadDose(NULL),
43 h_CrystalRadDoseSH(NULL),
45 h_LightYieldCrystal(NULL),
55 setDescription(
"Analyze simulations of CsI readings in BEAST. Requires HistoManager module.");
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);
61 string rootfilepath =
"output/waveforms.root";
62 addParam(
"waveformFilename", m_waveformFilename,
"Path to the root file to save waveforms", rootfilepath);
66 CsIStudyModule::~CsIStudyModule()
70 void CsIStudyModule::defineHisto()
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);
86 void CsIStudyModule::initialize()
91 m_simhits.isRequired();
92 m_digihits.isOptional();
94 fWF =
new TFile(m_waveformFilename.c_str(),
"recreate");
97 void CsIStudyModule::beginRun()
101 void CsIStudyModule::endRun()
108 void CsIStudyModule::event()
112 double eventWeight = eventMetaDataPtr->getGeneratedWeight();
114 h_Weight->Fill(eventWeight);
117 if (m_hits.getEntries() > 0) {
118 int hitNum = m_hits.getEntries();
120 double E_tmp[16] = {0};
123 double edeptodose = GeVtoJ / Mass * usInYr / Sampletime;
125 for (
int i = 0; i < hitNum; i++) {
127 CsiHit* aCsIHit = m_hits[i];
132 E_tmp[m_cellID] += edep;
135 h_CrystalSpectrum->Fill(edep, eventWeight);
136 h_CrystalEdep->Fill(m_cellID, edep * eventWeight);
137 h_CrystalRadDose->Fill(m_cellID, edep * edeptodose * eventWeight);
140 h_NhitCrystal->Fill(m_cellID, eventWeight * 1.0e9 / Sampletime);
145 if (m_simhits.getEntries() > 0) {
146 int hitNum = m_simhits.getEntries();
149 double edeptodose = GeVtoJ / Mass * usInYr / Sampletime;
151 for (
int i = 0; i < hitNum; i++) {
158 h_CrystalRadDoseSH->Fill(m_cellID, eventWeight * edep * edeptodose);
162 h_LightYieldCrystal->Fill(m_cellID, eventWeight);
168 if (m_digihits.getEntries() > 0) {
169 int hitNum = m_digihits.getEntries();
171 for (
int i = 0; i < hitNum; i++) {
177 double charge = (double) aCsIDigiHit->
getCharge();
178 double height = (double) aCsIDigiHit->
getMaxVal();
179 double baseline = (double) aCsIDigiHit->
getBaseline();
181 vector<uint16_t>* waveform = aCsIDigiHit->
getWaveform();
184 h_TrueEdep->Fill(log10(trueEdep), eventWeight);
185 h_Charge->Fill(log10(charge - baseline), eventWeight);
186 h_Height->Fill(log10(height), eventWeight);
193 if (waveform->size()) {
198 sprintf(histoTitle,
"Waveform Hit No. %i, Cell No %i", i, cellID);
199 h_Waveform->SetTitle(histoTitle);
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));
217 void CsIStudyModule::terminate()
223 bool* CsIStudyModule::readDPPStatusBits(
char data)
226 bool* bit =
new bool[8];
229 for (
int i = 0; i < 8; i++) {
230 bit[i] = ((data >> i) & 0x01);