Public Member Functions | Private Attributes | List of all members
beamlinereco::WCTrackAlg Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-12-03/BeamlineReco/WCTrackAlg.h"

Public Member Functions

 WCTrackAlg (const fhicl::ParameterSet &pset)
 
 ~WCTrackAlg ()
 
void reconfigure (const fhicl::ParameterSet &pset)
 
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. More...
 
bool shouldSkipTrigger (std::vector< std::vector< WCHitList > > &good_hits, int &WCMissed, TH1F *&WCDist)
 
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)
 
void findTheHitPositions (WCHitList &track, std::vector< TVector3 > &xyz, int &WCMissed)
 
std::vector< float > Regression (const std::vector< TVector3 > &xyz, int &WCMissed)
 
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)
 
float projectToMagnetFace (WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector3 > &mag_int_list, std::vector< double > &dist_to_mag_axis_list)
 
void projectToDetector (WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector2 > &face_list, std::vector< TVector3 > &dir_list)
 
void calculateTrackKink_Dists (const std::vector< TVector3 > &xyz, std::vector< float > &track_stats, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list)
 
float calculateRecoP (float theta_x_us, float theta_x_ds, float bestTrackSlope)
 
void setMagneticField (float field)
 

Private Attributes

int fInitialDefault
 
bool fVerbose
 
bool fPickyTracks
 
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
 
int fWCMissed
 
unsigned int fNHits
 
float fBField
 
float fMidplaneIntercept
 
float fMidplaneSlope
 
art::ServiceHandle< beamlinegeo::BeamlineGeometryfGeo
 

Detailed Description

Definition at line 39 of file WCTrackAlg.h.

Constructor & Destructor Documentation

beamlinereco::WCTrackAlg::WCTrackAlg ( const fhicl::ParameterSet pset)

Definition at line 20 of file WCTrackAlg.cxx.

References reconfigure().

20  {
21  this->reconfigure(pset);
22  }
void reconfigure(const fhicl::ParameterSet &pset)
Definition: WCTrackAlg.cxx:27
beamlinereco::WCTrackAlg::~WCTrackAlg ( )

Definition at line 24 of file WCTrackAlg.cxx.

24 { }

Member Function Documentation

float beamlinereco::WCTrackAlg::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 at line 194 of file WCTrackAlg.cxx.

References calculateTheMomentum(), calculateTrackKink_Dists(), stan::math::fabs(), findTheHitPositions(), fInitialDefault, hits(), beamlinereco::WCHitList::hits, MECModelEnuComparisons::i, projectToDetector(), Regression(), track, X, Y, and Z.

Referenced by reconstructTracks().

