Skip to content

Veto and US time alignment #214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
siilieva wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Veto and US time alignment #214

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5271002
Add Veto and US time alignment
Mar 13, 2024
260d4da
Add methods to apply the Veto and US time corrections
Mar 13, 2024
a545b24
Read the constants for the US and Veto time alignment
Mar 13, 2024
f615edf
Remove unused parameter for the speed of light
Mar 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 19 additions & 7 deletions shipLHC/MuFilter.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<float> 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();
Expand Down Expand Up @@ -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;
Expand Down
6 changes: 5 additions & 1 deletion shipLHC/MuFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,17 +35,20 @@ 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);

void InitEvent(SNDLHCEventHeader *e);
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<float> 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<Float_t> GetConfParVector(TString name){return conf_vectors[name];}

/** Initialization of the detector is done here */
virtual void Initialize();
Expand Down Expand Up @@ -101,6 +104,7 @@ class MuFilter : public FairDetector
std::map<TString,Float_t> conf_floats;
std::map<TString,Int_t> conf_ints;
std::map<TString,TString> conf_strings;
std::map<TString,std::vector<Float_t>> conf_vectors;
SNDLHCEventHeader *eventHeader;

// Vector to store runs covered in the geometry file.
Expand Down
48 changes: 38 additions & 10 deletions shipLHC/MuFilterHit.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "MuFilterHit.h"
#include "MuFilter.h"
#include "ShipUnit.h"
#include "TROOT.h"
#include "FairRunSim.h"
#include "TGeoNavigator.h"
Expand All @@ -8,7 +9,6 @@
#include <TRandom.h>
#include <iomanip>

Double_t speedOfLight = TMath::C() *100./1000000000.0 ; // from m/sec to cm/ns
// ----- Default constructor -------------------------------------------
MuFilterHit::MuFilterHit()
: SndlhcHit()
Expand Down Expand Up @@ -167,15 +167,16 @@ std::map<Int_t,Float_t> MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive)
}

// ----- Public method Get List of time measurements -------------------------------------------
std::map<Int_t,Float_t> MuFilterHit::GetAllTimes(Bool_t mask)
std::map<Int_t,Float_t> MuFilterHit::GetAllTimes(Bool_t mask, Bool_t apply_t_corr, Double_t SipmDistance)
{
OptForTimeCorrections(mask, apply_t_corr, SipmDistance);
std::map<Int_t,Float_t> allTimes;
for (unsigned int s=0; s<nSides; ++s){
for (unsigned int j=0; j<nSiPMs; ++j){
unsigned int channel = j+s*nSiPMs;
if (signals[channel]> 0){
if (!fMasked[channel] || !mask){
allTimes[channel] = times[channel];
allTimes[channel] = fTimesHelper[channel];
}
}
}
Expand All @@ -184,9 +185,10 @@ std::map<Int_t,Float_t> 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.;
Expand All @@ -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;
}
}
Expand All @@ -206,17 +208,18 @@ 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<nSides; ++s){
for (unsigned int j=0; j<nSiPMs; ++j){
unsigned int channel = j+s*nSiPMs;
if (signals[channel]> 0){
if (!fMasked[channel] || !mask){
if (times[channel]<first[s]) {first[s] = times[channel];}
if (times[channel]<first[s]) {first[s] = fTimesHelper[channel];}
}
}
}
Expand All @@ -229,8 +232,9 @@ Float_t MuFilterHit::GetFastDeltaT(Bool_t mask)


// ----- Public method Get mean time -----------------
Float_t MuFilterHit::GetImpactT(Bool_t mask)
Float_t MuFilterHit::GetImpactT(Bool_t mask, Bool_t apply_t_corr, Double_t SipmDistance)
{
OptForTimeCorrections(mask, apply_t_corr, SipmDistance);
Float_t mean[] = {0,0};
Int_t count[] = {0,0};
Float_t dT = -999.;
Expand All @@ -248,14 +252,14 @@ Float_t MuFilterHit::GetImpactT(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;
}
}
}
}
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;
}
Expand Down Expand Up @@ -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<MuFilter*> (gROOT->GetListOfGlobals()->FindObject("MuFilter"));
for (unsigned int s=0; s<nSides; ++s){
for (unsigned int j=0; j<nSiPMs; ++j){
unsigned int channel = j+s*nSiPMs;
if (signals[channel]> 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)
Expand Down
19 changes: 14 additions & 5 deletions shipLHC/MuFilterHit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,12 @@ class MuFilterHit : public SndlhcHit
Float_t SumOfSignals(char* opt,Bool_t mask=kTRUE);
std::map<TString,Float_t> SumOfSignals(Bool_t mask=kTRUE);
std::map<Int_t,Float_t> GetAllSignals(Bool_t mask=kTRUE,Bool_t positive=kTRUE);
std::map<Int_t,Float_t> 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<Int_t,Float_t> 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;}
Expand All @@ -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);


};
Expand Down
46 changes: 45 additions & 1 deletion shipLHC/modifyGeoFileDict.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#!/usr/bin/env python
import json
import ROOT,os
from rootpyPickler import Pickler
from rootpyPickler import Unpickler
Expand All @@ -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
Expand Down