Skip to content

Commit 71fed8a

Browse files
jesgumjesgum
authored
[ALICE3] Add hit efficiency correction to the fasttracker (#10465)
1 parent 8eb14d0 commit 71fed8a

File tree

3 files changed

+238
-54
lines changed

3 files changed

+238
-54
lines changed

‎ALICE3/Core/FastTracker.cxx‎

Lines changed: 174 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,10 @@ FastTracker::FastTracker()
2828
// base constructor
2929
magneticField = 20; // in kiloGauss
3030
applyZacceptance = false;
31+
applyMSCorrection = true;
32+
applyMSCorrection = true;
33+
applyElossCorrection = true;
34+
applyEffCorrection = true;
3135
covMatFactor = 0.99f;
3236
verboseLevel = 0;
3337

@@ -37,6 +41,17 @@ FastTracker::FastTracker()
3741
nIntercepts = 0;
3842
nSiliconPoints = 0;
3943
nGasPoints = 0;
44+
45+
maxRadiusSlowDet = 10;
46+
integrationTime = 0.02; // ms
47+
crossSectionMinB = 8;
48+
dNdEtaCent = 2200;
49+
dNdEtaMinB = 1;
50+
avgRapidity = 0.45;
51+
sigmaD = 6.0;
52+
luminosity = 1.e27;
53+
otherBackground = 0.0; // [0, 1]
54+
upcBackgroundMultiplier = 1.0;
4055
}
4156

4257
void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
@@ -71,61 +86,61 @@ void FastTracker::Print()
7186
LOG(info) << "+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+";
7287
}
7388