204  {
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) {
224  WCHitList track;
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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< float > Regression(const std::vector< TVector3 > &xyz, int &WCMissed)
Definition: WCTrackAlg.cxx:320
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
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
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
void projectToDetector(WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector2 > &face_list, std::vector< TVector3 > &dir_list)
Definition: WCTrackAlg.cxx:442
Float_t X
Definition: plot.C:38
Float_t track
Definition: plot.C:35
void findTheHitPositions(WCHitList &track, std::vector< TVector3 > &xyz, int &WCMissed)
Definition: WCTrackAlg.cxx:277
float beamlinereco::WCTrackAlg::calculateRecoP ( float  theta_x_us,
float  theta_x_ds,
float  bestTrackSlope 
)

Definition at line 183 of file WCTrackAlg.cxx.

References std::atan(), std::cos(), stan::math::fabs(), fBField, fGeo, beamlinegeo::BeamlineGeometry::MagnetAngle(), beamlinegeo::BeamlineGeometry::MagnetEffectiveLength(), num, and std::sin().

Referenced by calculateTheMomentum().

183  {
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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
float MagnetEffectiveLength() const
Effective length of the magnet.
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
T atan(T number)
Definition: d0nt_math.hpp:66
float MagnetAngle() const
Magnet angle.
T sin(T number)
Definition: d0nt_math.hpp:132
int num
Definition: f2_nu.C:119
T cos(T number)
Definition: d0nt_math.hpp:78
void beamlinereco::WCTrackAlg::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 at line 364 of file WCTrackAlg.cxx.

References a, std::atan(), b, plot_validation_datamc::c, calculateRecoP(), om::cout, d, e, allTimeWatchdog::endl, stan::math::exp(), MakeMiniprodValidationCuts::f, fBField, findTheHitPositions(), fWCMissed, projectToMagnetFace(), and submit_syst::x.

Referenced by buildFourPointTracks().

369  {
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  }
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
const double a
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Float_t d
Definition: plot.C:236
OStream cout
Definition: OStream.cxx:6
const hit & b
Definition: hits.cxx:21
Float_t e
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 beamlinereco::WCTrackAlg::calculateTrackKink_Dists ( const std::vector< TVector3 > &  xyz,
std::vector< float > &  track_stats,
std::vector< double > &  y_kink_list,
std::vector< TVector3 > &  dist_list 
)

Definition at line 468 of file WCTrackAlg.cxx.

References std::atan(), fGeo, beamlinegeo::BeamlineGeometry::MagnetAngle(), beamlinegeo::BeamlineGeometry::MagnetMidplaneIntercept(), and std::tan().

Referenced by buildFourPointTracks().

471  {
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  }
float MagnetMidplaneIntercept() const
Intercept of the magnet midplane?
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
T atan(T number)
Definition: d0nt_math.hpp:66
float MagnetAngle() const
Magnet angle.
T tan(T number)
Definition: d0nt_math.hpp:144
void beamlinereco::WCTrackAlg::findTheHitPositions ( WCHitList track,
std::vector< TVector3 > &  xyz,
int WCMissed 
)

Definition at line 277 of file WCTrackAlg.cxx.

References beamlinegeo::BeamlineGeometry::BeamlineComponentPos(), fBeamlineCoordSystem, fGeo, fInitialDefault, beamlinereco::WCHitList::hits, MECModelEnuComparisons::i, elec2geo::pos, PandAna.Demos.tute_pid_validation::var, beamlinegeo::BeamlineGeometry::WCAngle(), X, and Y.

Referenced by buildFourPointTracks(), and calculateTheMomentum().

279  {
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  }
float WCAngle(unsigned int wc) const
WC angle.
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
Definition: WCTrackAlg.h:119
Definition: event.h:19
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
Float_t Y
Definition: plot.C:38
TVector3 BeamlineComponentPos(BeamlineComponent component, BeamlineCoordSystem system) const
std::vector< hit > hits
Definition: event.h:30
Float_t X
Definition: plot.C:38
void beamlinereco::WCTrackAlg::projectToDetector ( WCHitList best_track,
const std::vector< TVector3 > &  xyz,
std::vector< TVector2 > &  face_list,
std::vector< TVector3 > &  dir_list 
)

Definition at line 442 of file WCTrackAlg.cxx.

References beamlinegeo::BeamlineGeometry::BeamlineComponentPos(), fBeamlineCoordSystem, fGeo, and febshutoff_auto::start.

Referenced by buildFourPointTracks().

445  {
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  }
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
Definition: WCTrackAlg.h:119
const char * p
Definition: xmltok.h:285
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
TVector3 BeamlineComponentPos(BeamlineComponent component, BeamlineCoordSystem system) const
Eigen::VectorXd vec
float beamlinereco::WCTrackAlg::projectToMagnetFace ( WCHitList best_track,
const std::vector< TVector3 > &  xyz,
std::vector< TVector3 > &  mag_int_list,
std::vector< double > &  dist_to_mag_axis_list 
)

Definition at line 404 of file WCTrackAlg.cxx.

References beamlinegeo::BeamlineGeometry::BeamlineComponentPos(), std::cos(), fBeamlineCoordSystem, fGeo, beamlinegeo::BeamlineGeometry::MagnetAngle(), beamlinegeo::BeamlineGeometry::MagnetEffectiveLength(), cet::pow(), std::sin(), std::sqrt(), std::tan(), Y, and Z.

Referenced by calculateTheMomentum().

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  }
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
Definition: WCTrackAlg.h:119
float MagnetEffectiveLength() const
Effective length of the magnet.
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
Definition: WCTrackAlg.h:130
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
TVector3 BeamlineComponentPos(BeamlineComponent component, BeamlineCoordSystem system) const
float MagnetAngle() const
Magnet angle.
T sin(T number)
Definition: d0nt_math.hpp:132
T cos(T number)
Definition: d0nt_math.hpp:78
T tan(T number)
Definition: d0nt_math.hpp:144
void beamlinereco::WCTrackAlg::reconfigure ( const fhicl::ParameterSet pset)

Definition at line 27 of file WCTrackAlg.cxx.

References om::cout, allTimeWatchdog::endl, fBeamlineCoordSystem, fInitialDefault, fPickyTracks, fVerbose, and fhicl::ParameterSet::get().

Referenced by WCTrackAlg().

27  {
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  }
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
Definition: WCTrackAlg.h:119
T get(std::string const &key) const
Definition: ParameterSet.h:231
OStream cout
Definition: OStream.cxx:6
void beamlinereco::WCTrackAlg::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 at line 63 of file WCTrackAlg.cxx.

References buildFourPointTracks(), fInitialDefault, fNHits, fWCMissed, shouldSkipTrigger(), and compare_h5_caf::skip.

Referenced by beamlinereco::WCTrackReco::produce().

76  {
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  }
bool shouldSkipTrigger(std::vector< std::vector< WCHitList > > &good_hits, int &WCMissed, TH1F *&WCDist)
Definition: WCTrackAlg.cxx:120
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
std::vector< float > beamlinereco::WCTrackAlg::Regression ( const std::vector< TVector3 > &  xyz,
int WCMissed 
)

Definition at line 320 of file WCTrackAlg.cxx.

References MECModelEnuComparisons::i, and std::sqrt().

Referenced by buildFourPointTracks().

321  {
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  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
void beamlinereco::WCTrackAlg::setMagneticField ( float  field)

Definition at line 517 of file WCTrackAlg.cxx.

References fBField.

Referenced by beamlinereco::WCTrackReco::produce().

517  {
518  fBField = b_field;
519  return;
520  }
bool beamlinereco::WCTrackAlg::shouldSkipTrigger ( std::vector< std::vector< WCHitList > > &  good_hits,
int WCMissed,
TH1F *&  WCDist 
)

Definition at line 120 of file WCTrackAlg.cxx.

References om::cout, allTimeWatchdog::endl, fNHits, fPickyTracks, fVerbose, hits(), and compare_h5_caf::skip.

Referenced by reconstructTracks().

122  {
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  }
void hits()
Definition: readHits.C:15
OStream cout
Definition: OStream.cxx:6

Member Data Documentation

beamlinegeo::BeamlineCoordSystem beamlinereco::WCTrackAlg::fBeamlineCoordSystem
private
float beamlinereco::WCTrackAlg::fBField
private

Definition at line 125 of file WCTrackAlg.h.

Referenced by calculateRecoP(), calculateTheMomentum(), and setMagneticField().

art::ServiceHandle<beamlinegeo::BeamlineGeometry> beamlinereco::WCTrackAlg::fGeo
private
int beamlinereco::WCTrackAlg::fInitialDefault
private
float beamlinereco::WCTrackAlg::fMidplaneIntercept
private

Definition at line 126 of file WCTrackAlg.h.

float beamlinereco::WCTrackAlg::fMidplaneSlope
private

Definition at line 127 of file WCTrackAlg.h.

unsigned int beamlinereco::WCTrackAlg::fNHits
private

Definition at line 123 of file WCTrackAlg.h.

Referenced by reconstructTracks(), and shouldSkipTrigger().

bool beamlinereco::WCTrackAlg::fPickyTracks
private

Definition at line 118 of file WCTrackAlg.h.

Referenced by reconfigure(), and shouldSkipTrigger().

bool beamlinereco::WCTrackAlg::fVerbose
private

Definition at line 117 of file WCTrackAlg.h.

Referenced by reconfigure(), and shouldSkipTrigger().

int beamlinereco::WCTrackAlg::fWCMissed
private

Definition at line 122 of file WCTrackAlg.h.

Referenced by calculateTheMomentum(), and reconstructTracks().


The documentation for this class was generated from the following files: