From 41df05898f81bba93de81712ce5991b5249eb2d3 Mon Sep 17 00:00:00 2001 From: LucasEhinger Date: Thu, 21 May 2026 16:05:09 -0400 Subject: [PATCH] Fix merge conflicts. Restructure Ladhodo and lad kin --Rather than copying the goodladhits from lad hodo and having almost an exact duplicate for ladkin, the ladkin variables are removed, and ladkin only modifies the ladhod variables. Added a chi-quare and is_proton variable to the goodladhit. Accidentally auto-formatted LADKine in this commit, hence the large diff. Will attempt to break into two commits later. TODO: --Add a proton cut to ladkin. --Clean up once functionality is verified (beyond just compiling) --Other todo's in Ching Him's tracking PR --- src/THcGoodLADHit.h | 51 ++-- src/THcLADHodoscope.cxx | 31 ++- src/THcLADKine.cxx | 515 +++++++++++++++++++++------------------- 3 files changed, 309 insertions(+), 288 deletions(-) diff --git a/src/THcGoodLADHit.h b/src/THcGoodLADHit.h index c5d3394..68865f6 100644 --- a/src/THcGoodLADHit.h +++ b/src/THcGoodLADHit.h @@ -7,10 +7,13 @@ class THcGoodLADHit : public TObject { public: THcGoodLADHit() { for (int i = 0; i < 2; ++i) { - plane[i] = paddle[i] = track_id[i] = -1; - hit_tof[i] = hit_tof_rfcorr[i] = dTrk_horiz[i] = dTrk_vert[i] = hit_time[i] = hit_theta[i] = hit_phi[i] = hit_edep[i] = - hit_edep_amp[i] = hit_yPos[i] = hit_alpha[i] = hit_beta[i] = 1e30; + plane[i] = paddle[i] = -1; + hit_tof[i] = hit_tof_rfcorr[i] = hit_time[i] = hit_theta[i] = hit_phi[i] = hit_edep[i] = hit_edep_amp[i] = + hit_yPos[i] = hit_alpha[i] = hit_beta[i] = 1e30; + is_proton[i] = false; } + trk_chiSqr = 1e30; + track_id = -1; }; virtual ~THcGoodLADHit() = default; @@ -28,17 +31,11 @@ class THcGoodLADHit : public TObject { CheckHitIndex(hit); paddle[hit] = value; } - void SetTrackID(Int_t hit, Int_t value) { + void SetTrackID(Int_t value) { track_id = value; } + void SetTrkChiSqr(Double_t value) { trk_chiSqr = value; } + void SetIsProton(Int_t hit, Bool_t value) { CheckHitIndex(hit); - track_id[hit] = value; - } - void SetdTrkHoriz(Int_t hit, Double_t value) { - CheckHitIndex(hit); - dTrk_horiz[hit] = value; - } - void SetdTrkVert(Int_t hit, Double_t value) { - CheckHitIndex(hit); - dTrk_vert[hit] = value; + is_proton[hit] = value; } void SetHitTime(Int_t hit, Double_t value) { CheckHitIndex(hit); @@ -86,9 +83,9 @@ class THcGoodLADHit : public TObject { if (copy_plane == 0) { SetPlane(this_plane, copyhit->GetPlaneHit0()); SetPaddle(this_plane, copyhit->GetPaddleHit0()); - SetTrackID(this_plane, copyhit->GetTrackIDHit0()); - SetdTrkHoriz(this_plane, copyhit->GetdTrkHorizHit0()); - SetdTrkVert(this_plane, copyhit->GetdTrkVertHit0()); + SetTrackID(copyhit->GetTrackID()); + SetTrkChiSqr(copyhit->GetTrkChiSqr()); + SetIsProton(this_plane, copyhit->GetIsProtonHit0()); SetHitTime(this_plane, copyhit->GetHitTimeHit0()); SetHitTheta(this_plane, copyhit->GetHitThetaHit0()); SetHitPhi(this_plane, copyhit->GetHitPhiHit0()); @@ -101,9 +98,9 @@ class THcGoodLADHit : public TObject { } else if (copy_plane == 1) { SetPlane(this_plane, copyhit->GetPlaneHit1()); SetPaddle(this_plane, copyhit->GetPaddleHit1()); - SetTrackID(this_plane, copyhit->GetTrackIDHit1()); - SetdTrkHoriz(this_plane, copyhit->GetdTrkHorizHit1()); - SetdTrkVert(this_plane, copyhit->GetdTrkVertHit1()); + SetTrackID(copyhit->GetTrackID()); + SetTrkChiSqr(copyhit->GetTrkChiSqr()); + SetIsProton(this_plane, copyhit->GetIsProtonHit1()); SetHitTime(this_plane, copyhit->GetHitTimeHit1()); SetHitTheta(this_plane, copyhit->GetHitThetaHit1()); SetHitPhi(this_plane, copyhit->GetHitPhiHit1()); @@ -127,17 +124,15 @@ class THcGoodLADHit : public TObject { Int_t GetPaddleHit0() const { return paddle[0]; } Int_t GetPaddleHit1() const { return paddle[1]; } - Int_t GetTrackIDHit0() const { return track_id[0]; } - Int_t GetTrackIDHit1() const { return track_id[1]; } + Int_t GetTrackID() const { return track_id; } Double_t GetBetaHit0() const { return hit_beta[0]; } Double_t GetBetaHit1() const { return hit_beta[1]; } - Double_t GetdTrkHorizHit0() const { return dTrk_horiz[0]; } - Double_t GetdTrkHorizHit1() const { return dTrk_horiz[1]; } + Double_t GetTrkChiSqr() const { return trk_chiSqr; } - Double_t GetdTrkVertHit0() const { return dTrk_vert[0]; } - Double_t GetdTrkVertHit1() const { return dTrk_vert[1]; } + Double_t GetIsProtonHit0() const { return is_proton[0]; } + Double_t GetIsProtonHit1() const { return is_proton[1]; } Double_t GetHitTimeHit0() const { return hit_time[0]; } Double_t GetHitTimeHit1() const { return hit_time[1]; } @@ -169,9 +164,9 @@ class THcGoodLADHit : public TObject { protected: Int_t plane[2]; Int_t paddle[2]; - Int_t track_id[2]; - Double_t dTrk_horiz[2]; - Double_t dTrk_vert[2]; + Int_t track_id; + Double_t trk_chiSqr; + Double_t is_proton[2]; Double_t hit_time[2]; Double_t hit_beta[2]; Double_t hit_theta[2]; diff --git a/src/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index a050e1c..b2dfa6c 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -217,30 +217,36 @@ Int_t THcLADHodoscope::DefineVariables(EMode mode) { {"goodhit_n", "Number of good hits", "goodhit_n"}, {"goodhit_plane_0", "Good hit plane", "fGoodLADHits.THcGoodLADHit.GetPlaneHit0()"}, {"goodhit_paddle_0", "Good hit paddle", "fGoodLADHits.THcGoodLADHit.GetPaddleHit0()"}, - {"goodhit_trackid_0", "Good hit track ID", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit0()"}, {"goodhit_beta_0", "Good hit beta", "fGoodLADHits.THcGoodLADHit.GetBetaHit0()"}, - {"goodhit_dTrkHoriz_0", "Good hit horizontal trk proj - hit position", - "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit0()"}, - {"goodhit_dTrkVert_0", "Good hit vertical trk proj - hit position", - "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit0()"}, + {"goodhit_isProton_0", "Good hit is proton", "fGoodLADHits.THcGoodLADHit.GetIsProtonHit0()"}, {"goodhit_hittime_0", "Good hit time", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit0()"}, {"goodhit_hittheta_0", "Good hit theta", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit0()"}, {"goodhit_hitphi_0", "Good hit phi", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit0()"}, {"goodhit_hitedep_0", "Good hit energy deposition", "fGoodLADHits.THcGoodLADHit.GetHitEdepHit0()"}, {"goodhit_hitedep_amp_0", "Good hit energy deposition (amplitude)", "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit0()"}, + {"goodhit_hit_alpha_0", "Good hit alpha", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit0()"}, + {"goodhit_hit_ypos_0", "Good hit y position", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit0()"}, + {"goodhit_hit_tof_0", "Good hit time-of-flight", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit0()"}, + {"goodhit_hit_tof_rfcorr_0", "Good hit time-of-flight (RF corrected)", + "fGoodLADHits.THcGoodLADHit.GetHitTOFRFcorrHit0()"}, {"goodhit_plane_1", "Good hit plane (second plane)", "fGoodLADHits.THcGoodLADHit.GetPlaneHit1()"}, {"goodhit_paddle_1", "Good hit paddle (second plane)", "fGoodLADHits.THcGoodLADHit.GetPaddleHit1()"}, - {"goodhit_trackid_1", "Good hit track ID (second plane)", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit1()"}, {"goodhit_beta_1", "Good hit beta (second plane)", "fGoodLADHits.THcGoodLADHit.GetBetaHit1()"}, - {"goodhit_dTrkHoriz_1", "Good hit horizontal trk proj - hit position (second plane)", - "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit1()"}, - {"goodhit_dTrkVert_1", "Good hit vertical trk proj - hit position (second plane)", - "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit1()"}, + {"goodhit_isProton_1", "Good hit is proton (second plane)", "fGoodLADHits.THcGoodLADHit.GetIsProtonHit1()"}, {"goodhit_hittime_1", "Good hit time (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit1()"}, {"goodhit_hittheta_1", "Good hit theta (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit1()"}, {"goodhit_hitphi_1", "Good hit phi (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit1()"}, {"goodhit_hitedep_1", "Good hit energy deposition (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitEdepHit1()"}, + {"goodhit_hitedep_amp_1", "Good hit energy deposition (amplitude, second plane)", + "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit1()"}, + {"goodhit_hit_alpha_1", "Good hit alpha (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit1()"}, + {"goodhit_hit_ypos_1", "Good hit y position (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit1()"}, + {"goodhit_hit_tof_1", "Good hit time-of-flight (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit1()"}, + {"goodhit_hit_tof_rfcorr_1", "Good hit time-of-flight (RF corrected, second plane)", + "fGoodLADHits.THcGoodLADHit.GetHitTOFRFcorrHit1()"}, + {"goodhit_trackid", "Good hit track ID", "fGoodLADHits.THcGoodLADHit.GetTrackID()"}, + {"goodhit_chiSquare", "Good hit chi-square", "fGoodLADHits.THcGoodLADHit.GetTrkChiSqr()"}, {"good_hit_n_unique", "Number of unique good hits", "num_unique_good_hits"}, {"all_hits_n_unique", "Number of all hits, not just with tracks", "num_unique_hits"}, {0}}; @@ -579,11 +585,10 @@ TVector3 THcLADHodoscope::GetHitPositionLab(int plane, int paddle, double ypos) 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]->GetZpos()); + Double_t offset = fPlanes[plane]->GetPosCenter(paddle); + TVector3 hit = TVector3(offset, ypos, fPlanes[plane]->GetZpos()); hit.RotateY(fPlanes[plane]->GetTheta()); return hit; } - ClassImp(THcLADHodoscope) diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 7a04ef8..8b9a53e 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -1,5 +1,8 @@ #include "THcLADKine.h" +#include "Math/Factory.h" +#include "Math/Functor.h" +#include "Math/Minimizer.h" #include "THcGlobals.h" #include "THcHallCSpectrometer.h" #include "THcLADGEM.h" @@ -10,10 +13,6 @@ #include "TVector3.h" #include "VarDef.h" #include "VarType.h" -#include "Math/Minimizer.h" -#include "Math/Factory.h" -#include "Math/Functor.h" - ClassImp(THcLADKine) //_____________________________________________________________________________ @@ -22,7 +21,7 @@ ClassImp(THcLADKine) : THcPrimaryKine(name, description, spectro, primary_kine, 0.938), fSpecName(spectro), fGEM(nullptr), fHodoscope(nullptr), fVertexModuleName(vertex_module), fVertexModule(nullptr), fTrack(nullptr) { // Constructor - fGoodLADHits = new TClonesArray("THcGoodLADHit", MAXGOODHITS); + // fGoodLADHits = nullptr; goodhit_n = 0; fTVertex = kBig; fRFTime = kBig; @@ -32,25 +31,23 @@ ClassImp(THcLADKine) n_rf_offsets = 0; rf_period = 4.00801; // Default RF period in ns - fZCellMin = -15.0; // cm - fZCellMax = +15.0; // cm - 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 - fSigma_GEM = 0.1; // default GEM resolution in cm - fSigma_Hodo = 10; // default Hodoscope resolution in cm - - + fZCellMin = -15.0; // cm + fZCellMax = +15.0; // cm + 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 + fSigma_GEM = 0.1; // default GEM resolution in cm + fSigma_Hodo = 10; // default Hodoscope resolution in cm } //_____________________________________________________________________________ THcLADKine::~THcLADKine() { // Destructor DefineVariables(kDelete); - fGoodLADHits->Delete(); - delete fGoodLADHits; + // fGoodLADHits->Delete(); + // delete fGoodLADHits; delete[] fFixed_z; delete[] rf_offset; } @@ -58,7 +55,7 @@ THcLADKine::~THcLADKine() { void THcLADKine::Clear(Option_t *opt) { // Clear the object THcPrimaryKine::Clear(opt); - fGoodLADHits->Delete(); + // fGoodLADHits->Delete(); goodhit_n = 0; for (Int_t i = 0; i < n_rf_offsets; i++) { rf_offset[i] = 0.0; @@ -127,8 +124,6 @@ THaAnalysisObject::EStatus THcLADKine::Init(const TDatime &run_time) { return fStatus; } - - return THaPhysicsModule::Init(run_time); } //_____________________________________________________________________________ @@ -181,30 +176,33 @@ 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 }, - { "ladkin.chisq_cut", fchisq_cut, kDouble, 2, 1 }, - { 0 } }; + 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}, + {"ladkin.chisq_cut", fchisq_cut, kDouble, 2, 1}, {0}}; gHcParms->LoadParmValues((DBRequest *)&track_constraints, prefix); if (fZCellMin > fZCellMax) { - Double_t tmp = fZCellMin; fZCellMin = fZCellMax; fZCellMax = tmp; + Double_t tmp = fZCellMin; + fZCellMin = fZCellMax; + fZCellMax = tmp; } if (fThetaMin > fThetaMax) { - Double_t tmp = fThetaMin; fThetaMin = fThetaMax; fThetaMax = tmp; + Double_t tmp = fThetaMin; + fThetaMin = fThetaMax; + fThetaMax = tmp; } if (fPhiMin > fPhiMax) { - Double_t tmp = fPhiMin; fPhiMin = fPhiMax; fPhiMax = tmp; + Double_t tmp = fPhiMin; + fPhiMin = fPhiMax; + fPhiMax = tmp; } - 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)); + 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; } @@ -214,10 +212,10 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { fTrack = fSpectro->GetGoldenTrack(); CalculateTVertex(); - //Get LADGoodHits for track matching + // Get LADGoodHits for track matching TClonesArray *LADHits_unfiltered = fHodoscope->GetLADGoodHits(); Int_t nHits = LADHits_unfiltered->GetLast() + 1; - goodhit_n = 0; + goodhit_n = 0; ////////////////////////////////////////////////////////////////////////////// // Check track projection to vertex fGEMTracks = fGEM->GetTracks(); @@ -238,7 +236,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; @@ -251,21 +249,23 @@ 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()) ; - 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[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 + 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[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 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); - usedHodoHit[0][0]=kTRUE; + 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); + usedHodoHit[0][0] = kTRUE; if (bestchisq[0][0] < 0) { track->SetIsGoodTrack(kFALSE); track->SetChisq(bestchisq[0][0]); @@ -276,12 +276,12 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { 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]; + 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 (goodHit == nullptr) continue; @@ -300,138 +300,150 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { v_hits.push_back(v_hit2); v_resolutions.push_back(fSigma_GEM); v_resolutions.push_back(fSigma_GEM); - - if(num_hodo_hits==1){ + + if (num_hodo_hits == 1) { if (goodHit->GetPlaneHit0() >= 0 && goodHit->GetPlaneHit0() < 999) { - TVector3 hodo_hit_pos = fHodoscope->GetHitPositionLab(goodHit->GetPlaneHit0(), goodHit->GetPaddleHit0(), goodHit->GetHitYPosHit0()); + 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()); + 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()); + } 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); v_resolutions.push_back(fSigma_Hodo); } Double_t dir[3]; - dir[0]=gemdir[0]; - dir[1]=gemdir[1]; - dir[2]=gemdir[2]; + dir[0] = gemdir[0]; + dir[1] = gemdir[1]; + dir[2] = gemdir[2]; Double_t chisq = FitTrack(vertex, v_hits, v_resolutions, dir); - if (num_hodo_hits==1){ - if (chisq>=0 && chisq= 0 && chisq < bestchisq[1][0]) { + bestchisq[1][0] = chisq; + best_gemdir[1][0][0] = dir[0]; + 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= 0 && chisq < bestchisq[2][0]) { + bestchisq[2][0] = chisq; + best_gemdir[2][0][0] = dir[0]; + best_gemdir[2][0][1] = dir[1]; + best_gemdir[2][0][2] = dir[2]; + bestHodoHitIndex[2][0] = iHodoHit; + usedHodoHit[2][0] = kTRUE; } v_hits.pop_back(); v_hits.pop_back(); v_resolutions.pop_back(); - dir[0]=gemdir[0]; - dir[1]=gemdir[1]; - dir[2]=gemdir[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()); + dir[0] = gemdir[0]; + dir[1] = gemdir[1]; + dir[2] = gemdir[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); Double_t chisq1 = FitTrack(vertex, v_hits, v_resolutions, dir); - if (chisq1>=0 && chisq1= 0 && chisq1 < bestchisq[2][1]) { + bestchisq[2][1] = chisq1; + best_gemdir[2][1][0] = dir[0]; + best_gemdir[2][1][1] = dir[1]; + best_gemdir[2][1][2] = dir[2]; + bestHodoHitIndex[2][1] = iHodoHit; + usedHodoHit[2][1] = kTRUE; } v_hits.pop_back(); v_hits.push_back(hodo_hit_pos2); - dir[0]=gemdir[0]; - dir[1]=gemdir[1]; - dir[2]=gemdir[2]; + dir[0] = gemdir[0]; + dir[1] = gemdir[1]; + dir[2] = gemdir[2]; Double_t chisq2 = FitTrack(vertex, v_hits, v_resolutions, dir); - if (chisq2>=0 && chisq2= 0 && chisq2 < bestchisq[2][2]) { + bestchisq[2][2] = chisq2; + best_gemdir[2][2][0] = dir[0]; + best_gemdir[2][2][1] = dir[1]; + best_gemdir[2][2][2] = dir[2]; + bestHodoHitIndex[2][2] = iHodoHit; + usedHodoHit[2][2] = kTRUE; } } v_hits.clear(); v_resolutions.clear(); } - //bestchisq should always be positive? - //check between f-only, b-only, and single, which has the best chisq - hit_index_for_best_chisq = 2*3+0;// default to 2 hodo hit - if (usedHodoHit[2][1] && bestchisq[2][1] < bestchisq[2][2] && bestchisq[2][1] < bestchisq[1][0]){ - hit_index_for_best_chisq = 2*3+1; - } else if (usedHodoHit[2][2] && bestchisq[2][2] < bestchisq[2][1] && bestchisq[2][2] < bestchisq[1][0]){ - hit_index_for_best_chisq = 2*3+2; + // bestchisq should always be positive? + // check between f-only, b-only, and single, which has the best chisq + hit_index_for_best_chisq = 2 * 3 + 0; // default to 2 hodo hit + if (usedHodoHit[2][1] && bestchisq[2][1] < bestchisq[2][2] && bestchisq[2][1] < bestchisq[1][0]) { + hit_index_for_best_chisq = 2 * 3 + 1; + } else if (usedHodoHit[2][2] && bestchisq[2][2] < bestchisq[2][1] && bestchisq[2][2] < bestchisq[1][0]) { + hit_index_for_best_chisq = 2 * 3 + 2; } else if (usedHodoHit[1][0]) { - hit_index_for_best_chisq = 1*3+0; + hit_index_for_best_chisq = 1 * 3 + 0; } - if (usedHodoHit[2][0] && ( bestchisq[2][0] - bestchisq[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3]<=fchisq_cut[0])) { - //check if the 2 hodo hit fit is better than the best single hodo hit fit by more than the chisq cut, if so, use the 2 hodo hit fit - hit_index_for_best_chisq = 2*3+0; - }else if (usedHodoHit[0][0] && (bestchisq[hit_index_for_best_chisq/3][hit_index_for_best_chisq%3] - bestchisq[0][0]) > 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; + if (usedHodoHit[2][0] && + (bestchisq[2][0] - bestchisq[hit_index_for_best_chisq / 3][hit_index_for_best_chisq % 3] <= fchisq_cut[0])) { + // check if the 2 hodo hit fit is better than the best single hodo hit fit by more than the chisq cut, if so, use + // the 2 hodo hit fit + hit_index_for_best_chisq = 2 * 3 + 0; + } else if (usedHodoHit[0][0] && (bestchisq[hit_index_for_best_chisq / 3][hit_index_for_best_chisq % 3] - + bestchisq[0][0]) > 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; } - track->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->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); + 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)); + 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++; - if (besthit == nullptr || newhit == nullptr) { + if (besthit == nullptr) { track->SetBestHodoHit(nullptr); 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++; + 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) { + besthit->SetTrackID(track->GetTrackID()); + besthit->SetTrkChiSqr(track->GetChisq()); + ntmp++; } - 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()); + if (hit_index_for_best_chisq % 3 == 2 || hit_index_for_best_chisq % 3 == 0) { + besthit->SetTrackID(track->GetTrackID()); + besthit->SetTrkChiSqr(track->GetChisq()); ntmp++; } track->SetHasHodoHit(ntmp); - }else if (hit_index_for_best_chisq/3==1){ - newhit->CopyHit(0,0,besthit); - newhit->SetTrackID(0,track->GetTrackID()); + } else if (hit_index_for_best_chisq / 3 == 1) { + besthit->SetTrackID(track->GetTrackID()); + besthit->SetTrkChiSqr(track->GetChisq()); track->SetHasHodoHit(1); } - if (hit_index_for_best_chisq/3!=0){ - track->SetBestHodoHit(newhit); + if (hit_index_for_best_chisq / 3 != 0) { + track->SetBestHodoHit(besthit); isGoodTrack[iTrack] = true; } } @@ -441,10 +453,8 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { isGoodTrack[iTrack] = true; } - - if (track->GetGoodD0()) { - if (track->GetD0()>0 && track->GetD0() < fD0Cut_wVertex) { + if (track->GetD0() > 0 && track->GetD0() < fD0Cut_wVertex) { isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; } } else { @@ -452,7 +462,7 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; } } - if(track->GetChisq()<0) { + if (track->GetChisq() < 0) { isGoodTrack[iTrack] = false; } if (track->GetdT() > fTrk_dtCut) { @@ -460,9 +470,12 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } track->SetIsGoodTrack(isGoodTrack[iTrack]); } - // Calculate beta, alpha, tof, etc. for the hits - for (Int_t i = 0; i < goodhit_n; i++) { - THcGoodLADHit *goodhit = static_cast(fGoodLADHits->At(i)); + + ////////////////////////////////////////////////////////////////////////////// + // Calculate ToF and RF-corrected ToF for good LAD hits + + for (Int_t i = 0; i < nHits; i++) { + THcGoodLADHit *goodhit = static_cast(LADHits_unfiltered->At(i)); if (goodhit == nullptr) continue; @@ -484,8 +497,9 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } ////////////////////////////////////////////////////////////////////////////// - // Calculate tof for raw planes - // Get the number of planes + // Calculate tof for raw planes (not good lad hits) + // Not sure this is necessary. TODO: revisit and possibly remove + Int_t nPlanes = fHodoscope->GetNPlanes(); // Loop through the planes and get LAD hits @@ -511,8 +525,6 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { } } - // LADHits_unfiltered->Clear(); // Clear the unfiltered hits - // fGEMTracks->Clear(); // Clear the tracks return kOK; } //_____________________________________________________________________________ @@ -618,42 +630,44 @@ Int_t THcLADKine::DefineVariables(EMode mode) { RVarDef vars[] = { {"t_vertex", "Calculated vertex time", "fTVertex"}, {"t_vertex_RFcorr", "Calculated vertex time (RF corrected)", "fTVertex_RFcorr"}, - {"n_goodhits", "Number of good hits", "goodhit_n"}, - {"goodhit_plane_0", "Good hit plane", "fGoodLADHits.THcGoodLADHit.GetPlaneHit0()"}, - {"goodhit_paddle_0", "Good hit paddle", "fGoodLADHits.THcGoodLADHit.GetPaddleHit0()"}, - {"goodhit_trackid_0", "Good hit track ID", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit0()"}, - {"goodhit_beta_0", "Good hit beta", "fGoodLADHits.THcGoodLADHit.GetBetaHit0()"}, - {"goodhit_dTrkHoriz_0", "Good hit horizontal trk proj - hit position", - "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit0()"}, - {"goodhit_dTrkVert_0", "Good hit vertical trk proj - hit position", - "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit0()"}, - {"goodhit_hittime_0", "Good hit time", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit0()"}, - {"goodhit_hittheta_0", "Good hit theta", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit0()"}, - {"goodhit_hitphi_0", "Good hit phi", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit0()"}, - {"goodhit_hitedep_0", "Good hit energy deposition", "fGoodLADHits.THcGoodLADHit.GetHitEdepHit0()"}, - {"goodhit_hitedep_amp_0", "Good hit energy deposition (amplitude)", - "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit0()"}, - {"goodhit_plane_1", "Good hit plane (second plane)", "fGoodLADHits.THcGoodLADHit.GetPlaneHit1()"}, - {"goodhit_paddle_1", "Good hit paddle (second plane)", "fGoodLADHits.THcGoodLADHit.GetPaddleHit1()"}, - {"goodhit_trackid_1", "Good hit track ID (second plane)", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit1()"}, - {"goodhit_beta_1", "Good hit beta (second plane)", "fGoodLADHits.THcGoodLADHit.GetBetaHit1()"}, - {"goodhit_dTrkHoriz_1", "Good hit horizontal trk proj - hit position (second plane)", - "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit1()"}, - {"goodhit_dTrkVert_1", "Good hit vertical trk proj - hit position (second plane)", - "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit1()"}, - {"goodhit_hittime_1", "Good hit time (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit1()"}, - {"goodhit_hittheta_1", "Good hit theta (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit1()"}, - {"goodhit_hitphi_1", "Good hit phi (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit1()"}, - {"goodhit_hitedep_1", "Good hit energy deposition (second plane)", - "fGoodLADHits.THcGoodLADHit.GetHitEdepHit1()"}, - {"goodhit_hitedep_amp_1", "Good hit energy deposition (amplitude) (second plane)", - "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit1()"}, - {"goodhit_ypos_0", "Good hit y position", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit0()"}, - {"goodhit_ypos_1", "Good hit y position (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit1()"}, - {"goodhit_tof_0", "Good hit time of flight", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit0()"}, - {"goodhit_tof_1", "Good hit time of flight (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit1()"}, - {"goodhit_alpha_0", "Good hit alpha", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit0()"}, - {"goodhit_alpha_1", "Good hit alpha (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit1()"}, + // All below variables are defined and identical to those in ladhodo. Just use H/P.ladhod.var + + // {"n_goodhits", "Number of good hits", "goodhit_n"}, + // {"goodhit_plane_0", "Good hit plane", "fGoodLADHits.THcGoodLADHit.GetPlaneHit0()"}, + // {"goodhit_paddle_0", "Good hit paddle", "fGoodLADHits.THcGoodLADHit.GetPaddleHit0()"}, + // {"goodhit_trackid_0", "Good hit track ID", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit0()"}, + // {"goodhit_beta_0", "Good hit beta", "fGoodLADHits.THcGoodLADHit.GetBetaHit0()"}, + // {"goodhit_dTrkHoriz_0", "Good hit horizontal trk proj - hit position", + // "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit0()"}, + // {"goodhit_dTrkVert_0", "Good hit vertical trk proj - hit position", + // "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit0()"}, + // {"goodhit_hittime_0", "Good hit time", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit0()"}, + // {"goodhit_hittheta_0", "Good hit theta", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit0()"}, + // {"goodhit_hitphi_0", "Good hit phi", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit0()"}, + // {"goodhit_hitedep_0", "Good hit energy deposition", "fGoodLADHits.THcGoodLADHit.GetHitEdepHit0()"}, + // {"goodhit_hitedep_amp_0", "Good hit energy deposition (amplitude)", + // "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit0()"}, + // {"goodhit_plane_1", "Good hit plane (second plane)", "fGoodLADHits.THcGoodLADHit.GetPlaneHit1()"}, + // {"goodhit_paddle_1", "Good hit paddle (second plane)", "fGoodLADHits.THcGoodLADHit.GetPaddleHit1()"}, + // {"goodhit_trackid_1", "Good hit track ID (second plane)", "fGoodLADHits.THcGoodLADHit.GetTrackIDHit1()"}, + // {"goodhit_beta_1", "Good hit beta (second plane)", "fGoodLADHits.THcGoodLADHit.GetBetaHit1()"}, + // {"goodhit_dTrkHoriz_1", "Good hit horizontal trk proj - hit position (second plane)", + // "fGoodLADHits.THcGoodLADHit.GetdTrkHorizHit1()"}, + // {"goodhit_dTrkVert_1", "Good hit vertical trk proj - hit position (second plane)", + // "fGoodLADHits.THcGoodLADHit.GetdTrkVertHit1()"}, + // {"goodhit_hittime_1", "Good hit time (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTimeHit1()"}, + // {"goodhit_hittheta_1", "Good hit theta (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitThetaHit1()"}, + // {"goodhit_hitphi_1", "Good hit phi (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitPhiHit1()"}, + // {"goodhit_hitedep_1", "Good hit energy deposition (second plane)", + // "fGoodLADHits.THcGoodLADHit.GetHitEdepHit1()"}, + // {"goodhit_hitedep_amp_1", "Good hit energy deposition (amplitude) (second plane)", + // "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit1()"}, + // {"goodhit_ypos_0", "Good hit y position", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit0()"}, + // {"goodhit_ypos_1", "Good hit y position (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitYPosHit1()"}, + // {"goodhit_tof_0", "Good hit time of flight", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit0()"}, + // {"goodhit_tof_1", "Good hit time of flight (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitTOFHit1()"}, + // {"goodhit_alpha_0", "Good hit alpha", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit0()"}, + // {"goodhit_alpha_1", "Good hit alpha (second plane)", "fGoodLADHits.THcGoodLADHit.GetHitAlphaHit1()"}, {0}}; return DefineVarsFromList(vars, mode); } @@ -661,7 +675,8 @@ Int_t THcLADKine::DefineVariables(EMode mode) { } //_____________________________________________________________________________ -Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_positions, std::vector sp_resolutions, double dir[3]){ +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 } @@ -672,7 +687,7 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position if (nPoints != sp_resolutions.size()) { return -2; // Mismatch in number of points and resolutions } - for(int i = 0; i < nPoints; i++) { + for (int i = 0; i < nPoints; i++) { if (sp_resolutions[i] <= 0) { return -3; // Invalid resolution value } @@ -680,85 +695,91 @@ Double_t THcLADKine::FitTrack(TVector3 vertex, std::vector sp_position 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){ + // check if the gem track is pointing to the hodoscope hits if there are any, return negative chisq if not + if (nPoints > 2) { TVector3 dir_vec = sp_positions[1] - sp_positions[0]; - dir_vec = dir_vec.Unit(); + 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()) * dir_vec.X() + + // closest approach of the track to the point + 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; - 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 + 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; + // 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]; // in radians - double phi = params[1]; // in 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()) * 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; - 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; + ROOT::Math::Functor f( + [=](const double *params) { + 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); + 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()) * 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; + 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]); // do we want to weight the vertical and horizontal + // residuals differently based on detector performance? } - } - 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); + return chi2_local; + }, + 3); minimizer->SetFunction(f); double initial_params[3] = {dir[0], dir[1], vertex.Z()}; - 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->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) { delete minimizer; return -6; // Fit did not converge } const double *best_params = minimizer->X(); - 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 + 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 delete minimizer; - //return the chi2 of the fit + // return the chi2 of the fit return chi2; } -