Skip to content

Commit d03da77

Browse files
njacazionjacazio
committed
Optimize the TOF inactive area
1 parent 5dec98a commit d03da77

File tree

2 files changed

+97
-106
lines changed

2 files changed

+97
-106
lines changed

‎ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx‎

Lines changed: 95 additions & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -95,9 +95,9 @@ struct OnTheFlyTofPid {
9595
Configurable<float> outerTOFRadius{"outerTOFRadius", 80, "barrel outer TOF radius (cm)"};
9696
Configurable<float> innerTOFLength{"innerTOFLength", 124, "barrel inner TOF length (cm)"};
9797
Configurable<float> outerTOFLength{"outerTOFLength", 250, "barrel outer TOF length (cm)"};
98-
Configurable<float> innerTOFPixelArea{"innerTOFPixelArea", 1.0, "barrel inner TOF pixel area (mm^2)"};
98+
Configurable<std::array<float, 2>> innerTOFPixelDimension{"innerTOFPixelDimension", {0.1, 0.1}, "barrel inner TOF pixel dimension in Z and RPhi (cm)"};
9999
Configurable<float> innerTOFFractionOfInactiveArea{"innerTOFFractionOfInactiveArea", 0.f, "barrel inner TOF fraction of inactive area"};
100-
Configurable<float> outerTOFPixelArea{"outerTOFPixelArea", 25.0, "barrel outer TOF pixel area (mm^2)"};
100+
Configurable<std::array<float, 2>> outerTOFPixelDimension{"outerTOFPixelDimension", {0.1, 0.1}, "barrel outer TOF pixel dimension in Z and RPhi (cm)"};
101101
Configurable<float> outerTOFFractionOfInactiveArea{"outerTOFFractionOfInactiveArea", 0.f, "barrel outer TOF fraction of inactive area"};
102102
Configurable<float> innerTOFTimeReso{"innerTOFTimeReso", 20, "barrel inner TOF time error (ps)"};
103103
Configurable<float> outerTOFTimeReso{"outerTOFTimeReso", 20, "barrel outer TOF time error (ps)"};
@@ -189,8 +189,6 @@ struct OnTheFlyTofPid {
189189
histos.add("h2dEventTimeres", "h2dEventTimeres", kTH2F, {axisdNdeta, {300, 0, 300, "resolution"}});
190190
listEfficiency.setObject(new THashList);
191191
listEfficiency->Add(new TEfficiency("effEventTime", "effEventTime", plotsConfig.nBinsMult, 0.0f, plotsConfig.maxMultRange));
192-
listEfficiency->Add(new TEfficiency("MapInnerTOF", "MapInnerTOF", 100, -simConfig.innerTOFLength / 2, simConfig.innerTOFLength / 2, 100, 0, simConfig.innerTOFRadius * M_PI_2));
193-
listEfficiency->Add(new TEfficiency("MapOuterTOF", "MapOuterTOF", 100, -simConfig.outerTOFLength / 2, simConfig.outerTOFLength / 2, 100, 0, simConfig.outerTOFRadius * M_PI_2));
194192

195193
const AxisSpec axisMomentum{static_cast<int>(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} (GeV/#it{c})"};
196194
const AxisSpec axisRigidity{static_cast<int>(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} / |#it{z}| (GeV/#it{c})"};
@@ -204,12 +202,14 @@ struct OnTheFlyTofPid {
204202
histos.add("iTOF/h2dTrackLengthInnerVsPt", "h2dTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner});
205203
histos.add("iTOF/h2dTrackLengthInnerRecoVsPt", "h2dTrackLengthInnerRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner});
206204
histos.add("iTOF/h2dDeltaTrackLengthInnerVsPt", "h2dDeltaTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength});
205+
histos.add("iTOF/h2HitMap", "h2HitMap", kTH2F, {{100, -simConfig.innerTOFLength / 2, simConfig.innerTOFLength / 2}, {100, 0, simConfig.innerTOFRadius * 2 * M_PI}});
207206

208207
histos.add("oTOF/h2dVelocityVsMomentumOuter", "h2dVelocityVsMomentumOuter", kTH2F, {axisMomentum, axisVelocity});
209208
histos.add("oTOF/h2dVelocityVsRigidityOuter", "h2dVelocityVsRigidityOuter", kTH2F, {axisRigidity, axisVelocity});
210209
histos.add("oTOF/h2dTrackLengthOuterVsPt", "h2dTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter});
211210
histos.add("oTOF/h2dTrackLengthOuterRecoVsPt", "h2dTrackLengthOuterRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter});
212211
histos.add("oTOF/h2dDeltaTrackLengthOuterVsPt", "h2dDeltaTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength});
212+
histos.add("oTOF/h2HitMap", "h2HitMap", kTH2F, {{100, -simConfig.outerTOFLength / 2, simConfig.outerTOFLength / 2}, {100, 0, simConfig.outerTOFRadius * 2 * M_PI}});
213213

214214
const AxisSpec axisPt{static_cast<int>(plotsConfig.nBinsP), 0.0f, +4.0f, "#it{p}_{T} (GeV/#it{c})"};
215215
const AxisSpec axisEta{static_cast<int>(plotsConfig.nBinsEta), -2.0f, +2.0f, "#eta"};
@@ -256,33 +256,39 @@ struct OnTheFlyTofPid {
256256
private:
257257
const float layerRadius;
258258
const float layerLength;
259-
const float pixelArea;
259+
const float pixelDimensionZ;
260+
const float pixelDimensionRPhi;
260261
const float fractionInactive;
261262
const float magField;
262263

263-
float pixelSize;
264-
float circumference;
265-
float zLength;
266-
int nPixelsPhi;
267-
int nPixelsZ;
268-
float inactiveBorderPhi;
269-
float inactiveBorderZ;
264+
// TH2F* hHitMapInPixel = nullptr;
265+
// TH2F* hHitMap = nullptr;
266+
TAxis* axisZ = nullptr;
267+
TAxis* axisRPhi = nullptr;
268+
TAxis* axisInPixelZ = nullptr;
269+
TAxis* axisInPixelRPhi = nullptr;
270270

271271
public:
272-
TOFLayerEfficiency(float r, float l, float pA, float fIA, float m)
273-
: layerRadius(r), layerLength(l), pixelArea(pA), fractionInactive(fIA), magField(m)
272+
TOFLayerEfficiency(float r, float l, std::array<float, 2> pDimensions, float fIA, float m)
273+
: layerRadius(r), layerLength(l), pixelDimensionZ(pDimensions[0]), pixelDimensionRPhi(pDimensions[1]), fractionInactive(fIA), magField(m)
274274
{
275275
// Assuming square pixels for simplicity
276-
pixelSize = std::sqrt(pA);
277-
circumference = 2.0f * M_PI * layerRadius * 10.0f; // convert cm to mm
278-
zLength = layerLength * 10.0f; // convert cm to mm
279-
nPixelsPhi = static_cast<int>(circumference / pixelSize);
280-
nPixelsZ = static_cast<int>(zLength / pixelSize);
281-
inactiveBorderPhi = pixelSize * std::sqrt(fractionInactive);
282-
inactiveBorderZ = pixelSize * std::sqrt(fractionInactive);
276+
const float circumference = 2.0f * M_PI * layerRadius;
277+
axisZ = new TAxis(static_cast<int>(layerLength / pixelDimensionZ), -layerLength / 2, layerLength);
278+
axisRPhi = new TAxis(static_cast<int>(circumference / pixelDimensionRPhi), 0.f, circumference);
279+
280+
const float inactiveBorderRPhi = pixelDimensionRPhi * std::sqrt(fractionInactive);
281+
const float inactiveBorderZ = pixelDimensionZ * std::sqrt(fractionInactive);
282+
const double arrayRPhi[4] = {0.f, inactiveBorderRPhi, pixelDimensionRPhi - inactiveBorderRPhi, pixelDimensionRPhi};
283+
axisInPixelRPhi = new TAxis(3, arrayRPhi);
284+
const double arrayZ[4] = {0.f, inactiveBorderZ, pixelDimensionZ - inactiveBorderZ, pixelDimensionZ};
285+
axisInPixelZ = new TAxis(3, arrayZ);
286+
287+
// hHitMap = new TH2F(Form("hHitMap_R%.0f", layerRadius), "HitMap;z (cm); r#phi (cm)", 1000, -1000, 1000, 1000, -1000, 1000);
288+
// hHitMapInPixel = new TH2F(Form("hHitMapInPixel_R%.0f", layerRadius), "HitMapInPixel;z (cm); r#phi (cm)", 1000, -10, 10, 1000, -10, 10);
283289
}
284290

285-
bool isInTOFActiveArea(std::array<float, 3> hitPosition = {0.0f, 0.0f, 0.0f})
291+
bool isInTOFActiveArea(std::array<float, 3> hitPosition)
286292
{
287293
if (fractionInactive <= 0.0f) {
288294
return true;
@@ -293,113 +299,92 @@ struct OnTheFlyTofPid {
293299

294300
// Convert 3D position to cylindrical coordinates for area calculation
295301
const float phi = std::atan2(hitPosition[1], hitPosition[0]);
302+
const float rphi = phi * layerRadius;
296303
const float z = hitPosition[2];
297304
const float r = std::sqrt(hitPosition[0] * hitPosition[0] + hitPosition[1] * hitPosition[1]);
298305

299306
// Check if hit is within layer geometric acceptance
300-
if (std::abs(z) > layerLength / 2.0f || r < layerRadius - 1.0f || r > layerRadius + 1.0f) {
307+
if (std::abs(layerRadius - r) > 10.f) {
308+
LOG(warning) << "Hit out of TOF layer acceptance: r=" << r << " cm with respect to the layer radius " << layerRadius;
309+
return false;
310+
}
311+
if (std::abs(z) > layerLength / 2.0f) {
312+
LOG(warning) << "Hit out of TOF layer acceptance: z=" << z << " cm with respect to the layer length " << layerLength;
301313
return false;
302314
}
303315

304-
// Convert hit position to pixel coordinates
305-
const float phiNorm = (phi + M_PI) / (2.0f * M_PI); // normalize phi to [0,1]
306-
const float zNorm = (z * 10.0f + zLength / 2.0f) / zLength; // normalize z to [0,1]
316+
const int pixelIndexPhi = axisRPhi->FindBin(rphi);
317+
const int pixelIndexZ = axisZ->FindBin(z);
307318

308-
const int pixelPhi = static_cast<int>(phiNorm * nPixelsPhi) % nPixelsPhi;
309-
const int pixelZ = static_cast<int>(zNorm * nPixelsZ);
319+
// LOG(info) << "Hit pixel " << pixelIndexPhi << "/" << nPixelsRPhi << " and " << pixelIndexZ << "/" << nPixelsZ;
310320

311-
// Check if pixel is valid
312-
if (pixelZ < 0 || pixelZ >= nPixelsZ) {
313-
LOG(fatal) << "Pixel Z index out of bounds: " << pixelZ << " (nPixelsZ: " << nPixelsZ << ")";
314-
}
315-
if (pixelPhi < 0 || pixelPhi >= nPixelsPhi) {
316-
LOG(fatal) << "Pixel Phi index out of bounds: " << pixelPhi << " (nPixelsPhi: " << nPixelsPhi << ")";
317-
}
318-
// Determine if the pixel is active based on the fraction of inactive area
319-
const float phiPosInPixel = std::abs(phiNorm - pixelPhi);
320-
const float zPosInPixel = std::abs(zNorm - pixelZ);
321-
if (phiPosInPixel > inactiveBorderPhi) {
322-
return false;
323-
}
324-
if (zPosInPixel > inactiveBorderZ) {
321+
if (pixelIndexPhi <= 0 || pixelIndexPhi > axisRPhi->GetNbins() || pixelIndexZ <= 0 || pixelIndexZ > axisZ->GetNbins()) {
322+
// LOG(info) << "Hit out of TOF layer pixel range: pixelIndexPhi=" << pixelIndexPhi << ", pixelIndexZ=" << pixelIndexZ;
325323
return false;
326324
}
325+
// Calculate local position within the pixel
326+
const float localRPhi = (rphi - axisRPhi->GetBinCenter(pixelIndexPhi));
327+
const float localZ = (z - axisZ->GetBinCenter(pixelIndexZ));
328+
329+
// The difference between the hit and the pixel position cannot be greater than the size of the pixel
330+
if (std::abs(localRPhi - axisRPhi->GetBinCenter(pixelIndexPhi)) > axisRPhi->GetBinWidth(pixelIndexPhi)) {
331+
// LOG(warning) << "Local hit difference in phi is bigger than the pixel size";
332+
}
333+
if (std::abs(localZ - axisZ->GetBinCenter(pixelIndexZ)) > axisZ->GetBinWidth(pixelIndexZ)) {
334+
// LOG(warning) << "Local hit difference in z is bigger than the pixel size";
335+
}
336+
switch (axisInPixelRPhi->FindBin(localRPhi)) {
337+
case 1:
338+
case 3:
339+
return false;
340+
}
341+
switch (axisInPixelZ->FindBin(localZ)) {
342+
case 1:
343+
case 3:
344+
return false;
345+
}
346+
// hHitMap->Fill(z, rphi);
347+
// hHitMapInPixel->Fill(localZ, localRPhi);
348+
// if (static_cast<int>(hHitMapInPixel->GetEntries()) % 100000 == 0) {
349+
// hHitMap->SaveAs(Form("/tmp/%s.png", hHitMap->GetName()));
350+
// hHitMapInPixel->SaveAs(Form("/tmp/%s.png", hHitMapInPixel->GetName()));
351+
// }
327352
return true;
328353
}
329354

330-
/// Utility function to estimate hit position on a TOF layer
331-
/// \param track the track parameters
332-
/// \param radius the radius of the TOF layer
333-
/// \return estimated hit position [x, y, z] in cm
334-
std::array<float, 3> estimateHitPosition(const o2::track::TrackParCov& track)
335-
{
336-
std::array<float, 3> hitPosition = {0.0f, 0.0f, 0.0f};
337-
// Get the track circle parameters
338-
o2::math_utils::CircleXYf_t trcCircle;
339-
float sna, csa;
340-
track.getCircleParams(magField, trcCircle, sna, csa);
341-
342-
// Calculate intersection points with the cylinder of given radius
343-
const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
344-
345-
if (centerDistance < trcCircle.rC + layerRadius && centerDistance > std::fabs(trcCircle.rC - layerRadius)) {
346-
// Calculate intersection point using the same logic as computeTrackLength
347-
const float ux = trcCircle.xC / centerDistance;
348-
const float uy = trcCircle.yC / centerDistance;
349-
const float vx = -uy;
350-
const float vy = +ux;
351-
const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + layerRadius * layerRadius) / (2.0f * centerDistance);
352-
const float displace = (0.5f / centerDistance) * std::sqrt(
353-
(-centerDistance + trcCircle.rC - layerRadius) *
354-
(-centerDistance - trcCircle.rC + layerRadius) *
355-
(-centerDistance + trcCircle.rC + layerRadius) *
356-
(centerDistance + trcCircle.rC + layerRadius));
357-
358-
const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
359-
const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
360-
361-
// Choose the correct intersection point based on momentum direction
362-
std::array<float, 3> mom;
363-
track.getPxPyPzGlo(mom);
364-
const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
365-
const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
366-
367-
if (scalarProduct1 > scalarProduct2) {
368-
hitPosition[0] = point1[0];
369-
hitPosition[1] = point1[1];
370-
} else {
371-
hitPosition[0] = point2[0];
372-
hitPosition[1] = point2[1];
373-
}
374-
375-
// Estimate Z position using the track slope
376-
float trackLength = computeTrackLength(track, layerRadius, magField);
377-
if (trackLength > 0) {
378-
hitPosition[2] = track.getZ() + trackLength * track.getTgl();
379-
}
380-
}
381-
return hitPosition;
382-
}
383-
384355
/// Check if a track hits the active area (convenience function)
385356
/// \param track the track parameters for automatic hit position calculation
386357
/// \return true if the hit is in the active area
387358
bool isTrackInActiveArea(const o2::track::TrackParCov& track)
388359
{
389-
auto hitPosition = estimateHitPosition(track);
390-
return isInTOFActiveArea(hitPosition);
360+
if (fractionInactive <= 0.f) {
361+
return true;
362+
}
363+
float x, y, z;
364+
if (!track.getXatLabR(layerRadius, x, magField)) {
365+
LOG(warning) << "Could not propagate track to TOF layer at radius " << layerRadius << " cm";
366+
return false;
367+
}
368+
bool b;
369+
ROOT::Math::PositionVector3D hit = track.getXYZGloAt(x, magField, b);
370+
if (!b) {
371+
LOG(warning) << "Could not get hit position at radius " << layerRadius << " cm";
372+
return false;
373+
}
374+
hit.GetCoordinates(x, y, z);
375+
return isInTOFActiveArea(std::array<float, 3>{x, y, z});
391376
}
392377
};
393378

394379
bool isInInnerTOFActiveArea(const o2::track::TrackParCov& track)
395380
{
396-
static TOFLayerEfficiency innerTOFLayer(simConfig.innerTOFRadius, simConfig.innerTOFLength, simConfig.innerTOFPixelArea, simConfig.innerTOFFractionOfInactiveArea, simConfig.magneticField);
381+
static TOFLayerEfficiency innerTOFLayer(simConfig.innerTOFRadius, simConfig.innerTOFLength, simConfig.innerTOFPixelDimension, simConfig.innerTOFFractionOfInactiveArea, simConfig.magneticField);
397382
return innerTOFLayer.isTrackInActiveArea(track);
398383
}
399384

400385
bool isInOuterTOFActiveArea(const o2::track::TrackParCov& track)
401386
{
402-
static TOFLayerEfficiency outerTOFLayer(simConfig.outerTOFRadius, simConfig.outerTOFLength, simConfig.outerTOFPixelArea, simConfig.outerTOFFractionOfInactiveArea, simConfig.magneticField);
387+
static TOFLayerEfficiency outerTOFLayer(simConfig.outerTOFRadius, simConfig.outerTOFLength, simConfig.outerTOFPixelDimension, simConfig.outerTOFFractionOfInactiveArea, simConfig.magneticField);
403388
return outerTOFLayer.isTrackInActiveArea(track);
404389
}
405390

@@ -656,18 +641,22 @@ struct OnTheFlyTofPid {
656641
}
657642

658643
// Check if the track hit a sensitive area of the TOF
659-
const bool activeInnerTOF = true; // isInInnerTOFActiveArea(o2track);
660-
float x = 0.f;
661-
o2track.getXatLabR(simConfig.innerTOFRadius, x, simConfig.magneticField.value); // to get phi at inner TOF radius
662-
static_cast<TEfficiency*>(listEfficiency->At(1))->Fill(o2track.getZAt(x, simConfig.magneticField.value), simConfig.innerTOFRadius * o2track.getPhiAt(x, simConfig.magneticField.value), activeInnerTOF);
644+
const bool activeInnerTOF = isInInnerTOFActiveArea(o2track);
663645
if (!activeInnerTOF) {
664646
trackLengthInnerTOF = -999.f;
647+
} else {
648+
float x = 0.f;
649+
o2track.getXatLabR(simConfig.innerTOFRadius, x, simConfig.magneticField);
650+
histos.fill(HIST("iTOF/h2HitMap"), o2track.getZAt(x, simConfig.magneticField), simConfig.innerTOFRadius * o2track.getPhiAt(x, simConfig.magneticField));
665651
}
652+
666653
const bool activeOuterTOF = isInOuterTOFActiveArea(o2track);
667-
o2track.getXatLabR(simConfig.outerTOFRadius, x, simConfig.magneticField.value); // to get phi at outer TOF radius
668-
static_cast<TEfficiency*>(listEfficiency->At(2))->Fill(o2track.getZAt(x, simConfig.magneticField.value), simConfig.outerTOFRadius * o2track.getPhiAt(x, simConfig.magneticField.value), activeOuterTOF);
669654
if (!activeOuterTOF) {
670655
trackLengthOuterTOF = -999.f;
656+
} else {
657+
float x = 0.f;
658+
o2track.getXatLabR(simConfig.outerTOFRadius, x, simConfig.magneticField);
659+
histos.fill(HIST("oTOF/h2HitMap"), o2track.getZAt(x, simConfig.magneticField), simConfig.outerTOFRadius * o2track.getPhiAt(x, simConfig.magneticField));
671660
}
672661

673662
// get mass to calculate velocity

‎ALICE3/TableProducer/OTF/onTheFlyTracker.cxx‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -304,6 +304,7 @@ struct OnTheFlyTracker {
304304
hCovMatOK->GetXaxis()->SetBinLabel(2, "OK");
305305

306306
histos.add("hPtGenerated", "hPtGenerated", kTH1F, {axes.axisMomentum});
307+
histos.add("hPhiGenerated", "hPhiGenerated", kTH1F, {{100, 0.0f, 2 * M_PI, "#phi (rad)"}});
307308
histos.add("hPtGeneratedEl", "hPtGeneratedEl", kTH1F, {axes.axisMomentum});
308309
histos.add("hPtGeneratedPi", "hPtGeneratedPi", kTH1F, {axes.axisMomentum});
309310
histos.add("hPtGeneratedKa", "hPtGeneratedKa", kTH1F, {axes.axisMomentum});
@@ -614,6 +615,7 @@ struct OnTheFlyTracker {
614615
}
615616

616617
histos.fill(HIST("hPtGenerated"), mcParticle.pt());
618+
histos.fill(HIST("hPhiGenerated"), mcParticle.phi());
617619
if (std::abs(mcParticle.pdgCode()) == kElectron)
618620
histos.fill(HIST("hPtGeneratedEl"), mcParticle.pt());
619621
if (std::abs(mcParticle.pdgCode()) == kPiPlus)

0 commit comments

Comments
 (0)