43 unsigned int coordSystem = pset.
get<
unsigned int>(
"CoordSystem", 0);
44 switch (coordSystem) {
52 std::cout <<
"Unknown coordinate system (" << coordSystem <<
"). " 53 <<
"Options are 0: Beamline; 1: NOvA Detector." <<
std::endl;
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,
75 TH1F* & WCdistribution,
104 dist_to_mag_axis_list,
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)
138 if (
fNHits > 2 && WCMissed == 1)
140 if (
fNHits > 2 && WCMissed == 2)
142 if (
fNHits > 2 && WCMissed == 3)
144 if (
fNHits > 2 && WCMissed == 4)
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) {
165 if (WCMissed==1 || WCMissed==4)
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();
186 float alpha_1 = theta_x_us + magnet_angle*TMath::DegToRad();
187 float alpha_2 = theta_x_ds + magnet_angle*TMath::DegToRad();
189 float denom = (3.33564*(
sin(alpha_2) -
sin(alpha_1)))*
cos(
atan(bestTrackSlope));
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,
206 std::vector<TVector3> xyz(4, TVector3(0, 0, 0));
209 std::vector<float> track_stats;
210 std::vector<float> bestRegressionStats;
212 for (
int i = 0;
i < 3;
i++)
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]);
236 if (track_stats[2]<
fabs(bestResSq)) {
238 bestRegressionStats=track_stats;
239 bestResSq=track_stats[2];
254 calculateTheMomentum(best_track,xyz,reco_p, bestRegressionStats, mag_int_list, dist_to_mag_axis_list);
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.,
262 wc_hit_pos_list.push_back(wc_hit_pos);
263 reco_p_list.push_back(reco_p);
268 event_final_tracks.push_back(best_track);
278 std::vector<TVector3>& xyz,
281 std::vector<float> x_wires;
282 std::vector<float> y_wires;
283 for (
int i=0;
i<4; ++
i) {
288 for (
size_t iHit = 0; iHit < track.
hits.size(); ++iHit) {
291 if (var == 0 and
int(jHit != 2*(WCMissed-1)+var))
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);
302 for (
int iWC = 0; iWC < 4; ++iWC) {
303 if (iWC != WCMissed-1) {
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)));
310 if (iWC == WCMissed-1)
324 std::vector<float> RegressionValues;
332 float residualsquare = 0;
334 for (
int i = 0;
i < 4; ++
i) {
335 if (
i != WCMissed-1) {
337 sum_zz += xyz[
i].Z()*xyz[
i].Z();
339 sum_yz += xyz[
i].Y()*xyz[
i].Z();
344 slope = (Npoints*sum_yz-sum_y*sum_z)/(Npoints*sum_zz-sum_z*sum_z);
345 intercept = (sum_y-slope*sum_z)/Npoints;
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);
353 float avgresidual =
std::sqrt(residualsquare)/(Npoints-2);
355 RegressionValues.push_back(slope);
356 RegressionValues.push_back(intercept);
357 RegressionValues.push_back(avgresidual);
359 return RegressionValues;
365 std::vector<TVector3>& xyz,
367 std::vector<float> & BestTrackStats,
368 std::vector<TVector3>& mag_int_list,
369 std::vector<double> & dist_to_mag_axis_list) {
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);
383 float unscaled_reco_p =
calculateRecoP(theta_x_us,theta_x_ds,BestTrackStats[0]);
386 float dist_to_axis =
projectToMagnetFace(best_track,xyz,mag_int_list,dist_to_mag_axis_list);
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;
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;
405 const std::vector<TVector3>& xyz,
406 std::vector<TVector3>& mag_int_list,
407 std::vector<double>& dist_to_mag_axis_list
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());
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;
429 int sign_of_dist = 1;
430 if (xint_front < mag_front_x)
432 float dist_to_mag_axis = sign_of_dist*
sqrt(
pow(xint_front-mag_front_x,2)+
pow(zint_front-mag_front_z,2));
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);
437 return dist_to_mag_axis;
443 const std::vector<TVector3>& xyz,
444 std::vector<TVector2>& face_list,
445 std::vector<TVector3>& incoming_dir_list) {
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());
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.));
469 std::vector<float>& track_stats,
470 std::vector<double>& y_kink_list,
471 std::vector<TVector3>& dist_list) {
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();
490 y_kink_list.push_back(
atan(y_ds_slope)-
atan(y_us_slope));
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;
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;
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.));
float WCAngle(unsigned int wc) const
WC angle.
beamlinegeo::BeamlineCoordSystem fBeamlineCoordSystem
float MagnetMidplaneIntercept() const
Intercept of the magnet midplane?
fvar< T > fabs(const fvar< T > &x)
bool shouldSkipTrigger(std::vector< std::vector< WCHitList > > &good_hits, int &WCMissed, TH1F *&WCDist)
float MagnetEffectiveLength() const
Effective length of the magnet.
void setMagneticField(float field)
std::vector< float > Regression(const std::vector< TVector3 > &xyz, int &WCMissed)
art::ServiceHandle< beamlinegeo::BeamlineGeometry > fGeo
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)
TVector3 BeamlineComponentPos(BeamlineComponent component, BeamlineCoordSystem system) const
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)
T get(std::string const &key) const
fvar< T > exp(const fvar< T > &x)
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 calculateTrackKink_Dists(const std::vector< TVector3 > &xyz, std::vector< float > &track_stats, std::vector< double > &y_kink_list, std::vector< TVector3 > &dist_list)
float MagnetAngle() const
Magnet angle.
WCTrackAlg(const fhicl::ParameterSet &pset)
void projectToDetector(WCHitList &best_track, const std::vector< TVector3 > &xyz, std::vector< TVector2 > &face_list, std::vector< TVector3 > &dir_list)
Algorithm implementation for wire chamber reconstruction. Part of the beamline reconstruction for NOv...
float calculateRecoP(float theta_x_us, float theta_x_ds, float bestTrackSlope)
void findTheHitPositions(WCHitList &track, std::vector< TVector3 > &xyz, int &WCMissed)
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.
void reconfigure(const fhicl::ParameterSet &pset)