REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
TRestRawReadoutAnalysisProcess.cxx
1 
17 #include "TRestRawReadoutAnalysisProcess.h"
18 
19 #include <TLegend.h>
20 #include <TPaveText.h>
21 
22 using namespace std;
23 
24 ClassImp(TRestRawReadoutAnalysisProcess);
25 
26 TRestRawReadoutAnalysisProcess::TRestRawReadoutAnalysisProcess() { Initialize(); }
27 
28 TRestRawReadoutAnalysisProcess::~TRestRawReadoutAnalysisProcess() {}
29 
30 void TRestRawReadoutAnalysisProcess::Initialize() {
31  SetSectionName(this->ClassName());
32  SetLibraryVersion(LIBRARY_VERSION);
33 
34  fSignalEvent = nullptr;
35 
36  fReadout = nullptr;
37 }
38 
39 void TRestRawReadoutAnalysisProcess::InitProcess() {
40  fReadout = GetMetadata<TRestDetectorReadout>();
41  if (fReadout != nullptr) {
42  auto iter = fModuleHitMaps.begin();
43  while (iter != fModuleHitMaps.end()) {
44  TRestDetectorReadoutModule* mod = fReadout->GetReadoutModuleWithID(iter->first);
45  if (mod == nullptr) {
46  RESTWarning << "REST Warning(TRestRawReadoutAnalysisProcess): readout "
47  "module with id "
48  << iter->first << " not found!" << RESTendl;
49  } else {
50  fModuleHitMaps[iter->first] =
51  new TH2D((TString) "Hitmap_M" + ToString(iter->first),
52  (TString) "FirstX/Y Hitmap of Module " + ToString(iter->first),
53  mod->GetNumberOfChannels() / 2, 0, mod->GetNumberOfChannels() / 2,
54  mod->GetNumberOfChannels() / 2, mod->GetNumberOfChannels() / 2,
55  mod->GetNumberOfChannels());
56 
57  fModuleActivityX[iter->first] =
58  new TH1D((TString) "ActivityX_M" + ToString(iter->first),
59  (TString) "X Channel Activity of Module " + ToString(iter->first),
60  mod->GetNumberOfChannels() / 2, 0, mod->GetNumberOfChannels() / 2);
61 
62  fModuleActivityY[iter->first] =
63  new TH1D((TString) "ActivityY_M" + ToString(iter->first),
64  (TString) "Y Channel Activity of Module " + ToString(iter->first),
65  mod->GetNumberOfChannels() / 2, mod->GetNumberOfChannels() / 2,
66  mod->GetNumberOfChannels());
67 
68  fModuleBSLSigmaX[iter->first] =
69  new TH2D((TString) "BSLSX_M" + ToString(iter->first),
70  (TString) "X Channel Baseline Sigma of Module " + ToString(iter->first),
71  mod->GetNumberOfChannels() / 2, 0, mod->GetNumberOfChannels() / 2, 100, 0, 100);
72 
73  fModuleBSLSigmaY[iter->first] =
74  new TH2D((TString) "BSLSY_M" + ToString(iter->first),
75  (TString) "Y Channel Baseline Sigma of Module " + ToString(iter->first), 100, 0,
76  100, mod->GetNumberOfChannels() / 2, mod->GetNumberOfChannels() / 2,
77  mod->GetNumberOfChannels());
78  }
79  iter++;
80  }
81  }
82 }
83 
84 TRestEvent* TRestRawReadoutAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
85  fSignalEvent = (TRestRawSignalEvent*)inputEvent;
86  if (fReadout != nullptr) {
87  Double_t firstX_id = -1.;
88  Double_t firstY_id = -1.;
89  Double_t firstX_t = 512.0;
90  Double_t firstY_t = 512.0;
91 
92  Double_t lastX_id = -1.;
93  Double_t lastY_id = -1.;
94  Double_t lastX_t = 0.0;
95  Double_t lastY_t = 0.0;
96 
97  double nan = numeric_limits<double>::quiet_NaN();
98  this->SetObservableValue("FirstX", nan);
99  this->SetObservableValue("FirstY", nan);
100  this->SetObservableValue("LastX", nan);
101  this->SetObservableValue("LastY", nan);
102 
103  set<int> TriggeredModuleId;
104 
105  double XEnergySum = 0;
106  double XEnergyPosSum = 0;
107  double YEnergySum = 0;
108  double YEnergyPosSum = 0;
109 
110  // calculate firstx, firsty in position coordinate
111  for (int i = 0; i < fSignalEvent->GetNumberOfSignals(); i++) {
112  TRestRawSignal* sgnl = fSignalEvent->GetSignal(i);
113 
114  int p, m, c;
115  fReadout->GetPlaneModuleChannel(sgnl->GetID(), p, m, c);
116  TriggeredModuleId.insert(m);
117 
118  if (!TMath::IsNaN(fReadout->GetX(sgnl->GetID()))) {
119  XEnergySum += sgnl->GetIntegral();
120  XEnergyPosSum += fReadout->GetX(sgnl->GetID()) * sgnl->GetIntegral();
121  }
122  if (!TMath::IsNaN(fReadout->GetY(sgnl->GetID()))) {
123  YEnergySum += sgnl->GetIntegral();
124  YEnergyPosSum += fReadout->GetY(sgnl->GetID()) * sgnl->GetIntegral();
125  }
126 
127  if (sgnl->GetMaxPeakBin() < firstX_t) {
128  if (!TMath::IsNaN(fReadout->GetX(sgnl->GetID()))) {
129  firstX_id = sgnl->GetID();
130  firstX_t = sgnl->GetMaxPeakBin();
131  }
132  }
133  if (sgnl->GetMaxPeakBin() < firstY_t) {
134  if (!TMath::IsNaN(fReadout->GetY(sgnl->GetID()))) {
135  firstY_id = sgnl->GetID();
136  firstY_t = sgnl->GetMaxPeakBin();
137  }
138  }
139 
140  if (sgnl->GetMaxPeakBin() > lastX_t) {
141  if (!TMath::IsNaN(fReadout->GetX(sgnl->GetID()))) {
142  lastX_id = sgnl->GetID();
143  lastX_t = sgnl->GetMaxPeakBin();
144  }
145  }
146  if (sgnl->GetMaxPeakBin() > lastY_t) {
147  if (!TMath::IsNaN(fReadout->GetY(sgnl->GetID()))) {
148  lastY_id = sgnl->GetID();
149  lastY_t = sgnl->GetMaxPeakBin();
150  }
151  }
152  }
153 
154  this->SetObservableValue("MeanX", XEnergySum == 0 ? nan : XEnergyPosSum / XEnergySum);
155  this->SetObservableValue("MeanY", YEnergySum == 0 ? nan : YEnergyPosSum / YEnergySum);
156  this->SetObservableValue("NmodulesTriggered", (int)TriggeredModuleId.size());
157  this->SetObservableValue("TriggeredModuleId", TriggeredModuleId);
158 
159  // fill firstx/y hitmap in channel coordinate
160  map<int, int> modulefirstxchannel; // moduleid, firstx channelid
161  map<int, int> modulefirstychannel; // moduleid, firsty channelid
162  if (firstX_id > -1 && firstY_id > -1) {
163  double firstx = fReadout->GetX(firstX_id);
164  double firsty = fReadout->GetY(firstY_id);
165  double lastx = fReadout->GetX(lastX_id);
166  double lasty = fReadout->GetY(lastY_id);
167  this->SetObservableValue("FirstX", firstx);
168  this->SetObservableValue("FirstY", firsty);
169  this->SetObservableValue("LastX", lastx);
170  this->SetObservableValue("LastY", lasty);
171 
172  int mod1 = -1, mod2 = -1;
173  int channel1 = -1, channel2 = -1;
174  int plane = -1;
175  fReadout->GetPlaneModuleChannel(firstX_id, plane, mod1, channel1);
176  fReadout->GetPlaneModuleChannel(firstY_id, plane, mod2, channel2);
177  if (mod1 == mod2 && mod1 > -1) {
178  // consider the rotation of readout module, firstX may be from the Y channel!
179  int x = -1, y = -1;
180  int n = fReadout->GetReadoutModuleWithID(mod1)->GetNumberOfChannels() / 2;
181  if (channel1 >= n && channel2 < n) {
182  x = channel2;
183  y = channel1;
184  } else if (channel2 >= n && channel1 < n) {
185  x = channel1;
186  y = channel2;
187  }
188  modulefirstxchannel[mod1] = x;
189  modulefirstychannel[mod1] = y;
190  if (fModuleHitMaps.count(mod1) > 0) {
191  if (fModuleHitMaps[mod1] != nullptr) fModuleHitMaps[mod1]->Fill(x, y);
192  }
193  // cout << n<<" "<<channel1 <<" "<< channel2 << endl;
194  // cout << x << " " << y << endl;
195  // cout << fReadout->GetX(firstX_id) << " " << fReadout->GetY(firstY_id)
196  // << endl; cout << endl;
197 
198  RESTDebug << "TRestRawReadoutAnalysisProcess. Adding point to hitmap of "
199  "module : "
200  << mod1 << RESTendl;
201  RESTDebug << "Position on module(X, Y) : (" << x << ", " << y << ")" << RESTendl;
202  RESTDebug << "Absolute position:(X, Y) : (" << firstx << ", " << firsty << ")" << RESTendl;
203  }
204  }
205  this->SetObservableValue("ModuleFirstX", modulefirstxchannel);
206  this->SetObservableValue("ModuleFirstY", modulefirstychannel);
207 
208  // for each channel
209  vector<int> moduleid;
210  vector<int> channelid;
211  vector<double> baselinesigma;
212  vector<double> baseline;
213  vector<double> thresholdint;
214  // map<int, map<int, double>> modulebaselinesigma; // moduleid, channelid, baselinesigma
215  // map<int, map<int, double>> modulebaseline; // moduleid, channelid, baseline
216  // map<int, map<int, double>> modulethresholdint; // moduleid, channelid, thresholdintergal
217 
218  for (int i = 0; i < fSignalEvent->GetNumberOfSignals(); i++) {
219  TRestRawSignal* sgn = fSignalEvent->GetSignal(i);
220 
221  // channel histo
222  int plane = -1, mod = -1, channel = -1;
223  fReadout->GetPlaneModuleChannel(sgn->GetID(), plane, mod, channel);
224  if (mod != -1 && channel != -1) {
225  moduleid.push_back(mod);
226  channelid.push_back(channel);
227  baselinesigma.push_back(sgn->GetBaseLineSigma());
228  baseline.push_back(sgn->GetBaseLine());
229  thresholdint.push_back(sgn->GetThresholdIntegral());
230 
231  if (fModuleHitMaps.count(mod) > 0) {
232  fModuleActivityX[mod]->Fill(channel);
233  fModuleActivityY[mod]->Fill(channel);
234  fModuleBSLSigmaX[mod]->Fill(channel, sgn->GetBaseLineSigma());
235  fModuleBSLSigmaY[mod]->Fill(sgn->GetBaseLineSigma(), channel);
236  }
237  }
238  }
239  this->SetObservableValue("Module", moduleid);
240  this->SetObservableValue("Channel", channelid);
241  this->SetObservableValue("BaselineSigma", baselinesigma);
242  this->SetObservableValue("Baseline", baseline);
243  this->SetObservableValue("ThresholdIntegral", thresholdint);
244  }
245  return fSignalEvent;
246 }
247 
248 void TRestRawReadoutAnalysisProcess::EndProcess() {
249  if (fReadout != nullptr) {
250  {
251  auto iter = fModuleHitMaps.begin();
252  while (iter != fModuleHitMaps.end()) {
253  if (fModuleHitMaps[iter->first] != nullptr) {
254  TH2D* histo = fModuleHitMaps[iter->first];
255  histo->GetXaxis()->SetTitle("X Channel");
256  histo->GetYaxis()->SetTitle("Y Channel");
257  histo->Write();
258  }
259 
260  if (fModuleCanvasSave != "none") {
261  if (fModuleHitMaps[iter->first] != nullptr) {
262  TH2D* histo = fModuleHitMaps[iter->first];
263  TCanvas* c1 =
264  new TCanvas((TString) "Can_ModuleHitMap" + ToString(iter->first),
265  (TString) "Hit Map of Module " + ToString(iter->first), 800, 600);
266  histo->Draw("colz");
267  c1->Write();
268  c1->Print((TString)fModuleCanvasSave + "/M" + ToString(iter->first) + "_HitMap.png");
269  delete c1;
270  }
271 
272  auto h0 = fModuleHitMaps[iter->first];
273  h0->SetStats(false);
274  h0->Reset();
275 
276  if (fModuleActivityX[iter->first] != nullptr &&
277  fModuleActivityY[iter->first] != nullptr) {
278  TH1D* h1 = fModuleActivityX[iter->first];
279  TH1D* h2 = fModuleActivityY[iter->first];
280 
281  TCanvas* c1 = new TCanvas(
282  (TString) "Can_ModuleActivity" + ToString(iter->first),
283  (TString) "Channel Activity of Module " + ToString(iter->first), 800, 600);
284  c1->Divide(2, 2, 0, 0);
285  c1->cd(1);
286  h1->SetFillColor(kBlue);
287  h1->SetStats(false);
288  h1->GetYaxis()->SetTitle("Counts");
289  h1->Draw("bar");
290 
291  c1->cd(4);
292  h2->SetFillColor(kBlue);
293  h2->SetStats(false);
294  h2->GetXaxis()->SetTitle("Counts");
295  h2->Draw("hbar");
296 
297  c1->cd(3);
298  h0->Draw("colz");
299  c1->Write();
300  c1->Print((TString)fModuleCanvasSave + "/M" + ToString(iter->first) +
301  "_Activity.png");
302  delete c1;
303  }
304 
305  if (fModuleBSLSigmaX[iter->first] != nullptr &&
306  fModuleBSLSigmaY[iter->first] != nullptr) {
307  TH2D* h1 = fModuleBSLSigmaX[iter->first];
308  TH2D* h2 = fModuleBSLSigmaY[iter->first];
309 
310  TCanvas* c1 = new TCanvas(
311  (TString) "Can_ModuleBSLS" + ToString(iter->first),
312  (TString) "Channel Baseline Sigma of Module " + ToString(iter->first), 800, 600);
313  c1->Divide(2, 2, 0, 0);
314  c1->cd(1);
315  h1->SetFillColor(kBlue);
316  h1->SetStats(false);
317  h1->GetYaxis()->SetTitle("Baseline Sigma");
318  h1->Draw("colz");
319 
320  c1->cd(4);
321  h2->SetFillColor(kBlue);
322  h2->SetStats(false);
323  h2->GetXaxis()->SetTitle("Baseline Sigma");
324  h2->Draw("colz");
325 
326  c1->cd(3);
327  h0->Draw("colz");
328  c1->Write();
329  c1->Print((TString)fModuleCanvasSave + "/M" + ToString(iter->first) + "_BSLS.png");
330  delete c1;
331  }
332  }
333  iter++;
334  }
335  }
336  }
337 }
338 
339 // setting amplification:
340 // <parameter name="modulesAmp" value = "2-1:5-1.2:6-0.8:8-0.9" />
341 // setting readout modules to draw:
342 // <parameter name="modulesHist" value="2:5:6:8"/>
343 void TRestRawReadoutAnalysisProcess::InitFromConfigFile() {
344  fModuleCanvasSave = GetParameter("outputPlotPath", "none");
345  if (fModuleCanvasSave != "none") {
346  fSingleThreadOnly = true;
347  TRestTools::Execute("mkdir -p " + fModuleCanvasSave);
348  }
349 
350  string moduleHist = GetParameter("modulesHist", "");
351  auto histdef = Split(moduleHist, ":");
352  for (unsigned int i = 0; i < histdef.size(); i++) {
353  fModuleHitMaps[StringToInteger(histdef[i])] = nullptr;
354  fModuleActivityX[StringToInteger(histdef[i])] = nullptr;
355  fModuleActivityY[StringToInteger(histdef[i])] = nullptr;
356  fModuleBSLSigmaX[StringToInteger(histdef[i])] = nullptr;
357  fModuleBSLSigmaY[StringToInteger(histdef[i])] = nullptr;
358  }
359 }
TRestDetectorReadoutModule::GetNumberOfChannels
size_t GetNumberOfChannels() const
Returns the total number of channels defined inside the module.
Definition: TRestDetectorReadoutModule.h:178
TRestEvent
A base class for any REST event.
Definition: TRestEvent.h:38
TRestRawReadoutAnalysisProcess::EndProcess
void EndProcess() override
To be executed at the end of the run (outside event loop)
Definition: TRestRawReadoutAnalysisProcess.cxx:248
TRestRawReadoutAnalysisProcess::ProcessEvent
TRestEvent * ProcessEvent(TRestEvent *inputEvent) override
Process one event.
Definition: TRestRawReadoutAnalysisProcess.cxx:84
TRestRawReadoutAnalysisProcess::Initialize
void Initialize() override
Making default settings.
Definition: TRestRawReadoutAnalysisProcess.cxx:30
TRestRawReadoutAnalysisProcess::InitFromConfigFile
void InitFromConfigFile() override
To make settings from rml file. This method must be implemented in the derived class.
Definition: TRestRawReadoutAnalysisProcess.cxx:343
TRestRawReadoutAnalysisProcess::InitProcess
void InitProcess() override
To be executed at the beginning of the run (outside event loop)
Definition: TRestRawReadoutAnalysisProcess.cxx:39
TRestRawSignalEvent
An event container for time rawdata signals with fixed length.
Definition: TRestRawSignalEvent.h:40
TRestRawSignal
It defines a Short_t array with a physical parameter that evolves in time using a fixed time bin.
Definition: TRestRawSignal.h:36
TRestRawSignal::GetBaseLineSigma
Double_t GetBaseLineSigma() const
Definition: TRestRawSignal.h:113
TRestRawSignal::GetThresholdIntegral
Double_t GetThresholdIntegral()
It returns the integral of points identified over threshold. InitializePointsOverThreshold should hav...
Definition: src/TRestRawSignal.cxx:329
TRestRawSignal::GetID
Int_t GetID() const
Returns the value of signal ID.
Definition: TRestRawSignal.h:87
TRestRawSignal::GetBaseLine
Double_t GetBaseLine() const
Definition: TRestRawSignal.h:109
TRestRawSignal::GetIntegral
Double_t GetIntegral()
It returns the integral of points found in the region defined by fRange. If fRange was not defined th...
Definition: src/TRestRawSignal.cxx:302
TRestRawSignal::GetMaxPeakBin
Int_t GetMaxPeakBin()
It returns the bin at which the maximum peak amplitude happens.
Definition: src/TRestRawSignal.cxx:489
TRestTools::Execute
static std::string Execute(std::string cmd)
Executes a shell command and returns its output in a string.
Definition: TRestTools.cxx:1061
REST_StringHelper::Split
std::vector< std::string > Split(std::string in, std::string separator, bool allowBlankString=false, bool removeWhiteSpaces=false, int startPos=-1)
Split the input string according to the given separator. Returning a vector of fragments.
Definition: TRestStringHelper.cxx:236
REST_StringHelper::StringToInteger
Int_t StringToInteger(std::string in)
Gets an integer from a string.
Definition: TRestStringHelper.cxx:624
Generated on Mon Nov 25 2024 10:53:04 for REST-for-Physics by doxygen 1.9.1