R3BROOT
R3B analysis software
Loading...
Searching...
No Matches
R3BNeulandHitCalibrationEngine.cxx
Go to the documentation of this file.
1/******************************************************************************
2 * Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH *
3 * Copyright (C) 2019-2025 Members of R3B Collaboration *
4 * *
5 * This software is distributed under the terms of the *
6 * GNU General Public Licence (GPL) version 3, *
7 * copied verbatim in the file "LICENSE". *
8 * *
9 * In applying this license GSI does not waive the privileges and immunities *
10 * granted to it by virtue of its status as an Intergovernmental Organization *
11 * or submit itself to any jurisdiction. *
12 ******************************************************************************/
13
14#include "R3BNeulandHitCalibrationEngine.h"
15#include "R3BNeulandHitModulePar.h"
16#include "R3BNeulandHitPar.h"
17
18#include "FairLogger.h"
19
20#include "TCanvas.h"
21#include "TDirectory.h"
22#include <cmath>
23
24using DPair = std::array<Double_t, 2>;
25using CalibrationStatus = R3B::Neuland::Calibration::HitCalibrationBar::CalibrationStatus;
26
27namespace R3B::Neuland
28{
29
30 // NaN2Value maps NaNs and Infs to a finite value.
31 // Usefull if you want to fill a possible nan variable into a histogram
32 inline double NaN2Value(const double val, const double valIfNaN = 0.)
33 {
34 return (std::isfinite(val) ? val : valIfNaN);
35 }
36
37 namespace Calibration
38 {
39 HitCalibrationEngine::HitCalibrationEngine() {}
40
41 void HitCalibrationEngine::Init(const R3BNeulandHitPar* hitpar)
42 {
43 const auto nPlanes = (!hitpar ? MaxNumberOfPlanes : hitpar->GetNumberOfPlanes());
44 const auto totalBars = nPlanes * BarsPerPlane;
45 fHitMask.resize(nPlanes, 0L);
46
47 fBars.reserve(totalBars);
48 for (auto plane = 0; plane < nPlanes; ++plane)
49 for (auto bar = 0; bar < BarsPerPlane; ++bar)
50 fBars.emplace_back(plane * BarsPerPlane + bar + 1);
51
52 const auto nParameter = hitpar->GetNumModulePar();
53
54 for (auto i = 0; i < nParameter; i++)
55 {
56 const auto modulePar = hitpar->GetModuleParAt(i);
57 const auto barID = modulePar->GetModuleId() - 1;
58 fBars[barID].Update(modulePar);
59 }
60
61 fBarDistribution = TH1F(
62 "HitCalibrationEngine::BarDistribution", "Bar Distribution; Bar ID", totalBars, 0.5, totalBars + 0.5);
63 fBarDistribution.SetStats(kFALSE);
64 fBarDistribution.SetDirectory(nullptr);
65
66 fStoppedDistribution = TH1F("HitCalibrationEngine::StoppedDistribution",
67 "Stopped Bar Distribution; Bar ID",
68 totalBars,
69 0.5,
70 totalBars + 0.5);
71 fStoppedDistribution.SetStats(kFALSE);
72 fStoppedDistribution.SetDirectory(nullptr);
73
74 fInteractionsDistribution = TH1F("HitCalibrationEngine::InteractionsDistribution",
75 "Interactions Distribution; Number of Interactions",
76 250,
77 0.5,
78 250 + 0.5);
79 fInteractionsDistribution.SetStats(kFALSE);
80 fInteractionsDistribution.SetDirectory(nullptr);
81
82 fStoppedInteractionsDistribution = TH1F("HitCalibrationEngine::StoppedInteractionsDistribution",
83 "Stopped Interactions Distribution; Number of Interactions",
84 250,
85 0.5,
86 250 + 0.5);
87 fStoppedInteractionsDistribution.SetStats(kFALSE);
88 fStoppedInteractionsDistribution.SetDirectory(nullptr);
89 fStoppedInteractionsDistribution.SetLineColor(kRed);
90
91 fTrackLengthDistribution = TH1F("HitCalibrationEngine::TrackLengthDistribution",
92 "Track Length Distribution; Track Length / cm",
93 5000,
94 0.,
95 100.);
96 fTrackLengthDistribution.SetStats(kFALSE);
97 fTrackLengthDistribution.SetDirectory(nullptr);
98
99 fTotalTrackLengthDistribution = TH1F("HitCalibrationEngine::TotalTrackLengthDistribution",
100 "Total Track Length Distribution; Track Length / cm",
101 1000,
102 0.,
103 500.);
104 fTotalTrackLengthDistribution.SetStats(kFALSE);
105 fTotalTrackLengthDistribution.SetDirectory(nullptr);
106
107 fTotalStoppedTrackLengthDistribution = TH1F("HitCalibrationEngine::fTotalStoppedTrackLengthDistribution",
108 "Total Stopped Track Length Distribution; Track Length / cm",
109 1000,
110 0.,
111 500.);
112 fTotalStoppedTrackLengthDistribution.SetStats(kFALSE);
113 fTotalStoppedTrackLengthDistribution.SetDirectory(nullptr);
114 fTotalStoppedTrackLengthDistribution.SetLineColor(kRed);
115
116 fCorrelationMatrix = TH2F("HitCalibrationEngine::CorrelationMatrix",
117 "Correlation Matrix;Bar ID;Bar ID",
118 totalBars,
119 0.5,
120 0.5 + totalBars,
121 totalBars,
122 0.5,
123 0.5 + totalBars);
124 fCorrelationMatrix.SetStats(kFALSE);
125 fCorrelationMatrix.SetDirectory(nullptr);
126
127 fTrackEntryPointDistribution = TH3F("HitCalibrationEngine::TrackEntryPointDistribution",
128 "Entry Point Distribution; Z / cm; X / cm; Y / cm",
129 150,
130 0,
131 300,
132 130,
133 -130,
134 130,
135 130,
136 -130,
137 130);
138 fTrackEntryPointDistribution.SetStats(kFALSE);
139 fTrackEntryPointDistribution.SetDirectory(nullptr);
140
141 fTrackDirectionDistribution = TH3F("HitCalibrationEngine::TrackDirectionDistribution",
142 "Track Direction Distribution; Z / cm/ns; X / cm/ns; Y / cm/ns",
143 100,
144 -50,
145 50,
146 100,
147 -50,
148 50,
149 100,
150 -50,
151 50);
152 fTrackDirectionDistribution.SetStats(kFALSE);
153 fTrackDirectionDistribution.SetDirectory(nullptr);
154 }
155
156 void HitCalibrationEngine::Set(const Int_t id, const Int_t side, const Double_t time, const Int_t qdc)
157 {
158 const auto plane = id / BarsPerPlane;
159 const auto bar = id % BarsPerPlane;
160
161 fBars[id].Set(side, time, qdc);
162
163 fHitMask[plane] |= 1UL << bar;
164 }
165
166 void HitCalibrationEngine::Add(const R3BNeulandCosmicTrack& track, const UInt_t eventNumber)
167 {
168 fTotalTrackLengthDistribution.Fill(track.TotalTrackLength);
169 if (track.Stopped)
170 fTotalStoppedTrackLengthDistribution.Fill(track.TotalTrackLength);
171
172 fTrackEntryPointDistribution.Fill(track.EntryPoint.Z(), track.EntryPoint.X(), track.EntryPoint.Y());
173 fTrackDirectionDistribution.Fill(track.Direction.Z(), track.Direction.X(), track.Direction.Y());
174
175 const auto nPlanes = fHitMask.size();
176
177 fInteractionsDistribution.Fill(track.Interactions.size());
178 if (track.Stopped)
179 {
180 fStoppedDistribution.Fill(track.Interactions.back().BarID + 1);
181 fStoppedInteractionsDistribution.Fill(track.Interactions.size());
182 }
183
184 // first add correlations
185 for (auto plane = 0; plane < nPlanes; ++plane)
186 {
187 if (!fHitMask[plane])
188 continue;
189
190 for (auto bar = 0; bar < BarsPerPlane; ++bar)
191 {
192 if (!(fHitMask[plane] & (1UL << bar)))
193 continue;
194
195 for (auto otherPlane = 0; otherPlane < nPlanes; ++otherPlane)
196 {
197 if (!fHitMask[otherPlane])
198 continue;
199
200 for (auto otherBar = 0; otherBar < BarsPerPlane; ++otherBar)
201 {
202 if ((fHitMask[otherPlane] & (1UL << otherBar)) && (plane != otherPlane || bar != otherBar))
203 fCorrelationMatrix.Fill(BarsPerPlane * plane + bar + 1,
204 BarsPerPlane * otherPlane + otherBar + 1);
205 }
206 }
207 }
208 }
209
210 // if the track stopped in NeuLAND, do not take the last interaction
211 const auto nInteractions = track.Interactions.size() - (track.Stopped ? 1 : 0);
212 for (auto i = 0; i < nInteractions; ++i)
213 {
214 const auto& interaction = track.Interactions[i];
215 const auto barID = interaction.BarID;
216 const auto bar = barID % BarsPerPlane;
217 const auto plane = ModuleID2PlaneID(barID);
218
219 fBarDistribution.Fill(barID + 1);
220 fTrackLengthDistribution.Fill(interaction.TrackLength);
221
222 if (fBars[barID].Add(interaction.EntryTime,
223 interaction.EntryPosition,
224 interaction.ExitPosition,
225 interaction.Energy,
226 eventNumber))
227 {
228 fTSyncer.AddBarData(barID, fBars[barID].GetTime());
229 }
230
231 fBars[barID].Reset();
232 fHitMask[plane] &= ~(1UL << bar);
233 }
234
235 fTSyncer.DoEvent();
236
237 for (auto plane = 0; plane < nPlanes; ++plane)
238 {
239 if (fHitMask[plane] == 0UL)
240 continue;
241
242 for (auto bar = 0; bar < BarsPerPlane; ++bar)
243 {
244 if (fHitMask[plane] & (1UL << bar))
245 fBars[plane * BarsPerPlane + bar].Reset();
246 }
247 fHitMask[plane] = 0UL;
248 }
249 }
250
251 void HitCalibrationEngine::Reset()
252 {
253 for (auto p = 0; p < fHitMask.size(); ++p)
254 {
255 if (!fHitMask[p])
256 continue;
257
258 for (auto b = 0; b < BarsPerPlane; ++b)
259 {
260 if (fHitMask[p] & (1UL << b))
261 fBars[p * BarsPerPlane + b].Reset();
262 }
263 fHitMask[p] = 0UL;
264 }
265 }
266
267 std::vector<R3BNeulandHitModulePar> HitCalibrationEngine::Calibrate(TDirectory* histoDir)
268 {
269 const auto nPlanes = fHitMask.size();
270 const auto nBars = nPlanes * BarsPerPlane;
271 std::vector<R3BNeulandHitModulePar> allParameters;
272 if (histoDir)
273 {
274 histoDir->cd();
275 TCanvas canvasTracks("Tracks", "Tracks");
276 canvasTracks.SetWindowSize(1920, 1200);
277 canvasTracks.Divide(1, 2);
278 auto pad = canvasTracks.cd(1);
279 pad->Divide(5, 1);
280
281 pad->cd(1);
282 fBarDistribution.Draw();
283 pad->cd(2);
284 fStoppedDistribution.Draw();
285 pad->cd(3);
286 fInteractionsDistribution.Draw();
287 fStoppedInteractionsDistribution.Draw("SAME");
288 pad->cd(4);
289 fTotalTrackLengthDistribution.Draw();
290 fTotalStoppedTrackLengthDistribution.Draw("SAME");
291 pad->cd(5);
292 fTrackLengthDistribution.Draw();
293
294 pad = canvasTracks.cd(2);
295 pad->Divide(3, 1);
296 pad->cd(1);
297 fCorrelationMatrix.Draw("COLZ");
298 pad->cd(2);
299 fTrackEntryPointDistribution.Draw();
300 pad->cd(3);
301 fTrackDirectionDistribution.Draw();
302 canvasTracks.Write("Tracking");
303 }
304
305 TDirectory* planeDir = nullptr;
306 for (Int_t plane = 0; plane < nPlanes; ++plane)
307 {
308 if (histoDir)
309 planeDir = histoDir->mkdir(TString::Format("Plane_%d", plane + 1));
310
311 for (Int_t bar = 0; bar < BarsPerPlane; ++bar)
312 {
313 std::cout << "Calibrating Bar " << BarsPerPlane * plane + bar << "/" << BarsPerPlane * nPlanes
314 << "\r" << std::flush;
315 fBars[BarsPerPlane * plane + bar].Calibrate(planeDir);
316 }
317 }
318 if (histoDir)
319 histoDir->cd();
320
321 std::cout << "Syncing NeuLAND Bars... \r" << std::flush;
322 const auto tsync = fTSyncer.GetTSync(nPlanes);
323
324 for (Int_t plane = 0; plane < nPlanes; ++plane)
325 {
326 for (Int_t bar = 0; bar < BarsPerPlane; ++bar)
327 {
328 const auto barID = plane * BarsPerPlane + bar;
329 fBars[barID].SetGlobalTSync(tsync[barID].Value, tsync[barID].Error);
330
331 const auto calibStatus = fBars[barID].GetCalibrationStatus();
332
333 if (calibStatus != CalibrationStatus::completeFail && calibStatus != CalibrationStatus::noData)
334 allParameters.push_back(fBars[barID].GetParameters());
335 }
336 }
337
338 TH1F hTDiff("R3BNeulandCal2HitPar::TDiff", "TDiff;BarID;TDiff / ns", nBars, 0.5, 0.5 + nBars);
339 hTDiff.SetStats(false);
340 TH1F hTSync("R3BNeulandCal2HitPar::TSync", "TSync;BarID;TSync / ns", nBars, 0.5, 0.5 + nBars);
341 hTSync.SetStats(false);
342 TH1F hVEff("R3BNeulandCal2HitPar::VEff", "VEff;BarID;VEff / cm/ns", nBars, 0.5, 0.5 + nBars);
343 hVEff.SetStats(false);
344 TH1F hAtt("R3BNeulandCal2HitPar::Att",
345 "Attenuation Length;BarID;Attenuation Length / cm",
346 nBars,
347 0.5,
348 0.5 + nBars);
349 hAtt.SetStats(false);
350 TH1F hPosErr(
351 "R3BNeulandCal2HitPar::PosError", "Position Error;BarID;\\Delta Pos / cm", nBars, 0.5, 0.5 + nBars);
352 hPosErr.SetStats(false);
353 TH1F hEPosErr("R3BNeulandCal2HitPar::EnergyPosError",
354 "Energy Position Error;BarID;\\Delta Pos / cm",
355 nBars,
356 0.5,
357 0.5 + nBars);
358 hEPosErr.SetStats(false);
359 std::array<TH1F, 2> hGain = {
360 TH1F("R3BNeulandCal2HitPar::Gain1", "Gain;BarID;Gain / ns/MeV", nBars, 0.5, 0.5 + nBars),
361
362 TH1F("R3BNeulandCal2HitPar::Gain2", "Gain;BarID;Gain / ns/MeV", nBars, 0.5, 0.5 + nBars)
363 };
364 hGain[0].SetLineColor(kRed);
365 hGain[1].SetLineColor(kGreen);
366 hGain[0].SetStats(false);
367 hGain[1].SetStats(false);
368
369 std::array<TH1F, 2> hSat = {
370 TH1F("R3BNeulandCal2HitPar::Sat1", "Saturation;BarID;Saturation / 1/MeV", nBars, 0.5, 0.5 + nBars),
371 TH1F("R3BNeulandCal2HitPar::Sat2", "Saturation;BarID;Saturation / 1/MeV", nBars, 0.5, 0.5 + nBars)
372 };
373 hSat[0].SetLineColor(kRed);
374 hSat[1].SetLineColor(kGreen);
375 hSat[0].SetStats(false);
376 hSat[1].SetStats(false);
377
378 std::array<TH1F, 2> hThr = {
379 TH1F("R3BNeulandCal2HitPar::Thr1", "Threshold;BarID;Threshold / MeV", nBars, 0.5, 0.5 + nBars),
380 TH1F("R3BNeulandCal2HitPar::Thr2", "Threshold;BarID;Threshold / MeV", nBars, 0.5, 0.5 + nBars)
381 };
382 hThr[0].SetLineColor(kRed);
383 hThr[1].SetLineColor(kGreen);
384 hThr[0].SetStats(false);
385 hThr[1].SetStats(false);
386
387 for (auto& parameter : allParameters)
388 {
389 const auto barID = parameter.GetModuleId();
390 const auto attLength = NaN2Value(parameter.GetLightAttenuationLength());
391
392 hTDiff.Fill(barID, parameter.GetTDiff());
393 hTSync.Fill(barID, NaN2Value(parameter.GetTSync()));
394 hVEff.Fill(barID, NaN2Value(parameter.GetEffectiveSpeed()));
395 hAtt.Fill(barID, attLength);
396
397 for (Int_t side = 0; side < 2; ++side)
398 {
399 hGain[side].Fill(barID, NaN2Value(parameter.GetEnergyGain(side + 1)));
400 hSat[side].Fill(barID, NaN2Value(parameter.GetPMTSaturation(side + 1)));
401 hThr[side].Fill(barID, NaN2Value(parameter.GetPMTThreshold(side + 1)));
402 }
403 }
404
405 draw();
406
407 if (histoDir)
408 {
409 TCanvas canvasOverview("Overview", "Neuland Hit Calibration");
410 canvasOverview.Divide(3, 2);
411
412 canvasOverview.cd(1);
413 hTDiff.Draw("HIST");
414 canvasOverview.cd(2);
415 hVEff.Draw("HIST");
416 canvasOverview.cd(3);
417 hTSync.Draw("HIST");
418 canvasOverview.cd(4);
419 hGain[0].Draw("HIST");
420 hGain[1].Draw("HIST SAME");
421 canvasOverview.cd(5);
422 hAtt.Draw("HIST");
423 canvasOverview.cd(6);
424 hThr[0].Draw("HIST");
425 hThr[1].Draw("HIST SAME");
426 canvasOverview.Write("Overview");
427 }
428 return allParameters;
429 }
430
431 void HitCalibrationEngine::draw() const
432 {
433 const auto nPlanes = fHitMask.size();
434 int stats = 0;
435
436 std::vector<CalibrationStatus> statuses;
437 statuses.reserve(16);
438
439 std::cout << std::endl << std::endl << " |";
440 for (auto plane = 0; plane < nPlanes; ++plane)
441 {
442 std::cout << TString::Format("%2d|", plane + 1);
443 }
444 std::cout << std::endl << " ";
445 for (auto plane = 0; plane < nPlanes; ++plane)
446 {
447 std::cout << " ꓕꓕ";
448 }
449 std::cout << std::endl;
450 for (auto bar = BarsPerPlane - 1; bar >= 0; --bar)
451 {
452 std::cout << TString::Format("%2d |", bar + 1);
453 for (auto plane = 0; plane < nPlanes; ++plane)
454 {
455 const auto cstatus = fBars[plane * BarsPerPlane + bar].GetCalibrationStatus();
456
457 if (cstatus != CalibrationStatus::success &&
458 std::find(statuses.begin(), statuses.end(), cstatus) == statuses.end())
459 {
460 statuses.push_back(cstatus);
461 }
462 std::cout << HitCalibrationBar::GetCalibrationStatusAbbreviation(cstatus);
463 std::cout << "|";
464 }
465 std::cout << std::endl;
466 }
467 std::cout << " ";
468 for (auto plane = 0; plane < nPlanes; ++plane)
469 {
470 std::cout << " TT";
471 }
472 std::cout << std::endl << " |";
473 for (auto plane = 0; plane < nPlanes; ++plane)
474 {
475 std::cout << TString::Format("%2d|", plane + 1);
476 }
477 std::cout << std::endl << std::endl;
478 if (statuses.size() > 0)
479 {
480 for (const auto status : statuses)
481 {
482 std::cout << " '" << HitCalibrationBar::GetCalibrationStatusAbbreviation(status)
483 << "' = " << HitCalibrationBar::GetCalibrationStatusDescription(status) << std::endl;
484 }
485 std::cout << std::endl;
486 }
487 }
488
489 } // namespace Calibration
490} // namespace R3B::Neuland
DPair
std::array< Double_t, 2 > DPair
Definition R3BNeulandHitCalibrationBar.cxx:31
CalibrationStatus
R3B::Neuland::Calibration::HitCalibrationBar::CalibrationStatus CalibrationStatus
Definition R3BNeulandHitCalibrationEngine.cxx:25
R3B::Neuland::Calibration::HitCalibrationBar::GetCalibrationStatusDescription
static const char *const GetCalibrationStatusDescription(const CalibrationStatus status)
Definition R3BNeulandHitCalibrationBar.h:76
R3B::Neuland::Calibration::HitCalibrationBar::GetCalibrationStatusAbbreviation
static const char *const GetCalibrationStatusAbbreviation(const CalibrationStatus status)
Definition R3BNeulandHitCalibrationBar.h:80
R3B::Neuland::Calibration::HitCalibrationEngine::Calibrate
std::vector< R3BNeulandHitModulePar > Calibrate(TDirectory *histoDir=nullptr)
Definition R3BNeulandHitCalibrationEngine.cxx:267
R3B::Neuland::Calibration::HitCalibrationEngine::Set
void Set(const Int_t id, const Int_t side, const Double_t time, const Int_t qdc)
Definition R3BNeulandHitCalibrationEngine.cxx:156
R3B::Neuland::Calibration::HitCalibrationEngine::Add
void Add(const R3BNeulandCosmicTrack &track, const UInt_t eventNumber)
Definition R3BNeulandHitCalibrationEngine.cxx:166
R3BNeulandHitPar::GetModuleParAt
auto GetModuleParAt(Int_t idx) const -> R3BNeulandHitModulePar *
Method to get single parameter container for a specific module.
Definition R3BNeulandHitPar.h:93
R3BNeulandHitPar::GetNumModulePar
Int_t GetNumModulePar() const
Method to get number of modules storred in array.
Definition R3BNeulandHitPar.h:86
R3BNeulandHitPar::GetNumberOfPlanes
Int_t GetNumberOfPlanes() const
Definition R3BNeulandHitPar.h:114
R3B::Neuland
Simulation of NeuLAND Bar/Paddle.
Definition R3BNeulandAnalysisApp.cxx:80
R3B::Neuland::BarsPerPlane
constexpr auto BarsPerPlane
Definition R3BNeulandCommon.h:84
R3B::Neuland::ModuleID2PlaneID
constexpr auto ModuleID2PlaneID(int moduleID) -> int
Definition R3BNeulandCommon.h:99
R3B::Neuland::NaN2Value
double NaN2Value(const double val, const double valIfNaN=0.)
Definition R3BNeulandHitCalibrationEngine.cxx:32
R3B::Neuland::MaxNumberOfPlanes
constexpr auto MaxNumberOfPlanes
Definition R3BNeulandCommon.h:85
R3BNeulandCosmicTrack::Interactions
std::vector< Interaction > Interactions
Definition R3BNeulandCosmicTrack.h:50
nPlanes
constexpr Int_t nPlanes
Definition testNeulandTSyncer.C:56