REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
Loading...
Searching...
No Matches
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.push_back(fObservables[i]);
232 fTrack_LE_Threshold.push_back(energy);
233 nTracks_LE.push_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.push_back(fObservables[i]);
241 fTrack_HE_Threshold.push_back(energy);
242 nTracks_HE.push_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.push_back(fObservables[i]);
250 fTrack_En_Threshold.push_back(energy);
251 nTracks_En.push_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.push_back(fObservables[i]);
260 fTwistLowTailPercentage.push_back(tailPercentage);
261 fTwistLowValue.push_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.push_back(fObservables[i]);
270 fTwistHighTailPercentage.push_back(tailPercentage);
271 fTwistHighValue.push_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.push_back(fObservables[i]);
280 fTwistBalanceTailPercentage.push_back(tailPercentage);
281 fTwistBalanceValue.push_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.push_back(fObservables[i]);
290 fTwistRatioTailPercentage.push_back(tailPercentage);
291 fTwistRatioValue.push_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.push_back(fObservables[i]);
300 fTwistWeightedLowTailPercentage.push_back(tailPercentage);
301 fTwistWeightedLowValue.push_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.push_back(fObservables[i]);
310 fTwistWeightedHighTailPercentage.push_back(tailPercentage);
311 fTwistWeightedHighValue.push_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.push_back(fObservables[i]);
320 fTwistLowTailPercentage_X.push_back(tailPercentage);
321 fTwistLowValue_X.push_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.push_back(fObservables[i]);
330 fTwistHighTailPercentage_X.push_back(tailPercentage);
331 fTwistHighValue_X.push_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.push_back(fObservables[i]);
340 fTwistBalanceTailPercentage_X.push_back(tailPercentage);
341 fTwistBalanceValue_X.push_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.push_back(fObservables[i]);
350 fTwistRatioTailPercentage_X.push_back(tailPercentage);
351 fTwistRatioValue_X.push_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.push_back(fObservables[i]);
360 fTwistWeightedLowTailPercentage_X.push_back(tailPercentage);
361 fTwistWeightedLowValue_X.push_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.push_back(fObservables[i]);
370 fTwistWeightedHighTailPercentage_X.push_back(tailPercentage);
371 fTwistWeightedHighValue_X.push_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.push_back(fObservables[i]);
380 fTwistLowTailPercentage_Y.push_back(tailPercentage);
381 fTwistLowValue_Y.push_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.push_back(fObservables[i]);
390 fTwistHighTailPercentage_Y.push_back(tailPercentage);
391 fTwistHighValue_Y.push_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.push_back(fObservables[i]);
400 fTwistBalanceTailPercentage_Y.push_back(tailPercentage);
401 fTwistBalanceValue_Y.push_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.push_back(fObservables[i]);
410 fTwistRatioTailPercentage_Y.push_back(tailPercentage);
411 fTwistRatioValue_Y.push_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.push_back(fObservables[i]);
420 fTwistWeightedLowTailPercentage_Y.push_back(tailPercentage);
421 fTwistWeightedLowValue_Y.push_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.push_back(fObservables[i]);
430 fTwistWeightedHighTailPercentage_Y.push_back(tailPercentage);
431 fTwistWeightedHighValue_Y.push_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
1041 SetObservableValue((string) "SecondMaxTrackEnergy", trackSecondMaxEnergy);
1042 SetObservableValue((string) "SecondMaxTrackSigmaX", tckSecondMaxXYZ_SigmaX);
1043 SetObservableValue((string) "SecondMaxTrackSigmaY", tckSecondMaxXYZ_SigmaY);
1044 SetObservableValue((string) "SecondMaxTrackGaussSigmaX", tckSecondMaxXYZ_gausSigmaX);
1045 SetObservableValue((string) "SecondMaxTrackGaussSigmaY", tckSecondMaxXYZ_gausSigmaY);
1046
1047 SetObservableValue((string) "SecondMaxTrackEnergy_X", tckSecondMaxEnergy_X);
1048 SetObservableValue((string) "SecondMaxTrack_XZ_SigmaX", tckSecondMaxXZ_SigmaX);
1049 SetObservableValue((string) "SecondMaxTrack_XZ_SigmaZ", tckSecondMaxXZ_SigmaZ);
1050 SetObservableValue((string) "SecondMaxTrack_XZ_GaussSigmaX", tckSecondMaxXZ_gausSigmaX);
1051 SetObservableValue((string) "SecondMaxTrack_XZ_GaussSigmaZ", tckSecondMaxXZ_gausSigmaZ_XZ);
1052
1053 SetObservableValue((string) "SecondMaxTrackEnergy_Y", tckSecondMaxEnergy_Y);
1054 SetObservableValue((string) "SecondMaxTrack_YZ_SigmaY", tckSecondMaxYZ_SigmaY);
1055 SetObservableValue((string) "SecondMaxTrack_YZ_SigmaZ", tckSecondMaxYZ_SigmaZ);
1056 SetObservableValue((string) "SecondMaxTrack_YZ_GaussSigmaY", tckSecondMaxYZ_gausSigmaY);
1057 SetObservableValue((string) "SecondMaxTrack_YZ_GaussSigmaZ", tckSecondMaxYZ_gausSigmaZ_YZ);
1058
1059 SetObservableValue("SecondMaxTrackxySigmaGausBalance",
1060 (tckSecondMaxXZ_gausSigmaX - tckSecondMaxYZ_gausSigmaY) /
1061 (tckSecondMaxXZ_gausSigmaX + tckSecondMaxYZ_gausSigmaY));
1062
1063 SetObservableValue("SecondMaxTrackxySigmaBalance", (tckSecondMaxXZ_SigmaX - tckSecondMaxYZ_SigmaY) /
1064 (tckSecondMaxXZ_SigmaX + tckSecondMaxYZ_SigmaY));
1065 SetObservableValue("SecondMaxTrackZSigmaBalance", (tckSecondMaxXZ_SigmaZ - tckSecondMaxYZ_SigmaZ) /
1066 (tckSecondMaxXZ_SigmaZ + tckSecondMaxYZ_SigmaZ));
1067 SetObservableValue("SecondMaxTrackZSigmaGausBalance",
1068 (tckSecondMaxXZ_gausSigmaZ_XZ - tckSecondMaxYZ_gausSigmaZ_YZ) /
1069 (tckSecondMaxXZ_gausSigmaZ_XZ + tckSecondMaxYZ_gausSigmaZ_YZ));
1070 SetObservableValue("SecondMaxTrackEnergyBalanceXY", (tckSecondMaxEnergy_X - tckSecondMaxEnergy_Y) /
1071 (tckSecondMaxEnergy_X + tckSecondMaxEnergy_Y));
1072
1073 /* }}} */
1074
1075 /* {{{ Track Length observables (MaxTrackLength_XX) */
1076 Double_t tckLenX = fInputTrackEvent->GetMaxEnergyTrackLength("X");
1077 Double_t tckLenY = fInputTrackEvent->GetMaxEnergyTrackLength("Y");
1078 Double_t tckLenXYZ = fInputTrackEvent->GetMaxEnergyTrackLength();
1079
1080 SetObservableValue((string) "MaxTrackLength_X", tckLenX);
1081 SetObservableValue((string) "MaxTrackLength_Y", tckLenY);
1082 SetObservableValue((string) "MaxTrackLength_XYZ", tckLenXYZ);
1083 /* }}} */
1084
1085 /* {{{ Track Volume observables (MaxTrackVolume_XX) */
1086 Double_t tckVolX = fInputTrackEvent->GetMaxEnergyTrackVolume("X");
1087 Double_t tckVolY = fInputTrackEvent->GetMaxEnergyTrackVolume("Y");
1088 Double_t tckVolXYZ = fInputTrackEvent->GetMaxEnergyTrackVolume();
1089
1090 SetObservableValue((string) "MaxTrackVolume_X", tckVolX);
1091 SetObservableValue((string) "MaxTrackVolume_Y", tckVolY);
1092 SetObservableValue((string) "MaxTrackVolume_XYZ", tckVolXYZ);
1093 /* }}} */
1094
1095 /* {{{ Setting mean position for max energy tracks (MaxTrack_{x,y,z}Mean_XXX)
1096 */
1097
1099 Double_t maxX = 0, maxY = 0, maxZ = 0;
1100 Double_t sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1101 Double_t dX = 0, dY = 0, dZ = 0;
1102 Double_t dXZ = 0, dYZ = 0, dXYZ = 0;
1103
1104 // Main max track
1105 TRestTrack* tMax = fInputTrackEvent->GetMaxEnergyTrack();
1106
1107 if (tMax != nullptr) {
1108 maxX = tMax->GetMeanPosition().X();
1109 maxY = tMax->GetMeanPosition().Y();
1110 maxZ = tMax->GetMeanPosition().Z();
1111 }
1112
1113 SetObservableValue((string) "MaxTrack_Xmean_XYZ", maxX);
1114 SetObservableValue((string) "MaxTrack_Ymean_XYZ", maxY);
1115 SetObservableValue((string) "MaxTrack_Zmean_XYZ", maxZ);
1116
1117 // Second max track
1118 TRestTrack* tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack();
1119
1120 if (tSecondMax != nullptr) {
1121 sMaxX = tSecondMax->GetMeanPosition().X();
1122 sMaxY = tSecondMax->GetMeanPosition().Y();
1123 sMaxZ = tSecondMax->GetMeanPosition().Z();
1124 }
1125
1126 SetObservableValue((string) "SecondMaxTrack_Xmean_XYZ", sMaxX);
1127 SetObservableValue((string) "SecondMaxTrack_Ymean_XYZ", sMaxY);
1128 SetObservableValue((string) "SecondMaxTrack_Zmean_XYZ", sMaxZ);
1129
1130 if (sMaxX != 0) dX = abs(maxX - sMaxX);
1131 if (sMaxY != 0) dY = abs(maxY - sMaxY);
1132 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1133
1134 if (dX != 0 || dY != 0 || dZ != 0) dXYZ = TMath::Sqrt(dX * dX + dY * dY + dZ * dZ);
1135 SetObservableValue((string) "TracksDistance_XYZ", dXYZ);
1136
1138
1139 maxX = 0, maxY = 0, maxZ = 0;
1140 dX = 0, dY = 0, dZ = 0;
1141 tMax = fInputTrackEvent->GetMaxEnergyTrack("X");
1142 if (tMax != nullptr) {
1143 maxX = tMax->GetMeanPosition().X();
1144 maxZ = tMax->GetMeanPosition().Z();
1145 }
1146
1147 SetObservableValue((string) "MaxTrack_Xmean_X", maxX);
1148 SetObservableValue((string) "MaxTrack_Zmean_X", maxZ);
1149
1150 sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1151 tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack("X");
1152 if (tSecondMax != nullptr) {
1153 sMaxX = tSecondMax->GetMeanPosition().X();
1154 sMaxZ = tSecondMax->GetMeanPosition().Z();
1155 }
1156
1157 SetObservableValue((string) "SecondMaxTrack_Xmean_X", sMaxX);
1158 SetObservableValue((string) "SecondMaxTrack_Zmean_X", sMaxZ);
1159
1160 if (sMaxX != 0) dX = abs(maxX - sMaxX);
1161 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1162
1163 if (dX != 0 || dY != 0 || dZ != 0) dXZ = TMath::Sqrt(dX * dX + dZ * dZ);
1164 SetObservableValue((string) "TracksDistance_XZ", dXZ);
1165
1167
1168 maxX = 0, maxY = 0, maxZ = 0;
1169 dX = 0, dY = 0, dZ = 0;
1170
1171 tMax = fInputTrackEvent->GetMaxEnergyTrack("Y");
1172 if (tMax != nullptr) {
1173 maxY = tMax->GetMeanPosition().Y();
1174 maxZ = tMax->GetMeanPosition().Z();
1175 }
1176
1177 SetObservableValue((string) "MaxTrack_Ymean_Y", maxY);
1178 SetObservableValue((string) "MaxTrack_Zmean_Y", maxZ);
1179
1180 sMaxX = 0, sMaxY = 0, sMaxZ = 0;
1181 tSecondMax = fInputTrackEvent->GetSecondMaxEnergyTrack("Y");
1182 if (tSecondMax != nullptr) {
1183 sMaxY = tSecondMax->GetMeanPosition().Y();
1184 sMaxZ = tSecondMax->GetMeanPosition().Z();
1185 }
1186
1187 SetObservableValue((string) "SecondMaxTrack_Ymean_Y", sMaxY);
1188 SetObservableValue((string) "SecondMaxTrack_Zmean_Y", sMaxZ);
1189
1190 if (sMaxY != 0) dY = abs(maxY - sMaxY);
1191 if (sMaxZ != 0) dZ = abs(maxZ - sMaxZ);
1192
1193 if (dX != 0 || dY != 0 || dZ != 0) dYZ = TMath::Sqrt(dY * dY + dZ * dZ);
1194 SetObservableValue((string) "TracksDistance_YZ", dYZ);
1195
1196 Double_t dXZ_YZ = 0.5 * TMath::Sqrt(dYZ * dYZ + dXZ * dXZ);
1197 SetObservableValue((string) "TracksDistance_XZ_YZ", dXZ_YZ);
1198
1200 Double_t x = 0, y = 0, z = 0;
1201
1202 if (tXYZ != nullptr) {
1203 x = tXYZ->GetMeanPosition().X();
1204 y = tXYZ->GetMeanPosition().Y();
1205 z = tXYZ->GetMeanPosition().Z();
1206 } else {
1207 double zxz = 0;
1208 int i = 0;
1209 if (tX != nullptr) {
1210 x = tX->GetMeanPosition().X();
1211 zxz += tX->GetMeanPosition().Z();
1212 i++;
1213 }
1214 if (tY != nullptr) {
1215 y = tY->GetMeanPosition().Y();
1216 zxz += tY->GetMeanPosition().Z();
1217 i++;
1218 }
1219
1220 z = i == 0 ? i : zxz / i;
1221 }
1222
1223 SetObservableValue((string) "xMean", x);
1224
1225 SetObservableValue((string) "yMean", y);
1226
1227 SetObservableValue((string) "zMean", z);
1228 /* }}} */
1229
1232 Double_t evTimeDelay = 0;
1233 if (fPreviousEventTime.size() > 0) evTimeDelay = fInputTrackEvent->GetTime() - fPreviousEventTime.back();
1234 SetObservableValue((string) "EventTimeDelay", evTimeDelay);
1235
1236 Double_t meanRate = 0;
1237 if (fPreviousEventTime.size() == 100)
1238 meanRate = 100. / (fInputTrackEvent->GetTime() - fPreviousEventTime.front());
1239 SetObservableValue((string) "MeanRate_InHz", meanRate);
1240
1241 if (fCutsEnabled) {
1242 if (nTracksX < fNTracksXCut.X() || nTracksX > fNTracksXCut.Y()) return nullptr;
1243 if (nTracksY < fNTracksYCut.X() || nTracksY > fNTracksYCut.Y()) return nullptr;
1244 }
1245
1246 if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
1247
1248 return fOutputTrackEvent;
1249}
1250
1251//
1252// void TRestTrackAnalysisProcess::EndOfEventProcess() {
1253// fPreviousEventTime.push_back(fInputTrackEvent->GetTimeStamp());
1254// if (fPreviousEventTime.size() > 100) fPreviousEventTime.erase(fPreviousEventTime.begin());
1255//}
1256
1257void TRestTrackAnalysisProcess::EndProcess() {
1258 // Function to be executed once at the end of the process
1259 // (after all events have been processed)
1260
1261 // Start by calling the EndProcess function of the abstract class.
1262 // Comment this if you don't want it.
1263 // TRestEventProcess::EndProcess();
1264}
1265
1266void TRestTrackAnalysisProcess::InitFromConfigFile() {
1267 fNTracksXCut = StringTo2DVector(GetParameter("nTracksXCut", "(1,10)"));
1268 fNTracksYCut = StringTo2DVector(GetParameter("nTracksYCut", "(1,10)"));
1269 fDeltaEnergy = GetDblParameterWithUnits("deltaEnergy", 1);
1270
1271 if (GetParameter("cutsEnabled", "false") == "true") fCutsEnabled = true;
1272
1273 if (GetParameter("enableTwistParameters", "false") == "true") fEnableTwistParameters = true;
1274}
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:1266
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:1257
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 Mon Jan 20 2025 12:40:39 for REST-for-Physics by doxygen 1.9.8