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