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"
191 using namespace std;
192 
193 ClassImp(TRestTrackAnalysisProcess);
194 
195 TRestTrackAnalysisProcess::TRestTrackAnalysisProcess() { Initialize(); }
196 
197 TRestTrackAnalysisProcess::TRestTrackAnalysisProcess(const char* configFilename) {
198  Initialize();
199 
200  if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
201 }
202 
203 TRestTrackAnalysisProcess::~TRestTrackAnalysisProcess() { delete fOutputTrackEvent; }
204 
205 void TRestTrackAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
206 
207 void 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 
219 void TRestTrackAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
220  if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
221 }
222 
223 void 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 
435 TRestEvent* 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 
1257 void 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 
1266 void 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 }
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::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::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
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 Nov 25 2024 10:53:04 for REST-for-Physics by doxygen 1.9.1