eventAnalysis  7.0-49-g0ac7482
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TTRExReconModule.cxx
Go to the documentation of this file.
1 #include <sstream>
2 
3 #include "TDatum.hxx"
4 #include "TGeomInfo.hxx"
5 #include "TPrincipal.h"
6 #include "TRealDatum.hxx"
7 
8 #include "TTRExReconModule.hxx"
9 
13 
14 #define CVSTAG \
15  "\
16  $Name: 7.0$"
17 #define CVSID \
18  "\
19  $Id: eventAnalysis TTRExReconModule.cxx,2024/03/20:09:46:11,Alexander_J_Finch,lapw.lancs.ac.uk $"
20 
21 ND::TTRExReconModule::TTRExReconModule(const char* name, const char* title) {
22  SetNameTitle(name, title);
23  // Enable this module by default:
24  fIsEnabled = kTRUE;
25  fDescription = "TREx output";
27  fCVSID = CVSID;
28  fIsMC = false;
29  // total = 0;
30  // Tree variables
31  fPatterns = new TClonesArray("ND::TTRExReconModule::TTPCAnaPattern", 200);
32 }
33 
35 
36 Bool_t ND::TTRExReconModule::ProcessFirstEvent(ND::TND280Event& event) {
37  return true;
38 }
39 
40 bool ND::TTRExReconModule::IsFullEventWorthSaving(ND::TND280Event& event) {
41  return (fNPatterns ? true : false);
42 }
43 
45 
47  fOutputTree->Branch("NPatterns", &fNPatterns, "NPatterns/I", fBufferSize);
48  // A TClonesArray of TTPCAnaPatterns
49  fOutputTree->Branch("Patterns", &fPatterns, fBufferSize, fSplitLevel);
50 }
51 
52 bool ND::TTRExReconModule::FillTree(ND::TND280Event& event) {
53  fNPatterns = 0;
54  fIsMC = false;
55  fPatterns->Clear("C");
56 
57  // Check for MC-osity
58  fIsMC = event.GetContext().IsMC();
59 
60  // Getting the gas interactions result from TREx
61  ND::THandle<ND::TAlgorithmResult> TRExResult = event.GetFit("TREXRecon");
62 
63  if (TRExResult) {
64  ND::THandle<ND::TReconObjectContainer> PatternReco =
65  TRExResult->GetResultsContainer("TRExUnmergedReco");
66  ND::THandle<ND::TReconObjectContainer> StandardReco =
67  TRExResult->GetResultsContainer("GasInteractionOutput");
68 
69  int pattern_count = 0;
70  if (StandardReco) {
71  //****************************
72  // HANDLING OAEVENT OUTPUT *
73  //****************************
74 
75  // Loop through patterns
76  //std::cout << " TREX patterns available " << std::endl;
77  for (ND::TReconObjectContainer::iterator patternIt =
78  StandardReco->begin();
79  patternIt < StandardReco->end(); patternIt++) {
80  ND::THandle<ND::TReconVertex> c_pattern =
81  *patternIt; // Connected pattern
82  ND::THandle<ND::TReconPID> d_pattern_pid =
83  *patternIt; // Disconnected pattern with PID
84  ND::THandle<ND::TReconTrack> d_pattern_no_pid =
85  *patternIt; // Disconnected pattern with no PID
86  ND::THandle<ND::TReconBase> generic_pattern = *patternIt;
87 
88  // Construct output pattern object
89  TTPCAnaPattern* outputPattern =
90  new ((*fPatterns)[pattern_count]) TTPCAnaPattern;
91  outputPattern->PatternID = 0;
92  if ((c_pattern) && (c_pattern->Get<ND::TIntegerDatum>("PatternId")))
93  outputPattern->PatternID =
94  c_pattern->Get<ND::TIntegerDatum>("PatternId")->GetValue();
95  if ((d_pattern_pid) &&
96  (d_pattern_pid->Get<ND::TIntegerDatum>("PatternId")))
97  outputPattern->PatternID =
98  d_pattern_pid->Get<ND::TIntegerDatum>("PatternId")->GetValue();
99  if ((d_pattern_no_pid) &&
100  (d_pattern_no_pid->Get<ND::TIntegerDatum>("PatternId")))
101  outputPattern->PatternID =
102  d_pattern_no_pid->Get<ND::TIntegerDatum>("PatternId")->GetValue();
103 
104  //std::cout << " TREX pattern available: " << outputPattern->PatternID
105  // << std::endl;
106 
107  // Set S1S success flag
108  //-----------------
109  int tpc = 0;
110  if (c_pattern) tpc = GetTPCFromPattern(c_pattern);
111  if (d_pattern_pid) tpc = GetTPCFromPattern(d_pattern_pid);
112  if (d_pattern_no_pid) tpc = GetTPCFromPattern(d_pattern_no_pid);
113 
114  outputPattern->S1Sflag = true;
115 
116  if (tpc == 1) {
117  // Check for P0D veto
118  ND::THandle<ND::THitSelection> p0dhits = event.GetHitSelection("p0d");
119  if (p0dhits){
120  ND::THandle<ND::TAlgorithmResult> p0dhitsResult(
121  new ND::TAlgorithmResult(*p0dhits));
122  ND::TP0DVetoForTPC vetoForTPC;
123  THandle<TAlgorithmResult> vetoResult =
124  vetoForTPC.Process(*p0dhitsResult);
125  if (vetoResult) outputPattern->S1Sflag = false;
126  }
127  }
128 
129  if (tpc == 2) {
130  // Check for FGD1 veto
131  if (FGDHitVeto(1, 5, event, 0)) outputPattern->S1Sflag = false;
132  }
133 
134  if (tpc == 3) {
135  // Check for FGD2 veto
136  if (FGDHitVeto(2, 5, event, 0)) outputPattern->S1Sflag = false;
137  }
138 
139  //-----------------------
140  // CONNECTED PATTERN CASE
141  //-----------------------
142  if (c_pattern) {
143  outputPattern->TPC = GetTPCFromPattern(c_pattern);
144 
145  //************************************************
146  // LOOP 1: construct lists of paths and junctions *
147  //************************************************
148  std::list<ND::THandle<ND::TReconBase> > patternpaths;
149  std::list<ND::THandle<ND::TReconVertex> > patternjuncs;
150 
151  ND::THandle<ND::TReconObjectContainer> patt_con =
152  c_pattern->GetConstituents();
153  if (patt_con) {
154  for (ND::TReconObjectContainer::iterator conIt = patt_con->begin();
155  conIt < patt_con->end(); conIt++) {
156  ND::THandle<ND::TReconVertex> junc = *conIt;
157  ND::THandle<ND::TReconTrack> track = *conIt;
158  ND::THandle<ND::TReconPID> pid = *conIt;
159  ND::THandle<ND::TReconBase> path = *conIt;
160 
161  if (junc) patternjuncs.push_back(junc);
162  if (track || pid) patternpaths.push_back(path);
163  }
164  }
165 
166  outputPattern->NJunctions = patternjuncs.size();
167  outputPattern->NPaths = patternpaths.size();
168 
169  std::map<int, std::vector<int> > path_junctions;
170 
171  //******************************************
172  // LOOP 2: link paths and junctions via IDs *
173  //******************************************
174  if (patt_con) {
175  int junc_id = 0;
176  int path_id = 0;
177 
178  for (ND::TReconObjectContainer::iterator conIt = patt_con->begin();
179  conIt < patt_con->end(); conIt++) {
180  // LOOP OVER JUNCTIONS
181  //===================
182  ND::THandle<ND::TReconVertex> junc = *conIt;
183  if (junc) {
184  TTPCAnaJunction* outputJunction =
185  new ((*(outputPattern->Junctions))[junc_id])
187 
188  outputJunction->JunctionID = junc_id;
189  outputJunction->JunctionMatchingID = -1;
190  if (junc->Get<ND::TIntegerDatum>("JunctionId"))
191  outputJunction->JunctionMatchingID =
192  junc->Get<ND::TIntegerDatum>("JunctionId")->GetValue();
193  outputJunction->Position = junc->GetPosition();
194 
195  // Examine hits associated to junction
196  ND::THandle<ND::THitSelection> junchits = junc->GetHits();
197  TVector3 maxima(-999999.9, -999999.9, -999999.9);
198  TVector3 minima(999999.9, 999999.9, 999999.9);
199  int hitcounter = 1;
200 
201  if (junchits){
202  for (ND::THitSelection::iterator hitIt = junchits->begin();
203  hitIt != junchits->end(); hitIt++) {
204 
205  hitcounter++;
206 
207  ND::THandle<ND::THit> hit = *hitIt;
208  if (!hit) continue;
209 
210  double t0 = 0;
211  if (junc->Get<ND::TRealDatum>("T0"))
212  t0 = junc->Get<ND::TRealDatum>("T0")->GetValue();
213 
214  // position of the pad that measured the charge, i.e. the read
215  // out plane X.
216  double RPx = hit->GetPosition().X();
217  double Sense = int(ND::TGeomInfo::Get().TPC().GetDriftSense(
218  hit->GetGeomId()));
219 
220  // How far are we from the read out plane ?
221  double TimeOffset = 259.;
222 
223  // TODO: data/mc
224  double DriftVelocity = 0.0785;
225 
226  double DriftDistance =
227  (hit->GetTime() - t0 - TimeOffset) * DriftVelocity;
228  double hitX = RPx - (Sense * DriftDistance);
229 
230  double hitY = hit->GetPosition().Y();
231  double hitZ = hit->GetPosition().Z();
232 
233  if (hitX > maxima.X()) {
234  maxima.SetX(hitX);
235  }
236  if (hitX < minima.X()) {
237  minima.SetX(hitX);
238  }
239 
240  if (hitY > maxima.Y()) maxima.SetY(hitY);
241  if (hitY < minima.Y()) minima.SetY(hitY);
242 
243  if (hitZ > maxima.Z()) maxima.SetZ(hitZ);
244  if (hitZ < minima.Z()) minima.SetZ(hitZ);
245  }
246  }
247  outputJunction->MaximumCoordinates = maxima;
248  outputJunction->MinimumCoordinates = minima;
249 
250  // Loop over paths associated to junction
251  ND::THandle<ND::TReconObjectContainer> juncpaths =
252  junc->GetConstituents();
253  if (juncpaths) {
254  int njuncpaths = juncpaths->size();
255  outputJunction->NPaths = njuncpaths;
256  outputJunction->PathIDs = new int[njuncpaths];
257 
258  int juncpath_count = 0;
259  for (ND::TReconObjectContainer::iterator pathIt =
260  juncpaths->begin();
261  pathIt != juncpaths->end(); pathIt++) {
262  // Get path
263  ND::THandle<ND::TReconBase> juncpath = *pathIt;
264 
265  if (juncpath) {
266  // Loop over path list, find corresponding path ID.
267  for (std::list<ND::THandle<ND::TReconBase> >::iterator
268  listIt = patternpaths.begin();
269  listIt != patternpaths.end(); listIt++) {
270  if (juncpath == *listIt) {
271  // Add path ID to junction
272  int i = std::distance(patternpaths.begin(), listIt);
273  outputJunction->PathIDs[juncpath_count] = i;
274  // Add junction ID to path
275  path_junctions[i].push_back(junc_id);
276  break;
277  }
278  }
279  juncpath_count = juncpath_count + 1;
280  }
281  }
282  }
283  junc_id = junc_id + 1;
284  }
285 
286  // LOOP OVER PATHS
287  //===============
288  ND::THandle<ND::TReconPID> pid = *conIt;
289  ND::THandle<ND::TReconTrack> track = *conIt;
290 
291  if (pid) // We have a TReconPID
292  {
293  TTPCAnaPath* outputPath =
294  new ((*(outputPattern->Paths))[path_id]) TTPCAnaPath;
295 
296  // Fill path information
297  outputPath->PathID = path_id;
298  FillPathInfo(outputPath, pid, patternpaths, StandardReco,
299  generic_pattern);
300 
301  // Identify associated junctions
302  int npathjuncs = path_junctions[path_id].size();
303  outputPath->NJunctions = npathjuncs;
304  outputPath->JunctionIDs = new int[npathjuncs];
305  for (int i = 0; i < npathjuncs; i++) {
306  outputPath->JunctionIDs[i] = path_junctions[path_id][i];
307  }
308 
309  // Increment path count
310  path_id = path_id + 1;
311  }
312 
313  else if (track) // We have a TReconTrack
314  {
315  TTPCAnaPath* outputPath =
316  new ((*(outputPattern->Paths))[path_id]) TTPCAnaPath;
317 
318  // Fill path information
319  outputPath->PathID = path_id;
320  FillPathInfo(outputPath, track, patternpaths, StandardReco,
321  generic_pattern);
322 
323  // Identify associated junctions
324  int npathjuncs = path_junctions[path_id].size();
325  outputPath->NJunctions = npathjuncs;
326  outputPath->JunctionIDs = new int[npathjuncs];
327  for (int i = 0; i < npathjuncs; i++) {
328  outputPath->JunctionIDs[i] = path_junctions[path_id][i];
329  }
330 
331  // Increment path count
332  path_id = path_id + 1;
333  }
334  }
335  }
336  pattern_count = pattern_count + 1;
337  }
338 
339  //-------------------------
340  // DISCONNECTED PATTERN CASE
341  //-------------------------
342 
343  // Make empty list for convenience
344  std::list<ND::THandle<ND::TReconBase> > patternpaths;
345 
346  if (d_pattern_pid) {
347  outputPattern->TPC = GetTPCFromPattern(d_pattern_pid);
348  outputPattern->NJunctions = 0;
349  outputPattern->NPaths = 1;
350 
351  // Get the solitary through-going path
352  TTPCAnaPath* outputPath =
353  new ((*(outputPattern->Paths))[0]) TTPCAnaPath;
354 
355  outputPath->PathID = 0;
356  outputPath->NJunctions = 0;
357  outputPath->JunctionIDs = 0;
358 
359  FillPathInfo(outputPath, d_pattern_pid, patternpaths, StandardReco,
360  generic_pattern);
361 
362  pattern_count = pattern_count + 1;
363  } else if (d_pattern_no_pid) {
364  outputPattern->TPC = GetTPCFromPattern(d_pattern_no_pid);
365  outputPattern->NJunctions = 0;
366  outputPattern->NPaths = 1;
367 
368  // Get the solitary through-going path
369  TTPCAnaPath* outputPath =
370  new ((*(outputPattern->Paths))[0]) TTPCAnaPath;
371 
372  outputPath->PathID = 0;
373  outputPath->NJunctions = 0;
374  outputPath->JunctionIDs = 0;
375 
376  FillPathInfo(outputPath, d_pattern_no_pid, patternpaths, StandardReco,
377  generic_pattern);
378 
379  pattern_count = pattern_count + 1;
380  }
381  }
382 
383  fNPatterns = pattern_count;
384  } else if (PatternReco) {
385  //****************************************************
386  // HANDLING TREX OUTPUT
387  //****************************************************
388 
389  //std::cout << "TTRExReconModule requires oaEvent output." << std::endl;
390  }
391  }
392  return true;
393 }
394 
396  ND::THandle<ND::TReconVertex> pattern) {
397  int TPC = 0;
398  if (pattern->UsesDetector(ND::TReconBase::kTPC1)) TPC = 1;
399  if (pattern->UsesDetector(ND::TReconBase::kTPC2)) TPC = 2;
400  if (pattern->UsesDetector(ND::TReconBase::kTPC3)) TPC = 3;
401  return TPC;
402 }
403 
405  ND::THandle<ND::TReconTrack> pattern) {
406  int TPC = 0;
407  if (pattern->UsesDetector(ND::TReconBase::kTPC1)) TPC = 1;
408  if (pattern->UsesDetector(ND::TReconBase::kTPC2)) TPC = 2;
409  if (pattern->UsesDetector(ND::TReconBase::kTPC3)) TPC = 3;
410  return TPC;
411 }
412 
414  ND::THandle<ND::TReconPID> pattern) {
415  int TPC = 0;
416  if (pattern->UsesDetector(ND::TReconBase::kTPC1)) TPC = 1;
417  if (pattern->UsesDetector(ND::TReconBase::kTPC2)) TPC = 2;
418  if (pattern->UsesDetector(ND::TReconBase::kTPC3)) TPC = 3;
419  return TPC;
420 }
421 
423  TTPCAnaPath* outputPath, ND::THandle<ND::TReconPID> path,
424  std::list<ND::THandle<ND::TReconBase> > patternpaths,
425  ND::THandle<ND::TReconObjectContainer> patterns,
426  ND::THandle<ND::TReconBase> parent_pattern) {
427 
428  outputPath->PathMatchingID =
429  path->Get<ND::TIntegerDatum>("PathId")->GetValue();
430  outputPath->FirstPosition = path->GetPosition().Vect();
431  outputPath->Time = path->GetPosition().T();
432 
433  ND::TReconNodeContainer pathnodes = path->GetNodes();
434  int node_counter = 0;
435  for (ND::TReconNodeContainer::iterator nodeIt = pathnodes.begin();
436  nodeIt != pathnodes.end(); nodeIt++) {
437  node_counter++;
438 
439  ND::THandle<ND::TReconState> recon_state = (*nodeIt)->GetState();
440  if (!recon_state) continue;
441  ND::THandle<ND::TPIDState> pid_state =
442  (ND::THandle<ND::TPIDState>)recon_state;
443  if (!pid_state) continue;
444 
445  TVector3 node_pos = pid_state->GetPosition().Vect();
446  if (node_counter == 2) outputPath->LastPosition = node_pos;
447  }
448  if (node_counter < 2) outputPath->LastPosition = TVector3(9999, 9999, 9999);
449 
450  ND::THandle<ND::TPIDState> state = path->GetState();
451 
452  outputPath->Length = path->Get<ND::TRealDatum>("Length")->GetValue();
453  outputPath->T0Source = path->Get<ND::TIntegerDatum>("T0Source")->GetValue();
454  outputPath->LikFit = path->CheckStatus(ND::TReconBase::kLikelihoodFit);
455  outputPath->Success = path->CheckStatus(ND::TReconBase::kSuccess);
456  if (outputPath->LikFit) outputPath->FitLikelihood = path->Get<ND::TRealDatum>("FitLogLikelihood")->GetValue();
457  else outputPath->FitLikelihood = -999;
458 
459  outputPath->PID = path->GetParticleId();
460  outputPath->G4ID = TrackTruthInfo::GetG4TrajIDHits(*(path->GetHits()));
461 
462  outputPath->Momentum = path->GetMomentum();
463  if (state->GetMomentumVariance() > 0.) {
464  outputPath->MomentumError = sqrt(state->GetMomentumVariance());
465  };
466  outputPath->Direction = path->GetDirection();
467  outputPath->Charge = path->GetCharge();
468 
469  ND::THandle<ND::THitSelection> hits = path->GetHits();
470  if (hits)
471  outputPath->NClusters = hits->size();
472  else
473  outputPath->NClusters = 0;
474 
475  //(POSSIBLY) TEMPORARY VARIABLES: PARTICLE PULLS
476  outputPath->ElePull = 9999;
477  outputPath->EleExp = 9999;
478  outputPath->EleSigma = 9999;
479  outputPath->MuonPull = 9999;
480  outputPath->MuonExp = 9999;
481  outputPath->MuonSigma = 9999;
482  outputPath->PionPull = 9999;
483  outputPath->PionExp = 9999;
484  outputPath->PionSigma = 9999;
485  outputPath->KaonPull = 9999;
486  outputPath->KaonExp = 9999;
487  outputPath->PionSigma = 9999;
488  outputPath->ProtonPull = 9999;
489  outputPath->ProtonExp = 9999;
490  outputPath->ProtonSigma = 9999;
491  outputPath->PDist = 999999.;
492  outputPath->PEField = 999999.;
493  outputPath->T0Range[0] = 999999.;
494  outputPath->T0Range[1] = 999999.;
495  outputPath->T0RangeDeltaX[0] = 999999.;
496  outputPath->T0RangeDeltaX[1] = 999999.;
497 
498  if (path->Get<ND::TRealDatum>("tpcPid_PullEle")) outputPath->ElePull = (Double_t)path->Get<ND::TRealDatum>("tpcPid_PullEle")->GetValue();
499  if (path->Get<ND::TRealDatum>("tpcPid_dEdxexpEle")) outputPath->EleExp = (Double_t)path->Get<ND::TRealDatum>("tpcPid_dEdxexpEle")->GetValue();
500  if (path->Get<ND::TRealDatum>("tpcPid_SigmaEle")) outputPath->EleSigma = (Double_t)path->Get<ND::TRealDatum>("tpcPid_SigmaEle")->GetValue();
501 
502  if (path->Get<ND::TRealDatum>("tpcPid_PullMuon")) outputPath->MuonPull = (Double_t)path->Get<ND::TRealDatum>("tpcPid_PullMuon")->GetValue();
503  if (path->Get<ND::TRealDatum>("tpcPid_dEdxexpMuon")) outputPath->MuonExp = (Double_t)path->Get<ND::TRealDatum>("tpcPid_dEdxexpMuon")->GetValue();
504  if (path->Get<ND::TRealDatum>("tpcPid_SigmaMuon")) outputPath->MuonSigma = (Double_t)path->Get<ND::TRealDatum>("tpcPid_SigmaMuon")->GetValue();
505 
506  if (path->Get<ND::TRealDatum>("tpcPid_PullPion")) outputPath->PionPull = (Double_t)path->Get<ND::TRealDatum>("tpcPid_PullPion")->GetValue();
507  if (path->Get<ND::TRealDatum>("tpcPid_dEdxexpPion")) outputPath->PionExp = (Double_t)path->Get<ND::TRealDatum>("tpcPid_dEdxexpPion")->GetValue();
508  if (path->Get<ND::TRealDatum>("tpcPid_SigmaPion")) outputPath->PionSigma = (Double_t)path->Get<ND::TRealDatum>("tpcPid_SigmaPion")->GetValue();
509 
510  if (path->Get<ND::TRealDatum>("tpcPid_PullKaon")) outputPath->KaonPull = (Double_t)path->Get<ND::TRealDatum>("tpcPid_PullKaon")->GetValue();
511  if (path->Get<ND::TRealDatum>("tpcPid_dEdxexpKaon")) outputPath->KaonExp = (Double_t)path->Get<ND::TRealDatum>("tpcPid_dEdxexpKaon")->GetValue();
512  if (path->Get<ND::TRealDatum>("tpcPid_SigmaKaon")) outputPath->KaonSigma = (Double_t)path->Get<ND::TRealDatum>("tpcPid_SigmaKaon")->GetValue();
513 
514  if (path->Get<ND::TRealDatum>("tpcPid_PullProton")) outputPath->ProtonPull = (Double_t)path->Get<ND::TRealDatum>("tpcPid_PullProton")->GetValue();
515  if (path->Get<ND::TRealDatum>("tpcPid_dEdxexpProton")) outputPath->ProtonExp = (Double_t)path->Get<ND::TRealDatum>("tpcPid_dEdxexpProton")->GetValue();
516  if (path->Get<ND::TRealDatum>("tpcPid_SigmaProton")) outputPath->ProtonSigma = (Double_t)path->Get<ND::TRealDatum>("tpcPid_SigmaProton")->GetValue();
517 
518  if (path->Get<ND::TRealDatum>("RefitDistortionCorrMom")) outputPath->PDist = (Double_t)path->Get<ND::TRealDatum>("RefitDistortionCorrMom")->GetValue();
519  if (path->Get<ND::TRealDatum>("RefitEFieldDistortionMom")) outputPath->PEField = (Double_t)path->Get<ND::TRealDatum>("RefitEFieldDistortionMom")->GetValue();
520 
521  if (path->Get<ND::TRealDatum>("T0Range")){
522  outputPath->T0Range[0] = (Double_t) path->Get<ND::TRealDatum>("T0Range")->at(0);
523  outputPath->T0Range[1] = (Double_t) path->Get<ND::TRealDatum>("T0Range")->at(1);
524  };
525  if (path->Get<ND::TRealDatum>("T0Range_DeltaX")){
526  outputPath->T0RangeDeltaX[0] = (Double_t) path->Get<ND::TRealDatum>("T0Range_DeltaX")->at(0);
527  outputPath->T0RangeDeltaX[1] = (Double_t) path->Get<ND::TRealDatum>("T0Range_DeltaX")->at(1);
528  };
529 
530  // PLACEHOLDERS
531  outputPath->IsContained = false;
532  outputPath->NetCharge = 0;
533 
534  // PATH MATCHING
535  FillPathMatchingInfo(outputPath, path, patternpaths, patterns,
536  parent_pattern);
537 }
538 
540  TTPCAnaPath* outputPath, ND::THandle<ND::TReconTrack> path,
541  std::list<ND::THandle<ND::TReconBase> > patternpaths,
542  ND::THandle<ND::TReconObjectContainer> patterns,
543  ND::THandle<ND::TReconBase> parent_pattern) {
544 
545  outputPath->PathMatchingID = path->Get<ND::TIntegerDatum>("PathId")->GetValue();
546  outputPath->FirstPosition = path->GetPosition().Vect();
547  outputPath->Time = path->GetPosition().T();
548  ND::TReconNodeContainer pathnodes = path->GetNodes();
549  int node_counter = 0;
550  for (ND::TReconNodeContainer::iterator nodeIt = pathnodes.begin();
551  nodeIt != pathnodes.end(); nodeIt++) {
552  node_counter++;
553 
554  ND::THandle<ND::TReconState> recon_state = (*nodeIt)->GetState();
555  if (!recon_state) continue;
556  ND::THandle<ND::TTrackState> track_state =
557  (ND::THandle<ND::TTrackState>)recon_state;
558  if (!track_state) continue;
559 
560  TVector3 node_pos = track_state->GetPosition().Vect();
561  if (node_counter == 2) outputPath->LastPosition = node_pos;
562  }
563  if (node_counter < 2) outputPath->LastPosition = TVector3(9999, 9999, 9999);
564 
565  outputPath->Length = path->Get<ND::TRealDatum>("Length")->GetValue();
566  outputPath->T0Source = path->Get<ND::TIntegerDatum>("T0Source")->GetValue();
567  outputPath->LikFit = path->CheckStatus(ND::TReconBase::kLikelihoodFit);
568  outputPath->Success = path->CheckStatus(ND::TReconBase::kSuccess);
569  if (outputPath->LikFit) outputPath->FitLikelihood = path->Get<ND::TRealDatum>("FitLogLikelihood")->GetValue();
570  else outputPath->FitLikelihood = -999;
571 
572  outputPath->PID = 0;
573  outputPath->G4ID = TrackTruthInfo::GetG4TrajIDHits(*(path->GetHits()));
574 
575  outputPath->Momentum = -9.99;
576  outputPath->MomentumError = 0.;
577  outputPath->Direction = path->GetDirection();
578 
579  ND::THandle<ND::THitSelection> hits = path->GetHits();
580  if (hits)
581  outputPath->NClusters = hits->size();
582  else
583  outputPath->NClusters = 0;
584 
585  outputPath->ElePull = -9999;
586  outputPath->EleExp = -9999;
587  outputPath->EleSigma = -9999;
588  outputPath->MuonPull = -9999;
589  outputPath->MuonExp = -9999;
590  outputPath->MuonSigma = -9999;
591  outputPath->PionPull = -9999;
592  outputPath->PionExp = -9999;
593  outputPath->PionSigma = -9999;
594  outputPath->KaonPull = -9999;
595  outputPath->KaonExp = -9999;
596  outputPath->PionSigma = -9999;
597  outputPath->ProtonPull = -9999;
598  outputPath->ProtonExp = -9999;
599  outputPath->ProtonSigma = -9999;
600 
601  outputPath->PDist = 999999.;
602  outputPath->PEField = 999999.;
603  outputPath->T0Range[0] = 999999.;
604  outputPath->T0Range[1] = 999999.;
605  outputPath->T0RangeDeltaX[0] = 999999.;
606  outputPath->T0RangeDeltaX[1] = 999999.;
607 
608  if (path->Get<ND::TRealDatum>("RefitDistortionCorrMom")) outputPath->PDist = (Double_t)path->Get<ND::TRealDatum>("RefitDistortionCorrMom")->GetValue();
609  if (path->Get<ND::TRealDatum>("RefitEFieldDistortionMom")) outputPath->PEField = (Double_t)path->Get<ND::TRealDatum>("RefitEFieldDistortionMom")->GetValue();
610 
611  if (path->Get<ND::TRealDatum>("T0Range")){
612  outputPath->T0Range[0] = (Double_t) path->Get<ND::TRealDatum>("T0Range")->at(0);
613  outputPath->T0Range[1] = (Double_t) path->Get<ND::TRealDatum>("T0Range")->at(1);
614  };
615  if (path->Get<ND::TRealDatum>("T0Range_DeltaX")){
616  outputPath->T0RangeDeltaX[0] = (Double_t) path->Get<ND::TRealDatum>("T0Range_DeltaX")->at(0);
617  outputPath->T0RangeDeltaX[1] = (Double_t) path->Get<ND::TRealDatum>("T0Range_DeltaX")->at(1);
618  };
619 
620  // PLACEHOLDERS
621  outputPath->IsContained = false;
622  outputPath->NetCharge = 0;
623 
624  // PATH MATCHING
625  FillPathMatchingInfo(outputPath, path, patternpaths, patterns,
626  parent_pattern);
627 }
628 
629 bool ND::TTRExReconModule::FGDHitVeto(int FGD, int nhits_threshold,
630  ND::TND280Event& event, int mode) {
631  // Get FGD limits in Z
632 
633  double fgd_limits[2] = {0, 0};
634  double frac = 0;
635  if (FGD == 1) {
636  fgd_limits[0] = TGeomInfo::Get().FGD().FGD1ActiveMin().Z();
637  fgd_limits[1] = TGeomInfo::Get().FGD().FGD1ActiveMax().Z();
638  frac = 0.5;
639  }
640  if (FGD == 2) {
641  fgd_limits[0] = TGeomInfo::Get().FGD().FGD2ActiveMin().Z();
642  fgd_limits[1] = TGeomInfo::Get().FGD().FGD2ActiveMax().Z();
643  frac = 0.6;
644  }
645  double backside = fgd_limits[1] - frac * (fgd_limits[1] - fgd_limits[0]);
646 
647  // Build hit selection of hits in defined veto volume
648  ND::THitSelection endhits;
649  ND::THandle<ND::THitSelection> fgd_hits = event.GetHitSelection("fgd");
650  if (fgd_hits) {
651  for (ND::THitSelection::iterator hitIt = fgd_hits->begin();
652  hitIt != fgd_hits->end(); hitIt++) {
653  ND::THandle<ND::THit> hit = *hitIt;
654  if (hit) {
655  if ((hit->GetPosition().Z() > backside) &&
656  (hit->GetPosition().Z() < fgd_limits[1]))
657  endhits.push_back(hit);
658  }
659  }
660  }
661 
662  // Order hit selection in time
663  std::sort(endhits.begin(), endhits.end(), timeLessThan());
664 
665  // Examine hit selection to see if veto conditions are met
666  if (endhits.size() > 0) {
667  for (ND::THitSelection::iterator targetHit = endhits.begin();
668  targetHit < (endhits.end() - 1); targetHit++) {
669  double target_hit_time = (*targetHit)->GetTime();
670  double total_charge = 0;
671  int hit_count = 0;
672  // Loop over other hits, count charge inside time window of target hit.
673  for (ND::THitSelection::iterator otherHit = targetHit + 1;
674  otherHit < endhits.end(); otherHit++) {
675  double other_hit_time = (*otherHit)->GetTime();
676 
677  if ((other_hit_time - target_hit_time) > 100.0 * unit::ns) break;
678  hit_count = hit_count + 1;
679  total_charge = total_charge + (*otherHit)->GetCharge();
680  }
681  if (mode == 0) // Vetoing on number of hits in veto layers
682  {
683  if (hit_count >= 5) // Veto condition is met
684  return true;
685  }
686  if (mode == 1) // Vetoing on total charge deposited in veto layers
687  {
688  if (total_charge >= 25.0) // Veto condition is met
689  return true;
690  }
691  }
692  }
693 
694  return false;
695 }
696 
698  TTPCAnaPath* outputPath, ND::THandle<ND::TReconBase> path,
699  std::list<ND::THandle<ND::TReconBase> > patternpaths,
700  ND::THandle<ND::TReconObjectContainer> patterns,
701  ND::THandle<ND::TReconBase> parent_pattern) {
702  // Do the intra-pattern matching
703  int temp_NMatchingPaths = 0;
704  outputPath->NMatchingPaths = 0;
705  std::vector<int> temp_MatchingPathIDs;
706  outputPath->MatchingPathIDs = 0;
707  std::vector<double> temp_PathMatchingLikelihood;
708  outputPath->PathMatchingLikelihood = 0;
709 
710  for (std::list<ND::THandle<ND::TReconBase> >::iterator listIt =
711  patternpaths.begin();
712  listIt != patternpaths.end(); listIt++) {
713  ND::THandle<ND::TReconBase> testpath = *listIt;
714  if (!testpath) continue;
715 
716  bool lklhd_fit = testpath->CheckStatus(ND::TReconBase::kLikelihoodFit);
717  if (!lklhd_fit) continue;
718 
719  bool fit_success = testpath->CheckStatus(ND::TReconBase::kSuccess);
720  if (!fit_success) continue;
721 
722  if (!(testpath->Get<ND::TIntegerDatum>("PathId"))) continue;
723  int test_PathMatchingID =
724  testpath->Get<ND::TIntegerDatum>("PathId")->GetValue();
725  if (!(testpath->Get<ND::TIntegerDatum>("PathIdMatch"))) continue;
726  if (!(testpath->Get<ND::TRealDatum>("PathMatchLklhd"))) continue;
727 
728  int test_nmatchingpaths =
729  (int)testpath->Get<ND::TIntegerDatum>("PathIdMatch")->size();
730  for (int j = 0; j < test_nmatchingpaths; j++) {
731  int test_matchingpathid =
732  (int)testpath->Get<ND::TIntegerDatum>("PathIdMatch")->at(j);
733  if (test_matchingpathid == outputPath->PathMatchingID) {
734  temp_NMatchingPaths++;
735  temp_MatchingPathIDs.push_back(test_PathMatchingID);
736  temp_PathMatchingLikelihood.push_back(
737  (double)testpath->Get<ND::TRealDatum>("PathMatchLklhd")->at(j));
738  break;
739  }
740  }
741  }
742 
743  if (temp_NMatchingPaths) {
744  outputPath->NMatchingPaths = temp_NMatchingPaths;
745  int n_matching_paths = outputPath->NMatchingPaths;
746  outputPath->MatchingPathIDs = new int[n_matching_paths];
747  outputPath->PathMatchingLikelihood = new double[n_matching_paths];
748  for (int i = 0; i < n_matching_paths; i++) {
749  outputPath->MatchingPathIDs[i] = temp_MatchingPathIDs[i];
750  outputPath->PathMatchingLikelihood[i] = temp_PathMatchingLikelihood[i];
751  }
752  }
753 
754  // Do the inter-pattern matching
755  int temp_NMatchingPatterns = 0;
756  outputPath->NMatchingPatterns = 0;
757  std::vector<int> temp_MatchingPatternIDs;
758  outputPath->MatchingPatternIDs = 0;
759  std::vector<int> temp_MatchingPatternPathIDs;
760  outputPath->MatchingPatternPathIDs = 0;
761  std::vector<double> temp_PatternMatchingLikelihood;
762  outputPath->PatternMatchingLikelihood = 0;
763 
764  if (parent_pattern->Get<ND::TIntegerDatum>("PatternId")) {
765  int parent_pattern_ID =
766  parent_pattern->Get<ND::TIntegerDatum>("PatternId")->GetValue();
767  int path_ID = outputPath->PathMatchingID;
768 
769  for (ND::TReconObjectContainer::iterator pattIt = patterns->begin();
770  pattIt < patterns->end(); pattIt++) {
771  ND::THandle<ND::TReconBase> testpattern = *pattIt;
772  if (!testpattern) continue;
773  if (!(testpattern->Get<ND::TIntegerDatum>("PatternId"))) continue;
774  int testpattern_ID =
775  testpattern->Get<ND::TIntegerDatum>("PatternId")->GetValue();
776  if (testpattern_ID == parent_pattern_ID)
777  continue; // We don't have to look into the same pattern.
778 
779  ND::THandle<ND::TReconObjectContainer> patt_cons =
780  testpattern->GetConstituents();
781  if (patt_cons && patt_cons->size() > 1) { // Pattern with multiple paths
782  for (ND::TReconObjectContainer::iterator conIt = patt_cons->begin();
783  conIt < patt_cons->end(); conIt++) {
784  ND::THandle<ND::TReconBase> testpath = *conIt;
785 
786  if (!testpath) continue;
787  bool lklhd_fit =
788  testpath->CheckStatus(ND::TReconBase::kLikelihoodFit);
789  if (!lklhd_fit) continue;
790  bool fit_success = testpath->CheckStatus(ND::TReconBase::kSuccess);
791  if (!fit_success) continue;
792  if (!(testpath->Get<ND::TIntegerDatum>("PatternIdMatch"))) continue;
793  if (!(testpath->Get<ND::TIntegerDatum>("PatternPathIdMatch")))
794  continue;
795  if (!(testpath->Get<ND::TRealDatum>("PatternMatchLklhd"))) continue;
796  if (!(testpath->Get<ND::TIntegerDatum>("PathId"))) continue;
797  int testpath_ID =
798  testpath->Get<ND::TIntegerDatum>("PathId")->GetValue();
799 
800  int N_matchingpatt_test =
801  testpath->Get<ND::TIntegerDatum>("PatternIdMatch")->size();
802  for (int j = 0; j < N_matchingpatt_test; j++) {
803  int matchingpattID_test =
804  (int)testpath->Get<ND::TIntegerDatum>("PatternIdMatch")->at(j);
805  int matchingpattpathID_test =
806  (int)testpath->Get<ND::TIntegerDatum>("PatternPathIdMatch")
807  ->at(j);
808  if ((matchingpattID_test == parent_pattern_ID) &&
809  (matchingpattpathID_test == path_ID)) {
810  temp_NMatchingPatterns++;
811  temp_MatchingPatternIDs.push_back(testpattern_ID);
812  temp_MatchingPatternPathIDs.push_back(testpath_ID);
813  temp_PatternMatchingLikelihood.push_back(
814  testpath->Get<ND::TRealDatum>("PatternMatchLklhd")->at(j));
815  break;
816  }
817  }
818  }
819  } else { // Single path pattern
820  ND::THandle<ND::TReconBase> testpath = testpattern;
821 
822  if (!testpath) continue;
823  bool lklhd_fit = testpath->CheckStatus(ND::TReconBase::kLikelihoodFit);
824  if (!lklhd_fit) continue;
825  bool fit_success = testpath->CheckStatus(ND::TReconBase::kSuccess);
826  if (!fit_success) continue;
827  if (!(testpath->Get<ND::TIntegerDatum>("PatternIdMatch"))) continue;
828  if (!(testpath->Get<ND::TIntegerDatum>("PatternPathIdMatch"))) continue;
829  if (!(testpath->Get<ND::TRealDatum>("PatternMatchLklhd"))) continue;
830  if (!(testpath->Get<ND::TIntegerDatum>("PathId"))) continue;
831  int testpath_ID =
832  testpath->Get<ND::TIntegerDatum>("PathId")->GetValue();
833 
834  int N_matchingpatt_test =
835  testpath->Get<ND::TIntegerDatum>("PatternIdMatch")->size();
836  for (int j = 0; j < N_matchingpatt_test; j++) {
837  int matchingpattID_test =
838  (int)testpath->Get<ND::TIntegerDatum>("PatternIdMatch")->at(j);
839  int matchingpattpathID_test =
840  (int)testpath->Get<ND::TIntegerDatum>("PatternPathIdMatch")
841  ->at(j);
842  if ((matchingpattID_test == parent_pattern_ID) &&
843  (matchingpattpathID_test == path_ID)) {
844  temp_NMatchingPatterns++;
845  temp_MatchingPatternIDs.push_back(testpattern_ID);
846  temp_MatchingPatternPathIDs.push_back(testpath_ID);
847  temp_PatternMatchingLikelihood.push_back(
848  testpath->Get<ND::TRealDatum>("PatternMatchLklhd")->at(j));
849  break;
850  }
851  }
852  }
853  }
854 
855  if (temp_NMatchingPatterns) {
856  outputPath->NMatchingPatterns = temp_NMatchingPatterns;
857  int n_matching_patterns = outputPath->NMatchingPatterns;
858  outputPath->MatchingPatternIDs = new int[n_matching_patterns];
859  outputPath->MatchingPatternPathIDs = new int[n_matching_patterns];
860  outputPath->PatternMatchingLikelihood = new double[n_matching_patterns];
861  for (int i = 0; i < n_matching_patterns; i++) {
862  outputPath->MatchingPatternIDs[i] = temp_MatchingPatternIDs[i];
863  outputPath->MatchingPatternPathIDs[i] = temp_MatchingPatternPathIDs[i];
864  outputPath->PatternMatchingLikelihood[i] =
865  temp_PatternMatchingLikelihood[i];
866  }
867  }
868  }
869 }
870 
Int_t NMatchingPaths
^ The unique identifiers of the associated junctions
virtual Bool_t ProcessFirstEvent(ND::TND280Event &event)
Is called after the first event is loaded in.
Double_t ElePull
Pull-related values.
Int_t PatternID
The ID of the pattern.
TTPCAnaJunction An object describing a meeting point of two or more paths, and containing all associa...
TTRExReconModule(const char *name="TREx", const char *title="TPC Recon Extension Module")
TTPCAnaPath An object describing a track-like grouping of hits in the TPC.
ClassImp(ND::TBeamSummaryDataModule::TBeamSummaryData)
Int_t PathID
A unique identifier for the path within the pattern.
std::string fDescription
A longish descrition of the analysis.
Double_t * PathMatchingLikelihood
^ The unique identifiers of the paths
bool FGDHitVeto(int FGD, int nhits_threshold, ND::TND280Event &event, int mode)
Double_t FitLikelihood
Value of log likelihood from fit.
Double_t Length
The length of the path.
virtual void InitializeModule()
Initialize Module, override if necessary.
Double_t MomentumError
The momentum error of the path.
TVector3 LastPosition
The end position of the path.
Int_t NMatchingPatterns
^ The likelihood of matching to that path.
TClonesArray * Paths
The constituent paths of the pattern.
Int_t G4ID
G4 ID of the primary true trajectory associated with this path.
Int_t TPC
The TPC in which the pattern resides.
Bool_t S1Sflag
A flag to say whether this pattern would pass the TPC gas Stage 1 Selection.
Int_t JunctionID
A unique identifier for the junction within the pattern.
std::string fCVSID
Defined if an official tagged version.
Double_t Time
The T0 for the path.
Int_t NClusters
Number of clusters in the path.
Bool_t LikFit
True if the likelihood fit was performed.
Int_t NJunctions
The number of junctions in the pattern.
void SetNameTitle(char const *name, char const *title)
void FillPathInfo(TTPCAnaPath *outputPath, ND::THandle< ND::TReconPID > path, std::list< ND::THandle< ND::TReconBase > > patternpaths, ND::THandle< ND::TReconObjectContainer > patterns, ND::THandle< ND::TReconBase > parent_pattern)
Fills a TTPCAnaPath with information from either a TReconTrack or a TReconPID.
Int_t NJunctions
The number of junctions to which the path connects.
Double_t Momentum
The momentum of the path.
TVector3 FirstPosition
The start position of the path.
virtual void InitializeBranches()
Initialize Branches. Don&#39;t do anything else in this function.
Int_t T0Source
Enumerator of T0 sources.
TTPCAnaPattern The master object for any connected system of hits in a TPC.
#define CVSID
virtual bool FillTree(ND::TND280Event &)
Fill all the stuff that goes in the output tree.
TVector3 MaximumCoordinates
^ The unique identifiers of the associated paths.
virtual bool IsFullEventWorthSaving(ND::TND280Event &event)
Returns true if there is at least one recon object.
std::string fCVSTagName
Defined if an official tagged version.
Bool_t IsContained
True if the path does not leave the TPC; false otherwise.
TVector3 MinimumCoordinates
The most extreme hit positions in the negative direction.
TVector3 Direction
The direction of the path.
int GetTPCFromPattern(ND::THandle< ND::TReconVertex >)
Returns TPC in which the pattern was found for all 3 possible input types.
Int_t NPaths
The number of paths connected to the junction.
#define CVSTAG
Double_t NetCharge
The total charge left by the path on the micromegas.
Int_t NPaths
The number of paths.
Double_t Charge
Charge of the path.
TLorentzVector Position
The spatial coordinates of the junction.
TClonesArray * Junctions
The constituent junctions of the pattern.
Int_t PID
The PID of the path.
Bool_t Success
True if the likelihood fit succeeded.
void FillPathMatchingInfo(TTPCAnaPath *outputPath, ND::THandle< ND::TReconBase > path, std::list< ND::THandle< ND::TReconBase > > patternpaths, ND::THandle< ND::TReconObjectContainer > patterns, ND::THandle< ND::TReconBase > parent_pattern)

Package Summary
Package Name: eventAnalysis
Package Version: 7.0-49-g0ac7482
Package Manager:

Generated on Mon Mar 25 2024 14:43:59 for eventAnalysis by doxygen 1.8.5