From 5271002e5fc8d89c518e9a89396c674982ca70e7 Mon Sep 17 00:00:00 2001 From: siilieva Date: Wed, 13 Mar 2024 11:18:30 +0100 Subject: [PATCH 1/4] Add Veto and US time alignment There is the time walk and the time alignment wrt DS. Needed also a new method to read vectors from the geometry, which can then be used for parameter structures of arbitary size. --- shipLHC/MuFilter.cxx | 26 +++++++++++++++++++------- shipLHC/MuFilter.h | 6 +++++- 2 files changed, 24 insertions(+), 8 deletions(-) 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. From 260d4daeba8309ea22e764529a2523a724d94c71 Mon Sep 17 00:00:00 2001 From: siilieva Date: Wed, 13 Mar 2024 11:20:43 +0100 Subject: [PATCH 2/4] Add methods to apply the Veto and US time corrections Now every time-related method is provided with 1. a flag whether to apply time alignment or not 2. the distance to the SiPM so to correct for the light propagation time. Defaults of these are respectively false and 0, so to have backward compatibility. --- shipLHC/MuFilterHit.cxx | 47 +++++++++++++++++++++++++++++++++-------- shipLHC/MuFilterHit.h | 19 ++++++++++++----- 2 files changed, 52 insertions(+), 14 deletions(-) diff --git a/shipLHC/MuFilterHit.cxx b/shipLHC/MuFilterHit.cxx index 13577f3f43..61109b537a 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" @@ -167,15 +168,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 +186,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 +198,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 +209,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 +324,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); }; From a545b2443dc001302cfbea9e0da56f4b7138f0f6 Mon Sep 17 00:00:00 2001 From: siilieva Date: Wed, 13 Mar 2024 11:23:40 +0100 Subject: [PATCH 3/4] Read the constants for the US and Veto time alignment --- shipLHC/modifyGeoFileDict.py | 46 +++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) 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 From f615edf0a85b89d8e82bce6fa8e5b869f70f5c56 Mon Sep 17 00:00:00 2001 From: siilieva Date: Thu, 21 Mar 2024 16:13:59 +0100 Subject: [PATCH 4/4] Remove unused parameter for the speed of light It is defined in ShipUnits, if needed read from there. --- shipLHC/MuFilterHit.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/shipLHC/MuFilterHit.cxx b/shipLHC/MuFilterHit.cxx index 61109b537a..c74129e0ad 100644 --- a/shipLHC/MuFilterHit.cxx +++ b/shipLHC/MuFilterHit.cxx @@ -9,7 +9,6 @@ #include #include -Double_t speedOfLight = TMath::C() *100./1000000000.0 ; // from m/sec to cm/ns // ----- Default constructor ------------------------------------------- MuFilterHit::MuFilterHit() : SndlhcHit()