WCTrackAlg.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////////
2 /// \file WCTrackAlg.cxx
3 /// \brief Algorithm implementation for wire chamber reconstruction.
4 /// Part of the beamline reconstruction for NOvA test beam.
5 /// \author Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
6 /// Fan Gao(University of Pittsburgh) <fag16@pitt.edu>
7 /// Based on LArIAT code lariatsoft/LArIATRecoAlg/WCTrackBuilderAlg.cxx
8 /// by Greg Pulliam gkpullia@syr.edu
9 /// \date December 2018
10 /////////////////////////////////////////////////////////////////////////////////
11 
12 // nova
14 
17 
18 namespace beamlinereco {
19 
21  this->reconfigure(pset);
22  }
23 
25 
26  // ------------------------------------------------------------------------------
28 
29  // Just a number to use to initialize things before they get filled correctly.
30  fInitialDefault = -999;
31 
32  // These four are parameters from histos I produced from picky-good tracks
33  //fCentralYKink = pset.get<float >("CentralYKink", -0.01 );
34  //fSigmaYKink = pset.get<float >("SigmaYKink", 0.03 );
35  //fCentralYDist = pset.get<float >("CentralYDist", 0.69 );
36  //fSigmaYDist = pset.get<float >("SigmaYDist", 18.0 );
37 
38  // Reconstruction options
39  fVerbose = pset.get<bool>("Verbose", false);
40  fPickyTracks = pset.get<bool>("PickyTracks", false);
41 
42  // Get coordinate system
43  unsigned int coordSystem = pset.get<unsigned int>("CoordSystem", 0);
44  switch (coordSystem) {
45  case 0:
46  fBeamlineCoordSystem = BeamlineCoordSystem::Beamline;
47  break;
48  case 1:
49  fBeamlineCoordSystem = BeamlineCoordSystem::NOvADetector;
50  break;
51  default:
52  std::cout << "Unknown coordinate system (" << coordSystem << "). "
53  << "Options are 0: Beamline; 1: NOvA Detector." << std::endl;
54  abort();
55  }
56 
57  return;
58 
59  }
60 
61  // ------------------------------------------------------------------------------
62  /// Main function called for each event
63  void WCTrackAlg::reconstructTracks(std::vector<double> & reco_p_list,
64  std::vector<TVector3>& mag_int_list,
65  std::vector<double>& dist_to_mag_axis_list,
66  std::vector<TVector2>& face_list,
67  std::vector<std::vector<TVector3> > & wc_hit_pos_list,
68  std::vector<TVector3> & incoming_dir_list,
69  std::vector<WCHitList> & event_final_tracks,
70  std::vector<std::vector<WCHitList> > & good_hits,
71  std::vector<double> & y_kink_list,
72  std::vector<TVector3>& dist_list,
73  int & WCMissed,
74  //std::vector<TH2F*> & Recodiff,
75  TH1F* & WCdistribution,
76  float & residual) {
77 
78  WCMissed = fInitialDefault;
79 
80  // Determine if one should skip this trigger based on whether there is exactly one good hit in each wire chamber and axis
81  // If there isn't, continue after adding a new empty vector to the reco_p_array contianer.
82  // If there is, then move on with the function.
83  // This can be modified to permit more than one good hit in each wire chamber axis - see comments in function
84  bool skip = shouldSkipTrigger(good_hits, WCMissed, WCdistribution);
85  if (skip == true)
86  return;
87 
88  // Assign value to class member variable
89  fWCMissed = WCMissed;
90 
91  // Depending on if an event has a hit in all 4 WC or whether it missed WC2 or WC3 (but not both),
92  // we reconstruct the momentum differently.
93  // This code doesn't change from before we allowed 3 point tracks.
94  if (fNHits == 4) {
95 
96  // At this point, we should have a list of good hits with at least one good hit in X,Y for each WC.
97  // Now find all possible combinations of these hits that could form a track, sans consideration
98  // of the kinks or end displacements (for now). For the "exactly one" condition set in the above
99  // step, this won't matter, but if you want to set the condition to "at least one hit in each WC axis,
100  // this will give many combinations.
101  residual = buildFourPointTracks(good_hits,
102  reco_p_list,
103  mag_int_list,
104  dist_to_mag_axis_list,
105  face_list,
106  wc_hit_pos_list,
107  incoming_dir_list,
108  event_final_tracks,
109  y_kink_list,
110  dist_list,
111  WCMissed);
112 
113  }
114 
115  return;
116 
117  }
118 
119  // ------------------------------------------------------------------------------
120  bool WCTrackAlg::shouldSkipTrigger(std::vector<std::vector<WCHitList> > & good_hits,
121  int & WCMissed,
122  TH1F* & WCDist) {
123 
124  //Now determine if we want to skip
125  bool skip = false;
126  fNHits = 0;
127  for (size_t iWC = 0; iWC < 4; ++iWC) {
128  if (good_hits[iWC][0].hits.size() > 0 && good_hits[iWC][1].hits.size() > 0)
129  ++fNHits;
130  else
131  WCMissed = iWC+1;
132  }
133 
134  // With WCDist histogram, we can simply estimate the efficiency of each WC
135  WCDist->Fill(0);
136  if (fNHits > 2)
137  WCDist->Fill(1);
138  if (fNHits > 2 && WCMissed == 1)
139  WCDist->Fill(2);
140  if (fNHits > 2 && WCMissed == 2)
141  WCDist->Fill(3);
142  if (fNHits > 2 && WCMissed == 3)
143  WCDist->Fill(4);
144  if (fNHits > 2 && WCMissed == 4)
145  WCDist->Fill(5);
146  if (fNHits == 4)
147  WCDist->Fill(6);
148 
149  // If we don't have 3 or 4 hits, skip.
150  // For commissioning, only require 4 hits first
151  // if (fNHits < 3)
152  if (fNHits < 4)
153  skip = true;
154 
155  if (fPickyTracks) {
156  for (size_t iWC=0; iWC<4; ++iWC)
157  for (size_t iAX=0; iAX<2; ++iAX)
158  if (good_hits[iWC][iAX].hits.size() != 1) {
159  // Only allow events with 1 and only 1 hit on each axis
160  skip = true;
161  break;
162  }
163  }
164  else // Skip events with less than 3 X/Y hits or is missing the first or last WC
165  if (WCMissed==1 || WCMissed==4)
166  skip = true;
167 
168  if (skip) {
169  if (fVerbose)
170  std::cout << "Skipping this event." << std::endl;
171 
172  // Clear the hit lists for each WC/axis
173  for (size_t iWC = 0; iWC < good_hits.size(); ++iWC)
174  for (size_t iAx = 0; iAx < good_hits.at(iWC).size(); ++iAx)
175  good_hits.at(iWC).at(iAx).hits.clear();
176  }
177 
178  return skip;
179 
180  }
181 
182  // ------------------------------------------------------------------------------
183  float WCTrackAlg::calculateRecoP(float theta_x_us, float theta_x_ds, float bestTrackSlope) {
184  float magnet_angle = fGeo->MagnetAngle();
185  float magnet_eff_length = fGeo->MagnetEffectiveLength();
186  float alpha_1 = theta_x_us + magnet_angle*TMath::DegToRad();
187  float alpha_2 = theta_x_ds + magnet_angle*TMath::DegToRad();
188  float num = (fabs(fBField) * magnet_eff_length * 0.01 /*cm_to_m*/ * 1000. /*GeV_to_MeV*/ );
189  float denom = (3.33564*(sin(alpha_2) - sin(alpha_1)))*cos(atan(bestTrackSlope));
190  return num / denom;
191  }
192 
193  // ------------------------------------------------------------------------------
194  float WCTrackAlg::buildFourPointTracks(std::vector<std::vector<WCHitList> > & good_hits,
195  std::vector<double> & reco_p_list,
196  std::vector<TVector3>& mag_int_list,
197  std::vector<double>& dist_to_mag_axis_list,
198  std::vector<TVector2>& face_list,
199  std::vector<std::vector<TVector3> > & wc_hit_pos_list,
200  std::vector<TVector3> & dir_list,
201  std::vector<WCHitList> & event_final_tracks,
202  std::vector<double> & y_kink_list,
203  std::vector<TVector3>& dist_list,
204  int & WCMissed) {
205 
206  std::vector<TVector3> xyz(4, TVector3(0, 0, 0));
207  float reco_p=0;
208  WCHitList best_track;
209  std::vector<float> track_stats;
210  std::vector<float> bestRegressionStats;
211  float bestResSq = fInitialDefault;
212  for (int i = 0; i < 3; i++)
213  bestRegressionStats.push_back(fInitialDefault);
214 
215  // Loop over all combinations of hits, and find the positions
216  for (size_t iHit0 = 0; iHit0 < good_hits[0][0].hits.size(); ++iHit0) {
217  for (size_t iHit1 = 0; iHit1 < good_hits[0][1].hits.size(); ++iHit1) {
218  for (size_t iHit2 = 0; iHit2 < good_hits[1][0].hits.size(); ++iHit2) {
219  for (size_t iHit3 = 0; iHit3 < good_hits[1][1].hits.size(); ++iHit3) {
220  for (size_t iHit4 = 0; iHit4 < good_hits[2][0].hits.size(); ++iHit4) {
221  for (size_t iHit5 = 0; iHit5 < good_hits[2][1].hits.size(); ++iHit5) {
222  for (size_t iHit6 = 0; iHit6 < good_hits[3][0].hits.size(); ++iHit6) {
223  for (size_t iHit7 = 0; iHit7 < good_hits[3][1].hits.size(); ++iHit7) {
225  track.hits.push_back(good_hits[0][0].hits[iHit0]);
226  track.hits.push_back(good_hits[0][1].hits[iHit1]);
227  track.hits.push_back(good_hits[1][0].hits[iHit2]);
228  track.hits.push_back(good_hits[1][1].hits[iHit3]);
229  track.hits.push_back(good_hits[2][0].hits[iHit4]);
230  track.hits.push_back(good_hits[2][1].hits[iHit5]);
231  track.hits.push_back(good_hits[3][0].hits[iHit6]);
232  track.hits.push_back(good_hits[3][1].hits[iHit7]);
233  findTheHitPositions(track,xyz,WCMissed);
234  // Do regression on the four points to find the straightest track in Y
235  track_stats = Regression(xyz,WCMissed);
236  if (track_stats[2]<fabs(bestResSq)) {
237  best_track=track;
238  bestRegressionStats=track_stats;
239  bestResSq=track_stats[2];
240  }
241  }
242  }
243  }
244  }
245  }
246  }
247  }
248  }
249 
250  if (bestResSq<12) {
251 
252  // Now we should have the straightest track in Y, which will be the track that goes to the event.
253  // Now we get the momentum and projections onto the detector
254  calculateTheMomentum(best_track,xyz,reco_p, bestRegressionStats, mag_int_list, dist_to_mag_axis_list);
255 
256  std::vector<TVector3> wc_hit_pos;
257  for (size_t i = 0; i < 4; ++i) {
258  wc_hit_pos.push_back(TVector3(xyz[i].X()/10.,
259  xyz[i].Y()/10.,
260  xyz[i].Z()/10.));
261  }
262  wc_hit_pos_list.push_back(wc_hit_pos);
263  reco_p_list.push_back(reco_p);
264 
265  // We should also have the x,y,z points of the best_track, so now find where it hits the NOvA detector
266  projectToDetector(best_track,xyz,face_list,dir_list);
267  calculateTrackKink_Dists(xyz,bestRegressionStats,y_kink_list,dist_list);
268  event_final_tracks.push_back(best_track);
269 
270  }
271 
272  return bestResSq;
273 
274  }
275 
276  // ------------------------------------------------------------------------------
278  std::vector<TVector3>& xyz,
279  int & WCMissed) {
280 
281  std::vector<float> x_wires;
282  std::vector<float> y_wires;
283  for (int i=0; i<4; ++i) {
284  x_wires.push_back(fInitialDefault);
285  y_wires.push_back(fInitialDefault);
286  }
287 
288  for (size_t iHit = 0; iHit < track.hits.size(); ++iHit) {
289  int var = iHit % 2;
290  int jHit= iHit;
291  if (var == 0 and int(jHit != 2*(WCMissed-1)+var)) // Skip iHit for the missed WC
292  x_wires[(iHit-var)/2]=(track.hits.at(iHit).wire);
293  if (var == 1 and int(jHit != 2*(WCMissed-1)+var))
294  y_wires[(iHit-var)/2]=(track.hits.at(iHit).wire);
295  }
296 
297  std::vector<TVector3> pos = {fGeo->BeamlineComponentPos(BeamlineComponent::WC1, fBeamlineCoordSystem),
298  fGeo->BeamlineComponentPos(BeamlineComponent::WC2, fBeamlineCoordSystem),
299  fGeo->BeamlineComponentPos(BeamlineComponent::WC3, fBeamlineCoordSystem),
300  fGeo->BeamlineComponentPos(BeamlineComponent::WC4, fBeamlineCoordSystem)};
301 
302  for (int iWC = 0; iWC < 4; ++iWC) {
303  if (iWC != WCMissed-1) {
304  // Find the position of the hit, correcting for the rotation of the WC
305  // Unit mm
306  xyz[iWC] = TVector3(pos[iWC].X()*10. + x_wires[iWC] * TMath::Cos(TMath::DegToRad()*fGeo->WCAngle(iWC)),
307  pos[iWC].Y()*10. + y_wires[iWC],
308  pos[iWC].Z()*10. + x_wires[iWC] * TMath::Sin(TMath::DegToRad()*fGeo->WCAngle(iWC)));
309  }
310  if (iWC == WCMissed-1)
311  // Just making sure a value is set, should never be used
312  xyz[iWC] = TVector3(fInitialDefault, fInitialDefault, fInitialDefault);
313  }
314 
315  return;
316 
317  }
318 
319  // ------------------------------------------------------------------------------
320  std::vector<float> WCTrackAlg::Regression(const std::vector<TVector3>& xyz,
321  int & WCMissed) {
322 
323  // Calculate the linear regression on yz plane, since the track deflects in x-axis
324  std::vector<float> RegressionValues;
325  int Npoints = 0;
326  float sum_z = 0;
327  float sum_zz = 0;
328  float sum_y = 0;
329  float sum_yz = 0;
330  float intercept = 0;
331  float slope = 0;
332  float residualsquare = 0;
333  float residual = 0;
334  for (int i = 0; i < 4; ++i) {
335  if (i != WCMissed-1) { //turning WC# to array index
336  sum_y += xyz[i].Y();
337  sum_zz += xyz[i].Z()*xyz[i].Z();
338  sum_z += xyz[i].Z();
339  sum_yz += xyz[i].Y()*xyz[i].Z();
340  ++Npoints; //Residual is a function of number of points used.
341  }
342  }
343 
344  slope = (Npoints*sum_yz-sum_y*sum_z)/(Npoints*sum_zz-sum_z*sum_z);
345  intercept = (sum_y-slope*sum_z)/Npoints;
346 
347  for (int i = 0; i < 4; ++i) {
348  if (i != WCMissed-1) {
349  residual = (xyz[i].Y()-slope*xyz[i].Z()-intercept)/std::sqrt(1+slope*slope);
350  residualsquare += (residual)*(residual);
351  }
352  }
353  float avgresidual = std::sqrt(residualsquare)/(Npoints-2);
354 
355  RegressionValues.push_back(slope);
356  RegressionValues.push_back(intercept);
357  RegressionValues.push_back(avgresidual);
358 
359  return RegressionValues;
360 
361  }
362 
363  // ------------------------------------------------------------------------------
365  std::vector<TVector3>& xyz,
366  float & reco_p,
367  std::vector<float> & BestTrackStats,
368  std::vector<TVector3>& mag_int_list,
369  std::vector<double> & dist_to_mag_axis_list) {
370 
371  // We need the x,y,z again, which would be for the last track combination tried.
372  // So we need to find the positions again
373  findTheHitPositions(best_track,xyz,fWCMissed);
374 
375  // Calculate the angle of the track, in the x,z and y,z planes, in upstream(us) and downstream(ds) ends of the WC
376  float dx_us = xyz[1].X() - xyz[0].X();
377  float dx_ds = xyz[3].X() - xyz[2].X();
378  float dz_us = xyz[1].Z() - xyz[0].Z();
379  float dz_ds = xyz[3].Z() - xyz[2].Z();
380  float theta_x_us = atan(dx_us/dz_us);
381  float theta_x_ds = atan(dx_ds/dz_ds);
382 
383  float unscaled_reco_p = calculateRecoP(theta_x_us,theta_x_ds,BestTrackStats[0]);
384 
385  // project track to magnet face to find where it enters relative to the central axis
386  float dist_to_axis = projectToMagnetFace(best_track,xyz,mag_int_list,dist_to_mag_axis_list);
387 
388  // scale momentum based on where particle enters magnet
389  float a = 1.00068819;
390  float b = -0.04635663;
391  float c = 0.11454936;
392  float d = 0.00971984;
393  float e = -0.11526033;
394  float f = 0.01645662;
395  float x = dist_to_axis/100; // convert cm to m
396  float scale_factor = (a+b*x)*(1/(1+exp((x-c)/d)) - 1/(1+exp((x-e)/f)));
397  reco_p = scale_factor * unscaled_reco_p;
398 
399  std::cout << "B: " << fBField << " momentum: " << reco_p << std::endl;
400 
401  }
402 
403  // ------------------------------------------------------------------------------
405  const std::vector<TVector3>& xyz,
406  std::vector<TVector3>& mag_int_list,
407  std::vector<double>& dist_to_mag_axis_list
408  ) {
409 
410  // Get magnet geometry - need center of magnet, length, and angle to find equation of front face
411  float mag_pos_x = fGeo->BeamlineComponentPos(BeamlineComponent::Magnet, fBeamlineCoordSystem).X();
412  float mag_pos_z = fGeo->BeamlineComponentPos(BeamlineComponent::Magnet, fBeamlineCoordSystem).Z();
413  float magnet_angle = fGeo->MagnetAngle();
414  float mag_length = fGeo->MagnetEffectiveLength();
415 
416  // center point of front face of magnet:
417  float mag_front_x = mag_pos_x+mag_length/2*sin(magnet_angle*TMath::DegToRad());
418  float mag_front_z = mag_pos_z-mag_length/2*cos(magnet_angle*TMath::DegToRad());
419 
420  // WC track intercept with front face of magnet
421  float dxdz_us = (xyz[1].X() - xyz[0].X())/(xyz[1].Z() - xyz[0].Z());
422  float dydz_us = (xyz[1].Y() - xyz[0].Y())/(xyz[1].Z() - xyz[0].Z());
423  float tanang = tan((90-magnet_angle)*TMath::DegToRad());
424  float zint_front = (xyz[1].X()/10 - dxdz_us*xyz[1].Z()/10 + tanang*mag_front_z-mag_front_x)/(tanang-dxdz_us);
425  float xint_front = tanang*zint_front - tanang*mag_front_z+mag_front_x;
426  float yint_front = dydz_us*(zint_front-xyz[1].Z()/10)+xyz[1].Y()/10;
427 
428  // Transverse distance from point of intercept to central axis of magnet
429  int sign_of_dist = 1;
430  if (xint_front < mag_front_x)
431  sign_of_dist = -1;
432  float dist_to_mag_axis = sign_of_dist*sqrt(pow(xint_front-mag_front_x,2)+pow(zint_front-mag_front_z,2));
433 
434  mag_int_list.push_back(TVector3(xint_front,yint_front,zint_front));
435  dist_to_mag_axis_list.push_back(dist_to_mag_axis);
436 
437  return dist_to_mag_axis;
438 
439  }
440 
441  // ------------------------------------------------------------------------------
443  const std::vector<TVector3>& xyz,
444  std::vector<TVector2>& face_list,
445  std::vector<TVector3>& incoming_dir_list) {
446 
447  // Direction of track
448  float dx_ds = xyz[3].X() - xyz[2].X();
449  float dy_ds = xyz[3].Y() - xyz[2].Y();
450  float dz_ds = xyz[3].Z() - xyz[2].Z();
451  TVector3 vec(dx_ds, dy_ds, dz_ds);
452  incoming_dir_list.push_back(vec.Unit());
453 
454  // Get interseption with the front face of the detector
455  TVector3 nova_face = fGeo->BeamlineComponentPos(BeamlineComponent::NOvA,
457  TVector3 normal(0, 0, 1), point(0, 0, nova_face.Z()*10.);
458  TVector3 start = xyz[3];
459  float p = (normal*(point-start))/(normal*vec);
460  TVector3 end_int = start + (p*vec);
461  face_list.push_back(TVector2(end_int.X()/10., end_int.Y()/10.));
462 
463  return;
464 
465  }
466 
467  // ------------------------------------------------------------------------------
468  void WCTrackAlg::calculateTrackKink_Dists(const std::vector<TVector3>& xyz,
469  std::vector<float>& track_stats,
470  std::vector<double>& y_kink_list,
471  std::vector<TVector3>& dist_list) {
472 
473  // Determine the equations of the lines, x=mz+b and y=mz+b, for the US and DS legs
474  float dx_us = xyz[1].X() - xyz[0].X();
475  float dy_us = xyz[1].Y() - xyz[0].Y();
476  float dz_us = xyz[1].Z() - xyz[0].Z();
477  float dx_ds = xyz[3].X() - xyz[2].X();
478  float dy_ds = xyz[3].Y() - xyz[2].Y();
479  float dz_ds = xyz[3].Z() - xyz[2].Z();
480  float x_us_slope = dx_us/dz_us;
481  float y_us_slope = dy_us/dz_us;
482  float x_ds_slope = dx_ds/dz_ds;
483  float y_ds_slope = dy_ds/dz_ds;
484  float x_us_int = xyz[1].X() - x_us_slope*xyz[1].Z();
485  float y_us_int = xyz[1].Y() - y_us_slope*xyz[1].Z();
486  float x_ds_int = xyz[3].X() - x_ds_slope*xyz[3].Z();
487  float y_ds_int = xyz[3].Y() - y_ds_slope*xyz[3].Z();
488 
489  // First Y_kink, the angle difference between the us and ds legs.
490  y_kink_list.push_back(atan(y_ds_slope)-atan(y_us_slope));
491 
492  // Because we need a way to compare tracks, regardless of current setting, we have to have a standard for the midplane.
493  // We will use the normal midplane (fMP_X=fMP_M=1)
494 
495  float magnet_angle = fGeo->MagnetAngle()*TMath::DegToRad();
496  float magnet_midplane = fGeo->MagnetMidplaneIntercept()*10.;
497 
498  // Get US intersection with midplane
499  float z_mp_us = (magnet_midplane-x_us_int)/(x_us_slope-1/(tan(magnet_angle)));
500  float x_mp_us = x_us_slope*z_mp_us+x_us_int;
501  float y_mp_us = y_us_slope*z_mp_us+y_us_int;
502 
503  // Get DS intersection with midplane
504  float z_mp_ds = (magnet_midplane-x_ds_int)/(x_ds_slope-1/(tan(magnet_angle)));
505  float x_mp_ds = x_ds_slope*z_mp_ds+x_ds_int;
506  float y_mp_ds = y_ds_slope*z_mp_ds+y_ds_int;
507 
508  dist_list.push_back(TVector3((x_mp_ds-x_mp_us)/10.,
509  (y_mp_ds-y_mp_us)/10.,
510  (z_mp_ds-z_mp_us)/10.));
511 
512  return;
513 
514  }
515 
516  // ------------------------------------------------------------------------------
517  void WCTrackAlg::setMagneticField(float b_field) {
518  fBField = b_field;
519  return;
520  }
521 
522 }
float WCAngle(unsigned int wc) const
WC angle.
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
Definition: WCTrackAlg.h:119
float MagnetMidplaneIntercept() const
Intercept of the magnet midplane?
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
bool shouldSkipTrigger(std::vector< std::vector< WCHitList > > &good_hits, int &WCMissed, TH1F *&WCDist)
Definition: WCTrackAlg.cxx:120
float MagnetEffectiveLength() const
Effective length of the magnet.
void setMagneticField(float field)
Definition: WCTrackAlg.cxx:517
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Definition: event.h:19
std::vector< float > Regression(const std::vector< TVector3 > &xyz, int &WCMissed)
Definition: WCTrackAlg.cxx:320
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
Float_t Y
Definition: plot.C:38
std::vector< WCHit > hits
float projectToMagnetFace(WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list)
Definition: WCTrackAlg.cxx:404
T atan(T number)
Definition: d0nt_math.hpp:66
Float_t Z
Definition: plot.C:38
TVector3 BeamlineComponentPos(BeamlineComponent component, BeamlineCoordSystem system) const
void hits()
Definition: readHits.C:15
void calculateTheMomentum(WCHitList &best_track, std::vector< TVector3 > &xyz, float &reco_p, std::vector< float > &BestTrackStats, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list)
Definition: WCTrackAlg.cxx:364
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
float buildFourPointTracks(std::vector< std::vector< WCHitList > > &good_hits, std::vector< double > &reco_p_list, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list, std::vector< TVector2 > &fact_list, std::vector< std::vector< TVector3 > > &wc_hit_pos_list, std::vector< TVector3 > &dir_list, std::vector< WCHitList > &event_final_tracks, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list, int &WCMissed)
Definition: WCTrackAlg.cxx:194
Float_t d
Definition: plot.C:236
void calculateTrackKink_Dists(const std::vector< TVector3 > &xyz, std::vector< float > &track_stats, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list)
Definition: WCTrackAlg.cxx:468
Eigen::VectorXd vec
float MagnetAngle() const
Magnet angle.
WCTrackAlg(const fhicl::ParameterSet &pset)
Definition: WCTrackAlg.cxx:20
OStream cout
Definition: OStream.cxx:6
void projectToDetector(WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector2 > &face_list, std::vector< TVector3 > &dir_list)
Definition: WCTrackAlg.cxx:442
T sin(T number)
Definition: d0nt_math.hpp:132
int num
Definition: f2_nu.C:119
Algorithm implementation for wire chamber reconstruction. Part of the beamline reconstruction for NOv...
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
T tan(T number)
Definition: d0nt_math.hpp:144
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
Float_t track
Definition: plot.C:35
float calculateRecoP(float theta_x_us, float theta_x_ds, float bestTrackSlope)
Definition: WCTrackAlg.cxx:183
void findTheHitPositions(WCHitList &track, std::vector< TVector3 > &xyz, int &WCMissed)
Definition: WCTrackAlg.cxx:277
void reconstructTracks(std::vector< double > &reco_p_list, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list, std::vector< TVector2 > &face_list, std::vector< std::vector< TVector3 > > &wc_hit_pos_list, std::vector< TVector3 > &dir_list, std::vector< WCHitList > &event_final_tracks, std::vector< std::vector< WCHitList > > &good_hits, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list, int &WCMissed, TH1F *&WCdistribution, float &residual)
Main function called for each event.
Definition: WCTrackAlg.cxx:63
void reconfigure(const fhicl::ParameterSet &pset)
Definition: WCTrackAlg.cxx:27