diff --git a/shipLHC/MuFilter.cxx b/shipLHC/MuFilter.cxx index bc28e2f97d..a7a132bc59 100644 --- a/shipLHC/MuFilter.cxx +++ b/shipLHC/MuFilter.cxx @@ -575,14 +575,20 @@ void MuFilter::GetLocalPosition(Int_t fDetectorID, TVector3& vLeft, TVector3& vR vRight.SetXYZ(locB[0],locB[1],locB[2]); } -Float_t MuFilter::GetCorrectedTime(Int_t fDetectorID, Int_t channel, Double_t rawTime, Double_t L){ +// Calculate the time walk correction +Float_t MuFilter::GetTimeWalk(Int_t fDetectorID, Int_t channel, Float_t qdc, TString tag){ + std::vector vec = conf_vectors["MuFilter/TW_DSa_"+std::to_string(fDetectorID) + +"_"+std::to_string(channel)+tag]; + /* A = vec[0]; B = vec[1]; C = vec[2]; + D = vec[3]; E = vec[4]; F = vec[5]; */ + Float_t Q = qdc-vec[0]; + return vec[3]*Q/(vec[1]+vec[2]*pow(Q,2)) + vec[4]*Q + vec[5]; +} + +Float_t MuFilter::GetCorrectedTime(Int_t fDetectorID, Int_t channel, Double_t rawTime, Double_t L, Float_t qdc){ /* expect time in u.ns and path length to sipm u.cm */ -/* calibration implemented only for DS! */ -/* channel 0 left or top, channel 1 right */ - if (fDetectorID<30000){ - LOG(ERROR) << "MuFilter::GetCorrectedTime: not yet implemented for Veto and DS, no correction applied" ; - return rawTime; - } +/* returns the time in u.ns */ +/* for DS: channel 0 left or top, channel 1 right */ TString tag = ""; if (eventHeader){ Int_t fRunNumber = eventHeader->GetRunId(); @@ -614,6 +620,12 @@ Float_t MuFilter::GetCorrectedTime(Int_t fDetectorID, Int_t channel, Double_t ra last_time_alignment_tag = tag; } } + if (fDetectorID<30000){ + return rawTime + GetTimeWalk(fDetectorID, channel, qdc, last_time_alignment_tag) + + conf_vectors["MuFilter/TW_DSa_" + +std::to_string(fDetectorID)+"_"+std::to_string(channel) + +last_time_alignment_tag].back(); + } Float_t cor = rawTime; int l = (fDetectorID-30000)/1000; int ichannel60 = fDetectorID%1000; diff --git a/shipLHC/MuFilter.h b/shipLHC/MuFilter.h index 63577c1eea..4d2f9398fc 100644 --- a/shipLHC/MuFilter.h +++ b/shipLHC/MuFilter.h @@ -35,7 +35,8 @@ class MuFilter : public FairDetector /** Getposition **/ void GetPosition(Int_t id, TVector3& vLeft, TVector3& vRight); // or top and bottom void GetLocalPosition(Int_t id, TVector3& vLeft, TVector3& vRight); - Float_t GetCorrectedTime(Int_t id, Int_t c, Double_t t, Double_t L); + Float_t GetTimeWalk(Int_t id, Int_t c, Float_t qdc, TString tag); + Float_t GetCorrectedTime(Int_t id, Int_t c, Double_t t, Double_t L=0, Float_t QDC=-1.); Int_t GetnSiPMs(Int_t detID); Int_t GetnSides(Int_t detID); @@ -43,9 +44,11 @@ class MuFilter : public FairDetector void SetConfPar(TString name, Float_t value){conf_floats[name]=value;} void SetConfPar(TString name, Int_t value){conf_ints[name]=value;} void SetConfPar(TString name, TString value){conf_strings[name]=value;} + void SetConfPar(TString name, std::vector values){conf_vectors[name]=values;} Float_t GetConfParF(TString name){return conf_floats[name];} Int_t GetConfParI(TString name){return conf_ints[name];} TString GetConfParS(TString name){return conf_strings[name];} + std::vector GetConfParVector(TString name){return conf_vectors[name];} /** Initialization of the detector is done here */ virtual void Initialize(); @@ -101,6 +104,7 @@ class MuFilter : public FairDetector std::map conf_floats; std::map conf_ints; std::map conf_strings; + std::map> conf_vectors; SNDLHCEventHeader *eventHeader; // Vector to store runs covered in the geometry file. diff --git a/shipLHC/MuFilterHit.cxx b/shipLHC/MuFilterHit.cxx index 13577f3f43..c74129e0ad 100644 --- a/shipLHC/MuFilterHit.cxx +++ b/shipLHC/MuFilterHit.cxx @@ -1,5 +1,6 @@ #include "MuFilterHit.h" #include "MuFilter.h" +#include "ShipUnit.h" #include "TROOT.h" #include "FairRunSim.h" #include "TGeoNavigator.h" @@ -8,7 +9,6 @@ #include #include -Double_t speedOfLight = TMath::C() *100./1000000000.0 ; // from m/sec to cm/ns // ----- Default constructor ------------------------------------------- MuFilterHit::MuFilterHit() : SndlhcHit() @@ -167,15 +167,16 @@ std::map MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive) } // ----- Public method Get List of time measurements ------------------------------------------- -std::map MuFilterHit::GetAllTimes(Bool_t mask) +std::map MuFilterHit::GetAllTimes(Bool_t mask, Bool_t apply_t_corr, Double_t SipmDistance) { + OptForTimeCorrections(mask, apply_t_corr, SipmDistance); std::map allTimes; for (unsigned int s=0; s 0){ if (!fMasked[channel] || !mask){ - allTimes[channel] = times[channel]; + allTimes[channel] = fTimesHelper[channel]; } } } @@ -184,9 +185,10 @@ std::map MuFilterHit::GetAllTimes(Bool_t mask) } // ----- Public method Get time difference mean Left - mean Right ----------------- -Float_t MuFilterHit::GetDeltaT(Bool_t mask) +Float_t MuFilterHit::GetDeltaT(Bool_t mask, Bool_t apply_t_corr, Double_t SipmDistance) // based on mean TDC measured on Left and Right { + OptForTimeCorrections(mask, apply_t_corr, SipmDistance); Float_t mean[] = {0,0}; Int_t count[] = {0,0}; Float_t dT = -999.; @@ -195,7 +197,7 @@ Float_t MuFilterHit::GetDeltaT(Bool_t mask) unsigned int channel = j+s*nSiPMs; if (signals[channel]> 0){ if (!fMasked[channel] || !mask){ - mean[s] += times[channel]; + mean[s] += fTimesHelper[channel]; count[s] += 1; } } @@ -206,9 +208,10 @@ Float_t MuFilterHit::GetDeltaT(Bool_t mask) } return dT; } -Float_t MuFilterHit::GetFastDeltaT(Bool_t mask) +Float_t MuFilterHit::GetFastDeltaT(Bool_t mask, Bool_t apply_t_corr, Double_t SipmDistance) // based on fastest (earliest) TDC measured on Left and Right { + OptForTimeCorrections(mask, apply_t_corr, SipmDistance); Float_t first[] = {1E20,1E20}; Float_t dT = -999.; for (unsigned int s=0; s 0){ if (!fMasked[channel] || !mask){ - if (times[channel] 0){ if (!fMasked[channel] || !mask){ - mean[s] += times[channel]; + mean[s] += fTimesHelper[channel]; count[s] += 1; } } } } if (count[0]>0 && count[1]>0) { - dT = (mean[0]/count[0] + mean[1]/count[1])/2.*6.25 - dL/2.; // TDC to ns = 6.25 + dT = (mean[0]/count[0] + mean[1]/count[1])/2.*ShipUnit::snd_TDC2ns - dL/2.; // TDC to ns = 6.25 } return dT; } @@ -319,6 +323,30 @@ void MuFilterHit::Print() const } std::cout << std::endl; } + +// ----- Private method to opt for time corrections ------------------------------------------- +void MuFilterHit::OptForTimeCorrections(Bool_t mask, Bool_t apply, Double_t SipmDistance) +{ + // simply copy the times if no time corrections is required + if (apply == kFALSE) memcpy(fTimesHelper, times, sizeof(fTimesHelper)); + else { // apply the time corrections + MuFilter* MuFilterDet = dynamic_cast (gROOT->GetListOfGlobals()->FindObject("MuFilter")); + for (unsigned int s=0; s 0){ + if (!fMasked[channel] || !mask){ + fTimesHelper[channel] = (MuFilterDet->GetCorrectedTime(fDetectorID,channel, + times[channel]*ShipUnit::snd_TDC2ns, + SipmDistance, + signals[channel])/ShipUnit::snd_TDC2ns); + } + } + else fTimesHelper[channel] = times[channel]; + } + } + } +} // ------------------------------------------------------------------------- ClassImp(MuFilterHit) diff --git a/shipLHC/MuFilterHit.h b/shipLHC/MuFilterHit.h index 835774dd13..260a474c79 100644 --- a/shipLHC/MuFilterHit.h +++ b/shipLHC/MuFilterHit.h @@ -29,10 +29,12 @@ class MuFilterHit : public SndlhcHit Float_t SumOfSignals(char* opt,Bool_t mask=kTRUE); std::map SumOfSignals(Bool_t mask=kTRUE); std::map GetAllSignals(Bool_t mask=kTRUE,Bool_t positive=kTRUE); - std::map GetAllTimes(Bool_t mask=kTRUE); - Float_t GetDeltaT(Bool_t mask=kTRUE); - Float_t GetFastDeltaT(Bool_t mask=kTRUE); - Float_t GetImpactT(Bool_t mask=kTRUE); + std::map GetAllTimes(Bool_t mask=kTRUE, + Bool_t apply_t_corr=kFALSE, + Double_t SipmDistance=0.); + Float_t GetDeltaT(Bool_t mask=kTRUE, Bool_t apply_t_corr=kFALSE, Double_t SipmDistance=0.); + Float_t GetFastDeltaT(Bool_t mask=kTRUE, Bool_t apply_t_corr=kFALSE, Double_t SipmDistance=0.); + Float_t GetImpactT(Bool_t mask=kTRUE, Bool_t apply_t_corr=kFALSE, Double_t SipmDistance=0.); bool isValid() const {return flag;} bool isMasked(Int_t i) const {return fMasked[i];} void SetMasked(Int_t i) {fMasked[i]=kTRUE;} @@ -44,11 +46,18 @@ class MuFilterHit : public SndlhcHit /** Copy constructor **/ MuFilterHit(const MuFilterHit& hit); MuFilterHit operator=(const MuFilterHit& hit); + /** Function to set the fTimesHelper array **/ + void OptForTimeCorrections(Bool_t mask, Bool_t apply, Double_t SipmDistance); Float_t flag; ///< flag Float_t fMasked[16]; /// masked signal + /* Helper container for time-corrected or raw times + depending if time alignment is requested. + The unit is clock cycles same as for the SndlhcHit's + times[16] data member */ + Float_t fTimesHelper[16]; - ClassDef(MuFilterHit,5); + ClassDef(MuFilterHit,6); }; diff --git a/shipLHC/modifyGeoFileDict.py b/shipLHC/modifyGeoFileDict.py index ec12c722f2..8827e09b8e 100644 --- a/shipLHC/modifyGeoFileDict.py +++ b/shipLHC/modifyGeoFileDict.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +import json import ROOT,os from rootpyPickler import Pickler from rootpyPickler import Unpickler @@ -21,7 +22,50 @@ def modifyDicts(): pkl = Unpickler(fg) sGeo = pkl.load('ShipGeo') fg.Close() - + + # Veto and US part + sample_runs = ['5097','5408','5999'] + tag_runs = ["4361","5117", "5478"] + for n, run in enumerate(sample_runs): + fa = open('/eos/experiment/sndlhc/calibration/MuFilter/VetoUS/alignmentparams_run'+run+'.json') + ftw = open('/eos/experiment/sndlhc/calibration/MuFilter/VetoUS/twparams_run'+run+'.json') + d_tw = json.load(ftw) + d_a = json.load(fa) + constants_twa= {} + for i in d_tw.items(): + detID = int(i[0].split('_')[0]) + chanID = int(i[0].split('_')[1]) + # safety net for missing tw constants + if len(i[1]) != 6: + print('Got', len(i[1]),'out of 6 time walk parameters for', detID, chanID, + '\nCheck json values!') + return + if detID not in constants_twa: constants_twa[detID]={} + constants_twa[detID][chanID]=[float(i[1][0]),float(i[1][1]),float(i[1][2]), + float(i[1][3]),float(i[1][4]),float(i[1][5])] + # the part relative to DS is put at the vector's end + for i in d_a.items(): + detID = int(i[0].split('_')[0]) + chanID = int(i[0].split('_')[1]) + constants_twa[detID][chanID].append(float(i[1][0])) + # put zeroes for uncalibrated channels + detID_list = [] + subsystem = {"Veto":1,"Upstream":2} + for sys in subsystem: + for nplanes in range(sGeo.MuFilter['N'+sys+'Planes']): + for nbar in range(sGeo.MuFilter['N'+sys+'Bars']): + detID_list.append(int(subsystem[sys]*1e4+nplanes*1e3+nbar)) + for detID in detID_list: + if detID not in constants_twa: constants_twa[detID]={} + for chanID in range(16): + if chanID not in constants_twa[detID]: + constants_twa[detID][chanID] = [0,1E-20,0,0,0,0,0] + # put these into the detector map before going to next run + for detID in constants_twa: + for chanID in constants_twa[detID]: + setattr(sGeo.MuFilter,'TW_DSa_'+str(detID)+'_'+str(chanID)+'t_'+tag_runs[n], constants_twa[detID][chanID]) + sGeo.MuFilter['TW_DSa_'+str(detID)+'_'+str(chanID)+'t_'+tag_runs[n]] = constants_twa[detID][chanID] + # DS part setattr(sGeo.MuFilter,'DsPropSpeed',14.9*u.cm/u.nanosecond) sGeo.MuFilter['DsPropSpeed'] = 14.9*u.cm/u.nanosecond