Skip to content

Commit 7856c3d

Browse files
preghenellapreghenella
authored
Improvements in TOF start-time determination (#50)
1 parent d10269d commit 7856c3d

File tree

4 files changed

+115
-34
lines changed

4 files changed

+115
-34
lines changed

‎examples/scripts/tof.sh‎

Lines changed: 14 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,18 @@
11
#! /usr/bin/env bash
22

3-
NJOBS=5 # number of max parallel runs
4-
NRUNS=10 # number of runs
5-
NEVENTS=10000 # number of events in a run
3+
NJOBS=5 # number of max parallel runs
4+
NRUNS=10 # number of runs
5+
NEVENTS=10000 # number of events in a run
66

7-
BFIELD=5. # magnetic field [kG]
8-
SIGMAT=0.020 # time resolution [ns]
9-
TAILLX=1.0 # tail on left [q]
10-
TAILRX=1.3 # tail on right [q]
11-
TOFRAD=100. # TOF radius [cm]
12-
TOFLEN=200. # TOF half length [cm]
13-
TOFETA=1.443 # TOF max pseudorapidity
7+
LUTTAG="werner" # LUT tag name
8+
BFIELD=5. # magnetic field [kG]
9+
SIGMAT=0.020 # time resolution [ns]
10+
SIGMA0=0.200 # vertex time spread [ns]
11+
TAILLX=1.0 # tail on left [q]
12+
TAILRX=1.3 # tail on right [q]
13+
TOFRAD=100. # TOF radius [cm]
14+
TOFLEN=200. # TOF half length [cm]
15+
TOFETA=1.443 # TOF max pseudorapidity
1416

1517
### calculate max eta from geometry
1618
TOFETA=`awk -v a=$TOFRAD -v b=$TOFLEN 'BEGIN {th=atan2(a,b)*0.5; sth=sin(th); cth=cos(th); print -log(sth/cth)}'`
@@ -36,10 +38,11 @@ sed -i -e "s/set barrel_TimeResolution .*$/set barrel_TimeResolution ${SIGMAT}e\
3638
sed -i -e "s/set barrel_TailRight .*$/set barrel_TailRight ${TAILRX}/" propagate.tcl
3739
sed -i -e "s/set barrel_TailLeft .*$/set barrel_TailLeft ${TAILLX}/" propagate.tcl
3840
sed -i -e "s/double tof_sigmat = .*$/double tof_sigmat = ${SIGMAT}\;/" tof.C
41+
sed -i -e "s/double tof_sigma0 = .*$/double tof_sigma0 = ${SIGMA0}\;/" tof.C
3942

4043
### create LUTs
4144
BFIELDT=`awk -v a=$BFIELD 'BEGIN {print a*0.1}'`
42-
$DELPHESO2_ROOT/examples/scripts/create_luts.sh werner $BFIELDT $TOFRAD
45+
$DELPHESO2_ROOT/examples/scripts/create_luts.sh $LUTTAG $BFIELDT $TOFRAD
4346

4447
### loop over runs
4548
rm -f .running.* delphes.*.root

‎examples/smearing/tof.C‎

Lines changed: 26 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ R__LOAD_LIBRARY(libDelphesO2)
44
double tof_radius = 100.; // [cm]
55
double tof_length = 200.; // [cm]
66
double tof_sigmat = 0.02; // [ns]
7+
double tof_sigma0 = 0.20; // [ns]
78

89
void
910
tof(const char *inputFile = "delphes.root",
@@ -25,7 +26,7 @@ tof(const char *inputFile = "delphes.root",
2526

2627
// TOF layer
2728
o2::delphes::TOFLayer toflayer;
28-
toflayer.setup(tof_radius, tof_length, tof_sigmat);
29+
toflayer.setup(tof_radius, tof_length, tof_sigmat, tof_sigma0);
2930

3031
// smearer
3132
o2::delphes::TrackSmearer smearer;
@@ -51,12 +52,18 @@ tof(const char *inputFile = "delphes.root",
5152
// histograms
5253
auto hTime0 = new TH1F("hTime0", ";t_{0} (ns)", 1000, -1., 1.);
5354
auto hBetaP = new TH2F("hBetaP", ";#it{p} (GeV/#it{c});#beta", nbins, xbins, 1000, 0.1, 1.1);
54-
TH2 *hNsigmaP[5];
55+
TH2 *hNsigmaPt[5], *hNsigmaPt_true[5][5];
5556
const char *pname[5] = {"el", "mu", "pi", "ka", "pr"};
5657
const char *plabel[5] = {"e", "#mu", "#pi", "K", "p"};
57-
for (int i = 0; i < 5; ++i)
58-
hNsigmaP[i] = new TH2F(Form("hNsigmaP_%s", pname[i]), Form("#it{p} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
59-
58+
for (int i = 0; i < 5; ++i) {
59+
hNsigmaPt[i] = new TH2F(Form("hNsigmaPt_%s", pname[i]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
60+
for (int j = 0; j < 5; ++j) {
61+
hNsigmaPt_true[i][j] = new TH2F(Form("hNsigmaPt_%s_true_%s", pname[i], pname[j]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
62+
}
63+
}
64+
65+
std::map<int, int> pidmap = { {11, 0}, {13, 1}, {211, 2}, {321, 3}, {2212, 4} };
66+
6067
for (Int_t ientry = 0; ientry < numberOfEntries; ++ientry) {
6168

6269
// Load selected branches with data from specified event
@@ -93,25 +100,35 @@ tof(const char *inputFile = "delphes.root",
93100
// loop over tracks and do PID
94101
for (auto track : tof_tracks) {
95102

103+
auto pdg = std::abs(track->PID);
104+
auto ipdg = pidmap[pdg];
105+
96106
// fill beta-p
97107
auto p = track->P;
108+
auto pt = track->PT;
98109
auto beta = toflayer.getBeta(*track);
99110
hBetaP->Fill(p, beta);
100111

101112
// fill nsigma
102113
std::array<float, 5> deltat, nsigma;
103114
toflayer.makePID(*track, deltat, nsigma);
104-
for (int i = 0; i < 5; ++i)
105-
hNsigmaP[i]->Fill(p, nsigma[i]);
115+
for (int i = 0; i < 5; ++i) {
116+
hNsigmaPt[i]->Fill(pt, nsigma[i]);
117+
hNsigmaPt_true[i][ipdg]->Fill(pt, nsigma[i]);
118+
}
106119

107120
}
108121
}
109122

110123
auto fout = TFile::Open(outputFile, "RECREATE");
111124
hTime0->Write();
112125
hBetaP->Write();
113-
for (int i = 0; i < 5; ++i)
114-
hNsigmaP[i]->Write();
126+
for (int i = 0; i < 5; ++i) {
127+
hNsigmaPt[i]->Write();
128+
for (int j = 0; j < 5; ++j) {
129+
hNsigmaPt_true[i][j]->Write();
130+
}
131+
}
115132
fout->Close();
116133

117134
}

‎src/TOFLayer.cc‎

Lines changed: 73 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -12,11 +12,12 @@ namespace delphes
1212
/*****************************************************************/
1313

1414
void
15-
TOFLayer::setup(float radius, float length, float sigmat)
15+
TOFLayer::setup(float radius, float length, float sigmat, float sigma0)
1616
{
1717
mRadius = radius;
1818
mLength = length;
1919
mSigmaT = sigmat;
20+
mSigma0 = sigma0;
2021
}
2122

2223
/*****************************************************************/
@@ -43,7 +44,7 @@ TOFLayer::hasTOF(const Track &track)
4344
float
4445
TOFLayer::getBeta(const Track &track)
4546
{
46-
double tof = track.TOuter * 1.e9 - mTime0; // [ns]
47+
double tof = track.TOuter * 1.e9; // [ns]
4748
double L = track.L * 0.1; // [cm]
4849
double c = 29.9792458; // [cm/ns]
4950
return (L / tof / c);
@@ -57,7 +58,8 @@ TOFLayer::makePID(const Track &track, std::array<float, 5> &deltat, std::array<f
5758
double pmass[5] = {0.00051099891, 0.10565800, 0.13957000, 0.49367700, 0.93827200};
5859

5960
/** get info **/
60-
double tof = track.TOuter * 1.e9 - mTime0; // [ns]
61+
double tof = track.TOuter * 1.e9; // [ns]
62+
double etof = track.ErrorT * 1.e9; // [ns]
6163
double L = track.L * 0.1; // [cm]
6264
double p = track.P;
6365
double p2 = p * p;
@@ -70,7 +72,7 @@ TOFLayer::makePID(const Track &track, std::array<float, 5> &deltat, std::array<f
7072
double mass2 = pmass[ipart] * pmass[ipart];
7173
double texp = Lc / p * TMath::Sqrt(mass2 + p2);
7274
double etexp = Lc * mass2 / p2 / TMath::Sqrt(mass2 + p2) * ep;
73-
double sigma = TMath::Sqrt(etexp * etexp + mSigmaT * mSigmaT + mSigma0 * mSigma0);
75+
double sigma = TMath::Sqrt(etexp * etexp + etof * etof);
7476
deltat[ipart] = tof - texp;
7577
nsigma[ipart] = deltat[ipart] / sigma;
7678
}
@@ -85,13 +87,14 @@ TOFLayer::eventTime(std::vector<Track *> &tracks, std::array<float, 2> &tzero)
8587

8688
double sum = 0.;
8789
double sumw = 0.;
88-
90+
8991
for (auto &track : tracks) {
9092

9193
int pid = track->PID;
9294
double mass = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
9395
double mass2 = mass * mass;
9496
double tof = track->TOuter * 1.e9; // [ns]
97+
double etof = track->ErrorT * 1.e9; // [ns]
9598
double L = track->L * 0.1; // [cm]
9699
double p = track->P; // [GeV/c]
97100
double ep = track->ErrorP;
@@ -100,7 +103,7 @@ TOFLayer::eventTime(std::vector<Track *> &tracks, std::array<float, 2> &tzero)
100103
double Lc = L / c;
101104
double texp = Lc / p * TMath::Sqrt(mass2 + p2);
102105
double etexp = Lc * mass2 / p2 / TMath::Sqrt(mass2 + p2) * ep;
103-
double sigma = TMath::Sqrt(etexp * etexp + mSigmaT * mSigmaT);
106+
double sigma = TMath::Sqrt(etexp * etexp + etof * etof);
104107
double deltat = tof - texp;
105108

106109
double w = 1. / (sigma * sigma);
@@ -110,14 +113,74 @@ TOFLayer::eventTime(std::vector<Track *> &tracks, std::array<float, 2> &tzero)
110113
}
111114

112115
if (sumw <= 0.) {
113-
mTime0 = 0.;
114-
mSigma0 = 1000.;
116+
tzero[0] = 0.;
117+
tzero[1] = mSigma0;
115118
return false;
116119
}
117120

118-
mTime0 = tzero[0] = sum / sumw;
119-
mSigma0 = tzero[1] = sqrt(1. / sumw);
121+
tzero[0] = sum / sumw;
122+
tzero[1] = sqrt(1. / sumw);
123+
124+
// if we have many tracks, we use the start-time computed with all tracks
125+
126+
if (tracks.size() > 4) {
127+
for (auto &track : tracks) {
128+
track->TOuter -= tzero[0] * 1.e-9; // [s]
129+
track->ErrorT = std::hypot(track->ErrorT, tzero[1] * 1.e-9);
130+
}
131+
return true;
132+
}
133+
134+
// if we have few tracks, we do the combinations excluding the track of interest
135+
136+
std::vector<double> time0, sigma0;
137+
time0.reserve(tracks.size());
138+
sigma0.reserve(tracks.size());
139+
for (int itrack = 0; itrack < tracks.size(); ++itrack) {
140+
141+
time0[itrack] = 0;
142+
sigma0[itrack] = mSigma0;
143+
144+
sum = 0.;
145+
sumw = 0.;
146+
147+
for (int jtrack = 0; jtrack < tracks.size(); ++jtrack) {
148+
if (itrack == jtrack) continue; // do not use self
149+
150+
auto &track = tracks[jtrack];
151+
int pid = track->PID;
152+
double mass = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
153+
double mass2 = mass * mass;
154+
double tof = track->TOuter * 1.e9; // [ns]
155+
double etof = track->ErrorT * 1.e9; // [ns]
156+
double L = track->L * 0.1; // [cm]
157+
double p = track->P; // [GeV/c]
158+
double ep = track->ErrorP;
159+
double p2 = p * p;
160+
double c = 29.9792458; // [cm/ns]
161+
double Lc = L / c;
162+
double texp = Lc / p * TMath::Sqrt(mass2 + p2);
163+
double etexp = Lc * mass2 / p2 / TMath::Sqrt(mass2 + p2) * ep;
164+
double sigma = TMath::Sqrt(etexp * etexp + etof * etof);
165+
double deltat = tof - texp;
166+
double w = 1. / (sigma * sigma);
167+
sum += w * deltat;
168+
sumw += w;
169+
}
170+
171+
if (sumw <= 0.) continue;
172+
173+
time0[itrack] = sum / sumw;
174+
sigma0[itrack] = std::sqrt(1. / sumw);
120175

176+
}
177+
178+
for (int itrack = 0; itrack < tracks.size(); ++itrack) {
179+
auto &track = tracks[itrack];
180+
track->TOuter -= time0[itrack] * 1.e-9; // [s]
181+
track->ErrorT = std::hypot(track->ErrorT, sigma0[itrack] * 1.e-9);
182+
}
183+
121184
return true;
122185
}
123186

‎src/TOFLayer.hh‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ public:
1919

2020
enum { kBarrel, kForward }; // type of TOF detector
2121

22-
void setup(float radius, float length, float sigmat);
22+
void setup(float radius, float length, float sigmat, float sigma0);
2323
bool hasTOF(const Track &track);
2424
float getBeta(const Track &track);
2525
void makePID(const Track &track, std::array<float, 5> &deltat, std::array<float, 5> &nsigma);
@@ -35,9 +35,7 @@ protected:
3535
float mRadiusIn = 10.; // [cm]
3636
float mLength = 200.; // [cm]
3737
float mSigmaT = 0.02; // [ns]
38-
39-
float mTime0 = 0.; // [ns]
40-
float mSigma0 = 0.; // [ns]
38+
float mSigma0 = 0.200; // [ns]
4139

4240
};
4341

0 commit comments

Comments
 (0)