eventAnalysis  7.0-49-g0ac7482
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TTrackerReconModule.cxx
Go to the documentation of this file.
1 #include "TTrackerReconModule.hxx"
2 
3 #include <algorithm>
4 #include <functional>
5 
6 #include <TDatum.hxx>
7 #include <TIntegerDatum.hxx>
8 #include <TRealDatum.hxx>
9 #include <TReconNode.hxx>
10 #include <TReconShower.hxx>
11 #include <TReconTrack.hxx>
12 
13 #include <TComboHit.hxx>
14 #include <TG4HitSegment.hxx>
15 #include <TGeomInfo.hxx>
16 #include <TMultiHit.hxx>
17 #include <retrieveHitTruthInfo.hxx>
18 
19 #include <TOARuntimeParams.hxx>
20 
21 #include "ReconObjectUtils.hxx"
22 #include "TEventFolder.hxx"
23 #include "TRecPackManager.hxx"
24 #include "TValidationUtils.hxx"
25 
26 #include <TFieldManager.hxx>
27 #include <TrackTruthInfo.hxx>
28 
29 ClassImp(ND::TTrackerReconModule::TUnusedHit);
30 ND::TTrackerReconModule::TUnusedHit::~TUnusedHit() {}
31 
32 ClassImp(ND::TTrackerReconModule::TTrackerVertex);
33 ND::TTrackerReconModule::TTrackerVertex::~TTrackerVertex() {}
34 
35 ClassImp(ND::TTrackerReconModule::TTPCTrack);
36 ND::TTrackerReconModule::TTPCTrack::~TTPCTrack() {}
37 
38 ClassImp(ND::TTrackerReconModule::TFGDTrack);
39 ND::TTrackerReconModule::TFGDTrack::~TFGDTrack() {}
40 
41 ClassImp(ND::TTrackerReconModule::TTrueParticle);
42 ND::TTrackerReconModule::TTrueParticle::~TTrueParticle() {}
43 
44 ClassImp(ND::TTrackerReconModule::TTrackerNode);
45 ND::TTrackerReconModule::TTrackerNode::~TTrackerNode() {}
46 
47 ClassImp(ND::TTrackerReconModule::TTrackerConstituent);
48 ND::TTrackerReconModule::TTrackerConstituent::~TTrackerConstituent() {}
49 ND::TTrackerReconModule::TTrackerConstituent::TTrackerConstituent() {
50  NConstituents = 0;
51  return;
52 }
53 
54 ClassImp(ND::TTrackerReconModule::TTrackerResult);
55 ND::TTrackerReconModule::TTrackerResult::~TTrackerResult() {}
56 ND::TTrackerReconModule::TTrackerResult::TTrackerResult() {
57  Constituents =
58  new TClonesArray("ND::TTrackerReconModule::TTrackerConstituent", 10);
59  TPC = new TClonesArray("ND::TTrackerReconModule::TTPCTrack", 10);
60  FGD = new TClonesArray("ND::TTrackerReconModule::TFGDTrack", 10);
61  Nodes = new TClonesArray("ND::TTrackerReconModule::TTrackerNode", 10);
62 
63  NConstituents = 0;
64  NTotalConstituents = 0;
65  Constituents->Clear();
66 
67  NTPCs = 0;
68  TPC->Clear();
69 
70  NFGDs = 0;
71  FGD->Clear();
72 
73  NNodes = 0;
74  Nodes->Clear();
75 
76  return;
77 }
78 
79 ClassImp(ND::TTrackerReconModule::TTrackOther);
80 ND::TTrackerReconModule::TTrackOther::~TTrackOther() {}
81 ND::TTrackerReconModule::TTrackOther::TTrackOther() {
82  NHits = 0;
83  Hits = new TClonesArray("ND::TTrackerReconModule::TUnusedHit", 200);
84  return;
85 }
86 
87 #define CVSTAG \
88  "\
89  $Name$"
90 #define CVSID \
91  "\
92 "
93 
94 //*************************************************************
95 ND::TTrackerReconModule::TTrackerReconModule(const char* name,
96  const char* title) {
97  //*************************************************************
98  SetNameTitle(name, title);
99  // Enable this module by default:
100  fIsEnabled = kTRUE;
101  fDescription = "Tracker Recon Output";
102  fCVSTagName = CVSTAG;
103  fCVSID = CVSID;
104  // Tree variables
105  fNVertices = 0;
106  fVertices = new TClonesArray("ND::TTrackerReconModule::TTrackerVertex", 200);
107  fNTracks = 0;
108  fTracks = new TClonesArray("ND::TTrackerReconModule::TTrackerResult", 200);
109  fNFGDOther = 0;
110  fFGDOther = new TClonesArray("ND::TTrackerReconModule::TTrackOther", 200);
111  fNTPCOther = 0;
112  fTPCOther = new TClonesArray("ND::TTrackerReconModule::TTrackOther", 200);
113  fNTPCIso = 0;
114  fTPCIso = new TClonesArray("ND::TTrackerReconModule::TTPCTrack", 200);
115  fNTPCUnused = 0;
116  fTPCUnused = new TClonesArray("ND::TTrackerReconModule::TUnusedHit", 200);
117  fNFGDUnused = 0;
118  fFGDUnused = new TClonesArray("ND::TTrackerReconModule::TUnusedHit", 200);
119  fNTPCExtrapolation = 0;
120  fTPCExtrapolation =
121  new TClonesArray("ND::TTrackerReconModule::TTPCTrack", 200);
122 }
123 
124 //*************************************************************
125 ND::TTrackerReconModule::~TTrackerReconModule() {}
126 //*************************************************************
127 
128 //*************************************************************
129 Bool_t ND::TTrackerReconModule::ProcessFirstEvent(ND::TND280Event& event) {
130  //*************************************************************
131 
132  fMaxDrift = ND::TGeomInfo::Get().TPC().GetMaxDriftDistance(0, 0);
133  const TGeoNode* cath = ND::TGeomInfo::Get().TPC().GetCathode(0);
134  const TGeoVolume* cathvol = cath->GetVolume();
135  TGeoShape* cathshape = cathvol->GetShape();
136  TGeoBBox* box = dynamic_cast<TGeoBBox*>(cathshape);
137  fCathodeOffset = box->GetDX();
138 
139  return true;
140 }
141 
142 //*************************************************************
143 void ND::TTrackerReconModule::InitializeModule() {
144  fDriftVelocity = ND::TOARuntimeParams::Get().GetParameterD(
145  "tpcRecon.TPCgas.DriftSpeed");
146 }
147 //*************************************************************
148 
149 //*************************************************************
150 void ND::TTrackerReconModule::InitializeBranches() {
151  //*************************************************************
152  SetSplitLevel(1);
153 
154  // The number of objects that have been saved.
155  fOutputTree->Branch("NVertices", &fNVertices, "NVertices/I", fBufferSize)
156  ->SetTitle("The number of added vertices ");
157  fOutputTree->Branch("NTracks", &fNTracks, "NTracks/I", fBufferSize)
158  ->SetTitle("The number of trackerRecon results ");
159  fOutputTree->Branch("NFGDOther", &fNFGDOther, "NFGDOther/I", fBufferSize)
160  ->SetTitle(" The number of FGD tracks with no fit ");
161  fOutputTree->Branch("NTPCOther", &fNTPCOther, "NTPCOther/I", fBufferSize)
162  ->SetTitle(" The number of TPC tracks with no fit ");
163  fOutputTree->Branch("NTPCIso", &fNTPCIso, "NTPCIso/I", fBufferSize)
164  ->SetTitle("The number of isolated TPC tracks with fits ");
165  fOutputTree->Branch("NFGDUnused", &fNFGDUnused, "NFGDUnused/I", fBufferSize)
166  ->SetTitle(" The number of unused FGD hits");
167  fOutputTree->Branch("NTPCUnused", &fNTPCUnused, "NTPCUnused/I", fBufferSize)
168  ->SetTitle(" The number of unused TPC hits ");
169  fOutputTree->Branch("NTPCExtrapolation", &fNTPCExtrapolation,
170  "NTPCExtrapolation/I", fBufferSize)
171  ->SetTitle(" ");
172 
173  // A branch for the clone array of Tracker object information.
174  fOutputTree->Branch("Vertices", &fVertices, fBufferSize, fSplitLevel)
175  ->SetTitle(" The vector of trackerRecon vertices");
176  fOutputTree->Branch("Tracks", &fTracks, fBufferSize, fSplitLevel)
177  ->SetTitle("The vector of overall trackerRecon results ");
178  fOutputTree->Branch("FGDOther", &fFGDOther, fBufferSize, fSplitLevel)
179  ->SetTitle("The vector of FGD tracks with no fit ");
180  fOutputTree->Branch("TPCOther", &fTPCOther, fBufferSize, fSplitLevel)
181  ->SetTitle("The vector of TPC tracks with no fit ");
182  fOutputTree->Branch("TPCIso", &fTPCIso, fBufferSize, fSplitLevel)
183  ->SetTitle(" The vector of isolated TPC tracks with fits");
184  fOutputTree->Branch("TPCUnused", &fTPCUnused, fBufferSize, fSplitLevel)
185  ->SetTitle(" The vector of unused TPC hits");
186  fOutputTree->Branch("FGDUnused", &fFGDUnused, fBufferSize, fSplitLevel)
187  ->SetTitle(" The number of unused FGD hits");
188  fOutputTree->Branch("TPCExtrapolation", &fTPCExtrapolation, fBufferSize,
189  fSplitLevel)
190  ->SetTitle("The vector of TPC tracks extrapolatedinto the FGD following Claudio's (2010a)method ");
191  return;
192 }
193 
194 void ND::TTrackerReconModule::FillConfigTree(TTree* configTree) {
195  if (fFilledConfigTree) {
196  ND280NamedWarn(
197  "TTrackerReconModule",
198  "Module: " << this->GetName()
199  << " has already written the config tree this run.");
200  }
201 
202  configTree->Branch("MaxDrift", &fMaxDrift, "MaxDrift/D");
203  configTree->Branch("CathodeOffset", &fCathodeOffset, "CathodeOffset/D");
204  configTree->Branch("DriftVelocity", &fDriftVelocity, "DriftVelocity/D");
205  TAnalysisModuleBase::FillConfigTree(configTree);
206 }
207 
208 //*************************************************************
209 void ND::TTrackerReconModule::FillVertex(ND::THandle<ND::TReconVertex> in,
210  TClonesArray* out, int idx) {
211  //*************************************************************
212 
213  TTrackerVertex* vertex = new ((*out)[idx]) TTrackerVertex;
214 
215  vertex->AlgorithmName = in->GetAlgorithmName();
216  vertex->Status = in->CheckStatus(in->kSuccess);
217  vertex->Quality = in->GetQuality();
218  vertex->NDOF = in->GetNDOF();
219 
220  ND::THandle<ND::TVertexState> state = in->GetState();
221  if (state) {
222  vertex->Position = state->GetPosition();
223  vertex->Variance = state->GetPositionVariance();
224  } else {
225  vertex->Position = TLorentzVector(0., 0., 0., 0.);
226  vertex->Variance = TLorentzVector(0., 0., 0., 0.);
227  }
228  return;
229 }
230 
231 //*************************************************************
232 void ND::TTrackerReconModule::FillFGD(ND::THandle<ND::TReconTrack> in,
233  TClonesArray* out, int idx) {
234  //*************************************************************
235 
236  ND::TReconBase::Status dets[2] = {ND::TReconBase::kFGD1,
237  ND::TReconBase::kFGD2};
238 
239  TFGDTrack* fgd = new ((*out)[idx]) TFGDTrack;
240 
241  for (int i = 0; i < 2; i++) {
242  if (in->UsesDetector(dets[i])) fgd->Detector = i + 1;
243  }
244  fgd->UniqueID = in->GetUniqueID();
245  fgd->Ndof = in->GetNDOF();
246  fgd->Chi2 = in->GetQuality();
247  ND::THandle<ND::TTrackState> state = in->GetState();
248  if (state) {
249  fgd->Position = state->GetPosition();
250  fgd->Direction = state->GetDirection();
251  fgd->EDeposit = state->GetEDeposit();
252  } else {
253  // set bogus values
254  fgd->Position = TLorentzVector(0., 0., 0., 0.);
255  fgd->Direction = TVector3(0., 0., 0.);
256  fgd->EDeposit = 0.;
257  }
258 
259  return;
260 }
261 
262 //*************************************************************
263 void ND::TTrackerReconModule::FillTPC(ND::THandle<ND::TReconPID> in,
264  TClonesArray* out, int idx) {
265  //*************************************************************
266 
267  ND::TReconBase::Status dets[3] = {
268  ND::TReconBase::kTPC1, ND::TReconBase::kTPC2, ND::TReconBase::kTPC3};
269 
270  TTPCTrack* tpc = new ((*out)[idx]) TTPCTrack;
271 
272  for (int i = 0; i < 3; i++) {
273  if (in->UsesDetector(dets[i])) tpc->Detector = i + 1;
274  }
275  tpc->Ndof = in->GetNDOF();
276  tpc->Chi2 = in->GetQuality();
277  tpc->NHits = 0;
278  tpc->UniqueID = in->GetUniqueID();
279  if (in->GetHits()) tpc->NHits = GetNumberOfHits(in);
280 
281  ND::THandle<ND::TPIDState> state = in->GetState();
282  if (state) {
283  tpc->Momentum = state->GetMomentum();
284  tpc->MomentumError = sqrt(fabs(state->GetMomentumVariance()));
285  tpc->Charge = state->GetCharge();
286  tpc->Position = state->GetPosition();
287  tpc->Direction = state->GetDirection();
288  tpc->DirectionVariance = state->GetDirectionVariance();
289  } else {
290  // set bogus values
291  tpc->Momentum = 0.;
292  tpc->MomentumError = 0.;
293  tpc->Charge = 0.;
294  tpc->Position = TLorentzVector(0., 0., 0., 0.);
295  tpc->Direction = TVector3(0., 0., 0.);
296  tpc->PositionVariance = TLorentzVector(0., 0., 0., 0.);
297  tpc->DirectionVariance = TVector3(0., 0., 0.);
298  }
299 
300  if (in->Get<ND::TRealDatum>("tpcPid_NTrun"))
301  tpc->NTrun = in->Get<ND::TRealDatum>("tpcPid_NTrun")->GetValue();
302  else
303  tpc->NTrun = -0xABCDEF;
304 
305  if (in->Get<ND::TRealDatum>("tpcPid_Ccorr"))
306  tpc->Ccorr = in->Get<ND::TRealDatum>("tpcPid_Ccorr")->GetValue();
307  else
308  tpc->Ccorr = -0xABCDEF;
309 
310  if (in->Get<ND::TRealDatum>("tpcPid_PullEle"))
311  tpc->PullEle = in->Get<ND::TRealDatum>("tpcPid_PullEle")->GetValue();
312  else
313  tpc->PullEle = -0xABCDEF;
314 
315  if (in->Get<ND::TRealDatum>("tpcPid_PullMuon"))
316  tpc->PullMuon = in->Get<ND::TRealDatum>("tpcPid_PullMuon")->GetValue();
317  else
318  tpc->PullMuon = -0xABCDEF;
319 
320  if (in->Get<ND::TRealDatum>("tpcPid_PullPion"))
321  tpc->PullPion = in->Get<ND::TRealDatum>("tpcPid_PullPion")->GetValue();
322  else
323  tpc->PullPion = -0xABCDEF;
324 
325  if (in->Get<ND::TRealDatum>("tpcPid_PullKaon"))
326  tpc->PullKaon = in->Get<ND::TRealDatum>("tpcPid_PullKaon")->GetValue();
327  else
328  tpc->PullKaon = -0xABCDEF;
329 
330  if (in->Get<ND::TRealDatum>("tpcPid_PullProton"))
331  tpc->PullProton = in->Get<ND::TRealDatum>("tpcPid_PullProton")->GetValue();
332  else
333  tpc->PullProton = -0xABCDEF;
334 
335  if (in->Get<ND::TRealDatum>("tpcPid_dEdxexpEle"))
336  tpc->dEdxexpEle = in->Get<ND::TRealDatum>("tpcPid_dEdxexpEle")->GetValue();
337  else
338  tpc->dEdxexpEle = -0xABCDEF;
339 
340  if (in->Get<ND::TRealDatum>("tpcPid_dEdxexpMuon"))
341  tpc->dEdxexpMuon =
342  in->Get<ND::TRealDatum>("tpcPid_dEdxexpMuon")->GetValue();
343  else
344  tpc->dEdxexpMuon = -0xABCDEF;
345 
346  if (in->Get<ND::TRealDatum>("tpcPid_dEdxexpPion"))
347  tpc->dEdxexpPion =
348  in->Get<ND::TRealDatum>("tpcPid_dEdxexpPion")->GetValue();
349  else
350  tpc->dEdxexpPion = -0xABCDEF;
351 
352  if (in->Get<ND::TRealDatum>("tpcPid_dEdxexpKaon"))
353  tpc->dEdxexpKaon =
354  in->Get<ND::TRealDatum>("tpcPid_dEdxexpKaon")->GetValue();
355  else
356  tpc->dEdxexpKaon = -0xABCDEF;
357 
358  if (in->Get<ND::TRealDatum>("tpcPid_dEdxexpProton"))
359  tpc->dEdxexpProton =
360  in->Get<ND::TRealDatum>("tpcPid_dEdxexpProton")->GetValue();
361  else
362  tpc->dEdxexpProton = -0xABCDEF;
363 
364  if (in->Get<ND::TRealDatum>("tpcPid_SigmaEle"))
365  tpc->SigmaEle = in->Get<ND::TRealDatum>("tpcPid_SigmaEle")->GetValue();
366  else
367  tpc->SigmaEle = -0xABCDEF;
368 
369  if (in->Get<ND::TRealDatum>("tpcPid_SigmaMuon"))
370  tpc->SigmaMuon = in->Get<ND::TRealDatum>("tpcPid_SigmaMuon")->GetValue();
371  else
372  tpc->SigmaMuon = -0xABCDEF;
373 
374  if (in->Get<ND::TRealDatum>("tpcPid_SigmaPion"))
375  tpc->SigmaPion = in->Get<ND::TRealDatum>("tpcPid_SigmaPion")->GetValue();
376  else
377  tpc->SigmaPion = -0xABCDEF;
378 
379  if (in->Get<ND::TRealDatum>("tpcPid_SigmaKaon"))
380  tpc->SigmaKaon = in->Get<ND::TRealDatum>("tpcPid_SigmaKaon")->GetValue();
381  else
382  tpc->SigmaKaon = -0xABCDEF;
383 
384  if (in->Get<ND::TRealDatum>("tpcPid_SigmaProton"))
385  tpc->SigmaProton =
386  in->Get<ND::TRealDatum>("tpcPid_SigmaProton")->GetValue();
387  else
388  tpc->SigmaProton = -0xABCDEF;
389 
390  // now look for constituent track that made this PID
391  tpc->NConstituents = 0;
392  tpc->Sigma0 = 0;
393  tpc->Sigma1 = 0;
394  tpc->Sigma2 = 0;
395  tpc->MeanDrift = 0;
396  tpc->TrDirection = TVector3(0, 0, 0);
397  tpc->TrDirectionVar = TVector3(0, 0, 0);
398  tpc->TrCurvature = 0.0;
399  tpc->TrCurvatureVar = 0.0;
400 
401  // check if this is a composite object... if so fill constituents
402  if (in->GetConstituents()) {
403  ND::THandle<ND::TReconTrack> track =
404  GetConstituentTrack(in, tpc->NConstituents);
405  if (track) {
406  // track->ls();
407  ND::THandle<ND::TTrackState> astate = track->GetState();
408  if (astate) {
409  tpc->TrDirection = astate->GetDirection();
410  tpc->TrDirectionVar = TVector3(
411  astate->GetCovarianceValue(astate->GetDirectionIndex(),
412  astate->GetDirectionIndex()),
413  astate->GetCovarianceValue(astate->GetDirectionIndex() + 1,
414  astate->GetDirectionIndex() + 1),
415  astate->GetCovarianceValue(astate->GetDirectionIndex() + 2,
416  astate->GetDirectionIndex() + 2));
417 
418  tpc->TrCurvature = astate->GetCurvature();
419  tpc->TrCurvatureVar = astate->GetCovarianceValue(
420  astate->GetCurvatureIndex(), astate->GetCurvatureIndex());
421  }
422  if (track->Get<ND::TRealDatum>("Sigma0"))
423  tpc->Sigma0 = track->Get<ND::TRealDatum>("Sigma0")->GetValue();
424  if (track->Get<ND::TRealDatum>("Sigma1"))
425  tpc->Sigma1 = track->Get<ND::TRealDatum>("Sigma1")->GetValue();
426  if (track->Get<ND::TRealDatum>("Sigma2"))
427  tpc->Sigma2 = track->Get<ND::TRealDatum>("Sigma2")->GetValue();
428  if (track->Get<ND::TRealDatum>("MeanDrift"))
429  tpc->MeanDrift = track->Get<ND::TRealDatum>("MeanDrift")->GetValue();
430  }
431  }
432  tpc->HasExtrapolation = false;
433 }
434 
435 // Utility functions from Claudio's code
436 void getTangents(ND::THandle<ND::TTrackState> state, double& tx, double& etx,
437  double& ty, double& ety) {
438  tx = state->GetDirection()[0] / state->GetDirection()[2];
439  ty = state->GetDirection()[1] / state->GetDirection()[2];
440 
441  CLHEP::HepMatrix Jacobian(2, 3, 0);
442  CLHEP::HepMatrix Covariance(3, 3, 0);
443  CLHEP::HepMatrix Cov(2, 2, 0);
444 
445  for (int i = 0; i < 3; i++)
446  for (int j = 0; j < 3; j++) {
447  Covariance[i][j] = state->GetCovarianceValue(
448  state->GetDirectionIndex() + i, state->GetDirectionIndex() + j);
449  }
450 
451  Jacobian[0][0] = 1. / state->GetDirection()[2];
452  Jacobian[0][1] = 0.;
453  Jacobian[0][2] = -tx / state->GetDirection()[2];
454 
455  Jacobian[1][0] = 0.;
456  Jacobian[1][1] = 1. / state->GetDirection()[2];
457  Jacobian[1][2] = -ty / state->GetDirection()[2];
458 
459  Cov = Jacobian * Covariance * Jacobian.T();
460 
461  etx = TMath::Sqrt(Cov[0][0]);
462  ety = TMath::Sqrt(Cov[1][1]);
463 }
464 
465 //=========================================================//
466 bool insideFGD1(double zx, double zy, double xx, double yy) {
467  return ((zx > 140) && (zy > 140) && (zx < 500) && (zy < 500) &&
468  (fabs(yy) < 800) && (fabs(xx) < 800));
469 }
470 //=========================================================//
471 bool insideFGD2(double zx, double zy, double xx, double yy) {
472  return ((zx > 1530) && (zy > 1530) && (zx < 1900) && (zy < 1900) &&
473  (fabs(yy) < 800) && (fabs(xx) < 800));
474 }
475 //=========================================================//
476 bool borderFGD1(TVector3 pos) {
477  bool ret = true;
478  if ((pos.Y() < 850) && (pos.Y() > -800) && (pos.X() < 850)) ret = false;
479  if (pos.Z() > 500) ret = false;
480  if (fabs(pos.Z() - 120.85) < 0.0001) {
481  ret = true;
482  }
483  return ret;
484 }
485 //=========================================================//
486 bool borderFGD2(TVector3 pos) {
487  bool ret = true;
488  if ((pos.Y() < 850) && (pos.Y() > -800) && (fabs(pos.X()) < 850)) {
489  ret = false;
490  }
491  if (pos.Z() < 500) ret = false;
492  if (fabs(pos.Z() - 1478.85) < 0.0001) ret = true;
493  return ret;
494 }
495 
496 // This is the distance computation function used in the extrapolation analysis
497 double computeDistFGD(TVector3 f, TVector3 in, double tx, double ty,
498  double curv, bool isX) {
499  // reject hits far away, tracks with zero curvature, and hits downstream from
500  // the track
501  if ((fabs(curv) < 1.e-20) || (fabs(in.Z() - f.Z()) > 600) || (f.Z() > in.Z()))
502  return -9999.0;
503 
504  double r1 = fabs(1. / curv);
505  double phi1 = atan(ty);
506  double charge1 = 1.;
507  if (curv > 0) charge1 = -1.;
508  double zc1 = in.Z() - charge1 * r1 * sin(phi1);
509  double yc1 = in.Y() + charge1 * r1 * cos(phi1);
510  double delta = pow(r1, 2) - pow((zc1 - f.Z()), 2);
511  if (delta < 0) return -9999.0;
512  double ySol1 = yc1 + sqrt(delta);
513  double ySol2 = yc1 - sqrt(delta);
514  double ySol = ySol1;
515  if (fabs(ySol1 - in.Y()) > fabs(ySol2 - in.Y())) ySol = ySol2;
516  // check with linear extrapo
517 
518  // yInts currentl not being used, so commented out to remove compiler warnings
519  // double yIntS = in.Y()+ty*(f.Z()- in.Z());
520 
521  if (!isX) {
522  return f.Y() - ySol;
523  } else {
524  return f.X() - in.X() - tx * (f.Z() - in.Z());
525  }
526 }
527 
528 //*************************************************************
529 bool ND::TTrackerReconModule::ExtrapolateTPC(
530  ND::THandle<ND::TReconPID> in, TClonesArray* out, int idx,
531  ND::THandle<ND::THitSelection> fgd) {
532  ND::THandle<ND::TReconTrack> track = in->GetConstituents()->front();
533  if (!track) return false;
534  ND::TReconNodeContainer& nodes = track->GetNodes();
535 
536  // Only do the extrapolation matching if there are enough TPC points (36 for
537  // now)
538  if (nodes.size() < 36) return false;
539  // ignore TPC1 tracks
540  if (in->UsesDetector(ND::TReconBase::kTPC1)) return false;
541 
542  ND::THandle<ND::TTrackState> frontstate = nodes.front()->GetState();
543  if (!frontstate) return false;
544  if (!fgd) return false;
545 
546  double txin, tyin, etxin, etyin;
547  getTangents(frontstate, txin, etxin, tyin, etyin);
548  TVector3 pos = frontstate->GetPosition().Vect();
549  double curvature = track->GetCurvature();
550 
551  // start with bogus numbers
552  double ExtrapolatedVertexXX = 9999;
553  double ExtrapolatedVertexZY = 9999;
554  double ExtrapolatedVertexZX = 9999;
555  double ExtrapolatedVertexYY = 9999;
556  // std::vector<TVector3*> ExtrapolatedHits;
557 
558  // calculate the vertex and find the set of associated FGD hits
559  // the distance calculation function will reject hits from the wrong fgd
560  // by discarding hits downstream of the track origin and more than 60 cm away
561  // in Z
562  for (ND::THitSelection::iterator hit = fgd->begin(); hit != fgd->end();
563  hit++) {
564  TVector3 fgdPos = (*hit)->GetPosition();
565 
566  double dist =
567  computeDistFGD(fgdPos, pos, txin, tyin, curvature, (*hit)->IsXHit());
568  if (fabs(dist) >= 30) // we reject this hit (30 mm is the threshold)
569  continue;
570 
571  // ExtrapolatedHits.push_back(new TVector3(fgdPos));
572 
573  if ((*hit)->IsYHit() && fgdPos.Z() < ExtrapolatedVertexZY) {
574  ExtrapolatedVertexZY = fgdPos.Z();
575  ExtrapolatedVertexYY = fgdPos.Y();
576  } // Y coord
577 
578  if ((*hit)->IsXHit() && fgdPos.Z() < ExtrapolatedVertexZX) {
579  ExtrapolatedVertexZX = fgdPos.Z();
580  ExtrapolatedVertexXX = fgdPos.X();
581  } // X coord
582  } // loop on FGD hits
583 
584  if (ExtrapolatedVertexZX == 9999 || ExtrapolatedVertexZY == 9999)
585  return false;
586 
587  bool EnteringFGD = false;
588  // Now check whether this qualifies as an 'entering' track
589  if (insideFGD1(ExtrapolatedVertexZX, ExtrapolatedVertexZY,
590  ExtrapolatedVertexXX, ExtrapolatedVertexYY)) {
591  for (ND::THitSelection::iterator hit = fgd->begin(); hit != fgd->end();
592  hit++) {
593  TVector3 fgdPos = (*hit)->GetPosition();
594 
595  if (!borderFGD1(fgdPos)) continue;
596 
597  double dist =
598  computeDistFGD(fgdPos, pos, txin, tyin, curvature, (*hit)->IsXHit());
599  if (fabs(dist) < 150) { // cutoff here is 150 mm
600  EnteringFGD = true;
601  break;
602  }
603  }
604  }
605 
606  if (insideFGD2(ExtrapolatedVertexZX, ExtrapolatedVertexZY,
607  ExtrapolatedVertexXX, ExtrapolatedVertexYY)) {
608  for (ND::THitSelection::iterator hit = fgd->begin(); hit != fgd->end();
609  hit++) {
610  TVector3 fgdPos = (*hit)->GetPosition();
611 
612  if (!borderFGD2(fgdPos)) continue;
613 
614  double dist =
615  computeDistFGD(fgdPos, pos, txin, tyin, curvature, (*hit)->IsXHit());
616  if (fabs(dist) < 150) {
617  EnteringFGD = true;
618  break;
619  }
620  }
621  }
622 
623  // now we fill the rest of the info for the object
624  FillTPC(in, out, idx);
625  TTPCTrack* tpc = dynamic_cast<TTPCTrack*>((*out)[idx]);
626  if (!tpc) return false;
627 
628  tpc->EnteringFGD = EnteringFGD;
629  tpc->ExtrapolatedVertexYY = ExtrapolatedVertexYY;
630  tpc->ExtrapolatedVertexXX = ExtrapolatedVertexXX;
631  tpc->ExtrapolatedVertexZY = ExtrapolatedVertexZY;
632  tpc->ExtrapolatedVertexZX = ExtrapolatedVertexZX;
633  // std::copy(ExtrapolatedHits.begin(), ExtrapolatedHits.end(),
634  // std::back_inserter(tpc->ExtrapolatedHits));
635  tpc->HasExtrapolation = true;
636  return true;
637 }
638 //*************************************************************
639 
640 /// called recursively to look for bottom level track that built a PID. In case
641 /// the PID is made of several PIDs, the one with the track with smallest
642 /// relative curvature
643 /// error is returned.
644 ND::THandle<ND::TReconTrack> ND::TTrackerReconModule::GetConstituentTrack(
645  ND::THandle<ND::TReconPID> in, int& nconstit) {
646  ND::THandle<ND::TReconTrack> retval = in;
647 
648  if (in->GetConstituents()) {
649  // std::cout<<"In GetConstituentTrack, nconstit="<<nconstit<<std::endl;
650  ND::TReconObjectContainer::const_iterator it;
651  for (it = in->GetConstituents()->begin();
652  it != in->GetConstituents()->end(); it++) {
653  nconstit++;
654  ND::THandle<ND::TReconTrack> track = (*it);
655  ND::THandle<ND::TReconPID> pid = (*it);
656  if (pid) {
657  // std::cout<<" Found pid, dig deeper"<<std::endl;
658  ND::THandle<ND::TReconTrack> atrack =
659  GetConstituentTrack(pid, nconstit);
660  // if( atrack){
661  // std::cout<<" deeper, found track"<<std::endl;
662  // atrack->ls();
663  //}
664  if (retval && atrack) {
665  // compare atrack to retval, save best to retval
666  ND::THandle<ND::TTrackState> astate = atrack->GetState();
667  ND::THandle<ND::TTrackState> rstate = retval->GetState();
668  if (astate && !rstate) retval = atrack;
669  if (rstate && astate) {
670  double amom, amomerr;
671  double rmom, rmomerr;
672  GetMomentum(astate, amom, amomerr);
673  GetMomentum(rstate, rmom, rmomerr);
674  if (amom > 0.0 && rmom <= 0.0) retval = atrack;
675  if (amom > 0.0 && rmom > 0.0 && amomerr / amom < rmomerr / rmom)
676  retval = atrack;
677  }
678  } else if (atrack) {
679  retval = atrack;
680  }
681  } else if (track) {
682  // std::cout<<" Found track"<<std::endl;
683  // track->ls();
684  if (retval) {
685  // compare track to retval
686  // compare atrack to retval, save best to retval
687  ND::THandle<ND::TTrackState> astate = track->GetState();
688  ND::THandle<ND::TTrackState> rstate = retval->GetState();
689  if (astate && !rstate) retval = track;
690  if (rstate && astate) {
691  double amom, amomerr;
692  double rmom, rmomerr;
693  GetMomentum(astate, amom, amomerr);
694  GetMomentum(rstate, rmom, rmomerr);
695  if (amom > 0.0 && rmom <= 0.0) retval = track;
696  if (amom > 0.0 && rmom > 0.0 && amomerr / amom < rmomerr / rmom)
697  retval = track;
698  }
699  } else {
700  retval = track;
701  }
702  }
703  }
704  }
705 
706  return retval;
707 }
708 
709 //**************************************************************
710 ND::THandle<ND::THitSelection> ND::TTrackerReconModule::GetHitsFromComboHits(
711  ND::THandle<ND::THitSelection> tpc) {
712  //**************************************************************
713 
714  if (!tpc) return ND::THandle<ND::THitSelection>();
715 
716  ND::THandle<ND::THitSelection> hits(new ND::THitSelection());
717 
718  for (ND::THitSelection::const_iterator hsel = tpc->begin();
719  hsel != tpc->end(); hsel++) {
720  ND::THandle<ND::TComboHit> combo(*hsel);
721  if (combo) {
722  for (ND::THitSelection::const_iterator cHit = combo->GetHits().begin();
723  cHit != combo->GetHits().end(); cHit++) {
724  if (*cHit) hits->push_back(*cHit);
725  }
726  }
727  }
728 
729  return hits;
730 }
731 
732 //*************************************************************
733 void ND::TTrackerReconModule::FillOther(ND::THandle<ND::TReconTrack> in,
734  TClonesArray* out, int idx) {
735  //*************************************************************
736 
737  ND::TReconBase::Status dets[5] = {
738  ND::TReconBase::kTPC1, ND::TReconBase::kTPC2, ND::TReconBase::kTPC3,
739  ND::TReconBase::kFGD1, ND::TReconBase::kFGD2};
740 
741  TTrackOther* other = new ((*out)[idx]) TTrackOther;
742 
743  other->AlgorithmName = in->GetAlgorithmName();
744  other->Detector = -1;
745  for (int i = 0; i < 5; i++) {
746  if (in->UsesDetector(dets[i])) other->Detector = i + 1;
747  }
748  ND::THandle<ND::THitSelection> hits = GetHitsFromComboHits(in->GetHits());
749  // get track t0:
750  double t0 = in->GetPosition().T();
751 
752  other->EDeposit = 0.0;
753 
754  // Initialise mint and maxt to remove compiler warnings
755 
756  TVector3 minpos(0, 0, 999999);
757  double mint = 0;
758  TVector3 maxpos(0, 0, -999999);
759  double maxt = 0;
760  for (ND::THitSelection::const_iterator hsel = hits->begin();
761  hsel != hits->end(); hsel++) {
762  TUnusedHit* hit = new ((*(other->Hits))[other->NHits++]) TUnusedHit;
763  // hit->TotalCharge = (*hsel)->GetCharge();
764  hit->Position = (*hsel)->GetPosition();
765  hit->Variance = (*hsel)->GetUncertainty();
766  // hit->Time = (*hsel)->GetTime();
767  hit->TimeUnc = (*hsel)->GetTimeUncertainty();
768 
769  double peakTime = 0;
770  double maxCharge = 0;
771 
772  // Find the pick charge and Time of the Wave form
773  ND::THandle<ND::TMultiHit> mh = *hsel;
774  std::vector<ND::THandle<ND::TSingleHit> >::const_iterator mhit;
775  for (mhit = mh->begin(); mhit != mh->end(); mhit++) {
776  if (maxCharge < (*mhit)->GetCharge()) {
777  maxCharge = (*mhit)->GetCharge();
778  peakTime = (*mhit)->GetTime();
779  }
780  }
781 
782  hit->TotalCharge = maxCharge;
783 
784  ND::TGeometryId geomId = (*hsel)->GetGeomId();
785  double t = peakTime - t0;
786  double x = (fMaxDrift + fCathodeOffset - t * fDriftVelocity) *
787  ND::TGeomInfo::Get().TPC().GetDriftSense(geomId);
788 
789  hit->Position[0] = x;
790  hit->Time = t;
791 
792  other->EDeposit += hit->TotalCharge;
793  if (minpos[2] > hit->Position[2]) {
794  minpos = hit->Position;
795  mint = hit->Time;
796  }
797  if (maxpos[2] < hit->Position[2]) {
798  maxpos = hit->Position;
799  maxt = hit->Time;
800  }
801  }
802 
803  other->FrontPosition = TLorentzVector(minpos[0], minpos[1], minpos[2], mint);
804  other->BackPosition = TLorentzVector(maxpos[0], maxpos[1], maxpos[2], maxt);
805 
806  return;
807 }
808 
809 //*************************************************************
810 // recursively called to fill all constituents, constituents of constituents,
811 // and so on
812 // if new constituent is created return true, otherwise return false
813 bool ND::TTrackerReconModule::FillConstituentPid(TTrackerResult* result,
814  ND::THandle<ND::TReconPID> in,
815  int& idx) {
816  //*************************************************************
817 
818  // First check if this is single detector pid (again only check TPC, since FGD
819  // is track)
820  std::string fDetName = std::string(in->ConvertDetector());
821  if (fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)") {
822  FillTPC(in, result->TPC, result->NTPCs);
823  result->NTPCs++;
824  return false;
825  }
826 
827  // otherwise it is indeed a composite object so create new constituent object
828  TTrackerConstituent* constituent =
829  new ((*(result->Constituents))[result->NTotalConstituents++])
830  TTrackerConstituent;
831  idx = result->NTotalConstituents - 1;
832  constituent->AlgorithmName = in->GetAlgorithmName();
833  constituent->Detectors = GetDetectorNumber(in);
834  constituent->Status = in->CheckStatus(in->kSuccess);
835  constituent->NDOF = in->GetNDOF();
836  constituent->Chi2 = in->GetQuality();
837  if (constituent->NDOF > 0)
838  constituent->Quality = TMath::Prob(constituent->Chi2, constituent->NDOF);
839  else
840  result->Quality = -1.0;
841 
842  // fill constituent state information
843  ND::THandle<ND::TPIDState> state = in->GetState();
844  if (state) {
845  constituent->Position = state->GetPosition();
846  constituent->Variance = state->GetPositionVariance();
847  constituent->Direction = state->GetDirection();
848  constituent->isForward = (constituent->Direction.Z() > 0);
849  constituent->Momentum = state->GetMomentum();
850  constituent->MomentumError = sqrt(fabs(state->GetMomentumVariance()));
851  constituent->Charge = state->GetCharge();
852  } else {
853  constituent->Position = TLorentzVector(0., 0., 0., 0.);
854  constituent->Variance = TLorentzVector(0., 0., 0., 0.);
855  constituent->Direction = TVector3(0., 0., 0.);
856  constituent->isForward = 1;
857  constituent->Momentum = 0.0;
858  constituent->MomentumError = 0.0;
859  constituent->Charge = 0.0;
860  }
861  constituent->EDeposit = 0.0; // No edeposit variable for pid recon objects
862 
863  // Get state at front and back of track
864  constituent->NNodes = in->GetNodes().size();
865  if (constituent->NNodes > 0) {
866  ND::THandle<ND::TPIDState> firstState = in->GetNodes()[0]->GetState();
867  ND::THandle<ND::TPIDState> lastState =
868  in->GetNodes()[in->GetNodes().size() - 1]->GetState();
869 
870  if (firstState) {
871  constituent->FrontPosition = firstState->GetPosition();
872  constituent->FrontVariance = firstState->GetPositionVariance();
873  constituent->FrontDirection = firstState->GetDirection();
874  constituent->FrontMomentum = firstState->GetMomentum();
875  } else {
876  constituent->FrontPosition = TLorentzVector(0., 0., 0., 0.);
877  constituent->FrontVariance = TLorentzVector(0., 0., 0., 0.);
878  constituent->FrontDirection = TVector3(0., 0., 0.);
879  constituent->FrontMomentum = 0.;
880  }
881 
882  if (lastState) {
883  constituent->BackPosition = lastState->GetPosition();
884  constituent->BackVariance = lastState->GetPositionVariance();
885  constituent->BackDirection = lastState->GetDirection();
886  constituent->BackMomentum = lastState->GetMomentum();
887  } else {
888  constituent->BackPosition = TLorentzVector(0., 0., 0., 0.);
889  constituent->BackVariance = TLorentzVector(0., 0., 0., 0.);
890  constituent->BackDirection = TVector3(0., 0., 0.);
891  constituent->BackMomentum = 0.;
892  }
893  }
894 
895  // get number of hits if available
896  constituent->NHits = 0;
897  if (in->GetHits()) constituent->NHits = GetNumberOfHits(in);
898 
899  // Check if this is a single detector object? If it is, then just save it,
900  // and don't worry about
901  // any constituents inside of single detector
902  constituent->NConstituents = 0;
903  // check if this is a composite object... if so fill constituents
904  if (in->GetConstituents()) {
905  ND::TReconObjectContainer::const_iterator it;
906  for (it = in->GetConstituents()->begin();
907  it != in->GetConstituents()->end(); it++) {
908  ND::THandle<ND::TReconTrack> track = (*it);
909  ND::THandle<ND::TReconPID> pid = (*it);
910  if (pid) {
911  if (constituent->NConstituents >= 2) {
912  // std::cout<<"TTrackerReconModule::FillConstituentPID>pid
913  // nconstituents already 2, not adding more"<<std::endl;
914  } else {
915  bool addedconstituent = FillConstituentPid(
916  result, pid, constituent->ConstitIdx[constituent->NConstituents]);
917  if (addedconstituent) constituent->NConstituents++;
918  }
919  } else if (track) {
920  if (constituent->NConstituents >= 2) {
921  // std::cout<<"TTrackerReconModule::FillConstituentPID>track
922  // nconstituents already 2, not adding more"<<std::endl;
923  } else {
924  bool addedconstituent = FillConstituentTrack(
925  result, track,
926  constituent->ConstitIdx[constituent->NConstituents]);
927  if (addedconstituent) constituent->NConstituents++;
928  }
929  }
930  }
931  }
932 
933  return true;
934 }
935 
936 //*************************************************************
937 // recursively called to fill all constituents, constituents of constituents,
938 // and so on
939 // if new constituent is created return true, otherwise return false
940 bool ND::TTrackerReconModule::FillConstituentTrack(
941  TTrackerResult* result, ND::THandle<ND::TReconTrack> in, int& idx) {
942  //*************************************************************
943 
944  // First check if this is single detector
945  std::string fDetName = std::string(in->ConvertDetector());
946  if (fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)") {
947  FillTPC(in, result->TPC, result->NTPCs);
948  result->NTPCs++;
949  return false;
950  }
951  if (fDetName == "(FGD1)" || fDetName == "(FGD2)") {
952  FillFGD(in, result->FGD, result->NFGDs);
953  result->NFGDs++;
954  return false;
955  }
956 
957  // otherwise it is indeed a composite object so create new constituent object
958  TTrackerConstituent* constituent =
959  new ((*(result->Constituents))[result->NTotalConstituents++])
960  TTrackerConstituent;
961  idx = result->NTotalConstituents - 1;
962  constituent->AlgorithmName = in->GetAlgorithmName();
963  constituent->Detectors = GetDetectorNumber(in);
964  constituent->Status = in->CheckStatus(in->kSuccess);
965  constituent->NDOF = in->GetNDOF();
966  constituent->Chi2 = in->GetQuality();
967  if (constituent->NDOF > 0)
968  constituent->Quality = TMath::Prob(constituent->Chi2, constituent->NDOF);
969  else
970  result->Quality = -1.0;
971 
972  // fill constituent state information
973  ND::THandle<ND::TTrackState> state = in->GetState();
974  if (state) {
975  constituent->Position = state->GetPosition();
976  constituent->Variance = state->GetPositionVariance();
977  constituent->Direction = state->GetDirection();
978  constituent->DirectionVariance = state->GetDirectionVariance();
979  constituent->isForward = (constituent->Direction.Z() > 0);
980  GetMomentum(state, constituent->Momentum, constituent->MomentumError);
981  constituent->Charge = GetCharge(state);
982  constituent->EDeposit = state->GetEDeposit();
983  } else {
984  constituent->Position = TLorentzVector(0., 0., 0., 0.);
985  constituent->Variance = TLorentzVector(0., 0., 0., 0.);
986  constituent->Direction = TVector3(0., 0., 0.);
987  constituent->DirectionVariance = TVector3(0., 0., 0.);
988  constituent->isForward = 1;
989  constituent->Momentum = 0.0;
990  constituent->MomentumError = 0.0;
991  constituent->Charge = 0.0;
992  }
993 
994  // Get state at front and back of track
995  constituent->NNodes = in->GetNodes().size();
996  if (constituent->NNodes > 0) {
997  ND::THandle<ND::TTrackState> firstState = in->GetNodes()[0]->GetState();
998  ND::THandle<ND::TTrackState> lastState =
999  in->GetNodes()[in->GetNodes().size() - 1]->GetState();
1000 
1001  if (firstState) {
1002  constituent->FrontPosition = firstState->GetPosition();
1003  constituent->FrontVariance = firstState->GetPositionVariance();
1004  constituent->FrontDirection = firstState->GetDirection();
1005  double tmperr;
1006  GetMomentum(firstState, constituent->FrontMomentum, tmperr);
1007  } else {
1008  constituent->FrontPosition = TLorentzVector(0., 0., 0., 0.);
1009  constituent->FrontVariance = TLorentzVector(0., 0., 0., 0.);
1010  constituent->FrontDirection = TVector3(0., 0., 0.);
1011  constituent->FrontMomentum = 0.;
1012  }
1013 
1014  if (lastState) {
1015  constituent->BackPosition = lastState->GetPosition();
1016  constituent->BackVariance = lastState->GetPositionVariance();
1017  constituent->BackDirection = lastState->GetDirection();
1018  double tmperr;
1019  GetMomentum(lastState, constituent->BackMomentum, tmperr);
1020  } else {
1021  constituent->BackPosition = TLorentzVector(0., 0., 0., 0.);
1022  constituent->BackVariance = TLorentzVector(0., 0., 0., 0.);
1023  constituent->BackDirection = TVector3(0., 0., 0.);
1024  constituent->BackMomentum = 0.;
1025  }
1026  }
1027 
1028  // get number of hits if available
1029  constituent->NHits = 0;
1030  if (in->GetHits()) constituent->NHits = GetNumberOfHits(in);
1031 
1032  // Check if this is a single detector object? If it is, then just save it,
1033  // and don't worry about
1034  // any constituents inside of single detector
1035  constituent->NConstituents = 0;
1036  // check if this is a composite object... if so fill constituents
1037  if (in->GetConstituents()) {
1038  ND::TReconObjectContainer::const_iterator it;
1039  for (it = in->GetConstituents()->begin();
1040  it != in->GetConstituents()->end(); it++) {
1041  ND::THandle<ND::TReconTrack> track = (*it);
1042  ND::THandle<ND::TReconPID> pid = (*it);
1043  if (pid) {
1044  if (constituent->NConstituents >= 2) {
1045  // std::cout<<"TTrackerReconModule::FillConstituentTrack>pid
1046  // nconstituents already 2, not adding more"<<std::endl;
1047  } else {
1048  bool addedconstituent = FillConstituentPid(
1049  result, pid, constituent->ConstitIdx[constituent->NConstituents]);
1050  if (addedconstituent) constituent->NConstituents++;
1051  }
1052  } else if (track) {
1053  if (constituent->NConstituents >= 2) {
1054  // std::cout<<"TTrackerReconModule::FillConstituentPID>track
1055  // nconstituents already 2, not adding more"<<std::endl;
1056  } else {
1057  bool addedconstituent = FillConstituentTrack(
1058  result, track,
1059  constituent->ConstitIdx[constituent->NConstituents]);
1060  if (addedconstituent) constituent->NConstituents++;
1061  }
1062  }
1063  }
1064  }
1065 
1066  return true;
1067 }
1068 
1069 //*************************************************************
1070 void ND::TTrackerReconModule::FillFinalPid(TTrackerResult* result,
1071  ND::THandle<ND::TReconPID> reco) {
1072  //*************************************************************
1073 
1074  // Get state at vertex
1075  ND::THandle<ND::TPIDState> state = reco->GetState();
1076  if (state) {
1077  result->Position = state->GetPosition();
1078  result->Variance = state->GetPositionVariance();
1079  result->Direction = state->GetDirection();
1080  result->DirectionVariance = state->GetDirectionVariance();
1081  result->isForward = (result->Direction.Z() > 0);
1082  result->Momentum = state->GetMomentum();
1083  result->MomentumError = sqrt(fabs(state->GetMomentumVariance()));
1084  result->Charge = state->GetCharge();
1085  } else {
1086  result->Position = TLorentzVector(0., 0., 0., 0.);
1087  result->Variance = TLorentzVector(0., 0., 0., 0.);
1088  result->Direction = TVector3(0., 0., 0.);
1089  result->DirectionVariance = TVector3(0., 0., 0.);
1090  result->isForward = 1;
1091  result->Momentum = 0.;
1092  result->MomentumError = 0.;
1093  result->Charge = 0.;
1094  }
1095  result->EDeposit = 0.; // No edeposit variable for pid recon objects
1096 
1097  // Get state at front and back of track
1098  result->NNodes = reco->GetNodes().size();
1099  if (result->NNodes > 0) {
1100  // Loop over nodes
1101  for (int inode = 0; inode < result->NNodes; inode++) {
1102  ND::THandle<ND::TPIDState> iState = reco->GetNodes()[inode]->GetState();
1103  TTrackerNode* node = new ((*(result->Nodes))[inode]) TTrackerNode;
1104  node->EDeposit = 0.0;
1105  if (iState) {
1106  node->Charge = iState->GetCharge();
1107  node->Position = iState->GetPosition();
1108  node->Variance = iState->GetPositionVariance();
1109  node->Direction = iState->GetDirection();
1110  node->DirectionVariance = iState->GetDirectionVariance();
1111  node->Momentum = iState->GetMomentum();
1112  node->MomentumError = sqrt(fabs(iState->GetMomentumVariance()));
1113  } else {
1114  node->Charge = 0.0;
1115  node->Position = TLorentzVector(0., 0., 0., 0.);
1116  node->Variance = TLorentzVector(0., 0., 0., 0.);
1117  node->Direction = TVector3(0., 0., 0.);
1118  node->DirectionVariance = TVector3(0., 0., 0.);
1119  node->Momentum = 0.;
1120  node->MomentumError = 0.;
1121  }
1122  }
1123  }
1124 
1125  // use recpack to find length of track
1126  Trajectory traj;
1127  ND::THandle<ND::TReconBase> object = reco;
1128  bool ok = ND::converter().TReconBase_to_Trajectory(object, traj);
1129  result->Length = 0;
1130  if (ok) {
1131  double length = 0;
1132  ok = ND::rpman().matching_svc().compute_length(traj, length);
1133  if (ok) result->Length = length;
1134  }
1135 
1136  // get number of hits if available
1137  result->NHits = 0;
1138  if (reco->GetHits()) result->NHits = GetNumberOfHits(reco);
1139 
1140  // get Pid information
1141  result->Likelihoods.push_back(reco->GetPIDWeight());
1142  result->Pids.push_back(reco->GetParticleId());
1143  for (ND::TReconObjectContainer::const_iterator it =
1144  reco->GetAlternates().begin();
1145  it != reco->GetAlternates().end(); it++) {
1146  ND::THandle<ND::TReconPID> alter = *it;
1147  result->Likelihoods.push_back(alter->GetPIDWeight());
1148  result->Pids.push_back(alter->GetParticleId());
1149  }
1150 
1151  // Check if this is a single detector object? If it is, then just save it,
1152  // and don't worry about
1153  // any constituents inside of single detector
1154  result->NConstituents = 0;
1155  std::string fDetName = std::string(reco->ConvertDetector());
1156  if (fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)") {
1157  FillTPC(reco, result->TPC, result->NTPCs);
1158  result->NTPCs++;
1159  } else {
1160  // Note that FGD only occurs as TReconTrack, so no need to check that here.
1161  // Instead this must be a composite object... so fill constituents
1162  if (reco->GetConstituents()) {
1163  ND::TReconObjectContainer::const_iterator it;
1164  int iconstit = 0;
1165  for (it = reco->GetConstituents()->begin();
1166  it != reco->GetConstituents()->end(); it++) {
1167  ND::THandle<ND::TReconTrack> track = (*it);
1168  ND::THandle<ND::TReconPID> pid = (*it);
1169  if (pid) {
1170  if (result->NConstituents >= 2) {
1171  // std::cout<<"TTrackerReconModule::FillFinalPid>pid nconstituents
1172  // already 2, not adding more"<<std::endl;
1173  } else {
1174  bool addedconstituent = FillConstituentPid(
1175  result, pid, result->ConstitIdx[result->NConstituents]);
1176  if (addedconstituent) result->NConstituents++;
1177  }
1178  } else if (track) {
1179  if (result->NConstituents >= 2) {
1180  // std::cout<<"TTrackerReconModule::FillFinalPid>track nconstituents
1181  // already 2, not adding more"<<std::endl;
1182  } else {
1183  bool addedconstituent = FillConstituentTrack(
1184  result, track, result->ConstitIdx[result->NConstituents]);
1185  if (addedconstituent) result->NConstituents++;
1186  }
1187  }
1188  }
1189  iconstit++;
1190  }
1191  }
1192 
1193  return;
1194 }
1195 
1196 //*************************************************************
1197 void ND::TTrackerReconModule::GetMomentum(ND::THandle<ND::TTrackState> state,
1198  double& p, double& ep) {
1199  //*************************************************************
1200 
1201  // need to read in magnet field...
1202  double B = ND::TFieldManager::GetFieldValue(state->GetPosition()).Mag() /
1203  unit::tesla;
1204 
1205  // should get speed of light constant ~0.3 from somewhere
1206  p = 0.3 * B /
1207  (state->GetCurvature() *
1208  TMath::Sqrt(1. - state->GetDirection()[0] * state->GetDirection()[0]));
1209 
1210  CLHEP::HepMatrix Jacobian(2, 2, 0);
1211  CLHEP::HepMatrix Covariance(2, 2, 0);
1212  CLHEP::HepMatrix Cov(2, 2, 0);
1213 
1214  Covariance[0][0] = state->GetCovarianceValue(state->GetDirectionIndex(),
1215  state->GetDirectionIndex());
1216  Covariance[0][1] = state->GetCovarianceValue(state->GetDirectionIndex(),
1217  state->GetCurvatureIndex());
1218 
1219  Covariance[1][0] = state->GetCovarianceValue(state->GetCurvatureIndex(),
1220  state->GetDirectionIndex());
1221  Covariance[1][1] = state->GetCovarianceValue(state->GetCurvatureIndex(),
1222  state->GetCurvatureIndex());
1223 
1224  Jacobian[0][0] = 1.;
1225  Jacobian[0][1] = 0.;
1226 
1227  Jacobian[1][0] = p /
1228  (1. - state->GetDirection()[0] * state->GetDirection()[0]) *
1229  TMath::Abs(state->GetDirection()[0]);
1230  Jacobian[1][1] = -p / state->GetCurvature();
1231 
1232  Cov = Jacobian * Covariance * Jacobian.T();
1233 
1234  ep = TMath::Sqrt(Cov[1][1]);
1235 
1236  p = TMath::Abs(p);
1237 
1238  return;
1239 }
1240 
1241 //*************************************************************
1242 double ND::TTrackerReconModule::GetCharge(ND::THandle<ND::TTrackState> state) {
1243  //*************************************************************
1244  double q = (state->GetCurvature() > 0.0 ? -1.0 : 1.0);
1245  if (state->GetDirection()[2] < 0.0) q *= -1.0;
1246  return q;
1247 }
1248 
1249 //*************************************************************
1250 void ND::TTrackerReconModule::FillFinalTrack(
1251  TTrackerResult* result, ND::THandle<ND::TReconTrack> reco) {
1252  //*************************************************************
1253 
1254  // Get state at vertex
1255  ND::THandle<ND::TTrackState> state = reco->GetState();
1256  if (state) {
1257  result->Position = state->GetPosition();
1258  result->Variance = state->GetPositionVariance();
1259  result->Direction = state->GetDirection();
1260  result->DirectionVariance = state->GetDirectionVariance();
1261  result->isForward = (result->Direction.Z() > 0);
1262  result->EDeposit = state->GetEDeposit();
1263  // need conversion from curvature to momentum
1264  GetMomentum(state, result->Momentum, result->MomentumError);
1265  // need conversion from curvature, and direction to charge
1266  result->Charge = GetCharge(state);
1267  } else {
1268  result->Position = TLorentzVector(0., 0., 0., 0.);
1269  result->Variance = TLorentzVector(0., 0., 0., 0.);
1270  result->Direction = TVector3(0., 0., 0.);
1271  result->DirectionVariance = TVector3(0., 0., 0.);
1272  result->isForward = 1;
1273  result->Momentum = 0.;
1274  result->MomentumError = 0.;
1275  result->Charge = 0.;
1276  }
1277 
1278  // Get state at front and back of track
1279  result->NNodes = reco->GetNodes().size();
1280  if (result->NNodes > 0) {
1281  // Loop over nodes
1282  for (int inode = 0; inode < result->NNodes; inode++) {
1283  ND::THandle<ND::TTrackState> iState = reco->GetNodes()[inode]->GetState();
1284  TTrackerNode* node = new ((*(result->Nodes))[inode]) TTrackerNode;
1285  if (iState) {
1286  node->EDeposit = iState->GetEDeposit();
1287  node->Position = iState->GetPosition();
1288  node->Variance = iState->GetPositionVariance();
1289  node->Direction = iState->GetDirection();
1290  node->DirectionVariance = iState->GetDirectionVariance();
1291  // need conversion from curvature, and direction to charge
1292  GetMomentum(iState, node->Momentum, node->MomentumError);
1293  node->Charge = GetCharge(iState);
1294  } else {
1295  node->EDeposit = 0.0;
1296  node->Charge = 0.0;
1297  node->Position = TLorentzVector(0., 0., 0., 0.);
1298  node->Variance = TLorentzVector(0., 0., 0., 0.);
1299  node->Direction = TVector3(0., 0., 0.);
1300  node->DirectionVariance = TVector3(0., 0., 0.);
1301  node->Momentum = 0.;
1302  node->MomentumError = 0.;
1303  }
1304  }
1305  }
1306 
1307  // use recpack to find length of track
1308  Trajectory traj;
1309  ND::THandle<ND::TReconBase> object = reco;
1310  bool ok = ND::converter().TReconBase_to_Trajectory(object, traj);
1311  result->Length = 0;
1312  if (ok) {
1313  double length = 0;
1314  ok = ND::rpman().matching_svc().compute_length(traj, length);
1315  if (ok) result->Length = length;
1316  }
1317 
1318  // get number of hits if available
1319  result->NHits = 0;
1320  if (reco->GetHits()) result->NHits = GetNumberOfHits(reco);
1321 
1322  // Check if this is a single detector object? If it is, then just save it,
1323  // and don't worry about
1324  // any constituents inside of single detector
1325  result->NConstituents = 0;
1326  std::string fDetName = std::string(reco->ConvertDetector());
1327  if (fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)") {
1328  // I don't think I should ever get here, since I should have found TReconPID
1329  // for the TPC instead
1330  FillTPC(reco, result->TPC, result->NTPCs);
1331  result->NTPCs++;
1332  } else if (fDetName == "(FGD1)" || fDetName == "(FGD2)") {
1333  FillFGD(reco, result->FGD, result->NFGDs);
1334  result->NFGDs++;
1335  } else {
1336  // Instead this must be a composite object... so fill constituents
1337  if (reco->GetConstituents()) {
1338  ND::TReconObjectContainer::const_iterator it;
1339  for (it = reco->GetConstituents()->begin();
1340  it != reco->GetConstituents()->end(); it++) {
1341  ND::THandle<ND::TReconTrack> track = (*it);
1342  ND::THandle<ND::TReconPID> pid = (*it);
1343  if (pid) {
1344  if (result->NConstituents >= 2) {
1345  // std::cout<<"TTrackerReconModule::FillFinalTrack>pid nconstituents
1346  // already 2, not adding more"<<std::endl;
1347  } else {
1348  bool addedconstituent = FillConstituentPid(
1349  result, pid, result->ConstitIdx[result->NConstituents]);
1350  if (addedconstituent) result->NConstituents++;
1351  }
1352  } else if (track) {
1353  if (result->NConstituents >= 2) {
1354  // std::cout<<"TTrackerReconModule::FillFinalPid>track nconstituents
1355  // already 2, not adding more"<<std::endl;
1356  } else {
1357  bool addedconstituent = FillConstituentTrack(
1358  result, track, result->ConstitIdx[result->NConstituents]);
1359  if (addedconstituent) result->NConstituents++;
1360  }
1361  }
1362  }
1363  }
1364  }
1365 
1366  return;
1367 }
1368 
1369 //*****************************************************
1370 int ND::TTrackerReconModule::GetDetectorNumber(
1371  ND::THandle<ND::TReconBase> object) {
1372  //*****************************************************
1373 
1374  int det = 0;
1375 
1376  if (object->UsesDetector(ND::TReconBase::kTPC1)) det = 1;
1377  if (object->UsesDetector(ND::TReconBase::kTPC2)) det = det * 10 + 2;
1378  if (object->UsesDetector(ND::TReconBase::kTPC3)) det = det * 10 + 3;
1379  if (object->UsesDetector(ND::TReconBase::kFGD1)) det = det * 10 + 4;
1380  if (object->UsesDetector(ND::TReconBase::kFGD2)) det = det * 10 + 5;
1381  if (object->UsesDetector(ND::TReconBase::kDSECal)) det = det * 10 + 6;
1382  if (object->UsesDetector(ND::TReconBase::kTECal) ||
1383  object->UsesDetector(ND::TReconBase::kPECal))
1384  det = det * 10 + 7;
1385  if (object->UsesDetector(ND::TReconBase::kP0D)) det = det * 10 + 8;
1386  if (object->UsesDetector(ND::TReconBase::kSMRD)) det = det * 10 + 9;
1387 
1388  return det;
1389 }
1390 
1391 //*************************************************************
1392 bool ND::TTrackerReconModule::FillTree(ND::TND280Event& event) {
1393  //*************************************************************
1394  // ND280Debug("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!filling the TRACKER tree ");
1395 
1396  fTracks->Clear();
1397 
1398  fNVertices = 0;
1399  fNTracks = 0;
1400  fNFGDOther = 0;
1401  fNTPCOther = 0;
1402  fNTPCIso = 0;
1403  fNFGDUnused = 0;
1404  fNTPCUnused = 0;
1405  fNTPCExtrapolation = 0;
1406 
1407  fVertices->Clear();
1408  fTracks->Clear();
1409  fFGDOther->Clear();
1410  fTPCOther->Clear();
1411  fTPCIso->Clear();
1412  fTPCUnused->Clear();
1413  fFGDUnused->Clear();
1414  fTPCExtrapolation->Clear();
1415 
1416  // Get the main result to save.
1417  ND::THandle<ND::TAlgorithmResult> Result = event.GetFit("trackerRecon");
1418  if (!Result) return true;
1419 
1420  // Now go through each of the different ReconObjects to save:
1421  ND::THandle<ND::TReconObjectContainer> final =
1422  Result->GetResultsContainer("final");
1423  ND::THandle<ND::TReconObjectContainer> fgdOther =
1424  Result->GetResultsContainer("FGD:other");
1425  ND::THandle<ND::TReconObjectContainer> tpcOther =
1426  Result->GetResultsContainer("TPC:other");
1427  ND::THandle<ND::TReconObjectContainer> tpcIso =
1428  Result->GetResultsContainer("TPC:IsoObjects");
1429  ND::THandle<ND::TReconObjectContainer> vertices =
1430  Result->GetResultsContainer("vertices");
1431  ND::THandle<ND::THitSelection> tpcUnused =
1432  Result->GetHitSelection("TPC:unused");
1433  ND::THandle<ND::THitSelection> fgdUnused =
1434  Result->GetHitSelection("FGD:unused");
1435 
1436  // fill fgd other branches
1437  if (fgdOther) {
1438  for (ND::TReconObjectContainer::iterator it = fgdOther->begin();
1439  it != fgdOther->end(); it++) {
1440  ND::THandle<ND::TReconTrack> track = (*it);
1441  if (track) {
1442  FillOther(track, fFGDOther, fNFGDOther);
1443  fNFGDOther++;
1444  }
1445  }
1446  }
1447 
1448  // fill tpc other branches
1449  if (tpcOther) {
1450  for (ND::TReconObjectContainer::iterator it = tpcOther->begin();
1451  it != tpcOther->end(); it++) {
1452  ND::THandle<ND::TReconTrack> track = (*it);
1453  if (track) {
1454  FillOther(track, fTPCOther, fNTPCOther);
1455  fNTPCOther++;
1456  }
1457  }
1458  }
1459 
1460  // fill tpc isolated track branches
1461  // Are there ever any of these?
1462  if (tpcIso) {
1463  for (ND::TReconObjectContainer::iterator it = tpcIso->begin();
1464  it != tpcIso->end(); it++) {
1465  ND::THandle<ND::TReconPID> pid = (*it);
1466  if (pid) {
1467  FillTPC(pid, fTPCIso, fNTPCIso);
1468  fNTPCIso++;
1469  }
1470  }
1471  }
1472 
1473  // fill vertices
1474  // Again, are there ever any of these?
1475  if (vertices) {
1476  for (ND::TReconObjectContainer::iterator it = vertices->begin();
1477  it != vertices->end(); it++) {
1478  ND::THandle<ND::TReconVertex> vert = (*it);
1479  if (vert) {
1480  FillVertex(vert, fVertices, fNVertices);
1481  fNVertices++;
1482  }
1483  }
1484  }
1485 
1486  // fill unused fgd hits
1487  if (fgdUnused) {
1488  for (ND::THitSelection::const_iterator hsel = fgdUnused->begin();
1489  hsel != fgdUnused->end(); hsel++) {
1490  TUnusedHit* hit = new ((*fFGDUnused)[fNFGDUnused++]) TUnusedHit;
1491  hit->TotalCharge = (*hsel)->GetCharge();
1492  hit->Position = (*hsel)->GetPosition();
1493  hit->Variance = (*hsel)->GetUncertainty();
1494  hit->Time = (*hsel)->GetTime();
1495  hit->TimeUnc = (*hsel)->GetTimeUncertainty();
1496  }
1497  }
1498 
1499  // Now loop over the final fitted tracks or pids
1500  double t0avg = 0.0;
1501  int NT0 = 0;
1502  if (final) {
1503  for (ND::TReconObjectContainer::iterator tt = final->begin();
1504  tt != final->end(); tt++) {
1505  ND::THandle<ND::TReconPID> pid = *tt;
1506  ND::THandle<ND::TReconTrack> track = *tt;
1507  ND::THandle<ND::TReconBase> reco = *tt;
1508  if (reco) {
1509  TTrackerResult* result = new ((*fTracks)[fNTracks++]) TTrackerResult;
1510  // Save the unique object ID to allow matching to global recon objects
1511  result->UniqueID = reco->GetUniqueID();
1512  result->AlgorithmName = reco->GetAlgorithmName();
1513  result->Detectors = GetDetectorNumber(reco);
1514  result->Status = reco->CheckStatus(reco->kSuccess);
1515  result->NDOF = reco->GetNDOF();
1516  result->Chi2 = reco->GetQuality();
1517 
1518  CheckMatchingFailure(result, reco);
1519 
1520  if (result->NDOF > 0)
1521  result->Quality = TMath::Prob(result->Chi2, result->NDOF);
1522  else
1523  result->Quality = -1.0;
1524 
1525  if (pid) {
1526  ND::THandle<ND::TPIDState> state = pid->GetState();
1527  if (state) {
1528  t0avg += state->GetPosition().T();
1529  NT0++;
1530  }
1531  FillFinalPid(result, pid);
1532  } else if (track) {
1533  FillFinalTrack(result, track);
1534  }
1535 
1536  // Using oaUtility tool to associate the track and the a G4track
1537  ND::TG4Trajectory* G4track;
1538  G4track = NULL;
1539  double pur = -0xABCDEF;
1540  double eff = -0xABCDEF;
1541  int G4TrkId = TrackTruthInfo::GetG4TrajectoryID(reco, eff, pur);
1542  ND::THandle<ND::TG4TrajectoryContainer> G4trajectories =
1543  event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
1544  if (G4trajectories) {
1545  ND::TG4TrajectoryContainer::iterator trajIter =
1546  G4trajectories->begin();
1547  ND::TG4TrajectoryContainer::iterator trajEnd = G4trajectories->end();
1548  for (; trajIter != trajEnd; ++trajIter) {
1549  ND::TG4Trajectory* trajectory = &(trajIter->second);
1550  if (trajectory->GetTrackId() == G4TrkId) {
1551  G4track = new ND::TG4Trajectory(*trajectory);
1552  break;
1553  }
1554  }
1555  }
1556  if (G4track != NULL) {
1557  FillTrueParticle(G4track, pur, eff, result->TrueParticle);
1558  }
1559  }
1560  }
1561  }
1562 
1563  if (NT0) t0avg /= double(NT0);
1564 
1565  // fill unused tpc hits... note we will make use of the average track time
1566  // to guess the t0 of the unused hits.
1567  if (tpcUnused) {
1568  ND::THandle<ND::THitSelection> hstpc = GetHitsFromComboHits(tpcUnused);
1569  for (ND::THitSelection::const_iterator hsel = hstpc->begin();
1570  hsel != hstpc->end(); hsel++) {
1571  TUnusedHit* hit = new ((*fTPCUnused)[fNTPCUnused++]) TUnusedHit;
1572  // hit->TotalCharge = (*hsel)->GetCharge();
1573  hit->Position = (*hsel)->GetPosition();
1574  hit->Variance = (*hsel)->GetUncertainty();
1575  // hit->Time = (*hsel)->GetTime();
1576  hit->TimeUnc = (*hsel)->GetTimeUncertainty();
1577 
1578  double peakTime = 0;
1579  double maxCharge = 0;
1580 
1581  // Find the pick charge and Time of the Wave form
1582  ND::THandle<ND::TMultiHit> mh = *hsel;
1583  std::vector<ND::THandle<ND::TSingleHit> >::const_iterator mhit;
1584  for (mhit = mh->begin(); mhit != mh->end(); mhit++) {
1585  if (maxCharge < (*mhit)->GetCharge()) {
1586  maxCharge = (*mhit)->GetCharge();
1587  peakTime = (*mhit)->GetTime();
1588  }
1589  }
1590 
1591  hit->TotalCharge = maxCharge;
1592 
1593  ND::TGeometryId geomId = (*hsel)->GetGeomId();
1594  double t = peakTime - t0avg;
1595  double x = (fMaxDrift + fCathodeOffset - t * fDriftVelocity) *
1596  ND::TGeomInfo::Get().TPC().GetDriftSense(geomId);
1597 
1598  hit->Position[0] = x;
1599  hit->Time = t;
1600  }
1601  }
1602 
1603  // Do the hit association with Claudio's algorithm; this uses the raw TPC
1604  // result (DR, Sept '10)
1605  ND::THandle<ND::TAlgorithmResult> TpcResult = event.GetFit("tpcRecon");
1606  if (!TpcResult) return true;
1607 
1608  ND::THandle<ND::TReconObjectContainer> TrackCont =
1609  TpcResult->GetResultsContainer("TPCPids");
1610  ND::THandle<ND::THitSelection> fgdHits = event.GetHitSelection("fgd");
1611  if (TrackCont && fgdHits) {
1612  for (ND::TReconObjectContainer::iterator it = TrackCont->begin();
1613  it != TrackCont->end(); it++) {
1614  ND::THandle<ND::TReconPID> pid = *it;
1615  if (pid) {
1616  if (ExtrapolateTPC(pid, fTPCExtrapolation, fNTPCExtrapolation, fgdHits))
1617  fNTPCExtrapolation++;
1618  }
1619  }
1620  }
1621 
1622  return true;
1623 }
1624 
1625 ND::TG4Trajectory ND::TTrackerReconModule::GetTrajectory(
1626  ND::TND280Event& event, ND::THandle<ND::TReconBase> reco) {
1627  ND::TG4Trajectory traj;
1628  ND::THandle<ND::THitSelection> hits = reco->GetHits();
1629 
1630  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
1631  event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
1632  if (!trajectories) return ND::TG4Trajectory();
1633  std::vector<int> search;
1634 
1635  search.resize(trajectories->size());
1636 
1637  for (ND::THitSelection::const_iterator ht = hits->begin(); ht != hits->end();
1638  ht++) { // Loop over reconstructed hits
1639  // Get Combo hit
1640  ND::THandle<ND::TComboHit> ch = *ht;
1641 
1642  // Loop over hits in Combo
1643  for (ND::THitSelection::const_iterator ch_member = ch->GetHits().begin();
1644  ch_member != ch->GetHits().end(); ++ch_member) {
1645  ND::THandle<ND::TMultiHit> mh = *ch_member;
1646 
1647  if (!mh) { // this is the old method with no waveforms.
1648  std::vector<ND::TG4VHit*> contributors =
1649  HitTruthInfo::GetHitTruthInfo(*ch_member);
1650 
1651  for (std::vector<ND::TG4VHit*>::const_iterator h = contributors.begin();
1652  h != contributors.end(); ++h) {
1653  if (!*h) continue;
1654 
1655  ND::TG4HitSegment* g4hitseg = dynamic_cast<ND::TG4HitSegment*>(*h);
1656  if (!g4hitseg) continue;
1657 
1658  for (unsigned int i = 0; i != g4hitseg->GetContributors().size();
1659  ++i) {
1660  unsigned int vv = g4hitseg->GetContributors()[i];
1661  if (vv >= search.size()) search.resize(vv + 1);
1662  search[vv] = search[vv] + 1;
1663  }
1664  }
1665  } else { // This is the new method with waveform
1666  for (ND::TMultiHit::iterator wfit = mh->begin(); wfit != mh->end();
1667  wfit++) {
1668  std::vector<ND::TG4VHit*> contributors =
1669  HitTruthInfo::GetHitTruthInfo(*wfit);
1670 
1671  for (std::vector<ND::TG4VHit*>::const_iterator h =
1672  contributors.begin();
1673  h != contributors.end(); ++h) {
1674  if (!*h) continue;
1675 
1676  ND::TG4HitSegment* g4hitseg = dynamic_cast<ND::TG4HitSegment*>(*h);
1677  if (!g4hitseg) continue;
1678 
1679  for (unsigned int i = 0; i != g4hitseg->GetContributors().size();
1680  ++i) {
1681  unsigned int vv = g4hitseg->GetContributors()[i];
1682  if (vv >= search.size()) search.resize(vv + 1);
1683  search[vv] = search[vv] + 1;
1684  }
1685  }
1686  }
1687  }
1688  }
1689  }
1690 
1691  int maxim = 0;
1692  int imax = -1;
1693 
1694  for (unsigned int i = 0; i < search.size(); i++) {
1695  if (search[i] > maxim) {
1696  maxim = search[i];
1697  imax = i;
1698  }
1699  }
1700 
1701  if (imax >= 0) {
1702  return (*trajectories)[imax];
1703  } else {
1704  std::cerr << " TRACK NOT FOUND " << hits->size() << " " << search.size()
1705  << std::endl;
1706  return ND::TG4Trajectory();
1707  }
1708 }
1709 
1710 //*************************************************************
1711 void ND::TTrackerReconModule::FillTrueParticle(
1712  ND::TG4Trajectory* G4track, double pur, double eff,
1713  ND::TTrackerReconModule::TTrueParticle& trueParticle) {
1714  //*************************************************************
1715 
1716  if (G4track) {
1717  trueParticle.ID = G4track->GetTrackId();
1718  trueParticle.Pur = pur;
1719  trueParticle.Eff = eff;
1720 
1721  ND::TG4PrimaryVertex G4vertex;
1722  bool ok = GetG4Vertex(G4track, G4vertex);
1723  FillTrueVertex(ok, G4vertex, 0, 0, trueParticle.Vertex);
1724 
1725  } else {
1726  // std::cout<<"True G4 track not found"<<std::endl;
1727  trueParticle.ID = -1;
1728 
1729  ND::TG4PrimaryVertex G4vertex;
1730  FillTrueVertex(false, G4vertex, 0, 0, trueParticle.Vertex);
1731  }
1732  return;
1733 }
1734 
1735 //*************************************************************
1736 ND::THandle<ND::TG4Trajectory> ND::TTrackerReconModule::GetParent(
1737  ND::TG4Trajectory* G4track) {
1738  //*************************************************************
1739 
1740  const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
1741  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
1742  event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
1743 
1744  if (G4track->GetParentId() != 0) {
1745  ND::THandle<ND::TG4Trajectory> traj(
1746  new ND::TG4Trajectory(((*trajectories)[G4track->GetParentId()])));
1747  return traj;
1748  } else
1749  return ND::THandle<ND::TG4Trajectory>();
1750 }
1751 
1752 //*************************************************************
1753 ND::THandle<ND::TG4Trajectory> ND::TTrackerReconModule::GetParent(
1754  ND::THandle<TG4Trajectory> G4track) {
1755  //*************************************************************
1756 
1757  const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
1758  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
1759  event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
1760 
1761  if (G4track->GetParentId() != 0) {
1762  ND::THandle<ND::TG4Trajectory> traj(
1763  new ND::TG4Trajectory(((*trajectories)[G4track->GetParentId()])));
1764  return traj;
1765  } else
1766  return ND::THandle<ND::TG4Trajectory>();
1767 }
1768 
1769 //*************************************************************
1770 void ND::TTrackerReconModule::FillTrueVertex(
1771  bool found, const ND::TG4PrimaryVertex& G4vertex, double pur, double eff,
1772  ND::TTrueVertex& vertex) {
1773  //*************************************************************
1774 
1775  if (found) {
1776  vertex.Position = G4vertex.GetPosition();
1777  vertex.ID = G4vertex.GetInteractionNumber();
1778  } else {
1779  vertex.Position = TLorentzVector(0, 0, 0, 0);
1780  vertex.ID = -1;
1781  }
1782 }
1783 
1784 //*************************************************************
1785 bool ND::TTrackerReconModule::GetG4Vertex(ND::TG4Trajectory* G4track,
1786  ND::TG4PrimaryVertex& G4vertex) {
1787  //*************************************************************
1788 
1789  const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
1790  ND::THandle<ND::TG4PrimaryVertexContainer> vertices =
1791  event.Get<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00");
1792 
1793  if (vertices) {
1794  std::vector<ND::TG4PrimaryVertex>::iterator it1;
1795  for (it1 = vertices->begin(); it1 != vertices->end(); it1++) {
1796  const ND::TG4PrimaryParticleContainer& particles =
1797  it1->GetPrimaryParticles();
1798 
1799  std::vector<ND::TG4PrimaryParticle>::const_iterator it2;
1800  for (it2 = particles.begin(); it2 != particles.end(); it2++) {
1801  if (it2->GetTrackId() == G4track->GetTrackId()) {
1802  G4vertex = *it1;
1803  return true;
1804  }
1805  }
1806  }
1807  }
1808 
1809  return false;
1810 }
1811 
1812 //*************************************************************
1813 bool ND::TTrackerReconModule::GetIncomingParticle(
1814  const ND::TG4PrimaryVertex& G4vertex, ND::TG4PrimaryParticle& incoming) {
1815  //*************************************************************
1816 
1817  for (std::vector<ND::TG4PrimaryVertex>::const_iterator infoVtxIter =
1818  G4vertex.GetInfoVertex().begin();
1819  infoVtxIter != G4vertex.GetInfoVertex().end(); ++infoVtxIter) {
1820  std::vector<ND::TG4PrimaryParticle>::const_iterator particleIter;
1821  for (particleIter = infoVtxIter->GetPrimaryParticles().begin();
1822  particleIter != infoVtxIter->GetPrimaryParticles().end();
1823  ++particleIter) {
1824  Int_t pdg = particleIter->GetPDGCode();
1825  if (particleIter->GetTrackId() == -1) {
1826  // is an incoming particle (including target)
1827  if (TMath::Abs(pdg) == 12 || TMath::Abs(pdg) == 14 ||
1828  TMath::Abs(pdg) == 16) {
1829  incoming = *particleIter;
1830  return true;
1831  }
1832  }
1833  }
1834  }
1835 
1836  return false;
1837 }
1838 
1839 //*************************************************************
1840 int ND::TTrackerReconModule::GetNumberOfHits(
1841  ND::THandle<ND::TReconBase> object) {
1842  //*************************************************************
1843 
1844  int nHits = 0;
1845  if (object->GetHits())
1846  nHits += object->GetHits()->size();
1847  else if (object->GetConstituents()) {
1848  ND::TReconObjectContainer::const_iterator it;
1849  for (it = object->GetConstituents()->begin();
1850  it != object->GetConstituents()->end(); it++) {
1851  ND::THandle<ND::TReconBase> obj = (*it);
1852  if (obj) {
1853  nHits += GetNumberOfHits(obj);
1854  }
1855  }
1856  }
1857 
1858  return nHits;
1859 }
1860 
1861 //*************************************************************
1862 void ND::TTrackerReconModule::CheckMatchingFailure(
1863  TTrackerResult* result, ND::THandle<ND::TReconBase> reco) {
1864  //*************************************************************
1865 
1866  bool useFGD1 = reco->UsesDetector(ND::TReconBase::kFGD1);
1867  bool useFGD2 = reco->UsesDetector(ND::TReconBase::kFGD2);
1868 
1869  ND::THandle<ND::TReconBase> fgdCons;
1870  ND::THandle<ND::TReconPID> fgdPid;
1871  ND::THandle<ND::TReconTrack> fgdTrack;
1872 
1873  if (useFGD1) {
1874  fgdCons =
1875  ReconObjectUtils::GetConstituentInDetector(reco, ND::TReconBase::kFGD1);
1876 
1877  if (fgdCons) {
1878  fgdPid = fgdCons;
1879  if (fgdPid) {
1880  ND::TReconObjectContainer::const_iterator it =
1881  fgdPid->GetConstituents()->begin();
1882  fgdTrack = *it;
1883  }
1884  }
1885 
1886  if (fgdTrack && fgdTrack->Has<ND::TIntegerDatum>("matchingFailure_flag"))
1887  result->matchingFailure_flag =
1888  fgdTrack->Get<ND::TIntegerDatum>("matchingFailure_flag")->GetValue();
1889  else
1890  result->matchingFailure_flag = 0;
1891  }
1892 
1893  else if (useFGD2) {
1894  fgdCons =
1895  ReconObjectUtils::GetConstituentInDetector(reco, ND::TReconBase::kFGD2);
1896 
1897  if (fgdCons) {
1898  fgdPid = fgdCons;
1899  if (fgdPid) {
1900  ND::TReconObjectContainer::const_iterator it =
1901  fgdPid->GetConstituents()->begin();
1902  fgdTrack = *it;
1903  }
1904  }
1905 
1906  if (fgdTrack && fgdTrack->Has<ND::TIntegerDatum>("matchingFailure_flag"))
1907  result->matchingFailure_flag =
1908  fgdTrack->Get<ND::TIntegerDatum>("matchingFailure_flag")->GetValue();
1909  else
1910  result->matchingFailure_flag = 0;
1911  }
1912 
1913  else
1914  result->matchingFailure_flag = 0;
1915 
1916  return;
1917 }
ND::TTrackerReconModule::TTPCTrack::ExtrapolatedVertexZY
double ExtrapolatedVertexZY
for ybar vertex, this is z coordinate in mm
Definition: TTrackerReconModule.hxx:155
ND::TTrackerReconModule::fVertices
TClonesArray * fVertices
The vector of trackerRecon vertices (none ever?).
Definition: TTrackerReconModule.hxx:507
ND::TTrackerReconModule::fTPCUnused
TClonesArray * fTPCUnused
The vector of unused TPC hits.
Definition: TTrackerReconModule.hxx:524
ND::TTrackerReconModule::ExtrapolateTPC
bool ExtrapolateTPC(ND::THandle< ND::TReconPID > in, TClonesArray *out, int idx, ND::THandle< ND::THitSelection > fgd)
Method to make a TPC object with extrapolated hit association information.
Definition: TTrackerReconModule.cxx:529
ND::TTrackerReconModule::TTrackerResult::TrueParticle
TTrueParticle TrueParticle
information about the true particleassociated with this track
Definition: TTrackerReconModule.hxx:316
ND::TTrackerReconModule::TTrackerConstituent::EDeposit
double EDeposit
The deposited charge for the constituent object(number of pe&#39;s).
Definition: TTrackerReconModule.hxx:240
ND::TTrackerReconModule::TTrackerResult
An object to describe a reconstructed trackerRecon result.
Definition: TTrackerReconModule.hxx:275
ND::TTrackerReconModule::TTrackerResult::AlgorithmName
std::string AlgorithmName
The name of the algorithm that created this object.
Definition: TTrackerReconModule.hxx:281
ND::TTrackerReconModule::fNTPCExtrapolation
Int_t fNTPCExtrapolation
The number of TPC tracks extrapolated into theFGD following Claudio&#39;s (2010a) method.
Definition: TTrackerReconModule.hxx:529
ND::TTrackerReconModule::TTrackerResult::Length
double Length
The total length of the object in mm.
Definition: TTrackerReconModule.hxx:296
ND::TTrackerReconModule::TTrackOther::Detector
int Detector
Detector used (1,2,3 for TPC, or 4,5 for FGD?)
Definition: TTrackerReconModule.hxx:201
ND::TTrackerReconModule::TTPCTrack::ExtrapolatedVertexXX
double ExtrapolatedVertexXX
for xbar vertex, this is x coordinate in mm
Definition: TTrackerReconModule.hxx:149
ND::TTrackerReconModule::GetDetectorNumber
int GetDetectorNumber(ND::THandle< ND::TReconBase > object)
Method to build an integer key that describes the detectors used to build the TRecon object...
Definition: TTrackerReconModule.cxx:1370
CVSTAG
#define CVSTAG
Definition: TTrackerReconModule.cxx:87
ND::TTrackerReconModule::TFGDTrack::Ndof
int Ndof
Number of degrees of freedom in FGD track fit.
Definition: TTrackerReconModule.hxx:168
ND::TTrackerReconModule::TFGDTrack::Direction
TVector3 Direction
the FGD track direction vector (dx,dy,dz)
Definition: TTrackerReconModule.hxx:172
ND::TTrackerReconModule::TTrackerConstituent::BackVariance
TLorentzVector BackVariance
The 4-vector position variance at theback of the track(var(x),var(y),var(z),var(t)) in mm^2...
Definition: TTrackerReconModule.hxx:244
ND::TTrackerReconModule::TTrackerResult::NTotalConstituents
Int_t NTotalConstituents
Number of all constituents, andconstituents-constituents...
Definition: TTrackerReconModule.hxx:290
ND::TTrackerReconModule::TTPCTrack::PullPion
double PullPion
Pull for TPC pid pion hypothesis.
Definition: TTrackerReconModule.hxx:107
ND::TTrackerReconModule::fCathodeOffset
double fCathodeOffset
TPC cathode offset used for tpc other hits.
Definition: TTrackerReconModule.hxx:534
ND::TTrackerReconModule::TTrackerResult::UniqueID
UInt_t UniqueID
Unique ID number to allow matching to Global Reconobject.
Definition: TTrackerReconModule.hxx:279
ND::TTrackerReconModule::TTrackerConstituent::Quality
double Quality
The quality of the fit.(probability)
Definition: TTrackerReconModule.hxx:231
ND::TTrackerReconModule::TTrackerConstituent::FrontPosition
TLorentzVector FrontPosition
The 4-vector position at the front ofthe track (x,y,z,t) in mm, ns.
Definition: TTrackerReconModule.hxx:241
ND::TTrackerReconModule::FillConfigTree
void FillConfigTree(TTree *configTree)
Definition: TTrackerReconModule.cxx:194
ND::TTrackerReconModule::TTPCTrack::TrCurvature
double TrCurvature
track curvature, units are 1/mm
Definition: TTrackerReconModule.hxx:142
ND::TTrackerReconModule::TUnusedHit
An object to describe a reconstructed track.
Definition: TTrackerReconModule.hxx:38
ND::TTrackerReconModule::TTPCTrack::UniqueID
UInt_t UniqueID
Unique ID number to allow matching to Global Reconobject.
Definition: TTrackerReconModule.hxx:86
ND::TTrackerReconModule::GetG4Vertex
bool GetG4Vertex(ND::TG4Trajectory *G4track, ND::TG4PrimaryVertex &G4vertex)
Method to find the best matching G4 vertex for a given G4 track.
Definition: TTrackerReconModule.cxx:1785
ND::TTrackerReconModule::TTrackerConstituent::MomentumError
double MomentumError
uncertainty in momentum MeV/c
Definition: TTrackerReconModule.hxx:255
ND::TTrackerReconModule::TTPCTrack::PullMuon
double PullMuon
Pull for TPC pid muon hypothesis.
Definition: TTrackerReconModule.hxx:106
ND::TTrackerReconModule::TTrackerResult::FGD
TClonesArray * FGD
FGD objects associated with track.
Definition: TTrackerReconModule.hxx:320
ND::TTrackerReconModule::TTrackerVertex::Quality
double Quality
The quality of the fit. Ie. the Prob(chi2,ndof)
Definition: TTrackerReconModule.hxx:65
ND::TTrackerReconModule::TTrackerResult::NDOF
int NDOF
The number of degrees of freedom.
Definition: TTrackerReconModule.hxx:285
ClassImp
ClassImp(ND::TBeamSummaryDataModule::TBeamSummaryData)
ND::TTrackerReconModule::TTrackerNode::Position
TLorentzVector Position
Position 4-vector (at node) x,y,z,t in mm, ns.
Definition: TTrackerReconModule.hxx:265
ND::TTrackerReconModule::fTPCIso
TClonesArray * fTPCIso
The vector of isolated TPC tracks with fits.
Definition: TTrackerReconModule.hxx:521
ND::TTrackerReconModule::TTrackerResult::NFGDs
Int_t NFGDs
Number of FGD Specific objects.
Definition: TTrackerReconModule.hxx:319
ND::TTrackerReconModule::InitializeBranches
virtual void InitializeBranches()
Initialize the branches in the trackerRecon eventAnalysis tree.
Definition: TTrackerReconModule.cxx:150
ND::TAnalysisModuleBase::fDescription
std::string fDescription
A longish descrition of the analysis.
Definition: TAnalysisModuleBase.hxx:188
insideFGD1
bool insideFGD1(double zx, double zy, double xx, double yy)
Definition: TTrackerReconModule.cxx:466
ND::TTrackerReconModule::FillFinalTrack
void FillFinalTrack(TTrackerResult *result, ND::THandle< ND::TReconTrack > reco)
Method to fill a trackerRecon result entry.
Definition: TTrackerReconModule.cxx:1250
ND::TTrueVertex::Position
TLorentzVector Position
The position of the vertex.
Definition: TTrueVertex.hxx:18
ND::TTrackerReconModule::TTrackerVertex
An object to describe a reconstructed vertex.
Definition: TTrackerReconModule.hxx:56
ND::TTrackerReconModule::GetConstituentTrack
ND::THandle< ND::TReconTrack > GetConstituentTrack(ND::THandle< ND::TReconPID > in, int &nconstit)
called recursively to look for bottom level track that built a PID.
Definition: TTrackerReconModule.cxx:644
ND::TTrackerReconModule::GetParent
ND::THandle< ND::TG4Trajectory > GetParent(ND::TG4Trajectory *G4track)
Method to get the parent G4 track for a given G4 track.
Definition: TTrackerReconModule.cxx:1736
ND::TTrackerReconModule::TTrackerResult::NNodes
int NNodes
The number of nodes (fgd hits + tpc tracks)
Definition: TTrackerReconModule.hxx:321
ND::TTrackerReconModule::TTPCTrack::Charge
double Charge
Charge from the TPC pid (+1, or -1)
Definition: TTrackerReconModule.hxx:94
ND::TTrackerReconModule::fNFGDUnused
Int_t fNFGDUnused
The number of unused FGD hits.
Definition: TTrackerReconModule.hxx:526
ND::TTrackerReconModule::FillTrueVertex
void FillTrueVertex(bool found, const ND::TG4PrimaryVertex &G4vertex, double pur, double eff, ND::TTrueVertex &vertex)
Method to fill the true vertex information from a given G4 vertex.
Definition: TTrackerReconModule.cxx:1770
ND::TTrackerReconModule::FillFGD
void FillFGD(ND::THandle< ND::TReconTrack > in, TClonesArray *out, int idx)
Method to fill trackerRecon tree with FGD information Input is a TReconPID of FGD fit information...
Definition: TTrackerReconModule.cxx:232
ND::TTrackerReconModule::FillVertex
void FillVertex(ND::THandle< ND::TReconVertex > in, TClonesArray *out, int idx)
Method to fill trackerRecon tree with vertex information Input is a TReconVertex. ...
Definition: TTrackerReconModule.cxx:209
ND::TTrackerReconModule::TTPCTrack::ExtrapolatedVertexZX
double ExtrapolatedVertexZX
for xbar vertex, this is z coordinate in mm
Definition: TTrackerReconModule.hxx:151
ND::TTrackerReconModule::fMaxDrift
double fMaxDrift
TPC maximum drift used for tpc other hits.
Definition: TTrackerReconModule.hxx:533
ND::TTrackerReconModule::TTPCTrack::SigmaPion
double SigmaPion
Sigma estimated width of TPC pid pion hypothesis.
Definition: TTrackerReconModule.hxx:117
ND::TTrackerReconModule::TTrackerConstituent::AlgorithmName
std::string AlgorithmName
algorithm that created this object.
Definition: TTrackerReconModule.hxx:228
ND::TTrackerReconModule::fNFGDOther
Int_t fNFGDOther
The number of FGD tracks with no fit (none ever?).
Definition: TTrackerReconModule.hxx:512
ND::TTrackerReconModule::TTrackerConstituent::Position
TLorentzVector Position
position 4-vector (x,y,z,t) in mm, ns
Definition: TTrackerReconModule.hxx:250
ND::TTrackerReconModule::TTPCTrack::MomentumError
double MomentumError
Uncertainty in the Momentum in MeV/c from theTPC pid.
Definition: TTrackerReconModule.hxx:93
ND::TTrackerReconModule::fNTPCOther
Int_t fNTPCOther
The number of TPC tracks with no fit.
Definition: TTrackerReconModule.hxx:515
borderFGD1
bool borderFGD1(TVector3 pos)
Definition: TTrackerReconModule.cxx:476
ND::TTrackerReconModule::TFGDTrack
An object to hold specific FGD track variables.
Definition: TTrackerReconModule.hxx:163
ND::TTrackerReconModule::TTrackOther::EDeposit
double EDeposit
The deposited charge for the object.
Definition: TTrackerReconModule.hxx:207
ND::TTrackerReconModule::TTrackerConstituent::FrontVariance
TLorentzVector FrontVariance
The 4-vector position variance at thefront of the track(var(x),var(y),var(z),var(t)) in mm^2...
Definition: TTrackerReconModule.hxx:243
ND::TTrackerReconModule::TTPCTrack::HasExtrapolation
bool HasExtrapolation
extrapolation method of vertex is calculated ornot
Definition: TTrackerReconModule.hxx:147
ND::TTrackerReconModule::fFGDUnused
TClonesArray * fFGDUnused
The vector of unused TPC hits.
Definition: TTrackerReconModule.hxx:527
ND::TTrackerReconModule::FillTrueParticle
void FillTrueParticle(ND::TG4Trajectory *G4track, double pur, double eff, ND::TTrackerReconModule::TTrueParticle &trueParticle)
Method to fill the true particle information given a G4track.
Definition: TTrackerReconModule.cxx:1711
ND::TTrackerReconModule::ProcessFirstEvent
virtual Bool_t ProcessFirstEvent(ND::TND280Event &event)
Method for special handling of first event.
Definition: TTrackerReconModule.cxx:129
ND::TTrackerReconModule::TTrueParticle
An object to describe the true G4 particle associated to the TTrackerResult.
Definition: TTrackerReconModule.hxx:181
ND::TTrackerReconModule::TTrackerVertex::Position
TLorentzVector Position
The position of the vertex components [0]=x [1]=y [2]=z [3]=t x,y,z are in mm, t is in ns...
Definition: TTrackerReconModule.hxx:72
ND::TTrackerReconModule::TTrackerNode::Variance
TLorentzVector Variance
Position variance 4-vector (at node)var(x),var(y),var(z),var(t) in mm^2, ns^2.
Definition: TTrackerReconModule.hxx:266
ND::TTrackerReconModule::TTPCTrack::PositionVariance
TLorentzVector PositionVariance
Variance in Position in mm^2, ns^2.
Definition: TTrackerReconModule.hxx:97
ND::TTrackerReconModule::TUnusedHit::Variance
TVector3 Variance
The position variance in mm.
Definition: TTrackerReconModule.hxx:46
ND::TTrackerReconModule::TTrackerResult::Pids
std::vector< int > Pids
the PID that goes with Likelihoods
Definition: TTrackerReconModule.hxx:309
ND::TTrackerReconModule::TTPCTrack::NConstituents
int NConstituents
Number of constituents (probably 1, for track!)
Definition: TTrackerReconModule.hxx:136
ND::TTrackerReconModule::TTrackerConstituent::ConstitIdx
Int_t ConstitIdx[2]
Index into Constituents inTTrackerResult::Constituents of this constituent&#39;sconstituents.
Definition: TTrackerReconModule.hxx:237
ND::TTrackerReconModule::TTrackerNode::Direction
TVector3 Direction
Direction vector (at node)
Definition: TTrackerReconModule.hxx:267
ND::TTrackerReconModule::TTPCTrack::Direction
TVector3 Direction
TPC pid direction vector in mm.
Definition: TTrackerReconModule.hxx:99
ND::TTrackerReconModule::TTrueParticle::Pur
double Pur
The purity for matching.
Definition: TTrackerReconModule.hxx:185
ND::TTrackerReconModule::FillConstituentTrack
bool FillConstituentTrack(TTrackerResult *result, ND::THandle< ND::TReconTrack > in, int &idx)
Method to fill information about constituent tracks.
Definition: TTrackerReconModule.cxx:940
ND::TTrackerReconModule::TTrackerResult::matchingFailure_flag
int matchingFailure_flag
Flag a object where the TPC-FGD matchingfailed.
Definition: TTrackerReconModule.hxx:297
ND::TTrackerReconModule::TTrackerConstituent::Charge
double Charge
The Charge of this constituent (+-1)
Definition: TTrackerReconModule.hxx:239
ND::TTrackerReconModule::TTPCTrack::Position
TLorentzVector Position
Position at which kinematics are reported in mm, ns.
Definition: TTrackerReconModule.hxx:96
ND::TAnalysisModuleBase::fBufferSize
Int_t fBufferSize
Buffer Size for TBranch.
Definition: TAnalysisModuleBase.hxx:184
ND::TTrackerReconModule::TTrackerConstituent::DirectionVariance
TVector3 DirectionVariance
direction variance vector
Definition: TTrackerReconModule.hxx:253
ND::TTrackerReconModule::TTPCTrack::DirectionVariance
TVector3 DirectionVariance
TPC pid variance in vector direction in mm^2.
Definition: TTrackerReconModule.hxx:101
ND::TTrackerReconModule::TFGDTrack::EDeposit
double EDeposit
the FGD track energy deposit (total number of pe&#39;s)
Definition: TTrackerReconModule.hxx:174
ND::TTrackerReconModule::TTrackerConstituent::Status
unsigned long Status
The status for the fit.
Definition: TTrackerReconModule.hxx:230
ND::TTrackerReconModule::TTrackerConstituent::FrontMomentum
double FrontMomentum
the momentum at the front of the track in MeV/c
Definition: TTrackerReconModule.hxx:248
ND::TTrackerReconModule::TTrackerResult::Likelihoods
std::vector< double > Likelihoods
the PID likelihoods for combined PID for now this is the product of the three TPC PID probabilities T...
Definition: TTrackerReconModule.hxx:308
ND::TTrackerReconModule::TTrackerVertex::NDOF
int NDOF
The number of degrees of freedom.
Definition: TTrackerReconModule.hxx:67
ND::TTrackerReconModule::fNVertices
Int_t fNVertices
Tree Entries.
Definition: TTrackerReconModule.hxx:505
ND::TTrackerReconModule::TTPCTrack::Sigma0
double Sigma0
TPC track diffusion parameters depend on diffusion model used Original TPC diffusion model was: sigma...
Definition: TTrackerReconModule.hxx:130
ND::TTrackerReconModule::GetNumberOfHits
int GetNumberOfHits(ND::THandle< ND::TReconBase > object)
Method to find the total number of hits used to build the reconstructed track/pid/etc.
Definition: TTrackerReconModule.cxx:1840
ND::TTrackerReconModule::TTrackerNode::MomentumError
double MomentumError
Track Momentum uncertainty (at node) in MeV/c.
Definition: TTrackerReconModule.hxx:270
ND::TTrackerReconModule::TTPCTrack::SigmaMuon
double SigmaMuon
Sigma estimated width of TPC pid muon hypothesis.
Definition: TTrackerReconModule.hxx:116
ND::TTrackerReconModule::TTPCTrack::dEdxexpKaon
double dEdxexpKaon
Estimated dE/dx for kaon hypothesis.
Definition: TTrackerReconModule.hxx:113
ND::TTrackerReconModule::TTrackerConstituent::BackDirection
TVector3 BackDirection
The direction vector at the back of the track.
Definition: TTrackerReconModule.hxx:247
ND::TTrackerReconModule::FillOther
void FillOther(ND::THandle< ND::TReconTrack > in, TClonesArray *out, int idx)
Method to fill trackerRecon tree with unfit track information Input is a TReconTrack that was not fit...
Definition: TTrackerReconModule.cxx:733
ND::TTrackerReconModule::TTrackerVertex::AlgorithmName
std::string AlgorithmName
The name of the algorithm that created this object.
Definition: TTrackerReconModule.hxx:61
computeDistFGD
double computeDistFGD(TVector3 f, TVector3 in, double tx, double ty, double curv, bool isX)
Definition: TTrackerReconModule.cxx:497
ND::TTrackerReconModule::TTrackerVertex::Variance
TLorentzVector Variance
The position of the vertex.
Definition: TTrackerReconModule.hxx:77
ND::TTrackerReconModule::InitializeModule
virtual void InitializeModule()
Method to initialize this module.
Definition: TTrackerReconModule.cxx:143
ND::TTrackerReconModule::TTrackerResult::Constituents
TClonesArray * Constituents
All constituents, and constituents-constituents...
Definition: TTrackerReconModule.hxx:292
ND::TTrackerReconModule::TTrackerReconModule
TTrackerReconModule(const char *name="Tracker", const char *title="Tracker Recon Module")
Constructor.
Definition: TTrackerReconModule.cxx:95
ND::TAnalysisModuleBase::fCVSID
std::string fCVSID
Defined if an official tagged version.
Definition: TAnalysisModuleBase.hxx:192
ND::TTrackerReconModule::~TTrackerReconModule
virtual ~TTrackerReconModule()
Destructor.
Definition: TTrackerReconModule.cxx:125
ND::TTrackerReconModule::TTPCTrack::Ccorr
double Ccorr
Corrected truncated mean charge deposit in PID.
Definition: TTrackerReconModule.hxx:104
ND::TTrackerReconModule::FillTree
virtual bool FillTree(ND::TND280Event &)
Main call to fill all of the branches of the trackerRecon eventAnalysis tree for a given event...
Definition: TTrackerReconModule.cxx:1392
ND::TTrackerReconModule::TFGDTrack::Chi2
double Chi2
Chi2 of the FGD track fit the FGD track starting position vector coordinates in order x...
Definition: TTrackerReconModule.hxx:169
ND::TTrackerReconModule::TTrackerNode
Object to describe the fit parameters at a track node.
Definition: TTrackerReconModule.hxx:260
ND::TTrackerReconModule::TTrackerNode::EDeposit
double EDeposit
The Energy Deposit (number of pe&#39;s)
Definition: TTrackerReconModule.hxx:264
ND::TTrackerReconModule::TTrackerResult::DirectionVariance
TVector3 DirectionVariance
track direction variance
Definition: TTrackerReconModule.hxx:313
ND::TTrackerReconModule::TTrackerResult::EDeposit
double EDeposit
The deposited charge for the object (number of pe&#39;s)
Definition: TTrackerReconModule.hxx:295
ND::TTrackerReconModule::TTPCTrack::SigmaKaon
double SigmaKaon
Sigma estimated width of TPC pid kaon hypothesis.
Definition: TTrackerReconModule.hxx:118
ND::TTrackerReconModule::GetTrajectory
ND::TG4Trajectory GetTrajectory(ND::TND280Event &event, ND::THandle< ND::TReconBase > reco)
Method to get the G4 trajectory that best matches the TRecon object.
Definition: TTrackerReconModule.cxx:1625
ND::TTrackerReconModule::TTrackerResult::Status
unsigned long Status
The status for the fit.
Definition: TTrackerReconModule.hxx:283
ND::TTrackerReconModule::TTPCTrack::dEdxexpMuon
double dEdxexpMuon
Estimated dE/dx for muon hypothesis.
Definition: TTrackerReconModule.hxx:111
ND::TTrackerReconModule::CheckMatchingFailure
void CheckMatchingFailure(TTrackerResult *result, ND::THandle< ND::TReconBase > reco)
Method to find FGD constituent track from a matching object and fill matchingFailure_flag.
Definition: TTrackerReconModule.cxx:1862
ND::TTrackerReconModule::TTPCTrack::PullEle
double PullEle
Pull for TPC pid electron hypothesis.
Definition: TTrackerReconModule.hxx:105
ND::TTrackerReconModule::TTrackerResult::Quality
double Quality
The quality of the fit.(probability(chi2,ndof))
Definition: TTrackerReconModule.hxx:284
ND::TAnalysisModuleBase::SetNameTitle
void SetNameTitle(char const *name, char const *title)
Definition: TAnalysisModuleBase.cxx:108
ND::TTrackerReconModule::TTrueParticle::Vertex
ND::TTrueVertex Vertex
True vertex associated to this TrueParticle.
Definition: TTrackerReconModule.hxx:187
ND::TTrackerReconModule::fFGDOther
TClonesArray * fFGDOther
The vector of FGD tracks with no fit.
Definition: TTrackerReconModule.hxx:513
ND::TTrackerReconModule::TTPCTrack::Momentum
double Momentum
Momentum of the TPC pid in MeV/c.
Definition: TTrackerReconModule.hxx:92
ND::TTrackerReconModule::GetIncomingParticle
bool GetIncomingParticle(const ND::TG4PrimaryVertex &G4vertex, ND::TG4PrimaryParticle &incoming)
Method to find the incident particle (neutrino) that gave rise to a given G4 track.
Definition: TTrackerReconModule.cxx:1813
ND::TTrackerReconModule::TTPCTrack::NHits
int NHits
number of clusters used in TPC fit
Definition: TTrackerReconModule.hxx:90
ND::TTrackerReconModule::TTrackerConstituent::BackPosition
TLorentzVector BackPosition
The 4-vector position of the back of thetrack (x,y,z,t) in mm, ns.
Definition: TTrackerReconModule.hxx:242
ND::TTrackerReconModule::TUnusedHit::Time
double Time
Time of the hit in ns.
Definition: TTrackerReconModule.hxx:48
ND::TTrackerReconModule::TTrackOther::AlgorithmName
std::string AlgorithmName
The name of the algorithm that created this object.
Definition: TTrackerReconModule.hxx:199
ND::TTrackerReconModule::TTPCTrack::SigmaEle
double SigmaEle
Sigma estimated width of TPC pid electron hypothesis.
Definition: TTrackerReconModule.hxx:115
ND::TTrackerReconModule::TTrackOther::BackPosition
TLorentzVector BackPosition
The position of the track at itsdownstream-most end (x,y,z,t) in mm,ns.
Definition: TTrackerReconModule.hxx:211
ND::TTrackerReconModule::TTrackerResult::Position
TLorentzVector Position
track position 4-vector (x,y,z,t) in mm, ns
Definition: TTrackerReconModule.hxx:310
ND::TTrackerReconModule::TTrackerNode::DirectionVariance
TVector3 DirectionVariance
Direction variance vector (at node)
Definition: TTrackerReconModule.hxx:268
ND::TTrackerReconModule::TTrackerConstituent::NConstituents
Int_t NConstituents
The number of constituents this constituent ismade of.
Definition: TTrackerReconModule.hxx:236
ND::TTrackerReconModule::FillConstituentPid
bool FillConstituentPid(TTrackerResult *result, ND::THandle< ND::TReconPID > in, int &idx)
Method to fill information about constituent pids.
Definition: TTrackerReconModule.cxx:813
ND::TTrackerReconModule::TTPCTrack::TrDirectionVar
TVector3 TrDirectionVar
variance in track direction vector
Definition: TTrackerReconModule.hxx:141
ND::TTrackerReconModule::TTPCTrack::Ndof
int Ndof
Number of degrees of freedom in TPC fit.
Definition: TTrackerReconModule.hxx:88
ND::TTrackerReconModule::TTPCTrack
TPC pid and track variables.
Definition: TTrackerReconModule.hxx:83
ND::TTrackerReconModule::TTPCTrack::ExtrapolatedVertexYY
double ExtrapolatedVertexYY
for ybar vertex, this is y coordinate in mm
Definition: TTrackerReconModule.hxx:153
CVSID
#define CVSID
Definition: TTrackerReconModule.cxx:90
ND::TTrackerReconModule::TTPCTrack::Chi2
double Chi2
TPC chi2 calculated after likelihood fit.
Definition: TTrackerReconModule.hxx:89
getTangents
void getTangents(ND::THandle< ND::TTrackState > state, double &tx, double &etx, double &ty, double &ety)
Definition: TTrackerReconModule.cxx:436
ND::TTrackerReconModule::TTPCTrack::dEdxexpPion
double dEdxexpPion
Estimated dE/dx for pion hypothesis.
Definition: TTrackerReconModule.hxx:112
ND::TTrackerReconModule::TTPCTrack::Sigma1
double Sigma1
TPC track diffusion sigma1 parameter.
Definition: TTrackerReconModule.hxx:131
ND::TTrackerReconModule::fTPCOther
TClonesArray * fTPCOther
The vector of TPC tracks with no fit.
Definition: TTrackerReconModule.hxx:516
ND::TAnalysisModuleBase::fSplitLevel
Int_t fSplitLevel
Split Level for TBranch.
Definition: TAnalysisModuleBase.hxx:186
ND::TTrackerReconModule::TTrackerConstituent::NDOF
int NDOF
The number of degrees of freedom.
Definition: TTrackerReconModule.hxx:232
ND::TTrackerReconModule::TTrackerConstituent::BackMomentum
double BackMomentum
the momentum at the back of the track iin MeV/c
Definition: TTrackerReconModule.hxx:249
ND::TTrueVertex::ID
Int_t ID
The vertex ID from G4.
Definition: TTrueVertex.hxx:21
ND::TTrackerReconModule::TUnusedHit::Position
TVector3 Position
The position of the hit component 0=x 1=y 2=z in mm.
Definition: TTrackerReconModule.hxx:44
ND::TTrackerReconModule::TTrackerResult::NTPCs
Int_t NTPCs
Number of TPC tracks used to build this track.
Definition: TTrackerReconModule.hxx:317
ND::TTrackerReconModule::TUnusedHit::TimeUnc
double TimeUnc
Time Uncertainty of hit in ns.
Definition: TTrackerReconModule.hxx:50
borderFGD2
bool borderFGD2(TVector3 pos)
Definition: TTrackerReconModule.cxx:486
ND::TTrackerReconModule::TTPCTrack::dEdxexpEle
double dEdxexpEle
Estimated dE/dx for electron hypothesis.
Definition: TTrackerReconModule.hxx:110
ND::TTrackerReconModule::TTrackOther
An object to describe a track that has no fit result.
Definition: TTrackerReconModule.hxx:194
ND::TTrackerReconModule::fNTPCIso
Int_t fNTPCIso
The number of isolated TPC tracks with fits (none ever?)
Definition: TTrackerReconModule.hxx:520
ND::TAnalysisModuleBase::fCVSTagName
std::string fCVSTagName
Defined if an official tagged version.
Definition: TAnalysisModuleBase.hxx:190
ND::TTrackerReconModule::fNTPCUnused
Int_t fNTPCUnused
The number of unused TPC hits.
Definition: TTrackerReconModule.hxx:523
ND::TTrueVertex
An object to describe the true G4 vertex associated to the TGlobalVertex.
Definition: TTrueVertex.hxx:14
ND::TTrackerReconModule::TFGDTrack::UniqueID
UInt_t UniqueID
Unique ID number to allow matching to Global Reconobject.
Definition: TTrackerReconModule.hxx:166
ND::TTrackerReconModule::fTracks
TClonesArray * fTracks
The vector of overall trackerRecon results.
Definition: TTrackerReconModule.hxx:510
ND::TTrackerReconModule::TTrackerResult::NConstituents
Int_t NConstituents
The number of constituents (tracks and pids) usedto build this track.
Definition: TTrackerReconModule.hxx:288
insideFGD2
bool insideFGD2(double zx, double zy, double xx, double yy)
Definition: TTrackerReconModule.cxx:471
ND::TTrackerReconModule::TTPCTrack::PullKaon
double PullKaon
Pull for TPC pid kaon hypothesis.
Definition: TTrackerReconModule.hxx:108
ND::TTrackerReconModule::TTPCTrack::TrCurvatureVar
double TrCurvatureVar
variance in track direction vector, units are(1mm)^2
Definition: TTrackerReconModule.hxx:143
ND::TTrackerReconModule::TTrueParticle::Eff
double Eff
The efficiency for matching.
Definition: TTrackerReconModule.hxx:186
ND::TTrackerReconModule::TTrackerResult::ConstitIdx
Int_t ConstitIdx[2]
Index into Constituents of the constituents usedto build this track.
Definition: TTrackerReconModule.hxx:289
ND::TTrackerReconModule::TTrackOther::FrontPosition
TLorentzVector FrontPosition
The position of the track at itsupstream-most end (x,y,z,t) in mm, ns.
Definition: TTrackerReconModule.hxx:209
ND::TTrackerReconModule::TTPCTrack::SigmaProton
double SigmaProton
Sigma estimated width of TPC pid proton hypothesis.
Definition: TTrackerReconModule.hxx:119
ND::TTrackerReconModule::TTrackerResult::TPC
TClonesArray * TPC
Information about the TPC pids/tracks used to buildthis track.
Definition: TTrackerReconModule.hxx:318
ND::TTrackerReconModule::fTPCExtrapolation
TClonesArray * fTPCExtrapolation
The vector of TPC tracks extrapolatedinto the FGD following Claudio&#39;s (2010a)method.
Definition: TTrackerReconModule.hxx:530
ND::TTrackerReconModule::TTPCTrack::dEdxexpProton
double dEdxexpProton
Estimated dE/dx for proton hypothesis.
Definition: TTrackerReconModule.hxx:114
ND::TTrackerReconModule::TTrackerConstituent
An object to describe a constituent of a trackerRecon result.
Definition: TTrackerReconModule.hxx:224
ND::TAnalysisModuleBase::SetSplitLevel
void SetSplitLevel(Int_t splitlevel)
ROOT output parameters, usually no need to touch.
Definition: TAnalysisModuleBase.hxx:132
ND::TTrackerReconModule::GetMomentum
void GetMomentum(ND::THandle< ND::TTrackState > state, double &p, double &ep)
Helper method to convert TTrackState curvature information into momentum and momentum uncertainty...
Definition: TTrackerReconModule.cxx:1197
ND::TTrackerReconModule::fDriftVelocity
double fDriftVelocity
TPC drift velocity used for tpc other hits.
Definition: TTrackerReconModule.hxx:535
ND::TTrackerReconModule::TTrackerResult::Nodes
TClonesArray * Nodes
Kinematics of the track at each node in the track fit.
Definition: TTrackerReconModule.hxx:323
ND::TTrackerReconModule::GetCharge
double GetCharge(ND::THandle< ND::TTrackState > state)
Helper method to convert TTrackState curvature information into track charge.
Definition: TTrackerReconModule.cxx:1242
ND::TTrackerReconModule::TTPCTrack::Sigma2
double Sigma2
TPC track diffusion sigma2 parameter.
Definition: TTrackerReconModule.hxx:132
ND::TTrackerReconModule::FillFinalPid
void FillFinalPid(TTrackerResult *result, ND::THandle< ND::TReconPID > reco)
Method to fill a trackerRecon result entry.
Definition: TTrackerReconModule.cxx:1070
ND::TTrackerReconModule::TTrackerResult::Direction
TVector3 Direction
track direction vector
Definition: TTrackerReconModule.hxx:312
ND::TTrackerReconModule::TTPCTrack::PullProton
double PullProton
Pull for TPC pid proton hypothesis.
Definition: TTrackerReconModule.hxx:109
ND::TTrackerReconModule::TTPCTrack::NTrun
double NTrun
70% of the number of clusters
Definition: TTrackerReconModule.hxx:103
ND::TTrackerReconModule::fNTracks
Int_t fNTracks
The number of trackerRecon results.
Definition: TTrackerReconModule.hxx:509
ND::TTrackerReconModule::TTrackerResult::Variance
TLorentzVector Variance
track position variance 4-vectorvar(x),var(y),var(z),var(t) in mm^2, ns^2
Definition: TTrackerReconModule.hxx:311
ND::TTrackerReconModule::GetHitsFromComboHits
ND::THandle< ND::THitSelection > GetHitsFromComboHits(ND::THandle< ND::THitSelection > tpc)
Method to get the tpc combo hits from the tpc hit selection (used for tpc other)
Definition: TTrackerReconModule.cxx:710
ND::TTrackerReconModule::TTPCTrack::TrDirection
TVector3 TrDirection
track direction vector
Definition: TTrackerReconModule.hxx:140
ND::TTrackerReconModule::TTrackerConstituent::Variance
TLorentzVector Variance
position variance 4-vector(var(x),var(y),var(z),var(t)) in mm^2, ns^2
Definition: TTrackerReconModule.hxx:251
ND::TTrackerReconModule::TTrackerConstituent::FrontDirection
TVector3 FrontDirection
The direction vector at the front of the track.
Definition: TTrackerReconModule.hxx:246
ND::TTrackerReconModule::TTPCTrack::MeanDrift
double MeanDrift
TPC track mean drift value used in diffusion model.
Definition: TTrackerReconModule.hxx:133
ND::TTrackerReconModule::TUnusedHit::TotalCharge
double TotalCharge
Deposited charge (the hit EDeposit)
Definition: TTrackerReconModule.hxx:42
ND::TTrackerReconModule::TTrackerNode::Momentum
double Momentum
Track Momentum (at node) in MeV/c.
Definition: TTrackerReconModule.hxx:269
ND::TTrackerReconModule::FillTPC
void FillTPC(ND::THandle< ND::TReconPID > in, TClonesArray *out, int idx)
Method to fill trackerRecon tree with TPC information Input is a TReconPID of TPC fit information...
Definition: TTrackerReconModule.cxx:263
ND::TAnalysisModuleBase::FillConfigTree
virtual void FillConfigTree(TTree *)
Definition: TAnalysisModuleBase.cxx:52
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