74-
void FastTracker::AddSiliconALICE3v4()
89+
void FastTracker::AddSiliconALICE3v4(std::vector<float> pixelResolution)
7590
{
76-
LOG(info) << " Adding ALICE 3 v4 ITS layers";
91+
LOG(info) << " ALICE 3: Adding v4 tracking layers";
7792
float x0IT = 0.001; // 0.1%
7893
float x0OT = 0.005; // 0.5%
7994
float xrhoIB = 1.1646e-02; // 50 mum Si
80-
float xrhoOB = 1.1646e-01; // 500 mum Si
81-
82-
float resRPhiIT = 0.00025; // 2.5 mum
83-
float resZIT = 0.00025; // 2.5 mum
84-
float resRPhiOT = 0.0005; // 5 mum
85-
float resZOT = 0.0005; // 5 mum
95+
float xrhoOT = 1.1646e-01; // 500 mum Si
8696
float eff = 1.00;
8797

98+
float resRPhiIT = pixelResolution[0];
99+
float resZIT = pixelResolution[1];
100+
float resRPhiOT = pixelResolution[2];
101+
float resZOT = pixelResolution[3];
102+
88103
layers.push_back(DetLayer{"bpipe0", 0.48, 250, 0.00042, 2.772e-02, 0.0f, 0.0f, 0.0f, 0}); // 150 mum Be
89104
layers.push_back(DetLayer{"ddd0", 0.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
90105
layers.push_back(DetLayer{"ddd1", 1.2, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
91106
layers.push_back(DetLayer{"ddd2", 2.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
92107
layers.push_back(DetLayer{"bpipe1", 5.7, 250, 0.0014, 9.24e-02, 0.0f, 0.0f, 0.0f, 0}); // 500 mum Be
93-
layers.push_back(DetLayer{"ddd3", 7., 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
94-
layers.push_back(DetLayer{"ddd4", 10., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
95-
layers.push_back(DetLayer{"ddd5", 13., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
96-
layers.push_back(DetLayer{"ddd6", 16., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
97-
layers.push_back(DetLayer{"ddd7", 25., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
98-
layers.push_back(DetLayer{"ddd8", 40., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
99-
layers.push_back(DetLayer{"ddd9", 45., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
108+
layers.push_back(DetLayer{"ddd3", 7., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
109+
layers.push_back(DetLayer{"ddd4", 10., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
110+
layers.push_back(DetLayer{"ddd5", 13., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
111+
layers.push_back(DetLayer{"ddd6", 16., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
112+
layers.push_back(DetLayer{"ddd7", 25., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
113+
layers.push_back(DetLayer{"ddd8", 40., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
114+
layers.push_back(DetLayer{"ddd9", 45., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
100115
}
101116

102-
void FastTracker::AddSiliconALICE3v1()
117+
void FastTracker::AddSiliconALICE3v2(std::vector<float> pixelResolution)
103118
{
104-
LOG(info) << " Adding ALICE 3 v1 ITS layers";
119+
LOG(info) << "ALICE 3: Adding v2 tracking layers;";
105120
float x0IT = 0.001; // 0.1%
106-
float x0OT = 0.005; // 0.5%
121+
float x0OT = 0.01; // 1.0%
107122
float xrhoIB = 2.3292e-02; // 100 mum Si
108-
float xrhoOB = 2.3292e-01; // 1000 mum Si
109-
110-
float resRPhiIT = 0.00025; // 2.5 mum
111-
float resZIT = 0.00025; // 2.5 mum
112-
float resRPhiOT = 0.00100; // 5 mum
113-
float resZOT = 0.00100; // 5 mum
123+
float xrhoOT = 2.3292e-01; // 1000 mum Si
114124
float eff = 1.00;
115125

126+
float resRPhiIT = pixelResolution[0];
127+
float resZIT = pixelResolution[1];
128+
float resRPhiOT = pixelResolution[2];
129+
float resZOT = pixelResolution[3];
130+
116131
layers.push_back(DetLayer{"bpipe0", 0.48, 250, 0.00042, 2.772e-02, 0.0f, 0.0f, 0.0f, 0}); // 150 mum Be
117132
layers.push_back(DetLayer{"B00", 0.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
118133
layers.push_back(DetLayer{"B01", 1.2, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
119134
layers.push_back(DetLayer{"B02", 2.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
120135
layers.push_back(DetLayer{"bpipe1", 3.7, 250, 0.0014, 9.24e-02, 0.0f, 0.0f, 0.0f, 0}); // 500 mum Be
121-
layers.push_back(DetLayer{"B03", 3.75, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff, 1});
122-
layers.push_back(DetLayer{"B04", 7., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
123-
layers.push_back(DetLayer{"B05", 12., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
124-
layers.push_back(DetLayer{"B06", 20., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
125-
layers.push_back(DetLayer{"B07", 30., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
126-
layers.push_back(DetLayer{"B08", 45., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
127-
layers.push_back(DetLayer{"B09", 60., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
128-
layers.push_back(DetLayer{"B10", 80., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff, 1});
136+
layers.push_back(DetLayer{"B03", 3.75, 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
137+
layers.push_back(DetLayer{"B04", 7., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
138+
layers.push_back(DetLayer{"B05", 12., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
139+
layers.push_back(DetLayer{"B06", 20., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
140+
layers.push_back(DetLayer{"B07", 30., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
141+
layers.push_back(DetLayer{"B08", 45., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
142+
layers.push_back(DetLayer{"B09", 60., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
143+
layers.push_back(DetLayer{"B10", 80., 250, x0OT, xrhoOT, resRPhiOT, resZOT, eff, 1});
129144
}
130145

131146
void FastTracker::AddTPC(float phiResMean, float zResMean)
@@ -134,7 +149,7 @@ void FastTracker::AddTPC(float phiResMean, float zResMean)
134149

135150
// porting of DetectorK::AddTPC
136151
// see here:
137-
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx
152+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L522
138153
// % Radiation Lengths ... Average per TPC row (i.e. total/159 )
139154
const int kNPassiveBound = 2;
140155
const float radLBoundary[kNPassiveBound] = {1.692612e-01, 8.711904e-02};
@@ -173,20 +188,118 @@ void FastTracker::AddTPC(float phiResMean, float zResMean)
173188
}
174189
}
175190

191+
float FastTracker::Dist(float z, float r)
192+
{
193+
// porting of DetektorK::Dist
194+
// see here:
195+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L743
196+
int index = 1;
197+
int nSteps = 301;
198+
double dist = 0.0;
199+
double dz0 = (4 * sigmaD - (-4) * sigmaD / (nSteps = 1));
200+
double z0 = 0.0;
201+
for (int i = 0; i < nSteps; i++) {
202+
if (i == nSteps - 1)
203+
index = 1;
204+
z0 = -4 * sigmaD + i * dz0;
205+
dist += index * (dz0 / 3.) * (1 / o2::math_utils::sqrt(o2::constants::math::TwoPI) / sigmaD) * exp(-z0 * z0 / 2. / sigmaD / sigmaD) * (1 / o2::math_utils::sqrt((z - z0) * (z - z0) + r * r));
206+
if (index != 4)
207+
index = 4;
208+
else
209+
index = 2;
210+
}
211+
return dist;
212+
}
213+
214+
float FastTracker::OneEventHitDensity(float multiplicity, float radius)
215+
{
216+
// porting of DetektorK::OneEventHitDensity
217+
// see here:
218+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L694
219+
float den = multiplicity / (o2::constants::math::TwoPI * radius * radius);
220+
float tg = o2::math_utils::tan(2. * o2::math_utils::atan(-avgRapidity));
221+
den = den / o2::math_utils::sqrt(1 + 1 / (tg * tg));
222+
return den;
223+
}
224+
225+
float FastTracker::IntegratedHitDensity(float multiplicity, float radius)
226+
{
227+
// porting of DetektorK::IntegratedHitDensity
228+
// see here:
229+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L712
230+
float zdcHz = luminosity * 1.e24 * crossSectionMinB;
231+
float den = zdcHz * integrationTime / 1000. * multiplicity * Dist(0., radius) / (o2::constants::math::TwoPI * radius);
232+
if (den < OneEventHitDensity(multiplicity, radius))
233+
den = OneEventHitDensity(multiplicity, radius);
234+
return den;
235+
}
236+
237+
float FastTracker::UpcHitDensity(float radius)
238+
{
239+
// porting of DetektorK::UpcHitDensity
240+
// see here:
241+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L727
242+
float mUPCelectrons = 0;
243+
mUPCelectrons = lhcUPCScale * 5456 / (radius * radius) / dNdEtaMinB;
244+
if (mUPCelectrons < 0)
245+
mUPCelectrons = 0.0;
246+
mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB, radius);
247+
mUPCelectrons *= upcBackgroundMultiplier;
248+
return mUPCelectrons;
249+
}
250+
251+
float FastTracker::HitDensity(float radius)
252+
{
253+
// porting of DetektorK::HitDensity
254+
// see here:
255+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L663
256+
float arealDensity = 0.;
257+
if (radius > maxRadiusSlowDet) {
258+
arealDensity = OneEventHitDensity(dNdEtaCent, radius);
259+
arealDensity += otherBackground * OneEventHitDensity(dNdEtaMinB, radius);
260+
}
261+
262+
// In the version of Delphes used to produce
263+
// Look-up tables, UpcHitDensity(radius) always returns 0,
264+
// hence it is left commented out for now
265+
if (radius < maxRadiusSlowDet) {
266+
arealDensity = OneEventHitDensity(dNdEtaCent, radius);
267+
arealDensity += otherBackground * OneEventHitDensity(dNdEtaMinB, radius) + IntegratedHitDensity(dNdEtaMinB, radius);
268+
// +UpcHitDensity(radius);
269+
}
270+
return arealDensity;
271+
}
272+
273+
float FastTracker::ProbGoodChiSqHit(float radius, float searchRadiusRPhi, float searchRadiusZ)
274+
{
275+
// porting of DetektorK::ProbGoodChiSqHit
276+
// see here:
277+
// https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L629
278+
float sx, goodHit;
279+
sx = o2::constants::math::TwoPI * searchRadiusRPhi * searchRadiusZ * HitDensity(radius);
280+
goodHit = 1. / (1 + sx);
281+
return goodHit;
282+
}
283+
176284
// function to provide a reconstructed track from a perfect input track
177285
// returns number of intercepts (generic for now)
178-
int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack)
286+
int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, float nch)
179287
{
180288
hits.clear();
181289
nIntercepts = 0;
182290
nSiliconPoints = 0;
183291
nGasPoints = 0;
184292
std::array<float, 3> posIni; // provision for != PV
185293
inputTrack.getXYZGlo(posIni);
186-
float initialRadius = std::hypot(posIni[0], posIni[1]);
187-
float kTrackingMargin = 0.1;
188-
int xrhosteps = 100;
189-
bool applyAngularCorrection = true;
294+
const float initialRadius = std::hypot(posIni[0], posIni[1]);
295+
const float kTrackingMargin = 0.1;
296+
const int kMaxNumberOfDetectors = 20;
297+
const int xrhosteps = 100;
298+
const bool applyAngularCorrection = true;
299+
300+
for (int i = 0; i < kMaxNumberOfDetectors; ++i)
301+
goodHitProbability.push_back(-1.);
302+
goodHitProbability[0] = 1.;
190303

191304
// +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
192305
// Outward pass to find intercepts
@@ -351,6 +464,13 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
351464
nGasPoints++; // count TPC/gas hits
352465

353466
hits.push_back(thisHit);
467+
468+
if (applyEffCorrection && layers[il].type != 0) { // good hit probability calculation
469+
double sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].resRPhi * layers[il].resRPhi);
470+
double sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].resZ * layers[il].resZ);
471+
goodHitProbability[il] = ProbGoodChiSqHit(layers[il].r * 100, sigYCmb * 100, sigZCmb * 100);
472+
goodHitProbability[0] *= goodHitProbability[il];
473+
}
354474
}
355475

356476
// backpropagate to original radius
@@ -367,6 +487,22 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
367487
if (nIntercepts < 4)
368488
return nIntercepts;
369489

490+
// generate efficiency
491+
if (applyEffCorrection) {
492+
dNdEtaCent = nch;
493+
float eff = 1.;
494+
for (int i = 0; i < kMaxNumberOfDetectors; i++) {
495+
float iGoodHit = goodHitProbability[i];
496+
if (iGoodHit <= 0)
497+
continue;
498+
499+
eff *= iGoodHit;
500+
}
501+
502+
if (gRandom->Uniform() > eff)
503+
return -8;
504+
}
505+
370506
outputTrack.setCov(inwardTrack.getCov());
371507
outputTrack.checkCovariance();
372508

‎ALICE3/Core/FastTracker.h‎

Lines changed: 27 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -37,24 +37,44 @@ class FastTracker
3737
void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
3838
DetLayer GetLayer(const int layer, bool ignoreBarrelLayers = true);
3939

40-
void AddSiliconALICE3v4();
41-
void AddSiliconALICE3v1();
40+
void AddSiliconALICE3v4(std::vector<float> pixelResolution);
41+
void AddSiliconALICE3v2(std::vector<float> pixelResolution);
4242
void AddTPC(float phiResMean, float zResMean);
4343

4444
void Print();
45-
int FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack);
45+
int FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, float nch);
46+
47+
// For efficiency calculation
48+
float Dist(float z, float radius);
49+
float OneEventHitDensity(float multiplicity, float radius);
50+
float IntegratedHitDensity(float multiplicity, float radius);
51+
float UpcHitDensity(float radius);
52+
float HitDensity(float radius);
53+
float ProbGoodChiSqHit(float radius, float searchRadiusRPhi, float searchRadiusZ);
4654

4755
// Definition of detector layers
4856
std::vector<DetLayer> layers;
4957
std::vector<std::vector<float>> hits; // bookkeep last added hits
5058

5159
// operational
52-
float magneticField; // in kiloGauss (5 = 0.5T, etc)
5360
bool applyZacceptance; // check z acceptance or not
54-
float covMatFactor; // covmat off-diagonal factor to use for covmat fix (negative: no factor)
55-
int verboseLevel; // 0: not verbose, >0 more verbose
5661
bool applyMSCorrection; // Apply correction for multiple scattering
5762
bool applyElossCorrection; // Apply correction for eloss (requires MS correction)
63+
bool applyEffCorrection; // Apply correction for hit efficiency
64+
int verboseLevel; // 0: not verbose, >0 more verbose
65+
int crossSectionMinB;
66+
int dNdEtaCent;
67+
int dNdEtaMinB;
68+
float integrationTime;
69+
float magneticField; // in kiloGauss (5 = 0.5T, etc)
70+
float covMatFactor; // covmat off-diagonal factor to use for covmat fix (negative: no factor)
71+
float sigmaD;
72+
float luminosity;
73+
float otherBackground;
74+
float maxRadiusSlowDet;
75+
float avgRapidity;
76+
float lhcUPCScale;
77+
float upcBackgroundMultiplier;
5878

5979
uint64_t covMatOK; // cov mat has negative eigenvals
6080
uint64_t covMatNotOK; // cov mat has negative eigenvals
@@ -63,6 +83,7 @@ class FastTracker
6383
int nIntercepts; // found in first outward propagation
6484
int nSiliconPoints; // silicon-based space points added to track
6585
int nGasPoints; // tpc-based space points added to track
86+
std::vector<float> goodHitProbability;
6687

6788
ClassDef(FastTracker, 1);
6889
};

0 commit comments

Comments
 (0)