From 1799bdc003e8d57caf5df040f63dc26ce082447c Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Thu, 26 Mar 2026 12:48:56 -0400 Subject: [PATCH 01/15] first iteration of tracking for GEMs Need to incoperate hodo hit in chisq --- src/THcLADGEM.cxx | 6 ++- src/THcLADGEMTrack.cxx | 4 ++ src/THcLADGEMTrack.h | 26 +++++++++- src/THcLADKine.cxx | 106 ++++++++++++++++++++++++++++++++++++++++- src/THcLADKine.h | 1 + 5 files changed, 139 insertions(+), 4 deletions(-) diff --git a/src/THcLADGEM.cxx b/src/THcLADGEM.cxx index f6049d8..a164365 100644 --- a/src/THcLADGEM.cxx +++ b/src/THcLADGEM.cxx @@ -230,10 +230,12 @@ Int_t THcLADGEM::DefineVariables(EMode mode) { {"trk.d0", "Track dist from vertex", "fGEMTracks.THcLADGEMTrack.GetD0()"}, {"trk.d0_good", "Track dist from true vertex", "fGEMTracks.THcLADGEMTrack.GetGoodD0()"}, {"trk.projz", "Projected z-vertex", "fGEMTracks.THcLADGEMTrack.GetProjVz()"}, + {"trk.projx", "Projected x-vertex", "fGEMTracks.THcLADGEMTrack.GetProjVx()"}, {"trk.projy", "Projected y-vertex", "fGEMTracks.THcLADGEMTrack.GetProjVy()"}, {"trk.has_hodo_hit", "Track has hodoscope hit", "fGEMTracks.THcLADGEMTrack.GetHasHodoHit()"}, - // {"trk.theta", "Track theta", "fGEMTracks.THcLADGEMTrack.GetTheta()"}, - // {"trk.phi", "Track phi", "fGEMTracks.THcLADGEMTrack.GetPhi()"}, + {"trk.chisq", "Track chi-squared", "fGEMTracks.THcLADGEMTrack.GetChisq()"}, + {"trk.theta", "Track theta", "fGEMTracks.THcLADGEMTrack.GetTheta()"}, + {"trk.phi", "Track phi", "fGEMTracks.THcLADGEMTrack.GetPhi()"}, {0}}; DefineVarsFromList(vars_trk, mode); } diff --git a/src/THcLADGEMTrack.cxx b/src/THcLADGEMTrack.cxx index 1a546e6..d0c0c25 100644 --- a/src/THcLADGEMTrack.cxx +++ b/src/THcLADGEMTrack.cxx @@ -5,10 +5,14 @@ THcLADGEMTrack::THcLADGEMTrack(Int_t nlayers) { fProjVz = -999.; fProjVy = -999.; + fProjVx = -999.; fNSp = 0; fD0 = -999; fhasGoodD0 = kFALSE; fSp = new GEM2DHits[nlayers]; + chisq = -999.; + ftheta = -999.; + fphi = -999.; } diff --git a/src/THcLADGEMTrack.h b/src/THcLADGEMTrack.h index 19e2424..fda593f 100644 --- a/src/THcLADGEMTrack.h +++ b/src/THcLADGEMTrack.h @@ -95,6 +95,7 @@ class THcLADGEMTrack : public TObject { // Track quantities Double_t GetProjVz() const { return fProjVz; } // projected z-vertex Double_t GetProjVy() const { return fProjVy; } // projected y-vertex + Double_t GetProjVx() const { return fProjVx; } // projected x-vertex Double_t GetD0() const { return fD0; } Bool_t GetGoodD0() const { return fhasGoodD0; } Double_t GetT() const { return fT; } @@ -107,6 +108,24 @@ class THcLADGEMTrack : public TObject { void SetGoodD0(Bool_t good) { fhasGoodD0 = good; } void SetZVertex(Double_t vz) { fProjVz = vz; } void SetYVertex(Double_t vy) { fProjVy = vy; } + void SetXVertex(Double_t vx) { fProjVx = vx; } + void SetProjVertex(Double_t vx, Double_t vy, Double_t vz) { + fProjVx = vx; + fProjVy = vy; + fProjVz = vz; + } + void SetAngles(Double_t theta, Double_t phi) { + ftheta = theta; + fphi = phi; + } + void GetAngles(Double_t &theta, Double_t &phi) const { + theta = ftheta; + phi = fphi; + } + Double_t GetTheta() const { return ftheta; } + Double_t GetPhi() const { return fphi; } + void SetChisq(Double_t c) { chisq = c; } + Double_t GetChisq() const { return chisq; } void SetTime(Double_t t, Double_t dt) { fT = t; fTdiff = dt; @@ -120,11 +139,13 @@ class THcLADGEMTrack : public TObject { int GetSpacePointID_0V() { return fSp[0].clusID[1]; } int GetSpacePointID_1U() { return fSp[1].clusID[0]; } int GetSpacePointID_1V() { return fSp[1].clusID[1]; } - + Double_t GetChi2() const { return chisq; } + void SetChi2(Double_t c) { chisq = c; } protected: Int_t fNSp; Double_t fProjVz; Double_t fProjVy; + Double_t fProjVx; Double_t fD0; Bool_t fhasGoodD0; Int_t fTrackID; @@ -133,6 +154,9 @@ class THcLADGEMTrack : public TObject { Double_t fT; Double_t fTdiff; bool fHasHodoHit; + Double_t chisq; + Double_t ftheta; // track angle between the track and z-azis in degrees + Double_t fphi; // track angle between the track and x-azis in degrees GEM2DHits *fSp; diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 6187c00..9ea8113 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -10,6 +10,10 @@ #include "TVector3.h" #include "VarDef.h" #include "VarType.h" +#include "Math/Minimizer.h" +#include "Math/Factory.h" +#include "Math/Functor.h" + ClassImp(THcLADKine) //_____________________________________________________________________________ @@ -208,7 +212,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } else { vertex.SetXYZ(0, 0, 0); } - + // Fix track vertex (for improved resolution on multifoils) if (fNfixed_z > 0 && track->GetGoodD0()) { std::vector distances; @@ -220,6 +224,27 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { vertex.SetZ(fFixed_z[min_index]); } + Double_t dir[3]; + dir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); + dir[1] = TMath::ATan((v_hit2.Y() - v_hit1.Y()) / (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); + dir[2] = vertex.Z(); + //Only fitting GEM hits for now, need to add LAD hits later + Double_t chisq=-kBig; + chisq= FitTrack(vertex, {v_hit1, v_hit2}, {0.1,0.1}, dir);//hardcoded resolution for now, should be parameterized based on detector performance + track->SetChisq(chisq);// set chisq even if the fit fails so we can see why + if(chisq<0) { + isGoodTrack[i] = false; + //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); + track->SetProjVertex(-kBig, -kBig, -kBig); + track->SetD0(-kBig); + continue; + } + track->SetAngles(dir[0], dir[1]); + track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); + double d0 = TMath::Abs(dir[2]-vertex.Z()); + track->SetD0(d0); + /* // Re-calculate d0 double numer = ((vertex - v_hit1).Cross((vertex - v_hit2))).Mag(); double denom = (v_hit2 - v_hit1).Mag(); @@ -234,6 +259,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { // Re-calculate vpy double vpy = v_hit1.Y() + (vertex.X() - v_hit1.X()) * (v_hit2.Y() - v_hit1.Y()) / (v_hit2.X() - v_hit1.X()); track->SetYVertex(vpy); + */ if (track->GetGoodD0()) { if (d0 < fD0Cut_wVertex) { @@ -550,3 +576,81 @@ Int_t THcLADKine::DefineVariables(EMode mode) { return kOK; } //_____________________________________________________________________________ + +Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_positions, std::vector sp_resolutions, double dir[3]){ + if (dir == nullptr || sp_positions.empty() || sp_resolutions.empty()) { + return -1; // Invalid input + } + int nPoints = sp_positions.size(); + if (nPoints < 2) { + return -1; // Not enough points to fit a track + } + if (nPoints != sp_resolutions.size()) { + return -2; // Mismatch in number of points and resolutions + } + for(int i = 0; i < nPoints; i++) { + if (sp_resolutions[i] <= 0) { + return -3; // Invalid resolution value + } + } + if (dir[0] < 0 || dir[0] > 180 || dir[1] < -180 || dir[1] > 180) { + return -4; // Invalid initial direction values + } + //requireing the track to originate from vertex (x,y), and only fitting for the track direction (theta, phi) and z vertex position. + double chi2 = -kBig; + ROOT::Math::Minimizer *minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); + minimizer->SetMaxFunctionCalls(1000); + minimizer->SetTolerance(1e-6); + ROOT::Math::Functor f([=](const double* params) { + double theta = params[0] * TMath::DegToRad(); // Convert to radians + double phi = params[1] * TMath::DegToRad(); // Convert to radians + double z = params[2]; + double chi2_local = 0; + for (int i = 0; i < nPoints; i++) { + //closest approach of the track to the point + double t = ((sp_positions[i].X() - vertex.X()) * TMath::Sin(theta) * TMath::Cos(phi) + + (sp_positions[i].Y() - vertex.Y()) * TMath::Sin(theta) * TMath::Sin(phi) + + (sp_positions[i].Z() - z) * TMath::Cos(theta)); + double x_closest = vertex.X() + t * TMath::Sin(theta) * TMath::Cos(phi); + double y_closest = vertex.Y() + t * TMath::Sin(theta) * TMath::Sin(phi); + double z_closest = z + t * TMath::Cos(theta); + double dx = sp_positions[i].X() - x_closest; + double dy = sp_positions[i].Y() - y_closest; + double dz = sp_positions[i].Z() - z_closest; + double dist2 = dx*dx +dz*dz; + if (i>2){ + // Assume these are hodoscope hit, so no chisq penalty if the hit is within the width of the paddle + double paddle_width = 22; // in cm, TODO: get actual paddle width from database + if (dist2 < (paddle_width/2.0)*(paddle_width/2.0)) { + dist2 = 0; + } + } + chi2_local += (dist2+ dy*dy) / (sp_resolutions[i]*sp_resolutions[i]); + } + return chi2_local; + },3); + minimizer->SetFunction(f); + + double initial_params[3] = {dir[0], dir[1], vertex.Z()}; + minimizer->SetVariable(0, "theta", initial_params[0], 0.1); + minimizer->SetVariableLimits(0, 45, 180);// Limit theta to be between 45 and 180 degrees to avoid unphysical solutions (tracks going backwards) + minimizer->SetVariable(1, "phi", initial_params[1], 0.1); + minimizer->SetVariableLimits(1, -90, 90);//Limit phi as the GEMs are only on the left side of the target in beam direction + minimizer->SetVariable(2, "z_vertex", initial_params[2], 0.1); + minimizer->SetVariableLimits(2, - 50, + 50); //the target length is 20cm, this should be more than wide enough + minimizer->Minimize(); + if (minimizer->Status() != 0) { + return -5; // Fit did not converge + } + const double *best_params = minimizer->X(); + dir[0] = best_params[0]; // theta in degree + dir[1] = best_params[1]; // phi in degree + dir[2] = best_params[2]; // z vertex position + chi2 = minimizer->MinValue(); + //clean up + delete minimizer; + + //return the chi2 of the fit + return chi2; +} + diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 2ed9ecb..ea19af4 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -60,6 +60,7 @@ class THcLADKine : public THcPrimaryKine { void CalculateTVertex(); Double_t CalculateToF(Double_t t_raw); Double_t CalculateTOFRFcorr(Double_t t_raw); + Double_t FitTrack(TVector3 vertex, std::vector sp_positions, std::vector sp_resolutions, double dir[3]); ClassDef(THcLADKine, 0) }; From bce9e733a497807b2eb08ea91d56c75e93a2c3fe Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Thu, 26 Mar 2026 15:55:22 -0400 Subject: [PATCH 02/15] use ATan2 to prevent divide by 0 --- src/THcLADGEMTrack.h | 2 -- src/THcLADKine.cxx | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/THcLADGEMTrack.h b/src/THcLADGEMTrack.h index fda593f..9d516f3 100644 --- a/src/THcLADGEMTrack.h +++ b/src/THcLADGEMTrack.h @@ -139,8 +139,6 @@ class THcLADGEMTrack : public TObject { int GetSpacePointID_0V() { return fSp[0].clusID[1]; } int GetSpacePointID_1U() { return fSp[1].clusID[0]; } int GetSpacePointID_1V() { return fSp[1].clusID[1]; } - Double_t GetChi2() const { return chisq; } - void SetChi2(Double_t c) { chisq = c; } protected: Int_t fNSp; Double_t fProjVz; diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 9ea8113..554fe82 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -582,7 +582,7 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position return -1; // Invalid input } int nPoints = sp_positions.size(); - if (nPoints < 2) { + if (nPoints < 2 || nPoints > 4) { return -1; // Not enough points to fit a track } if (nPoints != sp_resolutions.size()) { From 2fbcc9cd7a14ffc043c95db6c7538fcf4c4c99eb Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Thu, 26 Mar 2026 16:53:12 -0400 Subject: [PATCH 03/15] have the tracking limits in Db --- src/THcLADKine.cxx | 56 +++++++++++++++++++++++++++++++++++++--------- src/THcLADKine.h | 11 +++++++++ 2 files changed, 57 insertions(+), 10 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 554fe82..beb1736 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -31,6 +31,14 @@ ClassImp(THcLADKine) rf_offset = nullptr; n_rf_offsets = 0; rf_period = 4.00801; // Default RF period in ns + + fZCellMin = -15.0; // cm + fZCellMax = +15.0; // cm + fThetaMin = 60.0; // deg + fThetaMax = 170.0; // deg + fPhiMin = -50.0; // deg + fPhiMax = +50.0; // deg + } //_____________________________________________________________________________ THcLADKine::~THcLADKine() { @@ -183,6 +191,31 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { fFixed_z = nullptr; } + DBRequest track_constraints[] = {{ "ladkin.z_cell_min", &fZCellMin, kDouble, 0, 1 }, + { "ladkin.z_cell_max", &fZCellMax, kDouble, 0, 1 }, + { "ladkin.theta_min", &fThetaMin, kDouble, 0, 1 }, + { "ladkin.theta_max", &fThetaMax, kDouble, 0, 1 }, + { "ladkin.phi_min", &fPhiMin, kDouble, 0, 1 }, + { "ladkin.phi_max", &fPhiMax, kDouble, 0, 1 }, + { 0 } }; + gHcParms->LoadParmValues((DBRequest *)&track_constraints, prefix); + + if (fZCellMin > fZCellMax) { + Double_t tmp = fZCellMin; fZCellMin = fZCellMax; fZCellMax = tmp; + } + if (fThetaMin > fThetaMax) { + Double_t tmp = fThetaMin; fThetaMin = fThetaMax; fThetaMax = tmp; + } + if (fPhiMin > fPhiMax) { + Double_t tmp = fPhiMin; fPhiMin = fPhiMax; fPhiMax = tmp; + } + + fThetaMin = TMath::Max(0.0, TMath::Min(180.0, fThetaMin)); + fThetaMax = TMath::Max(0.0, TMath::Min(180.0, fThetaMax)); + fPhiMin = TMath::Max(-180.0,TMath::Min(180.0, fPhiMin)); + fPhiMax = TMath::Max(-180.0,TMath::Min(180.0, fPhiMax)); + + return err; } //_____________________________________________________________________________ @@ -226,7 +259,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { Double_t dir[3]; dir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); - dir[1] = TMath::ATan((v_hit2.Y() - v_hit1.Y()) / (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); + dir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); dir[2] = vertex.Z(); //Only fitting GEM hits for now, need to add LAD hits later Double_t chisq=-kBig; @@ -606,14 +639,17 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position double phi = params[1] * TMath::DegToRad(); // Convert to radians double z = params[2]; double chi2_local = 0; + double sx = TMath::Sin(theta) * TMath::Cos(phi); + double sy = TMath::Sin(theta) * TMath::Sin(phi); + double sz = TMath::Cos(theta); for (int i = 0; i < nPoints; i++) { //closest approach of the track to the point - double t = ((sp_positions[i].X() - vertex.X()) * TMath::Sin(theta) * TMath::Cos(phi) + - (sp_positions[i].Y() - vertex.Y()) * TMath::Sin(theta) * TMath::Sin(phi) + - (sp_positions[i].Z() - z) * TMath::Cos(theta)); - double x_closest = vertex.X() + t * TMath::Sin(theta) * TMath::Cos(phi); - double y_closest = vertex.Y() + t * TMath::Sin(theta) * TMath::Sin(phi); - double z_closest = z + t * TMath::Cos(theta); + double t = ((sp_positions[i].X() - vertex.X()) * sx + + (sp_positions[i].Y() - vertex.Y()) * sy + + (sp_positions[i].Z() - z) * sz); + double x_closest = vertex.X() + t * sx; + double y_closest = vertex.Y() + t * sy; + double z_closest = z + t * sz; double dx = sp_positions[i].X() - x_closest; double dy = sp_positions[i].Y() - y_closest; double dz = sp_positions[i].Z() - z_closest; @@ -633,11 +669,11 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position double initial_params[3] = {dir[0], dir[1], vertex.Z()}; minimizer->SetVariable(0, "theta", initial_params[0], 0.1); - minimizer->SetVariableLimits(0, 45, 180);// Limit theta to be between 45 and 180 degrees to avoid unphysical solutions (tracks going backwards) + minimizer->SetVariableLimits(0, fThetaMin, fThetaMax);// Limit theta to be between 45 and 180 degrees to avoid unphysical solutions (tracks going backwards) minimizer->SetVariable(1, "phi", initial_params[1], 0.1); - minimizer->SetVariableLimits(1, -90, 90);//Limit phi as the GEMs are only on the left side of the target in beam direction + minimizer->SetVariableLimits(1, fPhiMin, fPhiMax);//Limit phi as the GEMs are only on the left side of the target in beam direction minimizer->SetVariable(2, "z_vertex", initial_params[2], 0.1); - minimizer->SetVariableLimits(2, - 50, + 50); //the target length is 20cm, this should be more than wide enough + minimizer->SetVariableLimits(2, fZCellMin, fZCellMax); //the target length is 20cm, this should be more than wide enough minimizer->Minimize(); if (minimizer->Status() != 0) { return -5; // Fit did not converge diff --git a/src/THcLADKine.h b/src/THcLADKine.h index ea19af4..53a72ea 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -56,6 +56,17 @@ class THcLADKine : public THcPrimaryKine { Int_t n_rf_offsets; Double_t *rf_offset; Double_t rf_period; + + Double_t fZCellMin; // default -15.0 cm + Double_t fZCellMax; // default +15.0 cm + Double_t fThetaMin; // default 60.0 deg + Double_t fThetaMax; // default 170.0 deg + Double_t fPhiMin; // default -50.0 deg + Double_t fPhiMax; // default +50.0 deg + + + + virtual Int_t DefineVariables(EMode mode = kDefine); void CalculateTVertex(); Double_t CalculateToF(Double_t t_raw); From cc99d19110aef4825c35e9713ad78f783a650617 Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Thu, 26 Mar 2026 17:43:40 -0400 Subject: [PATCH 04/15] Add GetHitPositionLab in THcLADHodoscope --- src/THcLADHodoscope.cxx | 10 ++++++++++ src/THcLADHodoscope.h | 2 +- src/THcLADKine.cxx | 18 ++---------------- 3 files changed, 13 insertions(+), 17 deletions(-) diff --git a/src/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index 663c4f7..d5122fb 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -594,4 +594,14 @@ Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { return 0; } +TVector3 THcLADHodoscope::GetHitPosition(int plane, int paddle, double ypos=0) { + // Get the position of a hit in the hodoscope, given the plane and paddle number + // This is a placeholder function and should be implemented based on the actual geometry of the hodoscope + Double_t offset = fPlanes[plane]->GetPosCenter(paddle); + TVector3 hit = TVector3(offset, ypos + fPlanes[plane]->GetPosOffset(), fPlanes[plane]->GetZpos()); + hit->RotateY(fPlanes[plane]->GetTheta()); + return hit; +} + + ClassImp(THcLADHodoscope) diff --git a/src/THcLADHodoscope.h b/src/THcLADHodoscope.h index 3115abe..1adb3fa 100644 --- a/src/THcLADHodoscope.h +++ b/src/THcLADHodoscope.h @@ -61,7 +61,7 @@ class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { Double_t GetTDCThrs() const { return fTdc_Thrs; } Double_t GetEdep2MeV_int(Int_t iii) const { return fEdep2MeV_int[iii]; } Double_t GetEdep2MeV_amp(Int_t iii) const { return fEdep2MeV_amp[iii]; } - + TVector3 GetHitPositionLab(Int_t iplane, Int_t ipaddle, Double_t ypos=0); TClonesArray *GetLADGoodHits() { return fGoodLADHits; } TClonesArray *GetLADHits(Int_t plane) { return fPlanes[plane]->GetHits(); } diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index beb1736..c421793 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -277,22 +277,6 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); double d0 = TMath::Abs(dir[2]-vertex.Z()); track->SetD0(d0); - /* - // Re-calculate d0 - double numer = ((vertex - v_hit1).Cross((vertex - v_hit2))).Mag(); - double denom = (v_hit2 - v_hit1).Mag(); - // here we can put a range/fiducial cut on d0 taking into account the target size - double d0 = numer / denom; - track->SetD0(d0); - - // Re-calculate vpz - double vpz = v_hit1.Z() + (vertex.X() - v_hit1.X()) * (v_hit2.Z() - v_hit1.Z()) / (v_hit2.X() - v_hit1.X()); - track->SetZVertex(vpz); - - // Re-calculate vpy - double vpy = v_hit1.Y() + (vertex.X() - v_hit1.X()) * (v_hit2.Y() - v_hit1.Y()) / (v_hit2.X() - v_hit1.X()); - track->SetYVertex(vpy); - */ if (track->GetGoodD0()) { if (d0 < fD0Cut_wVertex) { @@ -327,6 +311,8 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { Int_t plane = hit->GetPlaneHit0(); Int_t paddle = hit->GetPaddleHit0(); + TVector3 hit_pos = fHodoscope->GetHitPositionLab(plane, paddle,hit->GetYposHit0()); + Int_t plane_loc; // Front vs back plane if (plane == 0 || plane == 2 || plane == 4) plane_loc = 0; From a911956371249200ab9531edcf2cf398043add58 Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 08:48:44 -0400 Subject: [PATCH 05/15] bug fixes --- src/THcLADHodoscope.cxx | 14 +++++++++++--- src/THcLADHodoscope.h | 2 +- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index d5122fb..b723c44 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -594,12 +594,20 @@ Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { return 0; } -TVector3 THcLADHodoscope::GetHitPosition(int plane, int paddle, double ypos=0) { +TVector3 THcLADHodoscope::GetHitPositionLab(int plane, int paddle, double ypos=0) { // Get the position of a hit in the hodoscope, given the plane and paddle number // This is a placeholder function and should be implemented based on the actual geometry of the hodoscope + if (plane < 0 || plane >= fNPlanes) { + cout << "[THcLADHodoscope] Error: Invalid plane number" << endl; + return TVector3(0, 0, 0); + } + if (paddle < 0 || paddle >= fNPaddle[plane]) { + cout << "[THcLADHodoscope] Error: Invalid paddle number" << endl; + return TVector3(0, 0, 0); + } Double_t offset = fPlanes[plane]->GetPosCenter(paddle); - TVector3 hit = TVector3(offset, ypos + fPlanes[plane]->GetPosOffset(), fPlanes[plane]->GetZpos()); - hit->RotateY(fPlanes[plane]->GetTheta()); + TVector3 hit = TVector3(offset, ypos , fPlanes[plane]->GetZpos()); + hit.RotateY(fPlanes[plane]->GetTheta()); return hit; } diff --git a/src/THcLADHodoscope.h b/src/THcLADHodoscope.h index 1adb3fa..10a2a9e 100644 --- a/src/THcLADHodoscope.h +++ b/src/THcLADHodoscope.h @@ -61,7 +61,7 @@ class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { Double_t GetTDCThrs() const { return fTdc_Thrs; } Double_t GetEdep2MeV_int(Int_t iii) const { return fEdep2MeV_int[iii]; } Double_t GetEdep2MeV_amp(Int_t iii) const { return fEdep2MeV_amp[iii]; } - TVector3 GetHitPositionLab(Int_t iplane, Int_t ipaddle, Double_t ypos=0); + TVector3 GetHitPositionLab(Int_t iplane, Int_t ipaddle, Double_t ypos = 0.0); TClonesArray *GetLADGoodHits() { return fGoodLADHits; } TClonesArray *GetLADHits(Int_t plane) { return fPlanes[plane]->GetHits(); } From 47006c3c58ae136e01478c12dcdf1b3da22dbc8b Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 12:35:44 -0400 Subject: [PATCH 06/15] re-write LADkine assuming hodo back/front matching would be done before --- src/THcLADGEMTrack.h | 4 ++ src/THcLADKine.cxx | 140 +++++++++++++++++++++++++++++++++---------- src/THcLADKine.h | 3 + 3 files changed, 114 insertions(+), 33 deletions(-) diff --git a/src/THcLADGEMTrack.h b/src/THcLADGEMTrack.h index 9d516f3..fe79666 100644 --- a/src/THcLADGEMTrack.h +++ b/src/THcLADGEMTrack.h @@ -6,6 +6,7 @@ #include "THcLADGEM.h" #include "TObject.h" #include "TVector3.h" +#include "THcGoodLADHit.h" class GEM2DHits { public: @@ -101,6 +102,8 @@ class THcLADGEMTrack : public TObject { Double_t GetT() const { return fT; } Double_t GetdT() const { return fTdiff; } bool GetHasHodoHit() const { return fHasHodoHit; } + THcGoodLADHit *GetBestHodoHit() const { return fBestHodoHit; } + void SetBestHodoHit(THcGoodLADHit *hit) { fBestHodoHit = hit; } void SetHasHodoHit(bool hasHodoHit) { fHasHodoHit = hasHodoHit; } void SetTrackID(int itrk) { fTrackID = itrk; } @@ -157,6 +160,7 @@ class THcLADGEMTrack : public TObject { Double_t fphi; // track angle between the track and x-azis in degrees GEM2DHits *fSp; + THcGoodLADHit *fBestHodoHit; // pointer to the best hodoscope hit associated with this track, if any private: THcLADGEMTrack(const THcLADGEMTrack &); diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index c421793..80dbeea 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -38,6 +38,8 @@ ClassImp(THcLADKine) fThetaMax = 170.0; // deg fPhiMin = -50.0; // deg fPhiMax = +50.0; // deg + fchisq_cut_withHodo = 30.0; // default chi2 cut for tracks with hodoscope hits + fchisq_cut_noHodo = 15.0; // default chi } //_____________________________________________________________________________ @@ -197,6 +199,8 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { { "ladkin.theta_max", &fThetaMax, kDouble, 0, 1 }, { "ladkin.phi_min", &fPhiMin, kDouble, 0, 1 }, { "ladkin.phi_max", &fPhiMax, kDouble, 0, 1 }, + { "ladkin.chisq_cut_withHodo", &fchisq_cut_withHodo, kDouble, 0, 1 }, + { "ladkin.chisq_cut_noHodo", &fchisq_cut_noHodo, kDouble, 0, 1 }, { 0 } }; gHcParms->LoadParmValues((DBRequest *)&track_constraints, prefix); @@ -224,6 +228,9 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { fTrack = fSpectro->GetGoldenTrack(); CalculateTVertex(); + //Get LADhodohits, + TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); + Int_t nHits = LADHits_unfiltered->GetLast() + 1; ////////////////////////////////////////////////////////////////////////////// // Check track projection to vertex fGEMTracks = fGEM->GetTracks(); @@ -257,26 +264,71 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { vertex.SetZ(fFixed_z[min_index]); } - Double_t dir[3]; - dir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); - dir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); - dir[2] = vertex.Z(); - //Only fitting GEM hits for now, need to add LAD hits later - Double_t chisq=-kBig; - chisq= FitTrack(vertex, {v_hit1, v_hit2}, {0.1,0.1}, dir);//hardcoded resolution for now, should be parameterized based on detector performance - track->SetChisq(chisq);// set chisq even if the fit fails so we can see why - if(chisq<0) { - isGoodTrack[i] = false; - //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); - track->SetProjVertex(-kBig, -kBig, -kBig); - track->SetD0(-kBig); - continue; + Double_t gemdir[3], best_gemdir[3]; + int bestHodoHitIndex=-1; + gemdir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); + gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); + gemdir[2] = vertex.Z(); + + Double_t bestchisq=kBig; + for (int iHodoHit=0; iHodoHit(LADHits_unfiltered->At(iHodoHit)); + if (hodo_hit == nullptr) + continue; + int plane = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetPlaneHit0() : hodo_hit->GetPlaneHit1(); + int paddle = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetPaddleHit0() : hodo_hit->GetPaddleHit1(); + double yPos = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetHitYPosHit0() : hodo_hit->GetHitYPosHit1(); + + TVector3 hodo_hitPos = fHodoscope->GetHitPosition(plane, paddle, yPos); + double res_hodo = 10; // cm, hardcoded for now, should be parameterized based on detector performance + Double_t dir[3]; + dir[0]=gemdir[0]; + dir[1]=gemdir[1]; + dir[2]=gemdir[2]; + double chisq = FitTrack(vertex, {v_hit1, v_hit2, v_hodo_hit}, {0.1, 0.1, res_hodo}, dir); + if (chisq>=0&&chisq=0&&bestchisqSetAngles(best_gemdir[0], best_gemdir[1]); + track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[2]); + track->SetChisq(bestchisq); + track->SetGoodD0(kTRUE); + track->SetD0(TMath::Abs(best_gemdir[2]-vertex.Z())*TMath::Sin(best_gemdir[0]*TMath::DegToRad())); + THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex)); + track->SetBestHodoHit(bestHodoHit); + track->SetHasBestHodoHit(kTRUE); + } + if(bestchisq<0||bestchisq>=fchisq_cut_withHodo) { + Double_t dir[3]; + dir[0]=gemdir[0]; + dir[1]=gemdir[1]; + dir[2]=gemdir[2]; + double chisq= FitTrack(vertex, {v_hit1, v_hit2}, {0.1,0.1}, dir);//hardcoded resolution for now, should be parameterized based on detector performance + track->SetChisq(chisq);// set chisq even if the fit fails so we can see why + if (chisq>=0&&chisqSetGoodD0(kTRUE); + isGoodTrack[i] = true; + track->SetAngles(dir[0], dir[1]); + track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); + track->SetD0(TMath::Abs(dir[2]-vertex.Z())*TMath::Sin(dir[0]*TMath::DegToRad())); + track->SetBestHodoHit(nullptr); + track->SetHasBestHodoHit(kFALSE); + }else { + isGoodTrack[i] = false; + //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); + track->SetProjVertex(-kBig, -kBig, -kBig); + track->SetD0(-kBig); + continue; + } } - track->SetAngles(dir[0], dir[1]); - track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); - double d0 = TMath::Abs(dir[2]-vertex.Z()); - track->SetD0(d0); if (track->GetGoodD0()) { if (d0 < fD0Cut_wVertex) { @@ -287,12 +339,14 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { isGoodTrack[i] = true; } } - + if(track->GetChisq()<0||track->GetChisq()>fchisq_cut_noHodo) { + isGoodTrack[i] = false; + } if (track->GetdT() > fTrk_dtCut) { isGoodTrack[i] = false; } } - +/* ////////////////////////////////////////////////////////////////////////////// // Remove hits with bad tracks and duplicate hits TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); @@ -304,14 +358,13 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { // Check if the hit has a track that points back to the origin. Int_t track_id = hit->GetTrackIDHit0(); - if (track_id < 0 || track_id >= ntracks || !isGoodTrack[track_id]) + if (track_id < 0 || track_id >= ntracks) continue; // Get the track parameters Int_t plane = hit->GetPlaneHit0(); Int_t paddle = hit->GetPaddleHit0(); - TVector3 hit_pos = fHodoscope->GetHitPositionLab(plane, paddle,hit->GetYposHit0()); Int_t plane_loc; // Front vs back plane if (plane == 0 || plane == 2 || plane == 4) @@ -329,7 +382,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { // Check if the hit is already in the list Int_t good_hit_plane = plane_loc ? goodhit->GetPlaneHit1() : goodhit->GetPlaneHit0(); Int_t good_hit_paddle = plane_loc ? goodhit->GetPaddleHit1() : goodhit->GetPaddleHit0(); - if (good_hit_plane == plane && good_hit_paddle == paddle) { + if (good_hit_plane == plane && good_hit_paddle == paddle) { //What if the same paddle is hit twice at different times? matching_hit_index = j; break; } @@ -357,7 +410,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { // Todo: Calculate beta, alpha, etc. for the hit // FIXME. One track can still have multiple (distinct) hodo hits that are counted as good } - +*/ ////////////////////////////////////////////////////////////////////////////// // Find matching front + back plane hits // Create a vector of pairs to store matching hits for each track @@ -615,6 +668,30 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position if (dir[0] < 0 || dir[0] > 180 || dir[1] < -180 || dir[1] > 180) { return -4; // Invalid initial direction values } + //check if the gem track is pointing to the hodoscope hits if there are any, return negative chisq if not + if (nPoints >2){ + double sx = TMath::Sin(dir[0] * TMath::DegToRad()) * TMath::Cos(dir[1] * TMath::DegToRad()); + double sy = TMath::Sin(dir[0] * TMath::DegToRad()) * TMath::Sin(dir[1] * TMath::DegToRad()); + double sz = TMath::Cos(dir[0] * TMath::DegToRad()); + for (int i = 2; i < nPoints; i++) { + //closest approach of the track to the point + double t = ((sp_positions[i].X() - sp_positions[0].X()) * sx + + (sp_positions[i].Y() - sp_positions[0].Y()) * sy + + (sp_positions[i].Z() - sp_positions[0].Z()) * sz); + double x_closest = sp_positions[0].X() + t * sx; + double y_closest = sp_positions[0].Y() + t * sy; + double z_closest = sp_positions[0].Z() + t * sz; + double dx = sp_positions[i].X() - x_closest; + double dy = sp_positions[i].Y() - y_closest; + double dz = sp_positions[i].Z() - z_closest; + double dist2 = dx*dx + dz*dz; + if (dist2 > (22*22)) { // if the track is more than 22cm away from the hodoscope hit, return negative chisq + return -5; // Track does not point to the hodoscope hit + } + } + } + + //requireing the track to originate from vertex (x,y), and only fitting for the track direction (theta, phi) and z vertex position. double chi2 = -kBig; ROOT::Math::Minimizer *minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); @@ -647,22 +724,19 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position dist2 = 0; } } - chi2_local += (dist2+ dy*dy) / (sp_resolutions[i]*sp_resolutions[i]); + chi2_local += (dist2+ dy*dy) / (sp_resolutions[i]*sp_resolutions[i]);//do we want to weight the vertical and horizontal residuals differently based on detector performance? } return chi2_local; },3); minimizer->SetFunction(f); double initial_params[3] = {dir[0], dir[1], vertex.Z()}; - minimizer->SetVariable(0, "theta", initial_params[0], 0.1); - minimizer->SetVariableLimits(0, fThetaMin, fThetaMax);// Limit theta to be between 45 and 180 degrees to avoid unphysical solutions (tracks going backwards) - minimizer->SetVariable(1, "phi", initial_params[1], 0.1); - minimizer->SetVariableLimits(1, fPhiMin, fPhiMax);//Limit phi as the GEMs are only on the left side of the target in beam direction - minimizer->SetVariable(2, "z_vertex", initial_params[2], 0.1); - minimizer->SetVariableLimits(2, fZCellMin, fZCellMax); //the target length is 20cm, this should be more than wide enough + minimizer->SetLimitedVariable(0, "theta", initial_params[0], 0.1, fThetaMin, fThetaMax);// Limit theta to avoid unphysical solutions (tracks going backwards) + minimizer->SetLimitedVariable(1, "phi", initial_params[1], 0.1, fPhiMin, fPhiMax);//Limit phi as the GEMs are only on the left side of the target in beam direction + minimizer->SetLimitedVariable(2, "z_vertex", initial_params[2], 0.1, fZCellMin, fZCellMax); //the target length is 20cm minimizer->Minimize(); if (minimizer->Status() != 0) { - return -5; // Fit did not converge + return -6; // Fit did not converge } const double *best_params = minimizer->X(); dir[0] = best_params[0]; // theta in degree diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 53a72ea..37b9c53 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -64,6 +64,9 @@ class THcLADKine : public THcPrimaryKine { Double_t fPhiMin; // default -50.0 deg Double_t fPhiMax; // default +50.0 deg + Double_t fchisq_cut_withHodo; + Double_t fchisq_cut_noHodo; + From 50a9394080f1ed4796b9ccbd665b850c7417abf3 Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 13:28:37 -0400 Subject: [PATCH 07/15] add lgic for multiple hodohit --- src/THcLADKine.cxx | 69 +++++++++++++++++++++++++++++----------------- src/THcLADKine.h | 3 ++ 2 files changed, 47 insertions(+), 25 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 80dbeea..1d50262 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -39,7 +39,10 @@ ClassImp(THcLADKine) fPhiMin = -50.0; // deg fPhiMax = +50.0; // deg fchisq_cut_withHodo = 30.0; // default chi2 cut for tracks with hodoscope hits - fchisq_cut_noHodo = 15.0; // default chi + fchisq_cut_noHodo = 15.0; // default chi2 cut for tracks without hodoscope hits + fSigma_GEM = 0.1; // default GEM resolution in cm + fSigma_Hodo = 10; // default Hodoscope resolution in cm + } //_____________________________________________________________________________ @@ -228,7 +231,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { fTrack = fSpectro->GetGoldenTrack(); CalculateTVertex(); - //Get LADhodohits, + //Get LADGoodHits for track matching TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); Int_t nHits = LADHits_unfiltered->GetLast() + 1; ////////////////////////////////////////////////////////////////////////////// @@ -268,31 +271,47 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { int bestHodoHitIndex=-1; gemdir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); - gemdir[2] = vertex.Z(); + gemdir[2] = vertex.Z(); //Hopefully the GEM track will be the same as elecctron vertex Z Double_t bestchisq=kBig; + //loop over hodoscope hits to find the best match to the track projection for (int iHodoHit=0; iHodoHit(LADHits_unfiltered->At(iHodoHit)); if (hodo_hit == nullptr) continue; - int plane = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetPlaneHit0() : hodo_hit->GetPlaneHit1(); - int paddle = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetPaddleHit0() : hodo_hit->GetPaddleHit1(); - double yPos = hodo_hit->GetPlaneHit0() >= 0 ? hodo_hit->GetHitYPosHit0() : hodo_hit->GetHitYPosHit1(); - - TVector3 hodo_hitPos = fHodoscope->GetHitPosition(plane, paddle, yPos); - double res_hodo = 10; // cm, hardcoded for now, should be parameterized based on detector performance - Double_t dir[3]; - dir[0]=gemdir[0]; - dir[1]=gemdir[1]; - dir[2]=gemdir[2]; - double chisq = FitTrack(vertex, {v_hit1, v_hit2, v_hodo_hit}, {0.1, 0.1, res_hodo}, dir); - if (chisq>=0&&chisq v_hits, v_resolutions; + v_hits.push_back(v_hit1); + v_hits.push_back(v_hit2); + v_resolutions.push_back(fSigma_GEM); + v_resolutions.push_back(fSigma_GEM); + + TVector3 v_hodo_hit; + if (hodo_hit->GetPlaneHit0() >= 0) { + v_hodo_hit = fHodoscope->GetHitPosition(hodo_hit->GetPlaneHit0(), hodo_hit->GetPaddleHit0(), hodo_hit->GetHitYPosHit0()); + v_hits.push_back(v_hodo_hit); + v_resolutions.push_back(fSigma_Hodo); + } + if (hodo_hit->GetPlaneHit1() >= 0) { + v_hodo_hit = fHodoscope->GetHitPosition(hodo_hit->GetPlaneHit1(), hodo_hit->GetPaddleHit1(), hodo_hit->GetHitYPosHit1()); + v_hits.push_back(v_hodo_hit); + v_resolutions.push_back(fSigma_Hodo); + } + Double_t dir[3]; + dir[0]=gemdir[0]; + dir[1]=gemdir[1]; + dir[2]=gemdir[2]; + double chisq = FitTrack(vertex, v_hits, v_resolutions, dir); + // Check if this is the best match so far, we only want to keep the best hodoscope hit for each track to avoid double counting hits + if (chisq>=0&&chisq=0&&bestchisq(LADHits_unfiltered->At(bestHodoHitIndex)); track->SetBestHodoHit(bestHodoHit); track->SetHasBestHodoHit(kTRUE); - } - if(bestchisq<0||bestchisq>=fchisq_cut_withHodo) { + }else { Double_t dir[3]; dir[0]=gemdir[0]; dir[1]=gemdir[1]; - dir[2]=gemdir[2]; - double chisq= FitTrack(vertex, {v_hit1, v_hit2}, {0.1,0.1}, dir);//hardcoded resolution for now, should be parameterized based on detector performance + dir[2]=gemdir[2];// if the track doesn't have a good hodoscope match, just fit to the GEM hits and vertex + double chisq= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); track->SetChisq(chisq);// set chisq even if the fit fails so we can see why if (chisq>=0&&chisqSetGoodD0(kTRUE); @@ -321,6 +339,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { track->SetBestHodoHit(nullptr); track->SetHasBestHodoHit(kFALSE); }else { + //Still bad track, mark as bad and set parameters to default values isGoodTrack[i] = false; //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 37b9c53..2a5de90 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -67,6 +67,9 @@ class THcLADKine : public THcPrimaryKine { Double_t fchisq_cut_withHodo; Double_t fchisq_cut_noHodo; + Double_t fSigma_GEM; // GEM resolution in cm, used for track fitting, should be set based on detector performance + Double_t fSigma_Hodo; // Hodoscope resolution in cm, used for track fitting, should be set based on detector performance + From f2635993e199051821c25f4de9f7e5c22f8d7e7b Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 14:55:05 -0400 Subject: [PATCH 08/15] bug fixes --- src/THcLADGEMTrack.h | 6 +- src/THcLADHodoscope.cxx | 2 +- src/THcLADHodoscope.h | 2 +- src/THcLADKine.cxx | 234 ++++++++++++---------------------------- src/THcLADKine.h | 3 +- 5 files changed, 75 insertions(+), 172 deletions(-) diff --git a/src/THcLADGEMTrack.h b/src/THcLADGEMTrack.h index fe79666..8caf042 100644 --- a/src/THcLADGEMTrack.h +++ b/src/THcLADGEMTrack.h @@ -101,11 +101,11 @@ class THcLADGEMTrack : public TObject { Bool_t GetGoodD0() const { return fhasGoodD0; } Double_t GetT() const { return fT; } Double_t GetdT() const { return fTdiff; } - bool GetHasHodoHit() const { return fHasHodoHit; } + short GetHasHodoHit() const { return fHasHodoHit; } THcGoodLADHit *GetBestHodoHit() const { return fBestHodoHit; } void SetBestHodoHit(THcGoodLADHit *hit) { fBestHodoHit = hit; } - void SetHasHodoHit(bool hasHodoHit) { fHasHodoHit = hasHodoHit; } + void SetHasHodoHit(short hasHodoHit) { fHasHodoHit = hasHodoHit; } void SetTrackID(int itrk) { fTrackID = itrk; } void SetD0(Double_t d0) { fD0 = d0; } void SetGoodD0(Bool_t good) { fhasGoodD0 = good; } @@ -154,7 +154,7 @@ class THcLADGEMTrack : public TObject { Int_t fClustID1; Double_t fT; Double_t fTdiff; - bool fHasHodoHit; + short fHasHodoHit; Double_t chisq; Double_t ftheta; // track angle between the track and z-azis in degrees Double_t fphi; // track angle between the track and x-azis in degrees diff --git a/src/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index b723c44..76561a3 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -594,7 +594,7 @@ Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { return 0; } -TVector3 THcLADHodoscope::GetHitPositionLab(int plane, int paddle, double ypos=0) { +TVector3 THcLADHodoscope::GetHitPositionLab(int plane, int paddle, double ypos) { // Get the position of a hit in the hodoscope, given the plane and paddle number // This is a placeholder function and should be implemented based on the actual geometry of the hodoscope if (plane < 0 || plane >= fNPlanes) { diff --git a/src/THcLADHodoscope.h b/src/THcLADHodoscope.h index 10a2a9e..7ff9a9d 100644 --- a/src/THcLADHodoscope.h +++ b/src/THcLADHodoscope.h @@ -61,7 +61,7 @@ class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { Double_t GetTDCThrs() const { return fTdc_Thrs; } Double_t GetEdep2MeV_int(Int_t iii) const { return fEdep2MeV_int[iii]; } Double_t GetEdep2MeV_amp(Int_t iii) const { return fEdep2MeV_amp[iii]; } - TVector3 GetHitPositionLab(Int_t iplane, Int_t ipaddle, Double_t ypos = 0.0); + TVector3 GetHitPositionLab(Int_t iplane, Int_t ipaddle, Double_t ypos ); TClonesArray *GetLADGoodHits() { return fGoodLADHits; } TClonesArray *GetLADHits(Int_t plane) { return fPlanes[plane]->GetHits(); } diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 1d50262..2b220c8 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -38,8 +38,9 @@ ClassImp(THcLADKine) fThetaMax = 170.0; // deg fPhiMin = -50.0; // deg fPhiMax = +50.0; // deg - fchisq_cut_withHodo = 30.0; // default chi2 cut for tracks with hodoscope hits - fchisq_cut_noHodo = 15.0; // default chi2 cut for tracks without hodoscope hits + fchisq_cut[0] = 20.0; // default chi2 cut for tracks with 2 hodo hits + fchisq_cut[1] = 15.0; // default chi2 cut for tracks with 1 hodo hit + fchisq_cut[2] = 10.0; // default chi2 cut for tracks with no hodo hits fSigma_GEM = 0.1; // default GEM resolution in cm fSigma_Hodo = 10; // default Hodoscope resolution in cm @@ -202,8 +203,7 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { { "ladkin.theta_max", &fThetaMax, kDouble, 0, 1 }, { "ladkin.phi_min", &fPhiMin, kDouble, 0, 1 }, { "ladkin.phi_max", &fPhiMax, kDouble, 0, 1 }, - { "ladkin.chisq_cut_withHodo", &fchisq_cut_withHodo, kDouble, 0, 1 }, - { "ladkin.chisq_cut_noHodo", &fchisq_cut_noHodo, kDouble, 0, 1 }, + { "ladkin.chisq_cut", fchisq_cut, kDouble, 3, 1 }, { 0 } }; gHcParms->LoadParmValues((DBRequest *)&track_constraints, prefix); @@ -234,13 +234,14 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { //Get LADGoodHits for track matching TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); Int_t nHits = LADHits_unfiltered->GetLast() + 1; + goodhit_n = 0; ////////////////////////////////////////////////////////////////////////////// // Check track projection to vertex fGEMTracks = fGEM->GetTracks(); Int_t ntracks = fGEMTracks->GetLast() + 1; std::vector isGoodTrack(ntracks, false); - for (Int_t i = 0; i < ntracks; i++) { - THcLADGEMTrack *track = static_cast(fGEMTracks->At(i)); + for (Int_t iTrack = 0; iTrack < ntracks; iTrack++) { + THcLADGEMTrack *track = static_cast(fGEMTracks->At(iTrack)); if (track == nullptr) continue; @@ -267,228 +268,131 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { vertex.SetZ(fFixed_z[min_index]); } - Double_t gemdir[3], best_gemdir[3]; - int bestHodoHitIndex=-1; + Double_t gemdir[3]; gemdir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); gemdir[2] = vertex.Z(); //Hopefully the GEM track will be the same as elecctron vertex Z - Double_t bestchisq=kBig; + Double_t bestchisq[2]={kBig,kBig};// 2 hodo hits, 1 hodo hit + Double_t best_gemdir[2][3]; + int bestHodoHitIndex[2]={-1,-1}; //loop over hodoscope hits to find the best match to the track projection for (int iHodoHit=0; iHodoHit(LADHits_unfiltered->At(iHodoHit)); if (hodo_hit == nullptr) continue; // There might be one or two hits with a signle goodhit, need to check both - std::vector v_hits, v_resolutions; + std::vector v_hits; + std::vector v_resolutions; v_hits.push_back(v_hit1); v_hits.push_back(v_hit2); v_resolutions.push_back(fSigma_GEM); v_resolutions.push_back(fSigma_GEM); - + int num_hodo_hits = 0; TVector3 v_hodo_hit; if (hodo_hit->GetPlaneHit0() >= 0) { - v_hodo_hit = fHodoscope->GetHitPosition(hodo_hit->GetPlaneHit0(), hodo_hit->GetPaddleHit0(), hodo_hit->GetHitYPosHit0()); + v_hodo_hit = fHodoscope->GetHitPositionLab(hodo_hit->GetPlaneHit0(), hodo_hit->GetPaddleHit0(), hodo_hit->GetHitYPosHit0()); v_hits.push_back(v_hodo_hit); v_resolutions.push_back(fSigma_Hodo); + num_hodo_hits++; } if (hodo_hit->GetPlaneHit1() >= 0) { - v_hodo_hit = fHodoscope->GetHitPosition(hodo_hit->GetPlaneHit1(), hodo_hit->GetPaddleHit1(), hodo_hit->GetHitYPosHit1()); + v_hodo_hit = fHodoscope->GetHitPositionLab(hodo_hit->GetPlaneHit1(), hodo_hit->GetPaddleHit1(), hodo_hit->GetHitYPosHit1()); v_hits.push_back(v_hodo_hit); v_resolutions.push_back(fSigma_Hodo); + num_hodo_hits++; } Double_t dir[3]; dir[0]=gemdir[0]; dir[1]=gemdir[1]; dir[2]=gemdir[2]; - double chisq = FitTrack(vertex, v_hits, v_resolutions, dir); + Double_t chisq = FitTrack(vertex, v_hits, v_resolutions, dir); // Check if this is the best match so far, we only want to keep the best hodoscope hit for each track to avoid double counting hits - if (chisq>=0&&chisq=0&&chisq=0&&bestchisqSetAngles(best_gemdir[0], best_gemdir[1]); - track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[2]); - track->SetChisq(bestchisq); + if(bestchisq[0]>=0&&bestchisq[0]SetAngles(best_gemdir[0][0], best_gemdir[0][1]); + track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[0][2]); + track->SetChisq(bestchisq[0]); track->SetGoodD0(kTRUE); - track->SetD0(TMath::Abs(best_gemdir[2]-vertex.Z())*TMath::Sin(best_gemdir[0]*TMath::DegToRad())); - THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex)); + track->SetD0(TMath::Abs(best_gemdir[0][2]-vertex.Z())*TMath::Sin(best_gemdir[0][0]*TMath::DegToRad())); + THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[0])); track->SetBestHodoHit(bestHodoHit); - track->SetHasBestHodoHit(kTRUE); + track->SetHasHodoHit(2); + THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); + goodhit->CopyHit(0,0,bestHodoHit); + goodhit->CopyHit(1,1,bestHodoHit); + goodhit->SetTrackID(0, iTrack); + goodhit->SetTrackID(1, iTrack); + goodhit_n++; + }else if(bestchisq[1]>=0&&bestchisq[1]SetAngles(best_gemdir[1][0], best_gemdir[1][1]); + track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[1][2]); + track->SetChisq(bestchisq[1]); + track->SetGoodD0(kTRUE); + track->SetD0(TMath::Abs(best_gemdir[1][2]-vertex.Z())*TMath::Sin(best_gemdir[1][0]*TMath::DegToRad())); + THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[1])); + track->SetBestHodoHit(bestHodoHit); + track->SetHasHodoHit(1); + THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); + goodhit->CopyHit(0,0,bestHodoHit); + goodhit->CopyHit(1,1,bestHodoHit); + goodhit->SetTrackID(0, iTrack); + goodhit->SetTrackID(1, iTrack); + goodhit_n++; }else { Double_t dir[3]; dir[0]=gemdir[0]; dir[1]=gemdir[1]; dir[2]=gemdir[2];// if the track doesn't have a good hodoscope match, just fit to the GEM hits and vertex - double chisq= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); + Double_t chisq= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); track->SetChisq(chisq);// set chisq even if the fit fails so we can see why - if (chisq>=0&&chisq=0&&chisqSetGoodD0(kTRUE); - isGoodTrack[i] = true; + isGoodTrack[iTrack] = true; track->SetAngles(dir[0], dir[1]); track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); track->SetD0(TMath::Abs(dir[2]-vertex.Z())*TMath::Sin(dir[0]*TMath::DegToRad())); track->SetBestHodoHit(nullptr); - track->SetHasBestHodoHit(kFALSE); + track->SetHasHodoHit(0); }else { //Still bad track, mark as bad and set parameters to default values - isGoodTrack[i] = false; - //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); track->SetProjVertex(-kBig, -kBig, -kBig); track->SetD0(-kBig); + track->SetGoodD0(kFALSE); continue; } } if (track->GetGoodD0()) { - if (d0 < fD0Cut_wVertex) { - isGoodTrack[i] = true; + if (track->GetD0()>0 && track->GetD0() < fD0Cut_wVertex) { + isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; } } else { - if (d0 < fD0Cut_noVertex) { - isGoodTrack[i] = true; + if (track->GetD0() > 0 && track->GetD0() < fD0Cut_noVertex) { + isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; } } - if(track->GetChisq()<0||track->GetChisq()>fchisq_cut_noHodo) { - isGoodTrack[i] = false; + if(track->GetChisq()<0) { + isGoodTrack[iTrack] = false; } if (track->GetdT() > fTrk_dtCut) { - isGoodTrack[i] = false; + isGoodTrack[iTrack] = false; } } -/* - ////////////////////////////////////////////////////////////////////////////// - // Remove hits with bad tracks and duplicate hits - TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); - Int_t nHits = LADHits_unfiltered->GetLast() + 1; - for (Int_t i = 0; i < nHits; i++) { - THcGoodLADHit *hit = static_cast(LADHits_unfiltered->At(i)); - if (hit == nullptr) - continue; - - // Check if the hit has a track that points back to the origin. - Int_t track_id = hit->GetTrackIDHit0(); - if (track_id < 0 || track_id >= ntracks) - continue; - - // Get the track parameters - Int_t plane = hit->GetPlaneHit0(); - Int_t paddle = hit->GetPaddleHit0(); - - - Int_t plane_loc; // Front vs back plane - if (plane == 0 || plane == 2 || plane == 4) - plane_loc = 0; - else - plane_loc = 1; - - Int_t matching_hit_index = -1; - Int_t partner_hit_index = -1; - for (Int_t j = 0; j < goodhit_n; j++) { - THcGoodLADHit *goodhit = static_cast(fGoodLADHits->At(j)); - if (goodhit == nullptr) - continue; - - // Check if the hit is already in the list - Int_t good_hit_plane = plane_loc ? goodhit->GetPlaneHit1() : goodhit->GetPlaneHit0(); - Int_t good_hit_paddle = plane_loc ? goodhit->GetPaddleHit1() : goodhit->GetPaddleHit0(); - if (good_hit_plane == plane && good_hit_paddle == paddle) { //What if the same paddle is hit twice at different times? - matching_hit_index = j; - break; - } - // Check if the hit has a partner hit - Int_t partner_hit_track_id = plane_loc ? goodhit->GetTrackIDHit1() : goodhit->GetTrackIDHit0(); - if (partner_hit_track_id == track_id) { - partner_hit_index = j; - } - } - if (matching_hit_index == -1 && partner_hit_index == -1) { - THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); - goodhit_n++; - goodhit->CopyHit(plane_loc, 0, hit); // Copy the new hit into the good hit - goodhit->SetTrackID(plane_loc, track_id); - } - if (matching_hit_index != -1) { - THcGoodLADHit *goodhit = static_cast(fGoodLADHits->At(matching_hit_index)); - if (abs(hit->GetdTrkHorizHit0()) < - (plane_loc ? abs(goodhit->GetdTrkHorizHit1()) : abs(goodhit->GetdTrkHorizHit0()))) { - - goodhit->CopyHit(plane_loc, 0, hit); // Copy the new hit into the good hit - goodhit->SetTrackID(plane_loc, track_id); - } - } - // Todo: Calculate beta, alpha, etc. for the hit - // FIXME. One track can still have multiple (distinct) hodo hits that are counted as good - } -*/ - ////////////////////////////////////////////////////////////////////////////// - // Find matching front + back plane hits - // Create a vector of pairs to store matching hits for each track - std::vector> matchingHits(ntracks, {-1, -1}); - std::vector dTrk_horiz_values; - - // TODO: This can probably be incorporated into the loop above - for (Int_t i = 0; i < goodhit_n; i++) { - THcGoodLADHit *goodhit = static_cast(fGoodLADHits->At(i)); - if (goodhit == nullptr) - continue; - - // Add dTrk_horiz values for both planes (if valid) to the vector - if (goodhit->GetdTrkHorizHit0() != 0) { - dTrk_horiz_values.push_back(goodhit->GetdTrkHorizHit0()); - } - if (goodhit->GetdTrkHorizHit1() != 0) { - dTrk_horiz_values.push_back(goodhit->GetdTrkHorizHit1()); - } - - Int_t trackID0 = goodhit->GetTrackIDHit0(); - Int_t trackID1 = goodhit->GetTrackIDHit1(); - - if (trackID0 >= 0 && trackID0 < ntracks) { - // Check if the track ID is already in the vector - if (matchingHits[trackID0].first == -1) { - // If it is, set the second hit index to the current index - matchingHits[trackID0].first = i; - } else if (dTrk_horiz_values[matchingHits[trackID0].first] > goodhit->GetdTrkHorizHit0()) { - matchingHits[trackID0].first = i; - } - } - if (trackID1 >= 0 && trackID1 < ntracks) { - // Check if the track ID is already in the vector - if (matchingHits[trackID1].second == -1) { - // If it is, set the second hit index to the current index - matchingHits[trackID1].second = i; - } else if (dTrk_horiz_values[matchingHits[trackID1].second] > goodhit->GetdTrkHorizHit1()) { - matchingHits[trackID1].second = i; - } - } - } - - // Loop over the matching hits and set the partner hit index - for (auto &match : matchingHits) { - if (match.first != -1 && match.second != -1) { - THcGoodLADHit *firstHit = static_cast(fGoodLADHits->At(match.first)); - THcGoodLADHit *secondHit = static_cast(fGoodLADHits->At(match.second)); - if (firstHit && secondHit) { - firstHit->CopyHit(1, 0, secondHit); // Copy the back (1) plane of the second hit into the first hit (0) - - // Remove the second hit from the good hits array - fGoodLADHits->RemoveAt(match.second); - goodhit_n--; - } - } - } - fGoodLADHits->Compress(); // Compress the array to remove null entries - // Calculate beta, alpha, tof, etc. for the hits for (Int_t i = 0; i < goodhit_n; i++) { THcGoodLADHit *goodhit = static_cast(fGoodLADHits->At(i)); diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 2a5de90..b2dbe49 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -64,8 +64,7 @@ class THcLADKine : public THcPrimaryKine { Double_t fPhiMin; // default -50.0 deg Double_t fPhiMax; // default +50.0 deg - Double_t fchisq_cut_withHodo; - Double_t fchisq_cut_noHodo; + Double_t fchisq_cut[3];//chisq cuts for track, 2 hodo hits, 1 hodo hit, and no hodo hits. Double_t fSigma_GEM; // GEM resolution in cm, used for track fitting, should be set based on detector performance Double_t fSigma_Hodo; // Hodoscope resolution in cm, used for track fitting, should be set based on detector performance From a620d88b339ef565ad5ddfe4d69f8f3df4ce1a07 Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 16:54:06 -0400 Subject: [PATCH 09/15] use use radians --- src/THcLADGEM.cxx | 1 + src/THcLADGEMTrack.h | 9 ++++++--- src/THcLADHodoPlane.cxx | 4 ++-- src/THcLADKine.cxx | 30 +++++++++++++++++------------- 4 files changed, 26 insertions(+), 18 deletions(-) diff --git a/src/THcLADGEM.cxx b/src/THcLADGEM.cxx index a164365..b926d02 100644 --- a/src/THcLADGEM.cxx +++ b/src/THcLADGEM.cxx @@ -236,6 +236,7 @@ Int_t THcLADGEM::DefineVariables(EMode mode) { {"trk.chisq", "Track chi-squared", "fGEMTracks.THcLADGEMTrack.GetChisq()"}, {"trk.theta", "Track theta", "fGEMTracks.THcLADGEMTrack.GetTheta()"}, {"trk.phi", "Track phi", "fGEMTracks.THcLADGEMTrack.GetPhi()"}, + {"trk.is_good", "Track is good", "fGEMTracks.THcLADGEMTrack.IsGoodTrack()"}, {0}}; DefineVarsFromList(vars_trk, mode); } diff --git a/src/THcLADGEMTrack.h b/src/THcLADGEMTrack.h index 8caf042..dd8c902 100644 --- a/src/THcLADGEMTrack.h +++ b/src/THcLADGEMTrack.h @@ -134,6 +134,9 @@ class THcLADGEMTrack : public TObject { fTdiff = dt; } + bool IsGoodTrack() const { return fIsGoodTrack; } + void SetIsGoodTrack(bool isGood) { fIsGoodTrack = isGood; } + GEM2DHits GetSpacePoint(int isp) { return fSp[isp]; } virtual void AddSpacePoint(GEM2DHits sp); Int_t GetSpID_0() const {return fSp[0].spID;}; @@ -156,9 +159,9 @@ class THcLADGEMTrack : public TObject { Double_t fTdiff; short fHasHodoHit; Double_t chisq; - Double_t ftheta; // track angle between the track and z-azis in degrees - Double_t fphi; // track angle between the track and x-azis in degrees - + Double_t ftheta; // track angle between the track and z-azis in radians + Double_t fphi; // track angle between the track and x-azis in radians + bool fIsGoodTrack; GEM2DHits *fSp; THcGoodLADHit *fBestHodoHit; // pointer to the best hodoscope hit associated with this track, if any diff --git a/src/THcLADHodoPlane.cxx b/src/THcLADHodoPlane.cxx index 7ea3440..3ab4e3c 100644 --- a/src/THcLADHodoPlane.cxx +++ b/src/THcLADHodoPlane.cxx @@ -1601,8 +1601,8 @@ Int_t THcLADHodoPlane::ProcessHits(TClonesArray *rawhits, Int_t nexthit) { ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetCorrectedTimes_FADC(adc_timec_top, adc_timec_btm, adc_time_corrected); ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetTopADCpeak(adcamp_top[i_good_top_adc_elem]); ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetBtmADCpeak(adcamp_btm[i_good_btm_adc_elem]); - ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetCalcPosition(fHitDistCorr); - ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetCalcPosition_FADC(fHitDistCorr_FADC); + ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetCalcPosition(fHitDistCorr+fPosOffset);//add yposition offset to match LAD geometry + ((THcLADHodoHit *)fHodoHits->At(fNScinHits))->SetCalcPosition_FADC(fHitDistCorr_FADC+fPosOffset); fGoodTopTdcTimeCorr.at(padnum - 1) = timec_top; fGoodBtmTdcTimeCorr.at(padnum - 1) = timec_btm; diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 2b220c8..c4d80a9 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -269,8 +269,8 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } Double_t gemdir[3]; - gemdir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) * TMath::RadToDeg(); - gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())) * TMath::RadToDeg(); + gemdir[0] = TMath::ACos((v_hit2.Z() - v_hit1.Z()) / (v_hit2 - v_hit1).Mag()) ; + gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())); gemdir[2] = vertex.Z(); //Hopefully the GEM track will be the same as elecctron vertex Z Double_t bestchisq[2]={kBig,kBig};// 2 hodo hits, 1 hodo hit @@ -324,7 +324,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[0][2]); track->SetChisq(bestchisq[0]); track->SetGoodD0(kTRUE); - track->SetD0(TMath::Abs(best_gemdir[0][2]-vertex.Z())*TMath::Sin(best_gemdir[0][0]*TMath::DegToRad())); + track->SetD0(TMath::Abs(best_gemdir[0][2]-vertex.Z())*TMath::Sin(best_gemdir[0][0])); THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[0])); track->SetBestHodoHit(bestHodoHit); track->SetHasHodoHit(2); @@ -340,7 +340,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[1][2]); track->SetChisq(bestchisq[1]); track->SetGoodD0(kTRUE); - track->SetD0(TMath::Abs(best_gemdir[1][2]-vertex.Z())*TMath::Sin(best_gemdir[1][0]*TMath::DegToRad())); + track->SetD0(TMath::Abs(best_gemdir[1][2]-vertex.Z())*TMath::Sin(best_gemdir[1][0])); THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[1])); track->SetBestHodoHit(bestHodoHit); track->SetHasHodoHit(1); @@ -362,7 +362,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { isGoodTrack[iTrack] = true; track->SetAngles(dir[0], dir[1]); track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); - track->SetD0(TMath::Abs(dir[2]-vertex.Z())*TMath::Sin(dir[0]*TMath::DegToRad())); + track->SetD0(TMath::Abs(dir[2]-vertex.Z())*TMath::Sin(dir[0])); track->SetBestHodoHit(nullptr); track->SetHasHodoHit(0); }else { @@ -373,6 +373,9 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { track->SetProjVertex(-kBig, -kBig, -kBig); track->SetD0(-kBig); track->SetGoodD0(kFALSE); + track->SetHasHodoHit(0); + track->SetBestHodoHit(nullptr); + track->SetIsGoodTrack(false); continue; } } @@ -392,6 +395,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { if (track->GetdT() > fTrk_dtCut) { isGoodTrack[iTrack] = false; } + track->SetIsGoodTrack(isGoodTrack[iTrack]); } // Calculate beta, alpha, tof, etc. for the hits for (Int_t i = 0; i < goodhit_n; i++) { @@ -588,14 +592,14 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position return -3; // Invalid resolution value } } - if (dir[0] < 0 || dir[0] > 180 || dir[1] < -180 || dir[1] > 180) { + if (dir[0] < 0 || dir[0] > TMath::Pi() || dir[1] < -TMath::Pi() || dir[1] > TMath::Pi()) { return -4; // Invalid initial direction values } //check if the gem track is pointing to the hodoscope hits if there are any, return negative chisq if not if (nPoints >2){ - double sx = TMath::Sin(dir[0] * TMath::DegToRad()) * TMath::Cos(dir[1] * TMath::DegToRad()); - double sy = TMath::Sin(dir[0] * TMath::DegToRad()) * TMath::Sin(dir[1] * TMath::DegToRad()); - double sz = TMath::Cos(dir[0] * TMath::DegToRad()); + double sx = TMath::Sin(dir[0]) * TMath::Cos(dir[1]); + double sy = TMath::Sin(dir[0]) * TMath::Sin(dir[1]); + double sz = TMath::Cos(dir[0]); for (int i = 2; i < nPoints; i++) { //closest approach of the track to the point double t = ((sp_positions[i].X() - sp_positions[0].X()) * sx + @@ -621,8 +625,8 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position minimizer->SetMaxFunctionCalls(1000); minimizer->SetTolerance(1e-6); ROOT::Math::Functor f([=](const double* params) { - double theta = params[0] * TMath::DegToRad(); // Convert to radians - double phi = params[1] * TMath::DegToRad(); // Convert to radians + double theta = params[0]; // in radians + double phi = params[1]; // in radians double z = params[2]; double chi2_local = 0; double sx = TMath::Sin(theta) * TMath::Cos(phi); @@ -662,8 +666,8 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position return -6; // Fit did not converge } const double *best_params = minimizer->X(); - dir[0] = best_params[0]; // theta in degree - dir[1] = best_params[1]; // phi in degree + dir[0] = best_params[0]; // theta in radians + dir[1] = best_params[1]; // phi in radians dir[2] = best_params[2]; // z vertex position chi2 = minimizer->MinValue(); //clean up From ea5eee954be3874332475ef39bc5066b2d15948a Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 16:55:19 -0400 Subject: [PATCH 10/15] initializing limits with radians as well --- src/THcLADKine.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index c4d80a9..c810950 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -34,10 +34,10 @@ ClassImp(THcLADKine) fZCellMin = -15.0; // cm fZCellMax = +15.0; // cm - fThetaMin = 60.0; // deg - fThetaMax = 170.0; // deg - fPhiMin = -50.0; // deg - fPhiMax = +50.0; // deg + fThetaMin = 60.0*TMath::DegToRad(); // rad + fThetaMax = 170.0*TMath::DegToRad(); // rad + fPhiMin = -50.0*TMath::DegToRad(); // rad + fPhiMax = +50.0*TMath::DegToRad(); // rad fchisq_cut[0] = 20.0; // default chi2 cut for tracks with 2 hodo hits fchisq_cut[1] = 15.0; // default chi2 cut for tracks with 1 hodo hit fchisq_cut[2] = 10.0; // default chi2 cut for tracks with no hodo hits From 73519076dc07dd9a1b370571e303b106651af8ab Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Fri, 27 Mar 2026 17:10:21 -0400 Subject: [PATCH 11/15] Also need to change the bounds checking to use radian --- src/THcLADKine.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index c810950..ad2be4e 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -217,10 +217,10 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { Double_t tmp = fPhiMin; fPhiMin = fPhiMax; fPhiMax = tmp; } - fThetaMin = TMath::Max(0.0, TMath::Min(180.0, fThetaMin)); - fThetaMax = TMath::Max(0.0, TMath::Min(180.0, fThetaMax)); - fPhiMin = TMath::Max(-180.0,TMath::Min(180.0, fPhiMin)); - fPhiMax = TMath::Max(-180.0,TMath::Min(180.0, fPhiMax)); + fThetaMin = TMath::Max(0.0, TMath::Min(TMath::Pi(), fThetaMin)); + fThetaMax = TMath::Max(0.0, TMath::Min(TMath::Pi(), fThetaMax)); + fPhiMin = TMath::Max(-TMath::Pi(),TMath::Min(TMath::Pi(), fPhiMin)); + fPhiMax = TMath::Max(-TMath::Pi(),TMath::Min(TMath::Pi(), fPhiMax)); return err; From ebf0657b2defcf3cdb7a104b3ecd808869b274ec Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Tue, 31 Mar 2026 13:25:16 -0400 Subject: [PATCH 12/15] chisq cut v2? --- src/THcLADKine.cxx | 237 +++++++++++++++++++++++++++++---------------- 1 file changed, 156 insertions(+), 81 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index ad2be4e..1ea9b7e 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -203,7 +203,7 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { { "ladkin.theta_max", &fThetaMax, kDouble, 0, 1 }, { "ladkin.phi_min", &fPhiMin, kDouble, 0, 1 }, { "ladkin.phi_max", &fPhiMax, kDouble, 0, 1 }, - { "ladkin.chisq_cut", fchisq_cut, kDouble, 3, 1 }, + { "ladkin.chisq_cut", fchisq_cut, kDouble, 2, 1 }, { 0 } }; gHcParms->LoadParmValues((DBRequest *)&track_constraints, prefix); @@ -273,113 +273,188 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { gemdir[1] = TMath::ATan2((v_hit2.Y() - v_hit1.Y()) , (v_hit2.X() - v_hit1.X())); gemdir[2] = vertex.Z(); //Hopefully the GEM track will be the same as elecctron vertex Z - Double_t bestchisq[2]={kBig,kBig};// 2 hodo hits, 1 hodo hit - Double_t best_gemdir[2][3]; - int bestHodoHitIndex[2]={-1,-1}; + Double_t bestchisq[3][3]={{kBig,kBig,kBig},{kBig,kBig,kBig},{kBig,kBig,kBig}};// no hodo hits, 1 hodo hit, 2 hodo hits+ (F-only, B-only, F&B) + Double_t best_gemdir[3][3][3]; + int bestHodoHitIndex[3][3]={{-1,-1,-1},{-1,-1,-1},{-1,-1,-1}}; + //fit GEM only first, if chisq is negative, skip the track and mark as bad + Double_t dir[3]; + dir[0]=gemdir[0]; + dir[1]=gemdir[1]; + dir[2]=gemdir[2];// if the track doesn't have a good hodoscope match, just fit to the GEM hits and vertex + bestchisq[0][0]= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); + if (bestchisq[0][0] < 0) { + track->SetIsGoodTrack(kFALSE); + track->SetChisq(bestchisq[0][0]); + track->SetAngles(-kBig, -kBig); + track->SetProjVertex(-kBig, -kBig, -kBig); + track->SetD0(-kBig); + track->SetHasHodoHit(kFALSE); + track->SetBestHodoHit(nullptr); + continue; + } + best_gemdir[0][0][0]=dir[0]; + best_gemdir[0][0][1]=dir[1]; + best_gemdir[0][0][2]=dir[2]; + int hit_index_for_best_chisq = -1; // (0-2)*3 for (no hodo hits, 1 hodo hit, 2 hodo hits) + (F-only, B-only)) //loop over hodoscope hits to find the best match to the track projection for (int iHodoHit=0; iHodoHit(LADHits_unfiltered->At(iHodoHit)); - if (hodo_hit == nullptr) + THcGoodLADHit *goodHit = static_cast(LADHits_unfiltered->At(iHodoHit)); + if (goodHit == nullptr) continue; - // There might be one or two hits with a signle goodhit, need to check both + // Check the number hits in goohit + int num_hodo_hits = 0; + if (goodHit->GetPlaneHit0() >= 0 && goodHit->GetPlaneHit0() < 999) { + num_hodo_hits++; + } + if (goodHit->GetPlaneHit1() >= 0 && goodHit->GetPlaneHit1() < 999) { + num_hodo_hits++; + } + std::vector v_hits; std::vector v_resolutions; v_hits.push_back(v_hit1); v_hits.push_back(v_hit2); v_resolutions.push_back(fSigma_GEM); v_resolutions.push_back(fSigma_GEM); - int num_hodo_hits = 0; - TVector3 v_hodo_hit; - if (hodo_hit->GetPlaneHit0() >= 0) { - v_hodo_hit = fHodoscope->GetHitPositionLab(hodo_hit->GetPlaneHit0(), hodo_hit->GetPaddleHit0(), hodo_hit->GetHitYPosHit0()); - v_hits.push_back(v_hodo_hit); + + if(num_hodo_hits==1){ + if (goodHit->GetPlaneHit0() >= 0 && goodHit->GetPlaneHit0() < 999) { + TVector3 hodo_hit_pos = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit0(), goodHit->GetPaddleHit0(), goodHit->GetHitYPosHit0()); + v_hits.push_back(hodo_hit_pos); + v_resolutions.push_back(fSigma_Hodo); + } else if (goodHit->GetPlaneHit1() >= 0 && goodHit->GetPlaneHit1() < 999) { + TVector3 hodo_hit_pos = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit1(), goodHit->GetPaddleHit1(), goodHit->GetHitYPosHit1()); + v_hits.push_back(hodo_hit_pos); + v_resolutions.push_back(fSigma_Hodo); + } + } else if(num_hodo_hits==2){ + TVector3 hodo_hit_pos1 = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit0(), goodHit->GetPaddleHit0(), goodHit->GetHitYPosHit0()); + TVector3 hodo_hit_pos2 = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit1(), goodHit->GetPaddleHit1(), goodHit->GetHitYPosHit1()); + v_hits.push_back(hodo_hit_pos1); + v_hits.push_back(hodo_hit_pos2); v_resolutions.push_back(fSigma_Hodo); - num_hodo_hits++; - } - if (hodo_hit->GetPlaneHit1() >= 0) { - v_hodo_hit = fHodoscope->GetHitPositionLab(hodo_hit->GetPlaneHit1(), hodo_hit->GetPaddleHit1(), hodo_hit->GetHitYPosHit1()); - v_hits.push_back(v_hodo_hit); v_resolutions.push_back(fSigma_Hodo); - num_hodo_hits++; } Double_t dir[3]; dir[0]=gemdir[0]; dir[1]=gemdir[1]; dir[2]=gemdir[2]; Double_t chisq = FitTrack(vertex, v_hits, v_resolutions, dir); - // Check if this is the best match so far, we only want to keep the best hodoscope hit for each track to avoid double counting hits - if (chisq>=0&&chisq=0 && chisq=0 && chisqGetHitPositionLab(goodHit->GetPlaneHit0(), goodHit->GetPaddleHit0(), goodHit->GetHitYPosHit0()); + TVector3 hodo_hit_pos2 = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit1(), goodHit->GetPaddleHit1(), goodHit->GetHitYPosHit1()); + v_hits.push_back(hodo_hit_pos1); + Double_t chisq1 = FitTrack(vertex, v_hits, v_resolutions, dir); + if (chisq1>=0 && chisq1=0 && chisq2=0&&bestchisq[0]SetAngles(best_gemdir[0][0], best_gemdir[0][1]); - track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[0][2]); - track->SetChisq(bestchisq[0]); - track->SetGoodD0(kTRUE); - track->SetD0(TMath::Abs(best_gemdir[0][2]-vertex.Z())*TMath::Sin(best_gemdir[0][0])); - THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[0])); - track->SetBestHodoHit(bestHodoHit); - track->SetHasHodoHit(2); - THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); - goodhit->CopyHit(0,0,bestHodoHit); - goodhit->CopyHit(1,1,bestHodoHit); - goodhit->SetTrackID(0, iTrack); - goodhit->SetTrackID(1, iTrack); - goodhit_n++; - }else if(bestchisq[1]>=0&&bestchisq[1]SetAngles(best_gemdir[1][0], best_gemdir[1][1]); - track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[1][2]); - track->SetChisq(bestchisq[1]); - track->SetGoodD0(kTRUE); - track->SetD0(TMath::Abs(best_gemdir[1][2]-vertex.Z())*TMath::Sin(best_gemdir[1][0])); - THcGoodLADHit *bestHodoHit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[1])); - track->SetBestHodoHit(bestHodoHit); - track->SetHasHodoHit(1); - THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); - goodhit->CopyHit(0,0,bestHodoHit); - goodhit->CopyHit(1,1,bestHodoHit); - goodhit->SetTrackID(0, iTrack); - goodhit->SetTrackID(1, iTrack); + //bestchisq should always be positive? + //check between f-only, b-only, and single, which has the best chisq + if (bestchisq[2][1] < bestchisq[2][1] && bestchisq[2][1] < bestchisq[1][0]){ + hit_index_for_best_chisq = 2*3+1; + } else if (bestchisq[2][2] < bestchisq[2][1] && bestchisq[2][2] < bestchisq[1][0]){ + hit_index_for_best_chisq = 2*3+2; + } else { + hit_index_for_best_chisq = 1*3+0; + } + cout<SetChisq(bestchisq[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3]); + track->SetAngles(best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][0], best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][1]); + track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][2]); + Double_t d0 = fabs(vertex.Z()-best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][2]); + track->SetD0(d0); + + track->SetHasHodoHit(hit_index_for_best_chisq/3); + if (hit_index_for_best_chisq >= 1) { + THcGoodLADHit *besthit = static_cast(LADHits_unfiltered->At(bestHodoHitIndex[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3])); + THcGoodLADHit *newhit = static_cast(fGoodLADHits->ConstructedAt(goodhit_n)); goodhit_n++; - }else { - Double_t dir[3]; - dir[0]=gemdir[0]; - dir[1]=gemdir[1]; - dir[2]=gemdir[2];// if the track doesn't have a good hodoscope match, just fit to the GEM hits and vertex - Double_t chisq= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); - track->SetChisq(chisq);// set chisq even if the fit fails so we can see why - if (chisq>=0&&chisqSetGoodD0(kTRUE); - isGoodTrack[iTrack] = true; - track->SetAngles(dir[0], dir[1]); - track->SetProjVertex(vertex.X(), vertex.Y(), dir[2]); - track->SetD0(TMath::Abs(dir[2]-vertex.Z())*TMath::Sin(dir[0])); + if (besthit == nullptr || newhit == nullptr) { track->SetBestHodoHit(nullptr); - track->SetHasHodoHit(0); - }else { - //Still bad track, mark as bad and set parameters to default values - isGoodTrack[iTrack] = false; - //cout<<"THcLADKine: Track fit failed for track "<SetAngles(-kBig, -kBig); - track->SetProjVertex(-kBig, -kBig, -kBig); - track->SetD0(-kBig); - track->SetGoodD0(kFALSE); - track->SetHasHodoHit(0); - track->SetBestHodoHit(nullptr); - track->SetIsGoodTrack(false); continue; } + if (goodhit_n <= MAXGOODHITS) { + if (hit_index_for_best_chisq/3==2){ + int ntmp=0; + if(hit_index_for_best_chisq%3==1 || hit_index_for_best_chisq%3==0){ + newhit->CopyHit(0,0,besthit); + newhit->SetTrackID(0,track->GetTrackID()); + ntmp++; + } else if(hit_index_for_best_chisq%3==2 || hit_index_for_best_chisq%3==0){ + newhit->CopyHit(1,1,besthit); + newhit->SetTrackID(1,track->GetTrackID()); + ntmp++; + } + track->SetHasHodoHit(ntmp); + }else if (hit_index_for_best_chisq/3==1){ + newhit->CopyHit(0,0,besthit); + newhit->SetTrackID(0,track->GetTrackID()); + track->SetHasHodoHit(1); + } + if (hit_index_for_best_chisq/3!=0){ + track->SetBestHodoHit(newhit); + isGoodTrack[iTrack] = true; + } + } + } else { + track->SetBestHodoHit(nullptr); + track->SetHasHodoHit(0); + isGoodTrack[iTrack] = true; } + + if (track->GetGoodD0()) { if (track->GetD0()>0 && track->GetD0() < fD0Cut_wVertex) { isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; From 4086f30c0c6f57c13c9e0b12adb1934603c46ddd Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Tue, 31 Mar 2026 13:34:02 -0400 Subject: [PATCH 13/15] remove debug print out --- src/THcLADKine.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 1ea9b7e..cb88c81 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -399,7 +399,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } else { hit_index_for_best_chisq = 1*3+0; } - cout< Date: Thu, 2 Apr 2026 16:00:02 -0400 Subject: [PATCH 14/15] fix bug in copying goodLADhit --- src/THcLADKine.cxx | 43 +++++++++++++++++++++++-------------------- src/THcLADKine.h | 2 +- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index cb88c81..623b0c9 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -40,7 +40,6 @@ ClassImp(THcLADKine) fPhiMax = +50.0*TMath::DegToRad(); // rad fchisq_cut[0] = 20.0; // default chi2 cut for tracks with 2 hodo hits fchisq_cut[1] = 15.0; // default chi2 cut for tracks with 1 hodo hit - fchisq_cut[2] = 10.0; // default chi2 cut for tracks with no hodo hits fSigma_GEM = 0.1; // default GEM resolution in cm fSigma_Hodo = 10; // default Hodoscope resolution in cm @@ -222,7 +221,6 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { fPhiMin = TMath::Max(-TMath::Pi(),TMath::Min(TMath::Pi(), fPhiMin)); fPhiMax = TMath::Max(-TMath::Pi(),TMath::Min(TMath::Pi(), fPhiMax)); - return err; } //_____________________________________________________________________________ @@ -274,6 +272,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { gemdir[2] = vertex.Z(); //Hopefully the GEM track will be the same as elecctron vertex Z Double_t bestchisq[3][3]={{kBig,kBig,kBig},{kBig,kBig,kBig},{kBig,kBig,kBig}};// no hodo hits, 1 hodo hit, 2 hodo hits+ (F-only, B-only, F&B) + bool usedHodoHit[3][3]={{kFALSE,kFALSE,kFALSE},{kFALSE,kFALSE,kFALSE},{kFALSE,kFALSE,kFALSE}}; Double_t best_gemdir[3][3][3]; int bestHodoHitIndex[3][3]={{-1,-1,-1},{-1,-1,-1},{-1,-1,-1}}; //fit GEM only first, if chisq is negative, skip the track and mark as bad @@ -282,6 +281,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { dir[1]=gemdir[1]; dir[2]=gemdir[2];// if the track doesn't have a good hodoscope match, just fit to the GEM hits and vertex bestchisq[0][0]= FitTrack(vertex, {v_hit1, v_hit2}, {fSigma_GEM, fSigma_GEM}, dir); + usedHodoHit[0][0]=kTRUE; if (bestchisq[0][0] < 0) { track->SetIsGoodTrack(kFALSE); track->SetChisq(bestchisq[0][0]); @@ -347,6 +347,8 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { best_gemdir[1][0][1]=dir[1]; best_gemdir[1][0][2]=dir[2]; bestHodoHitIndex[1][0]=iHodoHit; + usedHodoHit[1][0]=kTRUE; + } } else if(num_hodo_hits==2){ if (chisq>=0 && chisq fchisq_cut[1]) { //check if the GEM only fit is better than the best hodo hit fit by more than the chisq cut, if so, use the GEM only fit hit_index_for_best_chisq = 0*3+0; } - cout<<"using "<SetChisq(bestchisq[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3]); track->SetAngles(best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][0], best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][1]); track->SetProjVertex(vertex.X(), vertex.Y(), best_gemdir[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3][2]); @@ -431,7 +434,8 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { newhit->CopyHit(0,0,besthit); newhit->SetTrackID(0,track->GetTrackID()); ntmp++; - } else if(hit_index_for_best_chisq%3==2 || hit_index_for_best_chisq%3==0){ + } + if(hit_index_for_best_chisq%3==2 || hit_index_for_best_chisq%3==0){ newhit->CopyHit(1,1,besthit); newhit->SetTrackID(1,track->GetTrackID()); ntmp++; @@ -672,17 +676,16 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position } //check if the gem track is pointing to the hodoscope hits if there are any, return negative chisq if not if (nPoints >2){ - double sx = TMath::Sin(dir[0]) * TMath::Cos(dir[1]); - double sy = TMath::Sin(dir[0]) * TMath::Sin(dir[1]); - double sz = TMath::Cos(dir[0]); + TVector3 dir_vec = sp_positions[1] - sp_positions[0]; + dir_vec = dir_vec.Unit(); for (int i = 2; i < nPoints; i++) { //closest approach of the track to the point - double t = ((sp_positions[i].X() - sp_positions[0].X()) * sx + - (sp_positions[i].Y() - sp_positions[0].Y()) * sy + - (sp_positions[i].Z() - sp_positions[0].Z()) * sz); - double x_closest = sp_positions[0].X() + t * sx; - double y_closest = sp_positions[0].Y() + t * sy; - double z_closest = sp_positions[0].Z() + t * sz; + double t = ((sp_positions[i].X() - sp_positions[0].X()) * dir_vec.X() + + (sp_positions[i].Y() - sp_positions[0].Y()) * dir_vec.Y() + + (sp_positions[i].Z() - sp_positions[0].Z()) * dir_vec.Z()); + double x_closest = sp_positions[0].X() + t * dir_vec.X(); + double y_closest = sp_positions[0].Y() + t * dir_vec.Y(); + double z_closest = sp_positions[0].Z() + t * dir_vec.Z(); double dx = sp_positions[i].X() - x_closest; double dy = sp_positions[i].Y() - y_closest; double dz = sp_positions[i].Z() - z_closest; diff --git a/src/THcLADKine.h b/src/THcLADKine.h index b2dbe49..83aac0d 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -64,7 +64,7 @@ class THcLADKine : public THcPrimaryKine { Double_t fPhiMin; // default -50.0 deg Double_t fPhiMax; // default +50.0 deg - Double_t fchisq_cut[3];//chisq cuts for track, 2 hodo hits, 1 hodo hit, and no hodo hits. + Double_t fchisq_cut[2];//chisq difference between 1 and 2 hodo hit track fits, used to determine if we accept tracks with only 1 hodo hit (if chisq_2hit - chisq_1hit > fchisq_cut[0]), or if we have no hodo hits (if chisq_1hit - chisq_0hit < fchisq_cut[1]) Double_t fSigma_GEM; // GEM resolution in cm, used for track fitting, should be set based on detector performance Double_t fSigma_Hodo; // Hodoscope resolution in cm, used for track fitting, should be set based on detector performance From 60ad8c26a9d86ca35f06ee1ac2e667f50723e18a Mon Sep 17 00:00:00 2001 From: Ching Him Leung Date: Thu, 14 May 2026 12:02:49 -0400 Subject: [PATCH 15/15] delete minimizer even if fit fail --- src/THcLADKine.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 623b0c9..2ed247c 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -741,6 +741,7 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position minimizer->SetLimitedVariable(2, "z_vertex", initial_params[2], 0.1, fZCellMin, fZCellMax); //the target length is 20cm minimizer->Minimize(); if (minimizer->Status() != 0) { + delete minimizer; return -6; // Fit did not converge } const double *best_params = minimizer->X();