Skip to content

Commit

Permalink
Merge pull request #57 from DUNE/25-Track-Reco-With-Circunferences-fits
Browse files Browse the repository at this point in the history
25 track reco with circunferences fits
  • Loading branch information
valerpia authored Sep 26, 2024
2 parents b56a28f + 96745a4 commit 248cea2
Show file tree
Hide file tree
Showing 4 changed files with 65 additions and 97 deletions.
3 changes: 2 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -21,4 +21,5 @@ build/*
*.sh
*.png
*.gdml
# tests/*
tests/*
tests/python_tools/__pycache__/
29 changes: 7 additions & 22 deletions include/SANDRecoUtils.h
Original file line number Diff line number Diff line change
Expand Up @@ -253,23 +253,14 @@ class Circle
return {x_l(angle), y_l(angle)};
}

TVector2 GetDerivativeAt(double angle) const {
return {dx_derivative(angle), dy_derivative(angle)};
}

TVector2 GetDerivativeAt(double arg_x, double arg_y) const {
double angle = GetAngleFromPoint(arg_x, arg_y);
return {dx_derivative(angle), dy_derivative(angle)};
}
// TVector2 GetDerivativeAt(double angle) const {
// return {dx_derivative(angle), dy_derivative(angle)};
// }

double GetAngleFromPoint(double arg_x, double arg_y) const {
// return the angle [0, 2pi) corresponding to the point (x,y)
if(arg_y - center_y_ > 0){
return atan2(arg_y - center_y_, arg_x - center_x_);
}else{
return atan2(arg_y - center_y_, arg_x - center_x_) + 2 * TMath::Pi();
}
}
// TVector2 GetDerivativeAt(double arg_x, double arg_y) const {
// double angle = GetAngleFromPoint(arg_x, arg_y);
// return {dx_derivative(angle), dy_derivative(angle)};
// }

TF1* GetUpperSemiCircle() const {
// y as function of the x coordinate
Expand Down Expand Up @@ -659,12 +650,6 @@ double NLL(Helix& h,const std::vector<dg_wire>& digits); // negativ

double FunctorNLL(const double* p);

double GetDipAngleFromCircleLine(const Circle& circle,
const Line2D& line,
double Phi0,
int helicity,
TVector3& Momentum);

Helix GetHelixFromCircleLine(const Circle& circle,
const Line2D& line,
const Helix& true_helix,
Expand Down
46 changes: 14 additions & 32 deletions src/SANDRecoUtils.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -552,36 +552,6 @@ double RecoUtils::FunctorNLL(const double* p)
return nll;
}

double RecoUtils::GetDipAngleFromCircleLine(const Circle& circle,
const Line2D& line,
double Phi0,
int helicity,
TVector3& momentum){
/*
input:
circle : fitted trajectory in the ZY bending plane
line : fitted trajectory in the XZ plane
Phi0 : angle corresponding to the trajecotry starting point
Momentum : reference TVector3 that will be filled with particle
total momentum
NOTE: a guess for the initial vertex is needed because related to the dip angle
*/
double pt = circle.R() * 0.3 * 0.6;
// define theta as the angle [0, 2 pi) between the momentum vector and z axis
double theta = Phi0 - helicity * TMath::Pi() / 2.;
double py = pt * sin(theta);
double pz = pt * cos(theta);

// slope of the tangent in the zx plane
double pz_over_px = line.m();
double px = pz / pz_over_px;

// fill momentum
momentum = {px, py, pz};

return atan2(px, pt);
}

