diff --git a/src/THcLADGEM.cxx b/src/THcLADGEM.cxx index f6049d8..b926d02 100644 --- a/src/THcLADGEM.cxx +++ b/src/THcLADGEM.cxx @@ -230,10 +230,13 @@ 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()"}, + {"trk.is_good", "Track is good", "fGEMTracks.THcLADGEMTrack.IsGoodTrack()"}, {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..dd8c902 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: @@ -95,23 +96,47 @@ 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; } 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; } 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; } + 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;}; @@ -120,11 +145,11 @@ 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]; } - protected: Int_t fNSp; Double_t fProjVz; Double_t fProjVy; + Double_t fProjVx; Double_t fD0; Bool_t fhasGoodD0; Int_t fTrackID; @@ -132,9 +157,13 @@ 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 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 private: THcLADGEMTrack(const THcLADGEMTrack &); 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/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index 663c4f7..76561a3 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -594,4 +594,22 @@ Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { return 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) { + 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]->GetZpos()); + hit.RotateY(fPlanes[plane]->GetTheta()); + return hit; +} + + ClassImp(THcLADHodoscope) diff --git a/src/THcLADHodoscope.h b/src/THcLADHodoscope.h index 3115abe..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 ); TClonesArray *GetLADGoodHits() { return fGoodLADHits; } TClonesArray *GetLADHits(Int_t plane) { return fPlanes[plane]->GetHits(); } diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 6187c00..2ed247c 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) //_____________________________________________________________________________ @@ -27,6 +31,19 @@ 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*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() { @@ -179,6 +196,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 }, + { "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; + } + 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(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; } //_____________________________________________________________________________ @@ -187,13 +229,17 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { fTrack = fSpectro->GetGoldenTrack(); CalculateTVertex(); + //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; @@ -208,7 +254,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,158 +266,216 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { vertex.SetZ(fFixed_z[min_index]); } - // 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) { - isGoodTrack[i] = true; - } - } else { - if (d0 < fD0Cut_noVertex) { - isGoodTrack[i] = true; - } - } - - if (track->GetdT() > fTrk_dtCut) { - isGoodTrack[i] = 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 || !isGoodTrack[track_id]) + 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 + 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; + 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; - - // 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) + } + 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; - - // 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) { - matching_hit_index = j; - break; + // Check the number hits in goohit + int num_hodo_hits = 0; + if (goodHit->GetPlaneHit0() >= 0 && goodHit->GetPlaneHit0() < 999) { + num_hodo_hits++; } - // 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 (goodHit->GetPlaneHit1() >= 0 && goodHit->GetPlaneHit1() < 999) { + num_hodo_hits++; } - } - 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); + 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); + + 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); + v_resolutions.push_back(fSigma_Hodo); } + 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); + if (num_hodo_hits==1){ + if (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> 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()); + //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; } - if (goodhit->GetdTrkHorizHit1() != 0) { - dTrk_horiz_values.push_back(goodhit->GetdTrkHorizHit1()); + 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->SetD0(d0); - 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; + 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++; + if (besthit == nullptr || newhit == nullptr) { + track->SetBestHodoHit(nullptr); + continue; } - } - 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; + 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 || 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; } - } - // 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--; + + + if (track->GetGoodD0()) { + if (track->GetD0()>0 && track->GetD0() < fD0Cut_wVertex) { + isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; } + } else { + if (track->GetD0() > 0 && track->GetD0() < fD0Cut_noVertex) { + isGoodTrack[iTrack] = isGoodTrack[iTrack] && true; + } + } + if(track->GetChisq()<0) { + isGoodTrack[iTrack] = false; } + if (track->GetdT() > fTrk_dtCut) { + isGoodTrack[iTrack] = false; + } + track->SetIsGoodTrack(isGoodTrack[iTrack]); } - 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)); @@ -550,3 +654,105 @@ 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 || nPoints > 4) { + 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] > 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){ + 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()) * 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 + } + } + } + + + //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; + } + } + 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->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 + delete minimizer; + + //return the chi2 of the fit + return chi2; +} + diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 2ed9ecb..83aac0d 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -56,10 +56,27 @@ 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 + + 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 + + + + virtual Int_t DefineVariables(EMode mode = kDefine); 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) };