REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
TRestTrack2DAnalysisProcess.cxx
1/*************************************************************************
2 * This file is part of the REST software framework. *
3 * *
4 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) *
5 * For more information see http://gifna.unizar.es/trex *
6 * *
7 * REST is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * REST is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have a copy of the GNU General Public License along with *
18 * REST in $REST_PATH/LICENSE. *
19 * If not, see http://www.gnu.org/licenses/. *
20 * For the list of contributors see $REST_PATH/CREDITS. *
21 *************************************************************************/
22
99
100#include "TRestTrack2DAnalysisProcess.h"
101
102using namespace std;
103
104// Comparator function for sorting in descending order
105bool sortByValueDescending(const std::pair<int, Double_t>& a, const std::pair<int, Double_t>& b) {
106 return a.second > b.second; // Change to < for ascending order
107}
108
109ClassImp(TRestTrack2DAnalysisProcess);
110
111TRestTrack2DAnalysisProcess::TRestTrack2DAnalysisProcess() { Initialize(); }
112
113TRestTrack2DAnalysisProcess::TRestTrack2DAnalysisProcess(const char* configFilename) {
114 Initialize();
115 if (LoadConfigFromFile(configFilename)) {
116 LoadDefaultConfig();
117 }
118}
119
120TRestTrack2DAnalysisProcess::~TRestTrack2DAnalysisProcess() { delete fTrackEvent; }
121
122void TRestTrack2DAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
123
124void TRestTrack2DAnalysisProcess::Initialize() {
125 SetSectionName(this->ClassName());
126 SetLibraryVersion(LIBRARY_VERSION);
127
128 fTrackEvent = new TRestTrackEvent();
129}
130
131void TRestTrack2DAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
132 if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
133}
134
135void TRestTrack2DAnalysisProcess::InitProcess() {}
136
137TRestEvent* TRestTrack2DAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
138 fTrackEvent = (TRestTrackEvent*)inputEvent;
139
143
145 Double_t XZ_TotalEnergyX;
146 Double_t YZ_TotalEnergyY;
147
149 int NTracksX;
150 int NTracksY;
151
154 map<int, Double_t> XZ_EnergyX;
155 map<int, Double_t> YZ_EnergyY;
156
158 map<int, int> XZ_NHitsX;
159 map<int, int> YZ_NHitsY;
160
162 map<int, Double_t> XZ_SigmaX;
163 map<int, Double_t> XZ_SigmaZ;
164 map<int, Double_t> YZ_SigmaY;
165 map<int, Double_t> YZ_SigmaZ;
166
167 map<int, Double_t> XZ_YZ_SigmaXYBalance;
168 map<int, Double_t> XZ_YZ_SigmaZBalance;
169
170 map<int, Double_t> XZ_SkewZ;
171 map<int, Double_t> YZ_SkewZ;
172
174 map<int, Double_t> XZ_GaussSigmaX;
175 map<int, Double_t> XZ_GaussSigmaZ;
176 map<int, Double_t> YZ_GaussSigmaY;
177 map<int, Double_t> YZ_GaussSigmaZ;
178
179 map<int, Double_t> XZ_YZ_GaussSigmaXYBalance;
180 map<int, Double_t> XZ_YZ_GaussSigmaZBalance;
181
183 map<int, Double_t> XZ_LengthX;
184 map<int, Double_t> YZ_LengthY;
185
186 map<int, Double_t> XZ_VolumeX;
187 map<int, Double_t> YZ_VolumeY;
188
189 map<int, Double_t> XZ_MeanX;
190 map<int, Double_t> XZ_MeanZ;
191 map<int, Double_t> YZ_MeanY;
192 map<int, Double_t> YZ_MeanZ;
193
195 Double_t XZ_FirstSecondTracksDistanceXZ;
196 Double_t YZ_FirstSecondTracksDistanceYZ;
197
199 Double_t MaxTrack_XZ_YZ_SigmaZ;
200 Double_t MaxTrack_XZ_YZ_GaussSigmaZ;
201 Double_t MaxTrack_XZ_YZ_SkewXY;
202 Double_t MaxTrack_XZ_YZ_SkewZ;
203
207
208 NTracksX = fTrackEvent->GetNumberOfTracks("X");
209 NTracksY = fTrackEvent->GetNumberOfTracks("Y");
210
212 for (int tck = 0; tck < fTrackEvent->GetNumberOfTracks(); tck++) {
213 if (!fTrackEvent->isTopLevel(tck)) continue;
214
215 TRestTrack* t = fTrackEvent->GetTrack(tck);
216
217 if (t->isXZ()) {
218 XZ_NHitsX[t->GetTrackID()] = t->GetNumberOfHits();
219 XZ_EnergyX[t->GetTrackID()] = t->GetTrackEnergy();
220 XZ_SigmaX[t->GetTrackID()] = t->GetHits()->GetSigmaX();
221 XZ_SigmaZ[t->GetTrackID()] = t->GetHits()->GetSigmaZ2();
222 XZ_GaussSigmaX[t->GetTrackID()] = t->GetHits()->GetGaussSigmaX();
223 XZ_GaussSigmaZ[t->GetTrackID()] = t->GetHits()->GetGaussSigmaZ();
224 XZ_LengthX[t->GetTrackID()] = t->GetLength();
225 XZ_VolumeX[t->GetTrackID()] = t->GetVolume();
226 XZ_MeanX[t->GetTrackID()] = t->GetMeanPosition().X();
227 XZ_MeanZ[t->GetTrackID()] = t->GetMeanPosition().Z();
228 XZ_SkewZ[t->GetTrackID()] = t->GetHits()->GetSkewZ();
229
230 YZ_NHitsY[t->GetTrackID()] = 0;
231 YZ_EnergyY[t->GetTrackID()] = 0;
232 YZ_SigmaY[t->GetTrackID()] = 0;
233 YZ_SigmaZ[t->GetTrackID()] = 0;
234 YZ_GaussSigmaY[t->GetTrackID()] = 0;
235 YZ_GaussSigmaZ[t->GetTrackID()] = 0;
236 YZ_LengthY[t->GetTrackID()] = 0;
237 YZ_VolumeY[t->GetTrackID()] = 0;
238 YZ_MeanY[t->GetTrackID()] = 0;
239 YZ_MeanZ[t->GetTrackID()] = 0;
240 YZ_SkewZ[t->GetTrackID()] = 0;
241 } else if (t->isYZ()) {
242 XZ_NHitsX[t->GetTrackID()] = 0;
243 XZ_EnergyX[t->GetTrackID()] = 0;
244 XZ_SigmaX[t->GetTrackID()] = 0;
245 XZ_SigmaZ[t->GetTrackID()] = 0;
246 XZ_GaussSigmaX[t->GetTrackID()] = 0;
247 XZ_GaussSigmaZ[t->GetTrackID()] = 0;
248 XZ_LengthX[t->GetTrackID()] = 0;
249 XZ_VolumeX[t->GetTrackID()] = 0;
250 XZ_MeanX[t->GetTrackID()] = 0;
251 XZ_MeanZ[t->GetTrackID()] = 0;
252 XZ_SkewZ[t->GetTrackID()] = 0;
253
254 YZ_NHitsY[t->GetTrackID()] = t->GetNumberOfHits();
255 YZ_EnergyY[t->GetTrackID()] = t->GetTrackEnergy();
256 YZ_SigmaY[t->GetTrackID()] = t->GetHits()->GetSigmaY();
257 YZ_SigmaZ[t->GetTrackID()] = t->GetHits()->GetSigmaZ2();
258 YZ_GaussSigmaY[t->GetTrackID()] = t->GetHits()->GetGaussSigmaY();
259 YZ_GaussSigmaZ[t->GetTrackID()] = t->GetHits()->GetGaussSigmaZ();
260 YZ_LengthY[t->GetTrackID()] = t->GetLength();
261 YZ_VolumeY[t->GetTrackID()] = t->GetVolume();
262 YZ_MeanY[t->GetTrackID()] = t->GetMeanPosition().Y();
263 YZ_MeanZ[t->GetTrackID()] = t->GetMeanPosition().Z();
264 YZ_SkewZ[t->GetTrackID()] = t->GetHits()->GetSkewZ();
265 } else {
266 XZ_EnergyX[t->GetTrackID()] = 0;
267 XZ_SigmaX[t->GetTrackID()] = 0;
268 XZ_SigmaZ[t->GetTrackID()] = 0;
269 XZ_GaussSigmaX[t->GetTrackID()] = 0;
270 XZ_GaussSigmaZ[t->GetTrackID()] = 0;
271 XZ_LengthX[t->GetTrackID()] = 0;
272 XZ_VolumeX[t->GetTrackID()] = 0;
273 XZ_MeanX[t->GetTrackID()] = 0;
274 XZ_MeanZ[t->GetTrackID()] = 0;
275 XZ_SkewZ[t->GetTrackID()] = 0;
276
277 YZ_EnergyY[t->GetTrackID()] = 0;
278 YZ_SigmaY[t->GetTrackID()] = 0;
279 YZ_SigmaZ[t->GetTrackID()] = 0;
280 YZ_GaussSigmaY[t->GetTrackID()] = 0;
281 YZ_GaussSigmaZ[t->GetTrackID()] = 0;
282 YZ_LengthY[t->GetTrackID()] = 0;
283 YZ_VolumeY[t->GetTrackID()] = 0;
284 YZ_MeanY[t->GetTrackID()] = 0;
285 YZ_MeanZ[t->GetTrackID()] = 0;
286 YZ_SkewZ[t->GetTrackID()] = 0;
287 }
288 }
289
291 vector<pair<int, Double_t>> energiesX;
292 vector<pair<int, Double_t>> energiesY;
293
295 map<int, Double_t>::iterator it;
296 for (it = XZ_EnergyX.begin(); it != XZ_EnergyX.end(); it++) {
297 energiesX.emplace_back(it->first, it->second);
298 }
299 for (it = YZ_EnergyY.begin(); it != YZ_EnergyY.end(); it++) {
300 energiesY.emplace_back(it->first, it->second);
301 }
302
304 sort(energiesX.begin(), energiesX.end(), sortByValueDescending);
305 sort(energiesY.begin(), energiesY.end(), sortByValueDescending);
306
308 for (int i = 0; i < min(NTracksX, NTracksY); i++) {
309 XZ_YZ_SigmaXYBalance[i] = (XZ_SigmaX[energiesX[i].first] - YZ_SigmaY[energiesY[i].first]) /
310 (XZ_SigmaX[energiesX[i].first] + YZ_SigmaY[energiesY[i].first]);
311 XZ_YZ_SigmaZBalance[i] = (XZ_SigmaZ[energiesX[i].first] - YZ_SigmaZ[energiesY[i].first]) /
312 (XZ_SigmaZ[energiesX[i].first] + YZ_SigmaZ[energiesY[i].first]);
313 XZ_YZ_GaussSigmaXYBalance[i] =
314 (XZ_GaussSigmaX[energiesX[i].first] - YZ_GaussSigmaY[energiesY[i].first]) /
315 (XZ_GaussSigmaX[energiesX[i].first] + YZ_GaussSigmaY[energiesY[i].first]);
316 XZ_YZ_GaussSigmaZBalance[i] =
317 (XZ_GaussSigmaZ[energiesX[i].first] - YZ_GaussSigmaZ[energiesY[i].first]) /
318 (XZ_GaussSigmaZ[energiesX[i].first] + YZ_GaussSigmaZ[energiesY[i].first]);
319 }
320
322 if (fTrackEvent->GetNumberOfTracks() > 1) {
323 Double_t dXz = 0, dxZ = 0, dYz = 0, dyZ = 0;
324
325 dXz = abs(XZ_MeanX[energiesX[0].first] - XZ_MeanX[energiesX[1].first]);
326 dxZ = abs(XZ_MeanZ[energiesX[0].first] - XZ_MeanZ[energiesX[1].first]);
327
328 dYz = abs(YZ_MeanY[energiesY[0].first] - YZ_MeanY[energiesY[1].first]);
329 dyZ = abs(YZ_MeanZ[energiesY[0].first] - YZ_MeanZ[energiesY[1].first]);
330
331 XZ_FirstSecondTracksDistanceXZ = TMath::Sqrt(dXz * dXz + dxZ * dxZ);
332 YZ_FirstSecondTracksDistanceYZ = TMath::Sqrt(dYz * dYz + dyZ * dyZ);
333 } else {
334 XZ_FirstSecondTracksDistanceXZ = 0;
335 YZ_FirstSecondTracksDistanceYZ = 0;
336 }
337
339 XZ_TotalEnergyX = 0;
340 YZ_TotalEnergyY = 0;
341
342 for (auto pair : energiesX) {
343 XZ_TotalEnergyX += pair.second;
344 }
345 for (auto pair : energiesY) {
346 YZ_TotalEnergyY += pair.second;
347 }
348
350 TRestHits hits;
351 TRestHits* hitsXZ = nullptr;
352 TRestHits* hitsYZ = nullptr;
353 if (fTrackEvent->GetMaxEnergyTrack("X")) hitsXZ = fTrackEvent->GetMaxEnergyTrack("X")->GetHits();
354 if (fTrackEvent->GetMaxEnergyTrack("Y")) hitsYZ = fTrackEvent->GetMaxEnergyTrack("Y")->GetHits();
355
356 auto hitsBoth = {hitsXZ, hitsYZ};
357
358 for (auto arg : hitsBoth) {
359 if (arg == nullptr) {
360 continue;
361 }
362 for (int n = 0; n < int(arg->GetNumberOfHits()); n++) {
363 // your code in the existing loop, replacing `hits` by `arg`
364 Double_t eDep = arg->GetEnergy(n);
365 Double_t x = arg->GetX(n);
366 Double_t y = arg->GetY(n);
367 Double_t z = arg->GetZ(n);
368 auto time = arg->GetTime(n);
369 auto type = arg->GetType(n);
370
371 hits.AddHit({x, y, z}, eDep, time, type);
372 }
373 }
374
375 MaxTrack_XZ_YZ_SigmaZ = hits.GetSigmaZ2();
376 MaxTrack_XZ_YZ_GaussSigmaZ = hits.GetGaussSigmaZ();
377 MaxTrack_XZ_YZ_SkewXY = hits.GetSkewXY();
378 MaxTrack_XZ_YZ_SkewZ = hits.GetSkewZ();
379
383
384 // --- Number of tracks and energy --- //
385 SetObservableValue("NTracks", fTrackEvent->GetNumberOfTracks());
386 SetObservableValue("NTracksX", NTracksX);
387 SetObservableValue("NTracksY", NTracksY);
388 SetObservableValue("TotalEnergy", fTrackEvent->GetEnergy());
389 SetObservableValue("XZ_TotalEnergyX", XZ_TotalEnergyX);
390 SetObservableValue("YZ_TotalEnergyY", YZ_TotalEnergyY);
391
392 // --- Map observables --- //
393 SetObservableValue("Map_XZ_NHitsX", XZ_NHitsX);
394 SetObservableValue("Map_XZ_EnergyX", XZ_EnergyX);
395 SetObservableValue("Map_XZ_SigmaX", XZ_SigmaX);
396 SetObservableValue("Map_XZ_SigmaZ", XZ_SigmaZ);
397 SetObservableValue("Map_XZ_GaussSigmaX", XZ_GaussSigmaX);
398 SetObservableValue("Map_XZ_GaussSigmaZ", XZ_GaussSigmaZ);
399 SetObservableValue("Map_XZ_LengthX", XZ_LengthX);
400 SetObservableValue("Map_XZ_VolumeX", XZ_VolumeX);
401 SetObservableValue("Map_XZ_MeanX", XZ_MeanX);
402 SetObservableValue("Map_XZ_MeanZ", XZ_MeanZ);
403 SetObservableValue("Map_XZ_SkewZ", XZ_SkewZ);
404
405 SetObservableValue("Map_YZ_NHitsY", YZ_NHitsY);
406 SetObservableValue("Map_YZ_EnergyY", YZ_EnergyY);
407 SetObservableValue("Map_YZ_SigmaY", YZ_SigmaY);
408 SetObservableValue("Map_YZ_SigmaZ", YZ_SigmaZ);
409 SetObservableValue("Map_YZ_GaussSigmaY", YZ_GaussSigmaY);
410 SetObservableValue("Map_YZ_GaussSigmaZ", YZ_GaussSigmaZ);
411 SetObservableValue("Map_YZ_LengthY", YZ_LengthY);
412 SetObservableValue("Map_YZ_VolumeY", YZ_VolumeY);
413 SetObservableValue("Map_YZ_MeanY", YZ_MeanY);
414 SetObservableValue("Map_YZ_MeanZ", YZ_MeanZ);
415 SetObservableValue("Map_YZ_SkewZ", YZ_SkewZ);
416
417 SetObservableValue("Map_XZ_YZ_SigmaXYBalance", XZ_YZ_SigmaXYBalance);
418 SetObservableValue("Map_XZ_YZ_SigmaZBalance", XZ_YZ_SigmaZBalance);
419 SetObservableValue("Map_XZ_YZ_GaussSigmaXYBalance", XZ_YZ_GaussSigmaXYBalance);
420 SetObservableValue("Map_XZ_YZ_GaussSigmaZBalance", XZ_YZ_GaussSigmaZBalance);
421
422 // --- Max track observables --- //
423
424 // Check if the maximum (final) energy value in energiesX and Y is non-zero
425 bool hasNonZeroEnergyX = !energiesX.empty() && energiesX[0].second != 0;
426 bool hasNonZeroEnergyY = !energiesY.empty() && energiesY[0].second != 0;
427
428 SetObservableValue("MaxTrack_XZ_OK", hasNonZeroEnergyX);
429
430 if (hasNonZeroEnergyX) {
431 int energiesX0FirstKey =
432 energiesX[0].first; // Declare Keys outside to avoid error when accessing "energiesX[0].first"...
433
434 SetObservableValue("MaxTrack_XZ_NHitsX", XZ_NHitsX[energiesX0FirstKey]);
435 SetObservableValue("MaxTrack_XZ_EnergyX", XZ_EnergyX[energiesX0FirstKey]);
436 SetObservableValue("MaxTrack_XZ_SigmaX", XZ_SigmaX[energiesX0FirstKey]);
437 SetObservableValue("MaxTrack_XZ_SigmaZ", XZ_SigmaZ[energiesX0FirstKey]);
438 SetObservableValue("MaxTrack_XZ_GaussSigmaX", XZ_GaussSigmaX[energiesX0FirstKey]);
439 SetObservableValue("MaxTrack_XZ_GaussSigmaZ", XZ_GaussSigmaZ[energiesX0FirstKey]);
440 SetObservableValue("MaxTrack_XZ_LengthX", XZ_LengthX[energiesX0FirstKey]);
441 SetObservableValue("MaxTrack_XZ_VolumeX", XZ_VolumeX[energiesX0FirstKey]);
442 SetObservableValue("MaxTrack_XZ_MeanX", XZ_MeanX[energiesX0FirstKey]);
443 SetObservableValue("MaxTrack_XZ_MeanZ", XZ_MeanZ[energiesX0FirstKey]);
444 SetObservableValue("MaxTrack_XZ_SkewZ", XZ_SkewZ[energiesX0FirstKey]);
445 } else {
446 SetObservableValue("MaxTrack_XZ_NHitsX", 0);
447 SetObservableValue("MaxTrack_XZ_EnergyX", 0.);
448 SetObservableValue("MaxTrack_XZ_SigmaX", 0.);
449 SetObservableValue("MaxTrack_XZ_SigmaZ", 0.);
450 SetObservableValue("MaxTrack_XZ_GaussSigmaX", 0.);
451 SetObservableValue("MaxTrack_XZ_GaussSigmaZ", 0.);
452 SetObservableValue("MaxTrack_XZ_LengthX", 0.);
453 SetObservableValue("MaxTrack_XZ_VolumeX", 0.);
454 SetObservableValue("MaxTrack_XZ_MeanX", 0.);
455 SetObservableValue("MaxTrack_XZ_MeanZ", 0.);
456 SetObservableValue("MaxTrack_XZ_SkewZ", 0.);
457 }
458
459 SetObservableValue("MaxTrack_YZ_OK", hasNonZeroEnergyY);
460
461 if (hasNonZeroEnergyY) {
462 int energiesY0FirstKey = energiesY[0].first;
463
464 SetObservableValue("MaxTrack_YZ_NHitsY", YZ_NHitsY[energiesY0FirstKey]);
465 SetObservableValue("MaxTrack_YZ_EnergyY", YZ_EnergyY[energiesY0FirstKey]);
466 SetObservableValue("MaxTrack_YZ_SigmaY", YZ_SigmaY[energiesY0FirstKey]);
467 SetObservableValue("MaxTrack_YZ_SigmaZ", YZ_SigmaZ[energiesY0FirstKey]);
468 SetObservableValue("MaxTrack_YZ_GaussSigmaY", YZ_GaussSigmaY[energiesY0FirstKey]);
469 SetObservableValue("MaxTrack_YZ_GaussSigmaZ", YZ_GaussSigmaZ[energiesY0FirstKey]);
470 SetObservableValue("MaxTrack_YZ_LengthY", YZ_LengthY[energiesY0FirstKey]);
471 SetObservableValue("MaxTrack_YZ_VolumeY", YZ_VolumeY[energiesY0FirstKey]);
472 SetObservableValue("MaxTrack_YZ_MeanY", YZ_MeanY[energiesY0FirstKey]);
473 SetObservableValue("MaxTrack_YZ_MeanZ", YZ_MeanZ[energiesY0FirstKey]);
474 SetObservableValue("MaxTrack_YZ_SkewZ", YZ_SkewZ[energiesY0FirstKey]);
475 } else {
476 SetObservableValue("MaxTrack_YZ_NHitsY", 0);
477 SetObservableValue("MaxTrack_YZ_EnergyY", 0.);
478 SetObservableValue("MaxTrack_YZ_SigmaY", 0.);
479 SetObservableValue("MaxTrack_YZ_SigmaZ", 0.);
480 SetObservableValue("MaxTrack_YZ_GaussSigmaY", 0.);
481 SetObservableValue("MaxTrack_YZ_GaussSigmaZ", 0.);
482 SetObservableValue("MaxTrack_YZ_LengthY", 0.);
483 SetObservableValue("MaxTrack_YZ_VolumeY", 0.);
484 SetObservableValue("MaxTrack_YZ_MeanY", 0.);
485 SetObservableValue("MaxTrack_YZ_MeanZ", 0.);
486 SetObservableValue("MaxTrack_YZ_SkewZ", 0.);
487 }
488
489 SetObservableValue("MaxTrack_XZ_YZ_OK", hasNonZeroEnergyX && hasNonZeroEnergyY);
490
491 if (!energiesX.empty() && !energiesY.empty()) {
492 Double_t energiesX0SecondKey = energiesX[0].second;
493 Double_t energiesY0SecondKey = energiesY[0].second;
494
495 SetObservableValue("MaxTrack_XZ_YZ_Energy", energiesX0SecondKey + energiesY0SecondKey);
496 SetObservableValue("MaxTrack_XZ_YZ_MaxTrackEnergyPercentage",
497 (energiesX0SecondKey + energiesY0SecondKey) / fTrackEvent->GetEnergy());
498 SetObservableValue("MaxTrack_XZ_YZ_EnergyBalanceXY", (energiesX0SecondKey - energiesY0SecondKey) /
499 (energiesX0SecondKey + energiesY0SecondKey));
500
501 } else {
502 SetObservableValue("MaxTrack_XZ_YZ_Energy", 0.);
503 SetObservableValue("MaxTrack_XZ_YZ_MaxTrackEnergyPercentage", 0.);
504 SetObservableValue("MaxTrack_XZ_YZ_EnergyBalanceXY", 0.);
505 }
506
507 SetObservableValue("MaxTrack_XZ_YZ_SigmaXYBalance", XZ_YZ_SigmaXYBalance[0]);
508 SetObservableValue("MaxTrack_XZ_YZ_SigmaZBalance", XZ_YZ_SigmaZBalance[0]);
509 SetObservableValue("MaxTrack_XZ_YZ_GaussSigmaXYBalance", XZ_YZ_GaussSigmaXYBalance[0]);
510 SetObservableValue("MaxTrack_XZ_YZ_GaussSigmaZBalance", XZ_YZ_GaussSigmaZBalance[0]);
511
512 // --- Second max track observables --- //
513
514 // Check if the second maximum energy value in energiesX and Y is non-zero
515 bool hasNonZeroSecondMaxEnergyX = energiesX.size() > 1 && energiesX[1].second != 0;
516 bool hasNonZeroSecondMaxEnergyY = energiesY.size() > 1 && energiesY[1].second != 0;
517
518 SetObservableValue("SecondMaxTrack_XZ_OK", hasNonZeroSecondMaxEnergyX);
519
520 // Copy the SecondTrack keys immediately after checking the vector
521 if (hasNonZeroSecondMaxEnergyX) {
522 int energiesX1FirstKey =
523 energiesX[1].first; // Declare Keys outside to avoid error when accessing "energiesX[1].first"...
524
525 SetObservableValue("SecondMaxTrack_XZ_NHitsX", XZ_NHitsX[energiesX1FirstKey]);
526 SetObservableValue("SecondMaxTrack_XZ_EnergyX", XZ_EnergyX[energiesX1FirstKey]);
527 SetObservableValue("SecondMaxTrack_XZ_SigmaX", XZ_SigmaX[energiesX1FirstKey]);
528 SetObservableValue("SecondMaxTrack_XZ_SigmaZ", XZ_SigmaZ[energiesX1FirstKey]);
529 SetObservableValue("SecondMaxTrack_XZ_GaussSigmaX", XZ_GaussSigmaX[energiesX1FirstKey]);
530 SetObservableValue("SecondMaxTrack_XZ_GaussSigmaZ", XZ_GaussSigmaZ[energiesX1FirstKey]);
531 SetObservableValue("SecondMaxTrack_XZ_LengthX", XZ_LengthX[energiesX1FirstKey]);
532 SetObservableValue("SecondMaxTrack_XZ_VolumeX", XZ_VolumeX[energiesX1FirstKey]);
533 SetObservableValue("SecondMaxTrack_XZ_MeanX", XZ_MeanX[energiesX1FirstKey]);
534 SetObservableValue("SecondMaxTrack_XZ_MeanZ", XZ_MeanZ[energiesX1FirstKey]);
535 SetObservableValue("SecondMaxTrack_XZ_SkewZ", XZ_SkewZ[energiesX1FirstKey]);
536 } else {
537 SetObservableValue("SecondMaxTrack_XZ_NHitsX", 0);
538 SetObservableValue("SecondMaxTrack_XZ_EnergyX", 0.);
539 SetObservableValue("SecondMaxTrack_XZ_SigmaX", 0.);
540 SetObservableValue("SecondMaxTrack_XZ_SigmaZ", 0.);
541 SetObservableValue("SecondMaxTrack_XZ_GaussSigmaX", 0.);
542 SetObservableValue("SecondMaxTrack_XZ_GaussSigmaZ", 0.);
543 SetObservableValue("SecondMaxTrack_XZ_LengthX", 0.);
544 SetObservableValue("SecondMaxTrack_XZ_VolumeX", 0.);
545 SetObservableValue("SecondMaxTrack_XZ_MeanX", 0.);
546 SetObservableValue("SecondMaxTrack_XZ_MeanZ", 0.);
547 SetObservableValue("SecondMaxTrack_XZ_SkewZ", 0.);
548 }
549
550 SetObservableValue("SecondMaxTrack_YZ_OK", hasNonZeroSecondMaxEnergyY);
551
552 if (hasNonZeroSecondMaxEnergyY) {
553 int energiesY1FirstKey = energiesY[1].first;
554
555 SetObservableValue("SecondMaxTrack_YZ_NHitsY", YZ_NHitsY[energiesY1FirstKey]);
556 SetObservableValue("SecondMaxTrack_YZ_EnergyY", YZ_EnergyY[energiesY1FirstKey]);
557 SetObservableValue("SecondMaxTrack_YZ_SigmaY", YZ_SigmaY[energiesY1FirstKey]);
558 SetObservableValue("SecondMaxTrack_YZ_SigmaZ", YZ_SigmaZ[energiesY1FirstKey]);
559 SetObservableValue("SecondMaxTrack_YZ_GaussSigmaY", YZ_GaussSigmaY[energiesY1FirstKey]);
560 SetObservableValue("SecondMaxTrack_YZ_GaussSigmaZ", YZ_GaussSigmaZ[energiesY1FirstKey]);
561 SetObservableValue("SecondMaxTrack_YZ_LengthY", YZ_LengthY[energiesY1FirstKey]);
562 SetObservableValue("SecondMaxTrack_YZ_VolumeY", YZ_VolumeY[energiesY1FirstKey]);
563 SetObservableValue("SecondMaxTrack_YZ_MeanY", YZ_MeanY[energiesY1FirstKey]);
564 SetObservableValue("SecondMaxTrack_YZ_MeanZ", YZ_MeanZ[energiesY1FirstKey]);
565 SetObservableValue("SecondMaxTrack_YZ_SkewZ", YZ_SkewZ[energiesY1FirstKey]);
566 } else {
567 SetObservableValue("SecondMaxTrack_YZ_NHitsY", 0);
568 SetObservableValue("SecondMaxTrack_YZ_EnergyY", 0.);
569 SetObservableValue("SecondMaxTrack_YZ_SigmaY", 0.);
570 SetObservableValue("SecondMaxTrack_YZ_SigmaZ", 0.);
571 SetObservableValue("SecondMaxTrack_YZ_GaussSigmaY", 0.);
572 SetObservableValue("SecondMaxTrack_YZ_GaussSigmaZ", 0.);
573 SetObservableValue("SecondMaxTrack_YZ_LengthY", 0.);
574 SetObservableValue("SecondMaxTrack_YZ_VolumeY", 0.);
575 SetObservableValue("SecondMaxTrack_YZ_MeanY", 0.);
576 SetObservableValue("SecondMaxTrack_YZ_MeanZ", 0.);
577 SetObservableValue("SecondMaxTrack_YZ_SkewZ", 0.);
578 }
579
580 SetObservableValue("SecondMaxTrack_XZ_YZ_SigmaXYBalance", XZ_YZ_SigmaXYBalance[1]);
581 SetObservableValue("SecondMaxTrack_XZ_YZ_SigmaZBalance", XZ_YZ_SigmaZBalance[1]);
582 SetObservableValue("SecondMaxTrack_XZ_YZ_GaussSigmaXYBalance", XZ_YZ_GaussSigmaXYBalance[1]);
583 SetObservableValue("SecondMaxTrack_XZ_YZ_GaussSigmaZBalance", XZ_YZ_GaussSigmaZBalance[1]);
584
585 SetObservableValue("SecondMaxTrack_XZ_YZ_OK", hasNonZeroSecondMaxEnergyX && hasNonZeroSecondMaxEnergyY);
586
587 if (energiesY.size() > 1 && energiesX.size() > 1) {
588 Double_t energiesX1SecondKey = energiesX[1].second;
589 Double_t energiesY1SecondKey = energiesY[1].second;
590
591 if (fTrackEvent->GetNumberOfTracks() > 2) {
592 SetObservableValue("SecondMaxTrack_XZ_YZ_Energy", energiesX1SecondKey + energiesY1SecondKey);
593 SetObservableValue("SecondMaxTrack_XZ_YZ_EnergyPercentage",
594 (energiesX1SecondKey + energiesY1SecondKey) / fTrackEvent->GetEnergy());
595 SetObservableValue(
596 "SecondMaxTrack_XZ_YZ_EnergyBalanceXY",
597 (energiesX1SecondKey - energiesY1SecondKey) / (energiesX1SecondKey + energiesY1SecondKey));
598 } else {
599 SetObservableValue("SecondMaxTrack_XZ_YZ_Energy", 0.);
600 SetObservableValue("SecondMaxTrack_XZ_YZ_EnergyPercentage", 0.);
601 SetObservableValue("SecondMaxTrack_XZ_YZ_EnergyBalanceXY", 0.);
602 }
603 } else {
604 SetObservableValue("SecondMaxTrack_XZ_YZ_Energy", 0.);
605 SetObservableValue("SecondMaxTrack_XZ_YZ_EnergyPercentage", 0.);
606 SetObservableValue("SecondMaxTrack_XZ_YZ_EnergyBalanceXY", 0.);
607 }
608
609 // --- Distance observables between first two tracks --- //
610 SetObservableValue("XZ_FirstSecondTracksDistanceXZ", XZ_FirstSecondTracksDistanceXZ);
611 SetObservableValue("YZ_FirstSecondTracksDistanceYZ", YZ_FirstSecondTracksDistanceYZ);
612 SetObservableValue("XZ_YZ_FirstSecondTracksDistanceSum",
613 TMath::Sqrt(XZ_FirstSecondTracksDistanceXZ * XZ_FirstSecondTracksDistanceXZ +
614 YZ_FirstSecondTracksDistanceYZ * YZ_FirstSecondTracksDistanceYZ));
615
616 // --- Observables merging max tracks XZ and YZ --- //
617 SetObservableValue("MaxTrack_XZ_YZ_SigmaZ", MaxTrack_XZ_YZ_SigmaZ);
618 SetObservableValue("MaxTrack_XZ_YZ_GaussSigmaZ", MaxTrack_XZ_YZ_GaussSigmaZ);
619 SetObservableValue("MaxTrack_XZ_YZ_SkewXY", MaxTrack_XZ_YZ_SkewXY);
620 SetObservableValue("MaxTrack_XZ_YZ_SkewZ", MaxTrack_XZ_YZ_SkewZ);
621
622 return fTrackEvent;
623}
624
625void TRestTrack2DAnalysisProcess::EndProcess() {}
626
627void TRestTrack2DAnalysisProcess::InitFromConfigFile() {}
TRestEventProcess::SetObservableValue
void SetObservableValue(const std::string &name, const T &value)
Set observable value for AnalysisTree.
Definition: TRestEventProcess.h:151
TRestEvent
A base class for any REST event.
Definition: TRestEvent.h:38
TRestHits
It saves a 3-coordinate position and an energy for each punctual deposition.
Definition: TRestHits.h:39
TRestHits::GetSkewXY
Double_t GetSkewXY() const
It returns the 2-dimensional skewness on the XY-plane which is a measure of the hits distribution asy...
Definition: TRestHits.cxx:979
TRestHits::GetGaussSigmaX
Double_t GetGaussSigmaX(Double_t error=150.0, Int_t nHitsMin=100000)
It computes the gaussian sigma in the X-coordinate. It adds a hit to the right and a hit to the left,...
Definition: TRestHits.cxx:762
TRestHits::GetSigmaX
Double_t GetSigmaX() const
It calculates the hits standard deviation in the X-coordinate.
Definition: TRestHits.cxx:685
TRestHits::GetGaussSigmaZ
Double_t GetGaussSigmaZ(Double_t error=150.0, Int_t nHitsMin=100000)
It computes the gaussian sigma in the Z-coordinate. It adds a hit to the right and a hit to the left,...
Definition: TRestHits.cxx:907
TRestHits::AddHit
void AddHit(const TVector3 &position, Double_t energy, Double_t time=0, REST_HitType type=XYZ)
Adds a new hit to the list of hits using a TVector3.
Definition: TRestHits.cxx:345
TRestHits::GetGaussSigmaY
Double_t GetGaussSigmaY(Double_t error=150.0, Int_t nHitsMin=100000)
It computes the gaussian sigma in the Y-coordinate. It adds a hit to the right and a hit to the left,...
Definition: TRestHits.cxx:836
TRestHits::GetSigmaZ2
Double_t GetSigmaZ2() const
It returns the hits distribution variance on the Z-axis.
Definition: TRestHits.cxx:997
TRestHits::GetSkewZ
Double_t GetSkewZ() const
It returns the hits distribution skewness, or asymmetry on the Z-axis.
Definition: TRestHits.cxx:1013
TRestHits::GetSigmaY
Double_t GetSigmaY() const
It calculates the hits standard deviation in the Y-coordinate.
Definition: TRestHits.cxx:703
TRestMetadata::LoadConfigFromFile
Int_t LoadConfigFromFile(const std::string &configFilename, const std::string &sectionName="")
Give the file name, find out the corresponding section. Then call the main starter.
Definition: TRestMetadata.cxx:574
TRestMetadata::SetLibraryVersion
void SetLibraryVersion(TString version)
Set the library version of this metadata class.
Definition: TRestMetadata.cxx:2322
TRestMetadata::SetSectionName
void SetSectionName(std::string sName)
set the section name, clear the section content
Definition: TRestMetadata.h:329
TRestTrack2DAnalysisProcess
An analysis REST process to extract valuable information from Track type of data.
Definition: TRestTrack2DAnalysisProcess.h:31
TRestTrack2DAnalysisProcess::ProcessEvent
TRestEvent * ProcessEvent(TRestEvent *inputEvent) override
Process one event.
Definition: TRestTrack2DAnalysisProcess.cxx:137
TRestTrack2DAnalysisProcess::InitFromConfigFile
void InitFromConfigFile() override
To make settings from rml file. This method must be implemented in the derived class.
Definition: TRestTrack2DAnalysisProcess.cxx:627
TRestTrack2DAnalysisProcess::InitProcess
void InitProcess() override
To be executed at the beginning of the run (outside event loop)
Definition: TRestTrack2DAnalysisProcess.cxx:135
TRestTrack2DAnalysisProcess::EndProcess
void EndProcess() override
To be executed at the end of the run (outside event loop)
Definition: TRestTrack2DAnalysisProcess.cxx:625
TRestTrack2DAnalysisProcess::Initialize
void Initialize() override
Making default settings.
Definition: TRestTrack2DAnalysisProcess.cxx:124
Generated on Sun Dec 22 2024 05:10:16 for REST-for-Physics by doxygen 1.9.4