REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
TRestTrackAnalysisProcess.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
189
190#include "TRestTrackAnalysisProcess.h"
191using namespace std;
192
193ClassImp(TRestTrackAnalysisProcess);
194
195TRestTrackAnalysisProcess::TRestTrackAnalysisProcess() { Initialize(); }
196
197TRestTrackAnalysisProcess::TRestTrackAnalysisProcess(const char* configFilename) {
198 Initialize();
199
200 if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
201}
202
203TRestTrackAnalysisProcess::~TRestTrackAnalysisProcess() { delete fOutputTrackEvent; }
204
205void TRestTrackAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
206
207void TRestTrackAnalysisProcess::Initialize() {
208 SetSectionName(this->ClassName());
209 SetLibraryVersion(LIBRARY_VERSION);
210
211 fInputTrackEvent = nullptr;
212 fOutputTrackEvent = new TRestTrackEvent();
213
214 fCutsEnabled = false;
215
216 fEnableTwistParameters = false;
217}
218
219void TRestTrackAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
220 if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
221}
222
223void TRestTrackAnalysisProcess::InitProcess() {
224 std::vector<string> fObservables;
225 fObservables = TRestEventProcess::ReadObservables();
226
227 for (unsigned int i = 0; i < fObservables.size(); i++)
228 if (fObservables[i].find("nTracks_LE_") != string::npos) {
229 Double_t energy = StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
230
231 fTrack_LE_EnergyObservables.emplace_back(fObservables[i]);
232 fTrack_LE_Threshold.emplace_back(energy);
233 nTracks_LE.emplace_back(0);
234 }
235
236 for (unsigned int i = 0; i < fObservables.size(); i++)
237 if (fObservables[i].find("nTracks_HE_") != string::npos) {
238 Double_t energy = StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
239
240 fTrack_HE_EnergyObservables.emplace_back(fObservables[i]);
241 fTrack_HE_Threshold.emplace_back(energy);
242 nTracks_HE.emplace_back(0);
243 }
244
245 for (unsigned int i = 0; i < fObservables.size(); i++)
246 if (fObservables[i].find("nTracks_En_") != string::npos) {
247 Double_t energy = StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
248
249 fTrack_En_EnergyObservables.emplace_back(fObservables[i]);
250 fTrack_En_Threshold.emplace_back(energy);
251 nTracks_En.emplace_back(0);
252 }
253
254 for (unsigned int i = 0; i < fObservables.size(); i++)
255 if (fObservables[i].find("twistLow_") != string::npos) {
256 Double_t tailPercentage =
257 StringToDouble(fObservables[i].substr(9, fObservables[i].length()).c_str());
258
259 fTwistLowObservables.emplace_back(fObservables[i]);
260 fTwistLowTailPercentage.emplace_back(tailPercentage);
261 fTwistLowValue.emplace_back(0);
262 }
263
264 for (unsigned int i = 0; i < fObservables.size(); i++)
265 if (fObservables[i].find("twistHigh_") != string::npos) {
266 Double_t tailPercentage =
267 StringToDouble(fObservables[i].substr(10, fObservables[i].length()).c_str());
268
269 fTwistHighObservables.emplace_back(fObservables[i]);
270 fTwistHighTailPercentage.emplace_back(tailPercentage);
271 fTwistHighValue.emplace_back(0);
272 }
273
274 for (unsigned int i = 0; i < fObservables.size(); i++)
275 if (fObservables[i].find("twistBalance_") != string::npos) {
276 Double_t tailPercentage =
277 StringToDouble(fObservables[i].substr(13, fObservables[i].length()).c_str());
278
279 fTwistBalanceObservables.emplace_back(fObservables[i]);
280 fTwistBalanceTailPercentage.emplace_back(tailPercentage);
281 fTwistBalanceValue.emplace_back(0);
282 }
283
284 for (unsigned int i = 0; i < fObservables.size(); i++)
285 if (fObservables[i].find("twistRatio_") != string::npos) {
286 Double_t tailPercentage =
287 StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
288
289 fTwistRatioObservables.emplace_back(fObservables[i]);
290 fTwistRatioTailPercentage.emplace_back(tailPercentage);
291 fTwistRatioValue.emplace_back(0);
292 }
293
294 for (unsigned int i = 0; i < fObservables.size(); i++)
295 if (fObservables[i].find("twistWeightedLow_") != string::npos) {
296 Double_t tailPercentage =
297 StringToDouble(fObservables[i].substr(17, fObservables[i].length()).c_str());
298
299 fTwistWeightedLowObservables.emplace_back(fObservables[i]);
300 fTwistWeightedLowTailPercentage.emplace_back(tailPercentage);
301 fTwistWeightedLowValue.emplace_back(0);
302 }
303
304 for (unsigned int i = 0; i < fObservables.size(); i++)
305 if (fObservables[i].find("twistWeightedHigh_") != string::npos) {
306 Double_t tailPercentage =
307 StringToDouble(fObservables[i].substr(18, fObservables[i].length()).c_str());
308
309 fTwistWeightedHighObservables.emplace_back(fObservables[i]);
310 fTwistWeightedHighTailPercentage.emplace_back(tailPercentage);
311 fTwistWeightedHighValue.emplace_back(0);
312 }
313
314 for (unsigned int i = 0; i < fObservables.size(); i++)
315 if (fObservables[i].find("twistLow_X_") != string::npos) {
316 Double_t tailPercentage =
317 StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
318
319 fTwistLowObservables_X.emplace_back(fObservables[i]);
320 fTwistLowTailPercentage_X.emplace_back(tailPercentage);
321 fTwistLowValue_X.emplace_back(0);
322 }
323
324 for (unsigned int i = 0; i < fObservables.size(); i++)
325 if (fObservables[i].find("twistHigh_X_") != string::npos) {
326 Double_t tailPercentage =
327 StringToDouble(fObservables[i].substr(12, fObservables[i].length()).c_str());
328
329 fTwistHighObservables_X.emplace_back(fObservables[i]);
330 fTwistHighTailPercentage_X.emplace_back(tailPercentage);
331 fTwistHighValue_X.emplace_back(0);
332 }
333
334 for (unsigned int i = 0; i < fObservables.size(); i++)
335 if (fObservables[i].find("twistBalance_X_") != string::npos) {
336 Double_t tailPercentage =
337 StringToDouble(fObservables[i].substr(15, fObservables[i].length()).c_str());
338
339 fTwistBalanceObservables_X.emplace_back(fObservables[i]);
340 fTwistBalanceTailPercentage_X.emplace_back(tailPercentage);
341 fTwistBalanceValue_X.emplace_back(0);
342 }
343
344 for (unsigned int i = 0; i < fObservables.size(); i++)
345 if (fObservables[i].find("twistRatio_X_") != string::npos) {
346 Double_t tailPercentage =
347 StringToDouble(fObservables[i].substr(13, fObservables[i].length()).c_str());
348
349 fTwistRatioObservables_X.emplace_back(fObservables[i]);
350 fTwistRatioTailPercentage_X.emplace_back(tailPercentage);
351 fTwistRatioValue_X.emplace_back(0);
352 }
353
354 for (unsigned int i = 0; i < fObservables.size(); i++)
355 if (fObservables[i].find("twistWeightedLow_X_") != string::npos) {
356 Double_t tailPercentage =
357 StringToDouble(fObservables[i].substr(19, fObservables[i].length()).c_str());
358
359 fTwistWeightedLowObservables_X.emplace_back(fObservables[i]);
360 fTwistWeightedLowTailPercentage_X.emplace_back(tailPercentage);
361 fTwistWeightedLowValue_X.emplace_back(0);
362 }
363
364 for (unsigned int i = 0; i < fObservables.size(); i++)
365 if (fObservables[i].find("twistWeightedHigh_X_") != string::npos) {
366 Double_t tailPercentage =
367 StringToDouble(fObservables[i].substr(20, fObservables[i].length()).c_str());
368
369 fTwistWeightedHighObservables_X.emplace_back(fObservables[i]);
370 fTwistWeightedHighTailPercentage_X.emplace_back(tailPercentage);
371 fTwistWeightedHighValue_X.emplace_back(0);
372 }
373
374 for (unsigned int i = 0; i < fObservables.size(); i++)
375 if (fObservables[i].find("twistLow_Y_") != string::npos) {
376 Double_t tailPercentage =
377 StringToDouble(fObservables[i].substr(11, fObservables[i].length()).c_str());
378
379 fTwistLowObservables_Y.emplace_back(fObservables[i]);
380 fTwistLowTailPercentage_Y.emplace_back(tailPercentage);
381 fTwistLowValue_Y.emplace_back(0);
382 }
383
384 for (unsigned int i = 0; i < fObservables.size(); i++)
385 if (fObservables[i].find("twistHigh_Y_") != string::npos) {
386 Double_t tailPercentage =
387 StringToDouble(fObservables[i].substr(12, fObservables[i].length()).c_str());
388
389 fTwistHighObservables_Y.emplace_back(fObservables[i]);
390 fTwistHighTailPercentage_Y.emplace_back(tailPercentage);
391 fTwistHighValue_Y.emplace_back(0);
392 }
393
394 for (unsigned int i = 0; i < fObservables.size(); i++)
395 if (fObservables[i].find("twistBalance_Y_") != string::npos) {
396 Double_t tailPercentage =
397 StringToDouble(fObservables[i].substr(15, fObservables[i].length()).c_str());
398
399 fTwistBalanceObservables_Y.emplace_back(fObservables[i]);
400 fTwistBalanceTailPercentage_Y.emplace_back(tailPercentage);
401 fTwistBalanceValue_Y.emplace_back(0);
402 }
403
404 for (unsigned int i = 0; i < fObservables.size(); i++)
405 if (fObservables[i].find("twistRatio_Y_") != string::npos) {
406 Double_t tailPercentage =
407 StringToDouble(fObservables[i].substr(13, fObservables[i].length()).c_str());
408
409 fTwistRatioObservables_Y.emplace_back(fObservables[i]);
410 fTwistRatioTailPercentage_Y.emplace_back(tailPercentage);
411 fTwistRatioValue_Y.emplace_back(0);
412 }
413
414 for (unsigned int i = 0; i < fObservables.size(); i++)
415 if (fObservables[i].find("twistWeightedLow_Y_") != string::npos) {
416 Double_t tailPercentage =
417 StringToDouble(fObservables[i].substr(19, fObservables[i].length()).c_str());
418
419 fTwistWeightedLowObservables_Y.emplace_back(fObservables[i]);
420 fTwistWeightedLowTailPercentage_Y.emplace_back(tailPercentage);
421 fTwistWeightedLowValue_Y.emplace_back(0);
422 }
423
424 for (unsigned int i = 0; i < fObservables.size(); i++)
425 if (fObservables[i].find("twistWeightedHigh_Y_") != string::npos) {
426 Double_t tailPercentage =
427 StringToDouble(fObservables[i].substr(20, fObservables[i].length()).c_str());
428
429 fTwistWeightedHighObservables_Y.emplace_back(fObservables[i]);
430 fTwistWeightedHighTailPercentage_Y.emplace_back(tailPercentage);
431 fTwistWeightedHighValue_Y.emplace_back(0);
432 }
433}
434
435TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
436 fInputTrackEvent = (TRestTrackEvent*)inputEvent;
437
438 // Copying the input tracks to the output track
439 for (int tck = 0; tck < fInputTrackEvent->GetNumberOfTracks(); tck++)
440 fOutputTrackEvent->AddTrack(fInputTrackEvent->GetTrack(tck));
441
442 if (this->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
443 fInputTrackEvent->PrintOnlyTracks();
444
445 /* {{{ Number of tracks observables */
446 Int_t nTracksX = 0, nTracksY = 0, nTracksXYZ = 0;
447 nTracksX = fInputTrackEvent->GetNumberOfTracks("X");
448 nTracksY = fInputTrackEvent->GetNumberOfTracks("Y");
449 nTracksXYZ = fInputTrackEvent->GetNumberOfTracks("XYZ");
450
451 SetObservableValue((string) "trackEnergy", fInputTrackEvent->GetEnergy(""));
452
453 SetObservableValue((string) "nTracks_X", nTracksX);
454 SetObservableValue((string) "nTracks_Y", nTracksY);
455 SetObservableValue((string) "nTracks_XYZ", nTracksXYZ);
456 /* }}} */
457
458 /* {{{ Producing nTracks above/below threshold ( nTracks_LE/HE_XXX ) */
459 for (unsigned int n = 0; n < nTracks_HE.size(); n++) nTracks_HE[n] = 0;
460
461 for (unsigned int n = 0; n < nTracks_LE.size(); n++) nTracks_LE[n] = 0;
462
463 for (unsigned int n = 0; n < nTracks_En.size(); n++) nTracks_En[n] = 0;
464
465 for (int tck = 0; tck < fInputTrackEvent->GetNumberOfTracks(); tck++) {
466 if (!fInputTrackEvent->isTopLevel(tck)) continue;
467
468 TRestTrack* t = fInputTrackEvent->GetTrack(tck);
469 Double_t en = t->GetEnergy();
470
471 if ((t->isXYZ()) || (t->isXZ()) || (t->isYZ())) {
472 // cout<<"tracks "<<endl;
473 for (unsigned int n = 0; n < fTrack_HE_EnergyObservables.size(); n++)
474 if (en > fTrack_HE_Threshold[n]) nTracks_HE[n]++;
475
476 for (unsigned int n = 0; n < fTrack_LE_EnergyObservables.size(); n++)
477 if (en < fTrack_LE_Threshold[n]) nTracks_LE[n]++;
478
479 for (unsigned int n = 0; n < fTrack_En_EnergyObservables.size(); n++)
480 if (en > fTrack_En_Threshold[n] - fDeltaEnergy && en < fTrack_En_Threshold[n] + fDeltaEnergy)
481 nTracks_En[n]++;
482 }
483 }
484
485 for (unsigned int n = 0; n < fTrack_LE_EnergyObservables.size(); n++) {
486 SetObservableValue(fTrack_LE_EnergyObservables[n], nTracks_LE[n]);
487 }
488
489 for (unsigned int n = 0; n < fTrack_HE_EnergyObservables.size(); n++) {
490 SetObservableValue(fTrack_HE_EnergyObservables[n], nTracks_HE[n]);
491 }
492
493 for (unsigned int n = 0; n < fTrack_En_EnergyObservables.size(); n++) {
494 SetObservableValue(fTrack_En_EnergyObservables[n], nTracks_En[n]);
495 }
496 /* }}} */
497
498 TRestTrack* tXYZ = fInputTrackEvent->GetMaxEnergyTrack("XYZ");
499 TRestTrack* tX = fInputTrackEvent->GetMaxEnergyTrack("X");
500 TRestTrack* tY = fInputTrackEvent->GetMaxEnergyTrack("Y");
501
502 if (fEnableTwistParameters) {
503 /* {{{ Adding twist observables from XYZ track */
504
505 Double_t twist = -1, twistWeighted = -1;
506
507 for (unsigned int n = 0; n < fTwistWeightedHighValue.size(); n++) fTwistWeightedHighValue[n] = -1;
508
509 for (unsigned int n = 0; n < fTwistWeightedLowValue.size(); n++) fTwistWeightedLowValue[n] = -1;
510
511 for (unsigned int n = 0; n < fTwistLowValue.size(); n++) fTwistLowValue[n] = -1;
512
513 for (unsigned int n = 0; n < fTwistHighValue.size(); n++) fTwistHighValue[n] = -1;
514
515 for (unsigned int n = 0; n < fTwistBalanceValue.size(); n++) fTwistBalanceValue[n] = -1;
516
517 for (unsigned int n = 0; n < fTwistRatioValue.size(); n++) fTwistRatioValue[n] = -1;
518
519 if (tXYZ) {
520 TRestVolumeHits* hits = tXYZ->GetVolumeHits();
521 Int_t Nhits = hits->GetNumberOfHits();
522
523 twist = hits->GetHitsTwist(0, 0);
524 twistWeighted = hits->GetHitsTwistWeighted(0, 0);
525
526 for (unsigned int n = 0; n < fTwistLowObservables.size(); n++) {
527 Int_t Nend = fTwistLowTailPercentage[n] * Nhits / 100.;
528 Double_t twistStart = hits->GetHitsTwist(0, Nend);
529 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
530
531 if (twistStart < twistEnd)
532 fTwistLowValue[n] = twistStart;
533 else
534 fTwistLowValue[n] = twistEnd;
535 }
536
537 for (unsigned int n = 0; n < fTwistHighObservables.size(); n++) {
538 Int_t Nend = fTwistHighTailPercentage[n] * Nhits / 100.;
539 Double_t twistStart = hits->GetHitsTwist(0, Nend);
540 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
541
542 if (twistStart > twistEnd)
543 fTwistHighValue[n] = twistStart;
544 else
545 fTwistHighValue[n] = twistEnd;
546 }
547
548 for (unsigned int n = 0; n < fTwistBalanceObservables.size(); n++) {
549 Int_t Nend = fTwistBalanceTailPercentage[n] * Nhits / 100.;
550 Double_t twistStart = hits->GetHitsTwist(0, Nend);
551 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
552
553 if (twistStart < twistEnd)
554 fTwistBalanceValue[n] = (twistEnd - twistStart) / (twistEnd + twistStart);
555 else
556 fTwistBalanceValue[n] = (twistStart - twistEnd) / (twistEnd + twistStart);
557 }
558
559 for (unsigned int n = 0; n < fTwistRatioObservables.size(); n++) {
560 Int_t Nend = fTwistRatioTailPercentage[n] * Nhits / 100.;
561 Double_t twistStart = hits->GetHitsTwist(0, Nend);
562 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
563
564 if (twistStart > twistEnd) {
565 if (twistEnd <= 0) twistEnd = -1;
566 fTwistRatioValue[n] = twistStart / twistEnd;
567 } else {
568 if (twistStart <= 0) twistStart = -1;
569 fTwistRatioValue[n] = twistEnd / twistStart;
570 }
571 }
572
573 for (unsigned int n = 0; n < fTwistWeightedLowObservables.size(); n++) {
574 Int_t Nend = fTwistWeightedLowTailPercentage[n] * Nhits / 100.;
575 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
576 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
577
578 if (twistStart < twistEnd)
579 fTwistWeightedLowValue[n] = twistStart;
580 else
581 fTwistWeightedLowValue[n] = twistEnd;
582 }
583
584 for (unsigned int n = 0; n < fTwistWeightedHighObservables.size(); n++) {
585 Int_t Nend = fTwistWeightedHighTailPercentage[n] * Nhits / 100.;
586 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
587 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
588
589 if (twistStart > twistEnd)
590 fTwistWeightedHighValue[n] = twistStart;
591 else
592 fTwistWeightedHighValue[n] = twistEnd;
593 }
594 }
595
596 for (unsigned int n = 0; n < fTwistLowObservables.size(); n++) {
597 SetObservableValue((string)fTwistLowObservables[n], fTwistLowValue[n]);
598 }
599
600 for (unsigned int n = 0; n < fTwistHighObservables.size(); n++) {
601 SetObservableValue((string)fTwistHighObservables[n], fTwistHighValue[n]);
602 }
603
604 for (unsigned int n = 0; n < fTwistBalanceObservables.size(); n++) {
605 SetObservableValue((string)fTwistBalanceObservables[n], fTwistBalanceValue[n]);
606 }
607
608 for (unsigned int n = 0; n < fTwistRatioObservables.size(); n++) {
609 SetObservableValue((string)fTwistRatioObservables[n], fTwistRatioValue[n]);
610 }
611
612 for (unsigned int n = 0; n < fTwistWeightedLowObservables.size(); n++) {
613 SetObservableValue((string)fTwistWeightedLowObservables[n], fTwistWeightedLowValue[n]);
614 }
615
616 for (unsigned int n = 0; n < fTwistWeightedHighObservables.size(); n++) {
617 SetObservableValue((string)fTwistWeightedHighObservables[n], fTwistWeightedHighValue[n]);
618 }
619
620 SetObservableValue("twist", twist);
621
622 SetObservableValue("twistWeighted", twistWeighted);
623 /* }}} */
624
625 /* {{{ Adding twist observables from X track */
626 Double_t twist_X = -1, twistWeighted_X = -1;
627
628 for (unsigned int n = 0; n < fTwistWeightedHighValue_X.size(); n++) fTwistWeightedHighValue_X[n] = -1;
629
630 for (unsigned int n = 0; n < fTwistWeightedLowValue_X.size(); n++) fTwistWeightedLowValue_X[n] = -1;
631
632 for (unsigned int n = 0; n < fTwistLowValue_X.size(); n++) fTwistLowValue_X[n] = -1;
633
634 for (unsigned int n = 0; n < fTwistHighValue_X.size(); n++) fTwistHighValue_X[n] = -1;
635
636 for (unsigned int n = 0; n < fTwistBalanceValue_X.size(); n++) fTwistBalanceValue_X[n] = -1;
637
638 for (unsigned int n = 0; n < fTwistRatioValue_X.size(); n++) fTwistRatioValue_X[n] = -1;
639
640 if (tX) {
641 TRestVolumeHits* hits = tX->GetVolumeHits();
642 Int_t Nhits = hits->GetNumberOfHits();
643
644 twist_X = hits->GetHitsTwist(0, 0);
645 twistWeighted_X = hits->GetHitsTwistWeighted(0, 0);
646
647 for (unsigned int n = 0; n < fTwistLowObservables_X.size(); n++) {
648 Int_t Nend = fTwistLowTailPercentage_X[n] * Nhits / 100.;
649 Double_t twistStart = hits->GetHitsTwist(0, Nend);
650 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
651
652 if (twistStart < twistEnd)
653 fTwistLowValue_X[n] = twistStart;
654 else
655 fTwistLowValue_X[n] = twistEnd;
656 }
657
658 for (unsigned int n = 0; n < fTwistHighObservables_X.size(); n++) {
659 Int_t Nend = fTwistHighTailPercentage_X[n] * Nhits / 100.;
660 Double_t twistStart = hits->GetHitsTwist(0, Nend);
661 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
662
663 if (twistStart > twistEnd)
664 fTwistHighValue_X[n] = twistStart;
665 else
666 fTwistHighValue_X[n] = twistEnd;
667 }
668
669 for (unsigned int n = 0; n < fTwistBalanceObservables_X.size(); n++) {
670 Int_t Nend = fTwistBalanceTailPercentage_X[n] * Nhits / 100.;
671 Double_t twistStart = hits->GetHitsTwist(0, Nend);
672 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
673
674 if (twistStart < twistEnd)
675 fTwistBalanceValue_X[n] = (twistEnd - twistStart) / (twistEnd + twistStart);
676 else
677 fTwistBalanceValue_X[n] = (twistStart - twistEnd) / (twistEnd + twistStart);
678 }
679
680 for (unsigned int n = 0; n < fTwistRatioObservables_X.size(); n++) {
681 Int_t Nend = fTwistRatioTailPercentage_X[n] * Nhits / 100.;
682 Double_t twistStart = hits->GetHitsTwist(0, Nend);
683 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
684
685 if (twistStart > twistEnd) {
686 if (twistEnd <= 0) twistEnd = -1;
687 fTwistRatioValue[n] = twistStart / twistEnd;
688 } else {
689 if (twistStart <= 0) twistStart = -1;
690 fTwistRatioValue[n] = twistEnd / twistStart;
691 }
692 }
693
694 for (unsigned int n = 0; n < fTwistWeightedLowObservables_X.size(); n++) {
695 Int_t Nend = fTwistWeightedLowTailPercentage_X[n] * Nhits / 100.;
696 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
697 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
698
699 if (twistStart < twistEnd)
700 fTwistWeightedLowValue_X[n] = twistStart;
701 else
702 fTwistWeightedLowValue_X[n] = twistEnd;
703 }
704
705 for (unsigned int n = 0; n < fTwistWeightedHighObservables_X.size(); n++) {
706 Int_t Nend = fTwistWeightedHighTailPercentage_X[n] * Nhits / 100.;
707 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
708 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
709
710 if (twistStart > twistEnd)
711 fTwistWeightedHighValue_X[n] = twistStart;
712 else
713 fTwistWeightedHighValue_X[n] = twistEnd;
714 }
715 }
716
717 for (unsigned int n = 0; n < fTwistLowObservables_X.size(); n++) {
718 SetObservableValue((string)fTwistLowObservables_X[n], fTwistLowValue_X[n]);
719 }
720
721 for (unsigned int n = 0; n < fTwistHighObservables_X.size(); n++) {
722 SetObservableValue((string)fTwistHighObservables_X[n], fTwistHighValue_X[n]);
723 }
724
725 for (unsigned int n = 0; n < fTwistBalanceObservables_X.size(); n++) {
726 SetObservableValue((string)fTwistBalanceObservables_X[n], fTwistBalanceValue_X[n]);
727 }
728
729 for (unsigned int n = 0; n < fTwistRatioObservables_X.size(); n++) {
730 SetObservableValue((string)fTwistRatioObservables_X[n], fTwistRatioValue_X[n]);
731 }
732
733 for (unsigned int n = 0; n < fTwistWeightedLowObservables_X.size(); n++) {
734 SetObservableValue((string)fTwistWeightedLowObservables_X[n], fTwistWeightedLowValue_X[n]);
735 }
736
737 for (unsigned int n = 0; n < fTwistWeightedHighObservables_X.size(); n++) {
738 SetObservableValue((string)fTwistWeightedHighObservables_X[n], fTwistWeightedHighValue_X[n]);
739 }
740
741 SetObservableValue((string) "twist_X", twist_X);
742
743 SetObservableValue((string) "twistWeighted_X", twistWeighted_X);
744 /* }}} */
745
746 /* {{{ Adding twist observables from Y track */
747 Double_t twist_Y = -1, twistWeighted_Y = -1;
748
749 for (unsigned int n = 0; n < fTwistWeightedHighValue_Y.size(); n++) fTwistWeightedHighValue_Y[n] = -1;
750
751 for (unsigned int n = 0; n < fTwistWeightedLowValue_Y.size(); n++) fTwistWeightedLowValue_Y[n] = -1;
752
753 for (unsigned int n = 0; n < fTwistLowValue_Y.size(); n++) fTwistLowValue_Y[n] = -1;
754
755 for (unsigned int n = 0; n < fTwistHighValue_Y.size(); n++) fTwistHighValue_Y[n] = -1;
756
757 for (unsigned int n = 0; n < fTwistBalanceValue_Y.size(); n++) fTwistBalanceValue_Y[n] = -1;
758
759 for (unsigned int n = 0; n < fTwistRatioValue_Y.size(); n++) fTwistRatioValue_Y[n] = -1;
760
761 if (tY) {
762 TRestVolumeHits* hits = tY->GetVolumeHits();
763 Int_t Nhits = hits->GetNumberOfHits();
764
765 twist_Y = hits->GetHitsTwist(0, 0);
766 twistWeighted_Y = hits->GetHitsTwistWeighted(0, 0);
767
768 for (unsigned int n = 0; n < fTwistLowObservables_Y.size(); n++) {
769 Int_t Nend = fTwistLowTailPercentage_Y[n] * Nhits / 100.;
770 Double_t twistStart = hits->GetHitsTwist(0, Nend);
771 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
772
773 if (twistStart < twistEnd)
774 fTwistLowValue_Y[n] = twistStart;
775 else
776 fTwistLowValue_Y[n] = twistEnd;
777 }
778
779 for (unsigned int n = 0; n < fTwistHighObservables_Y.size(); n++) {
780 Int_t Nend = fTwistHighTailPercentage_Y[n] * Nhits / 100.;
781 Double_t twistStart = hits->GetHitsTwist(0, Nend);
782 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
783
784 if (twistStart > twistEnd)
785 fTwistHighValue_Y[n] = twistStart;
786 else
787 fTwistHighValue_Y[n] = twistEnd;
788 }
789
790 for (unsigned int n = 0; n < fTwistBalanceObservables_Y.size(); n++) {
791 Int_t Nend = fTwistBalanceTailPercentage_Y[n] * Nhits / 100.;
792 Double_t twistStart = hits->GetHitsTwist(0, Nend);
793 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
794
795 if (twistStart < twistEnd)
796 fTwistBalanceValue_Y[n] = (twistEnd - twistStart) / (twistEnd + twistStart);
797 else
798 fTwistBalanceValue_Y[n] = (twistStart - twistEnd) / (twistEnd + twistStart);
799 }
800
801 for (unsigned int n = 0; n < fTwistRatioObservables_Y.size(); n++) {
802 Int_t Nend = fTwistRatioTailPercentage_Y[n] * Nhits / 100.;
803 Double_t twistStart = hits->GetHitsTwist(0, Nend);
804 Double_t twistEnd = hits->GetHitsTwist(Nhits - Nend, Nhits);
805
806 if (twistStart > twistEnd) {
807 if (twistEnd <= 0) twistEnd = -1;
808 fTwistRatioValue[n] = twistStart / twistEnd;
809 } else {
810 if (twistStart <= 0) twistStart = -1;
811 fTwistRatioValue[n] = twistEnd / twistStart;
812 }
813 }
814
815 for (unsigned int n = 0; n < fTwistWeightedLowObservables_Y.size(); n++) {
816 Int_t Nend = fTwistWeightedLowTailPercentage_Y[n] * Nhits / 100.;
817 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
818 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
819
820 if (twistStart < twistEnd)
821 fTwistWeightedLowValue_Y[n] = twistStart;
822 else
823 fTwistWeightedLowValue_Y[n] = twistEnd;
824 }
825
826 for (unsigned int n = 0; n < fTwistWeightedHighObservables_Y.size(); n++) {
827 Int_t Nend = fTwistWeightedHighTailPercentage_Y[n] * Nhits / 100.;
828 Double_t twistStart = hits->GetHitsTwistWeighted(0, Nend);
829 Double_t twistEnd = hits->GetHitsTwistWeighted(Nhits - Nend, Nhits);
830
831 if (twistStart > twistEnd)
832 fTwistWeightedHighValue_Y[n] = twistStart;
833 else
834 fTwistWeightedHighValue_Y[n] = twistEnd;
835 }
836 }
837
838 for (unsigned int n = 0; n < fTwistLowObservables_Y.size(); n++) {
839 SetObservableValue((string)fTwistLowObservables_Y[n], fTwistLowValue_Y[n]);
840 }
841
842 for (unsigned int n = 0; n < fTwistHighObservables_Y.size(); n++) {
843 SetObservableValue((string)fTwistHighObservables_Y[n], fTwistHighValue_Y[n]);
844 }
845
846 for (unsigned int n = 0; n < fTwistBalanceObservables_Y.size(); n++) {
847 SetObservableValue((string)fTwistBalanceObservables_Y[n], fTwistBalanceValue_Y[n]);
848 }
849
850 for (unsigned int n = 0; n < fTwistRatioObservables_Y.size(); n++) {
851 SetObservableValue((string)fTwistRatioObservables_Y[n], fTwistRatioValue_Y[n]);
852 }
853
854 for (unsigned int n = 0; n < fTwistWeightedLowObservables_Y.size(); n++) {
855 SetObservableValue((string)fTwistWeightedLowObservables_Y[n], fTwistWeightedLowValue_Y[n]);
856 }
857
858 for (unsigned int n = 0; n < fTwistWeightedHighObservables_Y.size(); n++) {
859 SetObservableValue((string)fTwistWeightedHighObservables_Y[n], fTwistWeightedHighValue_Y[n]);
860 }
861
862 SetObservableValue("twist_Y", twist_Y);
863
864 SetObservableValue("twistWeighted_Y", twistWeighted_Y);
865 /* }}} */
866 }
867
868 /* {{{ Getting max track energies and track energy ratio */
869 Double_t tckMaxEnXYZ = 0, tckMaxEnX = 0, tckMaxEnY = 0;
870 Double_t tckMaxXYZ_SigmaX = 0, tckMaxXYZ_SigmaY = 0, tckMaxXYZ_SigmaZ = 0;
871 Double_t tckMaxXZ_SigmaX = 0;
872 Double_t tckMaxXZ_SigmaZ = 0;
873 Double_t tckMaxYZ_SigmaY = 0;
874 Double_t tckMaxYZ_SigmaZ = 0;
875 Double_t tckMaxXZ_nHits = 0;
876 Double_t tckMaxYZ_nHits = 0;
877 Double_t tckMaxXYZ_gausSigmaX = 0, tckMaxXYZ_gausSigmaY = 0, tckMaxXYZ_gausSigmaZ = 0;
878 Double_t tckMaxXZ_gausSigmaX = 0, tckMaxXZ_gausSigmaZ_XZ = 0;
879 Double_t tckMaxYZ_gausSigmaY = 0, tckMaxYZ_gausSigmaZ_YZ = 0;
880 Double_t tckMaxTrack_XYZ_GaussSigmaZ = 0;
881 Double_t tckMaxTrack_XYZ_SigmaZ2 = 0;
882 Double_t tckMaxTrack_XYZ_skewXY = 0, tckMaxTrack_XYZ_skewZ = 0;
883
884 if (fInputTrackEvent->GetMaxEnergyTrack()) {
885 tckMaxEnXYZ = fInputTrackEvent->GetMaxEnergyTrack()->GetEnergy();
886 tckMaxXYZ_SigmaX = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetSigmaX();
887 tckMaxXYZ_SigmaY = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetSigmaY();
888 tckMaxXYZ_SigmaZ = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetSigmaZ2();
889 RESTDebug << "id: " << fInputTrackEvent->GetID() << " " << fInputTrackEvent->GetSubEventTag()
890 << " tckMaxEnXYZ: " << tckMaxEnXYZ << RESTendl;
891 tckMaxXYZ_gausSigmaX = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetGaussSigmaX();
892 tckMaxXYZ_gausSigmaY = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetGaussSigmaY();
893 tckMaxXYZ_gausSigmaZ = fInputTrackEvent->GetMaxEnergyTrack()->GetHits()->GetGaussSigmaZ();
894 }
895
896 SetObservableValue((string) "MaxTrackSigmaX", tckMaxXYZ_SigmaX);
897 SetObservableValue((string) "MaxTrackSigmaY", tckMaxXYZ_SigmaY);
898 SetObservableValue((string) "MaxTrackSigmaZ", tckMaxXYZ_SigmaZ);
899 SetObservableValue((string) "MaxTrackGaussSigmaX", tckMaxXYZ_gausSigmaX);
900 SetObservableValue((string) "MaxTrackGaussSigmaY", tckMaxXYZ_gausSigmaY);
901 SetObservableValue((string) "MaxTrackGaussSigmaZ", tckMaxXYZ_gausSigmaZ);
902
903 if (fInputTrackEvent->GetMaxEnergyTrack("X")) {
904 tckMaxEnX = fInputTrackEvent->GetMaxEnergyTrack("X")->GetEnergy();
905 tckMaxXZ_SigmaX = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetSigmaX();
906 tckMaxXZ_SigmaZ = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetSigmaZ2();
907 tckMaxXZ_gausSigmaX = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetGaussSigmaX();
908 tckMaxXZ_gausSigmaZ_XZ = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetGaussSigmaZ();
909 tckMaxXZ_nHits = fInputTrackEvent->GetMaxEnergyTrack("X")->GetNumberOfHits();
910 RESTDebug << "id: " << fInputTrackEvent->GetID() << " " << fInputTrackEvent->GetSubEventTag()
911 << " tckMaxEnX: " << tckMaxEnX << RESTendl;
912 }
913
914 SetObservableValue((string) "MaxTrackEnergy_X", tckMaxEnX);
915 SetObservableValue((string) "MaxTrack_XZ_SigmaX", tckMaxXZ_SigmaX);
916 SetObservableValue((string) "MaxTrack_XZ_SigmaZ", tckMaxXZ_SigmaZ);
917 SetObservableValue((string) "MaxTrack_XZ_GaussSigmaX", tckMaxXZ_gausSigmaX);
918 SetObservableValue((string) "MaxTrack_XZ_GaussSigmaZ", tckMaxXZ_gausSigmaZ_XZ);
919 SetObservableValue((string) "MaxTrack_XZ_nHits", tckMaxXZ_nHits);
920
921 if (fInputTrackEvent->GetMaxEnergyTrack("Y")) {
922 tckMaxEnY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetEnergy();
923 tckMaxYZ_SigmaY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetSigmaY();
924 tckMaxYZ_SigmaZ = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetSigmaZ2();
925 tckMaxYZ_gausSigmaY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaY();
926 tckMaxYZ_gausSigmaZ_YZ = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaZ();
927 tckMaxYZ_nHits = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetNumberOfHits();
928 RESTDebug << "id: " << fInputTrackEvent->GetID() << " " << fInputTrackEvent->GetSubEventTag()
929 << " tckMaxEnY: " << tckMaxEnY << RESTendl;
930 }
931
932 SetObservableValue((string) "MaxTrackEnergy_Y", tckMaxEnY);
933 SetObservableValue((string) "MaxTrack_YZ_SigmaY", tckMaxYZ_SigmaY);
934 SetObservableValue((string) "MaxTrack_YZ_SigmaZ", tckMaxYZ_SigmaZ);
935 SetObservableValue((string) "MaxTrack_YZ_GaussSigmaY", tckMaxYZ_gausSigmaY);
936 SetObservableValue((string) "MaxTrack_YZ_GaussSigmaZ", tckMaxYZ_gausSigmaZ_YZ);
937 SetObservableValue((string) "MaxTrack_YZ_nHits", tckMaxYZ_nHits);
938
939 SetObservableValue("MaxTrackxySigmaGausBalance", (tckMaxXZ_gausSigmaX - tckMaxYZ_gausSigmaY) /
940 (tckMaxXZ_gausSigmaX + tckMaxYZ_gausSigmaY));
941
942 SetObservableValue("MaxTrackxySigmaBalance",
943 (tckMaxXZ_SigmaX - tckMaxYZ_SigmaY) / (tckMaxXZ_SigmaX + tckMaxYZ_SigmaY));
944
945 SetObservableValue("MaxTrackEnergyBalanceXY", (tckMaxEnX - tckMaxEnY) / (tckMaxEnX + tckMaxEnY));
946
947 Double_t tckMaxEnergy = tckMaxEnX + tckMaxEnY + tckMaxEnXYZ;
948
949 Double_t totalEnergy = fInputTrackEvent->GetEnergy();
950
951 Double_t trackEnergyRatio = (totalEnergy - tckMaxEnergy) / totalEnergy;
952
953 SetObservableValue((string) "MaxTrackEnergy", tckMaxEnergy);
954 SetObservableValue((string) "MaxTrackEnergyRatio", trackEnergyRatio);
955
956 SetObservableValue("MaxTrackZSigmaGausBalance", (tckMaxXZ_gausSigmaZ_XZ - tckMaxYZ_gausSigmaZ_YZ) /
957 (tckMaxXZ_gausSigmaZ_XZ + tckMaxYZ_gausSigmaZ_YZ));
958
959 SetObservableValue("MaxTrackZSigmaBalance",
960 (tckMaxXZ_SigmaZ - tckMaxYZ_SigmaZ) / (tckMaxXZ_SigmaZ + tckMaxYZ_SigmaZ));
961
962 SetObservableValue("MaxTrackEnergyBalanceXY", (tckMaxEnX - tckMaxEnY) / (tckMaxEnX + tckMaxEnY));
963
964 TRestHits hits;
965 TRestHits* hitsXZ = nullptr;
966 TRestHits* hitsYZ = nullptr;
967 if (fInputTrackEvent->GetMaxEnergyTrack("X"))
968 hitsXZ = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits();
969 if (fInputTrackEvent->GetMaxEnergyTrack("Y"))
970 hitsYZ = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits();
971
972 auto hitsBoth = {hitsXZ, hitsYZ};
973
974 for (auto arg : hitsBoth) {
975 if (arg == nullptr) continue;
976 for (unsigned int n = 0; n < arg->GetNumberOfHits(); n++) {
977 // your code in the existing loop, replacing `hits` by `arg`
978 Double_t eDep = arg->GetEnergy(n);
979
980 Double_t x = arg->GetX(n);
981 Double_t y = arg->GetY(n);
982 Double_t z = arg->GetZ(n);
983
984 auto time = arg->GetTime(n);
985 auto type = arg->GetType(n);
986
987 hits.AddHit({x, y, z}, eDep, time, type);
988 }
989 }
990 tckMaxTrack_XYZ_GaussSigmaZ = hits.GetGaussSigmaZ();
991 tckMaxTrack_XYZ_SigmaZ2 = hits.GetSigmaZ2();
992 tckMaxTrack_XYZ_skewZ = hits.GetSkewZ();
993 tckMaxTrack_XYZ_skewXY = hits.GetSkewXY();
994 SetObservableValue((string) "MaxTrack_XYZ_GaussSigmaZ", tckMaxTrack_XYZ_GaussSigmaZ);
995 SetObservableValue((string) "MaxTrack_XYZ_SigmaZ2", tckMaxTrack_XYZ_SigmaZ2);
996 SetObservableValue((string) "MaxTrack_XYZ_skewZ", tckMaxTrack_XYZ_skewZ);
997 SetObservableValue((string) "MaxTrack_XYZ_skewXY", tckMaxTrack_XYZ_skewXY);
998 /* }}} */
999
1000 /* {{{ Second Maximum Track Energy observable */
1001 Double_t tckSecondMaxXYZ_SigmaX = 0, tckSecondMaxXYZ_SigmaY = 0;
1002 Double_t tckSecondMaxXYZ_gausSigmaX = 0, tckSecondMaxXYZ_gausSigmaY = 0;
1003 Double_t tckSecondMaxXZ_gausSigmaX = 0, tckSecondMaxXZ_gausSigmaZ_XZ = 0;
1004 Double_t tckSecondMaxYZ_gausSigmaY = 0, tckSecondMaxYZ_gausSigmaZ_YZ = 0;
1005
1006 if (fInputTrackEvent->GetSecondMaxEnergyTrack() != nullptr) {
1007 tckSecondMaxXYZ_SigmaX = fInputTrackEvent->GetSecondMaxEnergyTrack()->GetHits()->GetSigmaX();
1008 tckSecondMaxXYZ_SigmaY = fInputTrackEvent->GetSecondMaxEnergyTrack()->GetHits()->GetSigmaY();
1009 tckSecondMaxXYZ_gausSigmaX = fInputTrackEvent->GetSecondMaxEnergyTrack()->GetHits()->GetGaussSigmaX();
1010 tckSecondMaxXYZ_gausSigmaY = fInputTrackEvent->GetSecondMaxEnergyTrack()->GetHits()->GetGaussSigmaY();
1011 }
1012
1013 Double_t tckSecondMaxEnergy_X = 0;
1014 Double_t tckSecondMaxXZ_SigmaX = 0;
1015 Double_t tckSecondMaxXZ_SigmaZ = 0;
1016 if (fInputTrackEvent->GetSecondMaxEnergyTrack("X") != nullptr) {
1017 tckSecondMaxXZ_SigmaX = fInputTrackEvent->GetSecondMaxEnergyTrack("X")->GetHits()->GetSigmaX();
1018 tckSecondMaxXZ_SigmaZ = fInputTrackEvent->GetSecondMaxEnergyTrack("X")->GetHits()->GetSigmaZ2();
1019 tckSecondMaxEnergy_X = fInputTrackEvent->GetSecondMaxEnergyTrack("X")->GetEnergy();
1020 tckSecondMaxXZ_gausSigmaX =
1021 fInputTrackEvent->GetSecondMaxEnergyTrack("X")->GetHits()->GetGaussSigmaX();
1022 tckSecondMaxXZ_gausSigmaZ_XZ =
1023 fInputTrackEvent->GetSecondMaxEnergyTrack("X")->GetHits()->GetGaussSigmaZ();
1024 }
1025
1026 Double_t tckSecondMaxEnergy_Y = 0;
1027 Double_t tckSecondMaxYZ_SigmaY = 0;
1028 Double_t tckSecondMaxYZ_SigmaZ = 0;
1029 if (fInputTrackEvent->GetSecondMaxEnergyTrack("Y") != nullptr) {
1030 tckSecondMaxYZ_SigmaY = fInputTrackEvent->GetSecondMaxEnergyTrack("Y")->GetHits()->GetSigmaY();
1031 tckSecondMaxYZ_SigmaZ = fInputTrackEvent->GetSecondMaxEnergyTrack("Y")->GetHits()->GetSigmaZ2();
1032 tckSecondMaxEnergy_Y = fInputTrackEvent->GetSecondMaxEnergyTrack("Y")->GetEnergy();
1033 tckSecondMaxYZ_gausSigmaY =
1034 fInputTrackEvent->GetSecondMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaY();
1035 tckSecondMaxYZ_gausSigmaZ_YZ =
1036 fInputTrackEvent->GetSecondMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaZ();
1037 }
1038
1039 Double_t trackSecondMaxEnergy = tckSecondMaxEnergy_X + tckSecondMaxEnergy_Y;
1040 if (fInputTrackEvent->GetSecondMaxEnergyTrack("XYZ") != nullptr) {
1041 trackSecondMaxEnergy = fInputTrackEvent->GetSecondMaxEnergyTrack("XYZ")->GetEnergy();
1042 }
1043
1044 SetObservableValue((string) "SecondMaxTrackEnergy", trackSecondMaxEnergy);
1045 SetObservableValue((string) "SecondMaxTrackSigmaX", tckSecondMaxXYZ_SigmaX);
1046 SetObservableValue((string) "SecondMaxTrackSigmaY", tckSecondMaxXYZ_SigmaY);
1047 SetObservableValue((string) "SecondMaxTrackGaussSigmaX", tckSecondMaxXYZ_gausSigmaX);
1048 SetObservableValue((string) "SecondMaxTrackGaussSigmaY", tckSecondMaxXYZ_gausSigmaY);
1049
1050 SetObservableValue((string) "SecondMaxTrackEnergy_X", tckSecondMaxEnergy_X);
1051 SetObservableValue((string) "SecondMaxTrack_XZ_SigmaX", tckSecondMaxXZ_SigmaX);
1052 SetObservableValue((string) "SecondMaxTrack_XZ_SigmaZ", tckSecondMaxXZ_SigmaZ);
1053 SetObservableValue((string) "SecondMaxTrack_XZ_GaussSigmaX", tckSecondMaxXZ_gausSigmaX);
1054 SetObservableValue((string) "SecondMaxTrack_XZ_GaussSigmaZ", tckSecondMaxXZ_gausSigmaZ_XZ);
1055
1056 SetObservableValue((string) "SecondMaxTrackEnergy_Y", tckSecondMaxEnergy_Y);
1057 SetObservableValue((string) "SecondMaxTrack_YZ_SigmaY", tckSecondMaxYZ_SigmaY);
1058 SetObservableValue((string) "SecondMaxTrack_YZ_SigmaZ", tckSecondMaxYZ_SigmaZ);
1059 SetObservableValue((string) "SecondMaxTrack_YZ_GaussSigmaY", tckSecondMaxYZ_gausSigmaY);
1060 SetObservableValue((string) "SecondMaxTrack_YZ_GaussSigmaZ", tckSecondMaxYZ_gausSigmaZ_YZ);
1061
1062 SetObservableValue("SecondMaxTrackxySigmaGausBalance",
1063 (tckSecondMaxXZ_gausSigmaX - tckSecondMaxYZ_gausSigmaY) /
1064 (tckSecondMaxXZ_gausSigmaX + tckSecondMaxYZ_gausSigmaY));
1065
1066 SetObservableValue("SecondMaxTrackxySigmaBalance", (tckSecondMaxXZ_SigmaX - tckSecondMaxYZ_SigmaY) /
1067 (tckSecondMaxXZ_SigmaX + tckSecondMaxYZ_SigmaY));
1068 SetObservableValue("SecondMaxTrackZSigmaBalance", (tckSecondMaxXZ_SigmaZ - tckSecondMaxYZ_SigmaZ) /
1069 (tckSecondMaxXZ_SigmaZ + tckSecondMaxYZ_SigmaZ));
1070 SetObservableValue("SecondMaxTrackZSigmaGausBalance",
1071 (tckSecondMaxXZ_gausSigmaZ_XZ - tckSecondMaxYZ_gausSigmaZ_YZ) /
1072 (tckSecondMaxXZ_gausSigmaZ_XZ + tckSecondMaxYZ_gausSigmaZ_YZ));
1073 SetObservableValue("SecondMaxTrackEnergyBalanceXY", (tckSecondMaxEnergy_X - tckSecondMaxEnergy_Y) /
1074 (tckSecondMaxEnergy_X + tckSecondMaxEnergy_Y));
1075
1076 /* }}} */
1077
1078 /* {{{ Track Length observables (MaxTrackLength_XX) */
1079 Double_t tckLenX = fInputTrackEvent->GetMaxEnergyTrackLength("X");
1080 Double_t tckLenY = fInputTrackEvent->GetMaxEnergyTrackLength("Y");
1081 Double_t tckLenXYZ = fInputTrackEvent->GetMaxEnergyTrackLength();
1082
1083 SetObservableValue((string) "MaxTrackLength_X", tckLenX);
1084 SetObservableValue((string) "MaxTrackLength_Y", tckLenY);
1085 SetObservableValue((string) "MaxTrackLength_XYZ", tckLenXYZ);
1086 /* }}} */
1087
1088 /* {{{ Track Volume observables (MaxTrackVolume_XX) */
1089 Double_t tckVolX = fInputTrackEvent->GetMaxEnergyTrackVolume("X");
1090 Double_t tckVolY = fInputTrackEvent->GetMaxEnergyTrackVolume("Y");
1091 Double_t tckVolXYZ = fInputTrackEvent->GetMaxEnergyTrackVolume();
1092
1093 SetObservableValue((string) "MaxTrackVolume_X", tckVolX);
1094 SetObservableValue((string) "MaxTrackVolume_Y", tckVolY);
1095 SetObservableValue((string) "MaxTrackVolume_XYZ", tckVolXYZ);
1096 /* }}} */
1097
1098 /* {{{ Setting mean position for max energy tracks (MaxTrack_{x,y,z}Mean_XXX)
1099 */
1100
1102 Double_t maxX = 0, maxY = 0, maxZ = 0;
1103 Double_t sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1104 Double_t dX = 0, dY = 0, dZ = 0;
1105 Double_t dXZ = 0, dYZ = 0, dXYZ = 0;
1106
1107 // Main max track
1108 TRestTrack* tMax = fInputTrackEvent->GetMaxEnergyTrack();
1109
1110 if (tMax != nullptr) {
1111 maxX = tMax->GetMeanPosition().X();
1112 maxY = tMax->GetMeanPosition().Y();
1113 maxZ = tMax->GetMeanPosition().Z();
1114 }
1115
1116 SetObservableValue((string) "MaxTrack_Xmean_XYZ", maxX);
1117 SetObservableValue((string) "MaxTrack_Ymean_XYZ", maxY);
1118 SetObservableValue((string) "MaxTrack_Zmean_XYZ", maxZ);
1119
1120 // Second max track
1121 TRestTrack* tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack();
1122
1123 if (tSecondMax != nullptr) {
1124 sMaxX = tSecondMax->GetMeanPosition().X();
1125 sMaxY = tSecondMax->GetMeanPosition().Y();
1126 sMaxZ = tSecondMax->GetMeanPosition().Z();
1127 }
1128
1129 SetObservableValue((string) "SecondMaxTrack_Xmean_XYZ", sMaxX);
1130 SetObservableValue((string) "SecondMaxTrack_Ymean_XYZ", sMaxY);
1131 SetObservableValue((string) "SecondMaxTrack_Zmean_XYZ", sMaxZ);
1132
1133 if (sMaxX != 0) dX = abs(maxX - sMaxX);
1134 if (sMaxY != 0) dY = abs(maxY - sMaxY);
1135 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1136
1137 if (dX != 0 || dY != 0 || dZ != 0) dXYZ = TMath::Sqrt(dX * dX + dY * dY + dZ * dZ);
1138 SetObservableValue((string) "TracksDistance_XYZ", dXYZ);
1139
1141
1142 maxX = 0, maxY = 0, maxZ = 0;
1143 dX = 0, dY = 0, dZ = 0;
1144 tMax = fInputTrackEvent->GetMaxEnergyTrack("X");
1145 if (tMax != nullptr) {
1146 maxX = tMax->GetMeanPosition().X();
1147 maxZ = tMax->GetMeanPosition().Z();
1148 }
1149
1150 SetObservableValue((string) "MaxTrack_Xmean_X", maxX);
1151 SetObservableValue((string) "MaxTrack_Zmean_X", maxZ);
1152
1153 sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1154 tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack("X");
1155 if (tSecondMax != nullptr) {
1156 sMaxX = tSecondMax->GetMeanPosition().X();
1157 sMaxZ = tSecondMax->GetMeanPosition().Z();
1158 }
1159
1160 SetObservableValue((string) "SecondMaxTrack_Xmean_X", sMaxX);
1161 SetObservableValue((string) "SecondMaxTrack_Zmean_X", sMaxZ);
1162
1163 if (sMaxX != 0) dX = abs(maxX - sMaxX);
1164 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1165
1166 if (dX != 0 || dY != 0 || dZ != 0) dXZ = TMath::Sqrt(dX * dX + dZ * dZ);
1167 SetObservableValue((string) "TracksDistance_XZ", dXZ);
1168
1170
1171 maxX = 0, maxY = 0, maxZ = 0;
1172 dX = 0, dY = 0, dZ = 0;
1173
1174 tMax = fInputTrackEvent->GetMaxEnergyTrack("Y");
1175 if (tMax != nullptr) {
1176 maxY = tMax->GetMeanPosition().Y();
1177 maxZ = tMax->GetMeanPosition().Z();
1178 }
1179
1180 SetObservableValue((string) "MaxTrack_Ymean_Y", maxY);
1181 SetObservableValue((string) "MaxTrack_Zmean_Y", maxZ);
1182
1183 sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1184 tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack("Y");
1185 if (tSecondMax != nullptr) {
1186 sMaxY = tSecondMax->GetMeanPosition().Y();
1187 sMaxZ = tSecondMax->GetMeanPosition().Z();
1188 }
1189
1190 SetObservableValue((string) "SecondMaxTrack_Ymean_Y", sMaxY);
1191 SetObservableValue((string) "SecondMaxTrack_Zmean_Y", sMaxZ);
1192
1193 if (sMaxY != 0) dY = abs(maxY - sMaxY);
1194 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1195
1196 if (dX != 0 || dY != 0 || dZ != 0) dYZ = TMath::Sqrt(dY * dY + dZ * dZ);
1197 SetObservableValue((string) "TracksDistance_YZ", dYZ);
1198
1199 Double_t dXZ_YZ = 0.5 * TMath::Sqrt(dYZ * dYZ + dXZ * dXZ);
1200 SetObservableValue((string) "TracksDistance_XZ_YZ", dXZ_YZ);
1201
1203 Double_t x = 0, y = 0, z = 0;
1204
1205 if (tXYZ != nullptr) {
1206 x = tXYZ->GetMeanPosition().X();
1207 y = tXYZ->GetMeanPosition().Y();
1208 z = tXYZ->GetMeanPosition().Z();
1209 } else {
1210 double zxz = 0;
1211 int i = 0;
1212 if (tX != nullptr) {
1213 x = tX->GetMeanPosition().X();
1214 zxz += tX->GetMeanPosition().Z();
1215 i++;
1216 }
1217 if (tY != nullptr) {
1218 y = tY->GetMeanPosition().Y();
1219 zxz += tY->GetMeanPosition().Z();
1220 i++;
1221 }
1222
1223 z = i == 0 ? i : zxz / i;
1224 }
1225
1226 SetObservableValue((string) "xMean", x);
1227
1228 SetObservableValue((string) "yMean", y);
1229
1230 SetObservableValue((string) "zMean", z);
1231 /* }}} */
1232
1235 Double_t evTimeDelay = 0;
1236 if (fPreviousEventTime.size() > 0) evTimeDelay = fInputTrackEvent->GetTime() - fPreviousEventTime.back();
1237 SetObservableValue((string) "EventTimeDelay", evTimeDelay);
1238
1239 Double_t meanRate = 0;
1240 if (fPreviousEventTime.size() == 100)
1241 meanRate = 100. / (fInputTrackEvent->GetTime() - fPreviousEventTime.front());
1242 SetObservableValue((string) "MeanRate_InHz", meanRate);
1243
1244 if (fCutsEnabled) {
1245 if (nTracksX < fNTracksXCut.X() || nTracksX > fNTracksXCut.Y()) return nullptr;
1246 if (nTracksY < fNTracksYCut.X() || nTracksY > fNTracksYCut.Y()) return nullptr;
1247 }
1248
1249 if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
1250
1251 return fOutputTrackEvent;
1252}
1253
1254//
1255// void TRestTrackAnalysisProcess::EndOfEventProcess() {
1256// fPreviousEventTime.emplace_back(fInputTrackEvent->GetTimeStamp());
1257// if (fPreviousEventTime.size() > 100) fPreviousEventTime.erase(fPreviousEventTime.begin());
1258//}
1259
1260void TRestTrackAnalysisProcess::EndProcess() {
1261 // Function to be executed once at the end of the process
1262 // (after all events have been processed)
1263
1264 // Start by calling the EndProcess function of the abstract class.
1265 // Comment this if you don't want it.
1266 // TRestEventProcess::EndProcess();
1267}
1268
1269void TRestTrackAnalysisProcess::InitFromConfigFile() {
1270 fNTracksXCut = StringTo2DVector(GetParameter("nTracksXCut", "(1,10)"));
1271 fNTracksYCut = StringTo2DVector(GetParameter("nTracksYCut", "(1,10)"));
1272 fDeltaEnergy = GetDblParameterWithUnits("deltaEnergy", 1);
1273
1274 if (GetParameter("cutsEnabled", "false") == "true") fCutsEnabled = true;
1275
1276 if (GetParameter("enableTwistParameters", "false") == "true") fEnableTwistParameters = true;
1277}
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::GetHitsTwistWeighted
Double_t GetHitsTwistWeighted(Int_t n, Int_t m) const
Same as GetHitsTwist but weighting with the energy.
Definition: TRestHits.cxx:1325
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::GetHitsTwist
Double_t GetHitsTwist(Int_t n, Int_t m) const
A parameter measuring how straight is a given sequence of hits. If the value is close to zero,...
Definition: TRestHits.cxx:1301
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::RESTendl
endl_t RESTendl
Termination flag object for TRestStringOutput.
Definition: TRestMetadata.h:173
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::GetVerboseLevel
TRestStringOutput::REST_Verbose_Level GetVerboseLevel()
returns the verboselevel in type of REST_Verbose_Level enumerator
Definition: TRestMetadata.h:308
TRestMetadata::SetSectionName
void SetSectionName(std::string sName)
set the section name, clear the section content
Definition: TRestMetadata.h:329
TRestMetadata::GetParameter
std::string GetParameter(std::string parName, TiXmlElement *e, TString defaultValue=PARAMETER_NOT_FOUND_STR)
Returns the value for the parameter named parName in the given section.
Definition: TRestMetadata.cxx:1511
TRestStringOutput::REST_Verbose_Level::REST_Debug
@ REST_Debug
+show the defined debug messages
TRestTrackAnalysisProcess
An analysis REST process to extract valuable information from Track type of data.
Definition: TRestTrackAnalysisProcess.h:31
TRestTrackAnalysisProcess::InitProcess
void InitProcess() override
To be executed at the beginning of the run (outside event loop)
Definition: TRestTrackAnalysisProcess.cxx:223
TRestTrackAnalysisProcess::InitFromConfigFile
void InitFromConfigFile() override
To make settings from rml file. This method must be implemented in the derived class.
Definition: TRestTrackAnalysisProcess.cxx:1269
TRestTrackAnalysisProcess::ProcessEvent
TRestEvent * ProcessEvent(TRestEvent *inputEvent) override
Process one event.
Definition: TRestTrackAnalysisProcess.cxx:435
TRestTrackAnalysisProcess::Initialize
void Initialize() override
Making default settings.
Definition: TRestTrackAnalysisProcess.cxx:207
TRestTrackAnalysisProcess::EndProcess
void EndProcess() override
To be executed at the end of the run (outside event loop)
Definition: TRestTrackAnalysisProcess.cxx:1260
REST_StringHelper::GetChar
Int_t GetChar(std::string hint="Press a KEY to continue ...")
Helps to pause the program, printing a message before pausing.
Definition: TRestStringHelper.cxx:198
REST_StringHelper::StringToDouble
Double_t StringToDouble(std::string in)
Gets a double from a string.
Definition: TRestStringHelper.cxx:602
REST_StringHelper::StringTo2DVector
TVector2 StringTo2DVector(std::string in)
Gets a 2D-vector from a string.
Definition: TRestStringHelper.cxx:723
Generated on Sun Nov 24 2024 05:10:15 for REST-for-Physics by doxygen 1.9.4