Helix RecoUtils::GetHelixFromCircleLine(const Circle& circle,
const Line2D& line,
const Helix& true_helix,
Expand All @@ -599,10 +569,22 @@ Helix RecoUtils::GetHelixFromCircleLine(const Circle& circle,

auto helicity = true_helix.h();

double Phi0 = circle.GetAngleFromPoint(vertex.Z() ,vertex.Y());
double pt = circle.R() * 0.29979 * 0.6;

double pz = pt * (helicity * (vertex.Y() - circle.center_y()) / circle.R());

double py = pt * (-helicity * (vertex.Z() - circle.center_x()) / circle.R());

double pz_over_px = line.m();

double dip_angle = RecoUtils::GetDipAngleFromCircleLine(circle, line, Phi0, helicity, momentum);
double px = pz / pz_over_px;

double Phi0 = TMath::ATan2(py, pz) + helicity * TMath::Pi()*0.5;;

momentum = {px, py, pz};

double dip_angle = TMath::ATan2(px, pt);

return Helix(circle.R(), dip_angle, Phi0, helicity, vertex);
}

Expand Down
84 changes: 42 additions & 42 deletions src/reconstructionNLLmethod.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1015,55 +1015,55 @@ TF1* GetRecoXZTrack(TF1* first_guess,
return SinFinalFit;
}

Helix Reconstruct(Circle FittedCircle,
Line2D FittedLine,
const std::vector<dg_wire*>& hor_wires,
const std::vector<dg_wire*>& ver_wires
){
/*
Define the reconstructed helix from
the separate reconstruction of the
track in the XZ plane and ZY plane
*/
double zc = FittedCircle.center_x();
double yc = FittedCircle.center_y();
double R = FittedCircle.R();
double m = FittedLine.m();
double q = FittedLine.q();
// Helix Reconstruct(Circle FittedCircle,
// Line2D FittedLine,
// const std::vector<dg_wire*>& hor_wires,
// const std::vector<dg_wire*>& ver_wires
// ){
// /*
// Define the reconstructed helix from
// the separate reconstruction of the
// track in the XZ plane and ZY plane
// */
// double zc = FittedCircle.center_x();
// double yc = FittedCircle.center_y();
// double R = FittedCircle.R();
// double m = FittedLine.m();
// double q = FittedLine.q();

bool forward_track = (ver_wires[0]->z < ver_wires[1]->z) ? 1 : 0;

TVector2 vertex_XZ = {ver_wires[0]->x, ver_wires[0]->x * m + q};

TVector2 vertex_ZY;
double upper_y = FittedCircle.GetUpperSemiCircle()->Eval(hor_wires[0]->z);
double lower_y = FittedCircle.GetLowerSemiCircle()->Eval(hor_wires[0]->z);
if(fabs(upper_y - hor_wires[0]->y) < fabs(lower_y - hor_wires[0]->y)){
vertex_ZY = {hor_wires[0]->z, upper_y};
}else{
vertex_ZY = {hor_wires[0]->z, lower_y};
}
// bool forward_track = (ver_wires[0]->z < ver_wires[1]->z) ? 1 : 0;

// TVector2 vertex_XZ = {ver_wires[0]->x, ver_wires[0]->x * m + q};

// TVector2 vertex_ZY;
// double upper_y = FittedCircle.GetUpperSemiCircle()->Eval(hor_wires[0]->z);
// double lower_y = FittedCircle.GetLowerSemiCircle()->Eval(hor_wires[0]->z);
// if(fabs(upper_y - hor_wires[0]->y) < fabs(lower_y - hor_wires[0]->y)){
// vertex_ZY = {hor_wires[0]->z, upper_y};
// }else{
// vertex_ZY = {hor_wires[0]->z, lower_y};
// }


double Phi0 = FittedCircle.GetAngleFromPoint(vertex_ZY.X(), vertex_ZY.Y());
auto track_direction_ZY = FittedCircle.GetDerivativeAt(vertex_ZY.X(), vertex_ZY.Y()).Unit();
TVector2 pt = R*0.3*0.6 * track_direction_ZY;
/*
reconstructed pt is the value of the tangent0
to the circle in the first point of the curve
*/
double pz = (forward_track) ? -1. * fabs(pt.X()) : fabs(pt.X());
double py = pt.Y();
double px = pz / FittedLine.m();
// double Phi0 = FittedCircle.GetAngleFromPoint(vertex_ZY.X(), vertex_ZY.Y());
// auto track_direction_ZY = FittedCircle.GetDerivativeAt(vertex_ZY.X(), vertex_ZY.Y()).Unit();
// TVector2 pt = R*0.3*0.6 * track_direction_ZY;
// /*
// reconstructed pt is the value of the tangent0
// to the circle in the first point of the curve
// */
// double pz = (forward_track) ? -1. * fabs(pt.X()) : fabs(pt.X());
// double py = pt.Y();
// double px = pz / FittedLine.m();

LOG("R", TString::Format("Reconstructed muon (px, py, pz) : (%f, %f, %f)",px,py,pz).Data());
// LOG("R", TString::Format("Reconstructed muon (px, py, pz) : (%f, %f, %f)",px,py,pz).Data());

double dip = atan2(px, pt.Mod());
// double dip = atan2(px, pt.Mod());

Helix h(R, dip, Phi0, 1, {vertex_XZ.X(), vertex_ZY.Y(), vertex_XZ.Y()});
// Helix h(R, dip, Phi0, 1, {vertex_XZ.X(), vertex_ZY.Y(), vertex_XZ.Y()});

return h;
}
// return h;
// }

// MAIN________________________________________________________________________

Expand Down

0 comments on commit 248cea2

Please sign in to comment.