diff --git a/src/THcLADHodoscope.cxx b/src/THcLADHodoscope.cxx index 663c4f7..ca9678b 100644 --- a/src/THcLADHodoscope.cxx +++ b/src/THcLADHodoscope.cxx @@ -229,20 +229,18 @@ Int_t THcLADHodoscope::DefineVariables(EMode mode) { {"goodhit_hitedep_0", "Good hit energy deposition", "fGoodLADHits.THcGoodLADHit.GetHitEdepHit0()"}, {"goodhit_hitedep_amp_0", "Good hit energy deposition (amplitude)", "fGoodLADHits.THcGoodLADHit.GetHitEdepAmpHit0()"}, - // Plane 1 not used at this stage. All hits are in plane 0, until matching can occur in THcLADKine - // {"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_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()"}, {"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}}; @@ -315,24 +313,25 @@ Int_t THcLADHodoscope::ReadDatabase(const TDatime &date) { fEdep2MeV_amp[ii] = 1.0; // Default value } - DBRequest list3[] = {{"ladcosmicflag", &fCosmicFlag, kInt, 0, 1}, - {"ladhodo_tdc_min", &fScinTdcMin, kDouble}, - {"ladhodo_tdc_max", &fScinTdcMax, kDouble}, - {"ladhodo_tdc_to_time", &fScinTdcToTime, kDouble}, - {"ladhodo_TopAdcTimeWindowMin", fHodoTopAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1}, - {"ladhodo_TopAdcTimeWindowMax", fHodoTopAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1}, - {"ladhodo_BtmAdcTimeWindowMin", fHodoBtmAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1}, - {"ladhodo_BtmAdcTimeWindowMax", fHodoBtmAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1}, - {"ladhodo_track_tolerance_vert", &fTrackToleranceVert, kDouble}, - {"ladhodo_track_tolerance_horiz", &fTrackToleranceHoriz, kDouble}, - {0}}; - - fTrackToleranceVert = 0.0; - fTrackToleranceHoriz = 0.0; - fCosmicFlag = 0; - fScinTdcMin = 0; - fScinTdcMax = 0; - fScinTdcToTime = 0; + DBRequest list3[] = {{"ladcosmicflag", &fCosmicFlag, kInt, 0, 1}, + {"ladhodo_tdc_min", &fScinTdcMin, kDouble}, + {"ladhodo_tdc_max", &fScinTdcMax, kDouble}, + {"ladhodo_tdc_to_time", &fScinTdcToTime, kDouble}, + {"ladhodo_TopAdcTimeWindowMin", fHodoTopAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1}, + {"ladhodo_TopAdcTimeWindowMax", fHodoTopAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1}, + {"ladhodo_BtmAdcTimeWindowMin", fHodoBtmAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1}, + {"ladhodo_BtmAdcTimeWindowMax", fHodoBtmAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1}, + {"ladhodo_matching_paddle_tol", &fMatchingPaddleTol, kInt, 0, 1}, + {"ladhodo_matching_time_tol", &fMatchingTimeTol, kDouble, 0, 1}, + {"ladhodo_matching_dy_tol", &fMatchingDyTol, kDouble, 0, 1}, + {0}}; + fMatchingPaddleTol = 0; + fMatchingTimeTol = 10.0; + fMatchingDyTol = 20.0; + fCosmicFlag = 0; + fScinTdcMin = 0; + fScinTdcMax = 0; + fScinTdcToTime = 0; gHcParms->LoadParmValues((DBRequest *)&list3, prefix_lad); @@ -451,7 +450,7 @@ Int_t THcLADHodoscope::CoarseProcess(TClonesArray &tracks) { return 0; } Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { fGoodLADHits->Delete(); - // Hodo tracking coes in FineProcess, to ensure it comes after GEM Coarse Process tracking, no matter what. + // Hodo tracking goes in FineProcess, to ensure it comes after GEM Coarse Process tracking, no matter what. fSpectro = dynamic_cast(GetApparatus()); fGEM = dynamic_cast(fSpectro->GetDetector("gem")); @@ -464,132 +463,107 @@ Int_t THcLADHodoscope::FineProcess(TClonesArray &tracks) { return 0; } - // TClonesArray* gemTracks = spectrometer->GetGEM()->GetTracks(); - Int_t ntracks = fGEMTracks->GetLast() + 1; - if (ntracks > 0) { - vector nPmtHit(ntracks); - vector timeAtFP(ntracks); - // Loop over all tracks and get corrected time, tof, beta... - num_unique_hits = 0; - for (Int_t ip = 0; ip < fNPlanes; ip++) { - if (strcmp(fPlanes[ip]->GetName(), "REFBAR") == 0) { - continue; - } - num_unique_hits += fPlanes[ip]->GetNScinHits(); - } - num_unique_good_hits = 0; - std::map, bool> hitUsedMap; + vector nPmtHit(ntracks); + vector timeAtFP(ntracks); - // Initialize the map for all hits - for (Int_t ip = 0; ip < fNPlanes; ip++) { - for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) { - THcLADHodoHit *hit = (THcLADHodoHit *)fPlanes[ip]->GetHits()->At(iphit); - int paddle = hit->GetPaddleNumber() - 1; - hitUsedMap[std::make_pair(ip, paddle)] = false; - } + // Loop over all tracks and get corrected time, tof, beta... + num_unique_hits = 0; + for (Int_t ip = 0; ip < fNPlanes; ip++) { + if (strcmp(fPlanes[ip]->GetName(), "REFBAR") == 0) { + continue; } - // Loop over all tracks - for (Int_t itrack = 0; itrack < ntracks; itrack++) { // Line 133 - nPmtHit[itrack] = 0; - timeAtFP[itrack] = 0; - - THcLADGEMTrack *theTrack = dynamic_cast(fGEMTracks->At(itrack)); - if (!theTrack) - return -1; - - // Calculate Theta, Phi, X and Y (intercepts with z=0 plane) from THcLADGEMTrack - TVector3 sp1(theTrack->GetX1(), theTrack->GetY1(), theTrack->GetZ1()); - TVector3 sp2(theTrack->GetX2(), theTrack->GetY2(), theTrack->GetZ2()); - - Double_t track_dz = sp2.Z() - sp1.Z(); - if (track_dz == 0) { - return 0; - } - Double_t track_dx = sp2.X() - sp1.X(); - Double_t track_dy = sp2.Y() - sp1.Y(); - - Double_t track_theta = TMath::ATan2(TMath::Sqrt(track_dx * track_dx + track_dy * track_dy), track_dz); - Double_t track_phi = TMath::ATan2(track_dy, track_dx); - - Double_t track_x0 = sp1.X() - sp1.Z() * track_dx / track_dz; - Double_t track_y0 = sp1.Y() - sp1.Z() * track_dy / track_dz; - - // Loop over all scintillator planes - for (Int_t ip = 0; ip < fNPlanes; ip++) { - if (strcmp(fPlanes[ip]->GetName(), "REFBAR") == 0) { - continue; - } - - Int_t fNScinHits = fPlanes[ip]->GetNScinHits(); - TClonesArray *hodoHits = fPlanes[ip]->GetHits(); - - Double_t zPos = fPlanes[ip]->GetZpos(); - Double_t dzPos = fPlanes[ip]->GetDzpos(); - Double_t planeTheta = fPlanes[ip]->GetTheta(); - - // Loop over hits with in a single plane - for (Int_t iphit = 0; iphit < fNScinHits; iphit++) { - // iphit is hit # within a plane - THcLADHodoHit *hit = (THcLADHodoHit *)hodoHits->At(iphit); - - Int_t paddle = hit->GetPaddleNumber() - 1; - Double_t zposition = zPos; // TODO: fix this (lad won't have this offset) - - Double_t track_TrnsCoord, track_LongCoord; - // track_TrnsCoord = - track_x0 * TMath::Sin(planeTheta - TMath::Pi() / 2); - track_TrnsCoord = track_x0 * TMath::Sin(track_theta - TMath::Pi() / 2) / - TMath::Sin(TMath::Pi() / 2 - track_theta + planeTheta); - track_LongCoord = track_y0; - - Double_t scinTrnsCoord, scinLongCoord; - // x & y cooridnates are with respect to the central angle pointing to the scintillator plane - scinTrnsCoord = -track_TrnsCoord + TMath::Tan(track_theta - planeTheta) * (zposition); // Line 183 - - scinLongCoord = (-track_LongCoord + - TMath::Tan(track_phi) / TMath::Cos(track_theta - planeTheta) * (zposition)); // Line 184 + num_unique_hits += fPlanes[ip]->GetNScinHits(); + } + num_unique_good_hits = 0; - Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset(); + // Key by (plane, hit_index) instead of (plane, paddle) so that multiple + // hits sharing the same paddle in the same plane are each tracked + // independently. + std::map, bool> hitUsedMap; - if ((TMath::Abs(scinCenter - scinTrnsCoord) < (fPlanes[ip]->GetSize() * 0.5 + fTrackToleranceHoriz)) && - (TMath::Abs(scinLongCoord - hit->GetCalcPosition()) < fTrackToleranceVert)) { - if (goodhit_n >= MAXGOODHITS) { - // cout << "Error: Too many \"good hits\"" << endl; - return -1; - } - theTrack->SetHasHodoHit(true); - // Check if the hit has already been used - if ((TMath::Abs(scinCenter - scinTrnsCoord) < 100)) // Hardcoded tolerance. Fix later - hitUsedMap[std::make_pair(ip, paddle)] = true; // Mark this hit as used + // Initialize the map for all hits + for (Int_t ip = 0; ip < fNPlanes; ip++) { + for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) { + hitUsedMap[std::make_pair(ip, iphit)] = false; + } + } - THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); + // Pair front and back hits to make good hits + for (Int_t ip = 0; ip < fNPlanes; ip += 2) { + for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) { + THcLADHodoHit *hit = (THcLADHodoHit *)fPlanes[ip]->GetHits()->At(iphit); + if (hitUsedMap[std::make_pair(ip, iphit)]) { + continue; + } + if (ip + 1 < fNPlanes) { + for (Int_t iphit2 = 0; iphit2 < fPlanes[ip + 1]->GetNScinHits(); iphit2++) { + THcLADHodoHit *hit2 = (THcLADHodoHit *)fPlanes[ip + 1]->GetHits()->At(iphit2); + if (hitUsedMap[std::make_pair(ip + 1, iphit2)]) { + continue; + } + + bool same_paddle = fabs(hit2->GetPaddleNumber() - hit->GetPaddleNumber()) <= fMatchingPaddleTol; + bool time_match = fabs(hit2->GetScinCorrectedTime() - hit->GetScinCorrectedTime()) < fMatchingTimeTol; + bool dy_match = fabs(hit2->GetCalcPosition() - hit->GetCalcPosition()) < fMatchingDyTol; + if (same_paddle && time_match && dy_match) { + hitUsedMap[std::make_pair(ip, iphit)] = true; + hitUsedMap[std::make_pair(ip + 1, iphit2)] = true; + THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); goodhit_n++; goodhit->SetPlane(0, ip); - goodhit->SetPaddle(0, paddle); - goodhit->SetTrackID(0, itrack); - goodhit->SetdTrkHoriz(0, scinCenter - scinTrnsCoord); - goodhit->SetdTrkVert(0, scinLongCoord - hit->GetCalcPosition()); + goodhit->SetPaddle(0, hit->GetPaddleNumber() - 1); goodhit->SetHitTime(0, hit->GetScinCorrectedTime()); - goodhit->SetHitTheta(0, track_theta); - goodhit->SetHitPhi(0, track_phi); goodhit->SetHitEdep(0, hit->GetPaddleADC()); goodhit->SetHitEdepAmp(0, hit->GetPaddleADCpeak()); goodhit->SetHitYPos(0, hit->GetCalcPosition()); - - } // condition for cenetr on a paddle - // ihhit++; - } // First loop over hits in a plane <--------- + goodhit->SetPlane(1, ip + 1); + goodhit->SetPaddle(1, hit2->GetPaddleNumber() - 1); + goodhit->SetHitTime(1, hit2->GetScinCorrectedTime()); + goodhit->SetHitEdep(1, hit2->GetPaddleADC()); + goodhit->SetHitEdepAmp(1, hit2->GetPaddleADCpeak()); + goodhit->SetHitYPos(1, hit2->GetCalcPosition()); + break; // front hit is now used; stop looking for more back matches for it + } + } } + } + } - } // Main loop over tracks ends here. - num_unique_good_hits = 0; - for (const auto &hit : hitUsedMap) { - if (hit.second) { - num_unique_good_hits++; - } + for (Int_t ip = 0; ip < fNPlanes; ip++) { + // skip reference bar plane + if (strcmp(fPlanes[ip]->GetName(), "REFBAR") == 0) + continue; + + for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++) { + THcLADHodoHit *hit = (THcLADHodoHit *)fPlanes[ip]->GetHits()->At(iphit); + auto key = std::make_pair(ip, iphit); + if (hitUsedMap[key]) + continue; + + // mark as used and create a single-plane good hit + hitUsedMap[key] = true; + THcGoodLADHit *goodhit = new ((*fGoodLADHits)[goodhit_n]) THcGoodLADHit(); + goodhit_n++; + + // Fill front planes (0,2,4) into index 0, back planes (1,3) into index 1 + Int_t hitIndex = (ip % 2 == 0) ? 0 : 1; + goodhit->SetPlane(hitIndex, ip); + goodhit->SetPaddle(hitIndex, hit->GetPaddleNumber() - 1); + goodhit->SetHitTime(hitIndex, hit->GetScinCorrectedTime()); + goodhit->SetHitEdep(hitIndex, hit->GetPaddleADC()); + goodhit->SetHitEdepAmp(hitIndex, hit->GetPaddleADCpeak()); + goodhit->SetHitYPos(hitIndex, hit->GetCalcPosition()); } - } // If condition for at least one track + } + + num_unique_good_hits = 0; + for (const auto &hit : hitUsedMap) { + if (hit.second) { + num_unique_good_hits++; + } + } return 0; } diff --git a/src/THcLADHodoscope.h b/src/THcLADHodoscope.h index 3115abe..3d3d654 100644 --- a/src/THcLADHodoscope.h +++ b/src/THcLADHodoscope.h @@ -5,12 +5,12 @@ #include "TH1F.h" #include "THaNonTrackingDetector.h" #include "THaSpectrometer.h" +#include "THcGoodLADHit.h" #include "THcHitList.h" #include "THcLADGEM.h" #include "THcLADHodoHit.h" #include "THcLADHodoPlane.h" #include "THcRawHodoHit.h" -#include "THcGoodLADHit.h" class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { public: @@ -77,8 +77,9 @@ class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { Double_t *fAdcTdcOffset; Double_t *fHodoSlop; Int_t fIsMC; - Double_t fTrackToleranceVert; - Double_t fTrackToleranceHoriz; + Int_t fMatchingPaddleTol; + Double_t fMatchingTimeTol; + Double_t fMatchingDyTol; // Output variables static const Int_t MAXGOODHITS = 500; @@ -91,8 +92,7 @@ class THcLADHodoscope : public THaNonTrackingDetector, public THcHitList { THaApparatus *fSpectro; THcLADGEM *fGEM; - - Int_t *fNScinHits; // [fNPlanes] + Int_t *fNScinHits; // [fNPlanes] // Time walk Double_t *fHodoVelFit; diff --git a/src/THcLADKine.cxx b/src/THcLADKine.cxx index 6187c00..3cd102b 100644 --- a/src/THcLADKine.cxx +++ b/src/THcLADKine.cxx @@ -110,6 +110,8 @@ THaAnalysisObject::EStatus THcLADKine::Init(const TDatime &run_time) { return fStatus; } + + return THaPhysicsModule::Init(run_time); } //_____________________________________________________________________________ @@ -131,16 +133,7 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { fglobal_time_offset = 0.0; fTrk_dtCut = 10.0; - // delete[] rf_offset; - rf_offset = new Double_t[2]; // shouldn't be length 2 FIXME - cout << "Reading LAD Kinematics parameters from database..." << endl; - THaVar *varptr; - varptr = gHcParms->Find("gen_run_number"); - Int_t *runnum; - if (varptr) { - runnum = (Int_t *)varptr->GetValuePointer(); // Assume correct type - } DBRequest list[] = {{"d0_cut_wVertex", &fD0Cut_wVertex, kDouble, 0, 1}, {"d0_cut_noVertex", &fD0Cut_noVertex, kDouble, 0, 1}, @@ -149,27 +142,19 @@ Int_t THcLADKine::ReadDatabase(const TDatime &date) { {"nfixed_z", &fNfixed_z, kInt, 0, 1}, {"global_time_offset", &fglobal_time_offset, kDouble, 0, 1}, {"trk_dt_cut", &fTrk_dtCut, kDouble, 0, 1}, - // {"_rf_offset", rf_offset, kDouble, 2}, {"_rf_period", &rf_period, kDouble, 0, 1}, {0}}; gHcParms->LoadParmValues((DBRequest *)&list, prefix); - DBRequest list_rf[] = {{"l_nr_rf_offsets", &n_rf_offsets, kInt, 0, 1}, {0}}; + DBRequest list_rf[] = {{"_nr_rf_offsets", &n_rf_offsets, kInt, 0, 1}, {0}}; gHcParms->LoadParmValues((DBRequest *)&list_rf, prefix); - rf_offset = new Double_t[n_rf_offsets]; - for (Int_t i = 0; i < n_rf_offsets; i++) { + rf_offset = new Double_t[n_rf_offsets * 3]; + for (Int_t i = 0; i < n_rf_offsets * 3; i++) { rf_offset[i] = 0.0; } - DBRequest list_rf_offsets[] = {{"_rf_offset", rf_offset, kDouble, n_rf_offsets}, {0}}; + DBRequest list_rf_offsets[] = {{"_rf_offset", &rf_offset[0], kDouble, (UInt_t)(n_rf_offsets * 3)}, {0}}; gHcParms->LoadParmValues((DBRequest *)&list_rf_offsets, prefix); - gHcParms->LoadParmValues((DBRequest *)&list_rf, prefix); - - // else { - // runnum = new Int_t[1]; - // gHcParms->Define("gen_run_number","Run Number", *runnum); - // } - delete[] fFixed_z; if (fNfixed_z > 0) { fFixed_z = new Double_t[fNfixed_z]; @@ -201,7 +186,6 @@ Int_t THcLADKine::Process(const THaEvData &evdata) { TVector3 v_hit1(track->GetX1(), track->GetY1(), track->GetZ1()); TVector3 v_hit2(track->GetX2(), track->GetY2(), track->GetZ2()); - TVector3 vertex; if (fVertexModule->HasVertex()) { vertex = fVertexModule->GetVertex(); track->SetGoodD0(kTRUE); @@ -475,13 +459,33 @@ void THcLADKine::CalculateTVertex() { fTVertex = fPtime - PathLengthCorr; - int fRFTimeIndex = (aparatus_prefix[0] == 'p') ? 0 : 1; - fRFTime = fTrigDet->Get_RF_TrigTime(fRFTimeIndex) + rf_offset[fRFTimeIndex]; - double tmp = fmod(fTVertex - fRFTime, rf_period); + THaVar *varptr; + varptr = gHcParms->Find("gen_run_number"); + Int_t *runnum; + if (varptr) { + runnum = (Int_t *)varptr->GetValuePointer(); // Assume correct type + } + int fRFTimeIndex = -1; + int run_idx = 0; + while (fRFTimeIndex < 0) { + if (run_idx >= 3 * n_rf_offsets || rf_offset[run_idx] > *runnum) { + fRFTimeIndex = run_idx - 3; // Go back 1 step (3 indices) to get the correct RF offset for the run + break; + } + run_idx += 3; + } + int fRFSpecIndex = (aparatus_prefix[0] == 'p') ? 0 : 1; + fRFTimeIndex += (1 + fRFSpecIndex); + // (aparatus_prefix[0] == 'p') ? 0 : 1; + + fRFTime = fTrigDet->Get_RF_TrigTime(fRFSpecIndex); + double tmp = fmod(fTVertex - fRFTime + rf_offset[fRFTimeIndex], rf_period); if (tmp > 2) tmp -= rf_period; fTVertex_RFcorr = fTVertex - tmp; + // cout << "Offset" << rf_offset[fRFTimeIndex] << " Offset Index " << fRFTimeIndex << " RF Time " << fRFTime + // << " TVertex before RF corr " << fTVertex << " TVertex after RF corr " << fTVertex_RFcorr << endl; return; } //_____________________________________________________________________________ @@ -495,6 +499,8 @@ Double_t THcLADKine::CalculateTOFRFcorr(Double_t t_raw) { Double_t tof = t_raw - fTVertex_RFcorr + fglobal_time_offset; + tof -= vertex.Z() / lightSpeed; + return tof; } //_____________________________________________________________________________ diff --git a/src/THcLADKine.h b/src/THcLADKine.h index 2ed9ecb..38a338a 100644 --- a/src/THcLADKine.h +++ b/src/THcLADKine.h @@ -56,6 +56,8 @@ class THcLADKine : public THcPrimaryKine { Int_t n_rf_offsets; Double_t *rf_offset; Double_t rf_period; + TVector3 vertex; + Double_t lightSpeed = 29.9792458; // Speed of light in cm/ns virtual Int_t DefineVariables(EMode mode = kDefine); void CalculateTVertex(); Double_t CalculateToF(Double_t t_raw);