Skip to content

Commit

Permalink
TPC Cluster error: should average over 1/sqrt(charge) but use non-sqr…
Browse files Browse the repository at this point in the history 
…t charge
  • Loading branch information
@davidrohr
davidrohr committed Mar 13, 2024
1 parent 94a4838 commit f1c87d2
Show file tree
Hide file tree
Showing 7 changed files with 41 additions and 37 deletions.
8 changes: 4 additions & 4 deletions GPU/GPUTracking/Base/GPUParam.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,12 +96,12 @@ struct GPUParam : public internal::GPUParam_t<GPUSettingsRec, GPUSettingsParam>
}
return 0.174533f + par.dAlpha * iSlice;
}
GPUd() float GetClusterErrorSeeding(int yz, int type, float z, float angle2) const;
GPUd() void GetClusterErrorsSeeding2(char sector, int row, float z, float sinPhi, float DzDs, float time, float avgCharge, float charge, float& ErrY2, float& ErrZ2) const;
GPUd() float GetClusterErrorSeeding(int yz, int type, float zDiff, float angle2) const;
GPUd() void GetClusterErrorsSeeding2(char sector, int row, float z, float sinPhi, float DzDs, float time, float avgInvCharge, float invCharge, float& ErrY2, float& ErrZ2) const;
GPUd() float GetSystematicClusterErrorIFC2(float x, float y, float z, bool sideC) const;

GPUd() float GetClusterError2(int yz, int type, float z, float angle2, float scaledMult, float scaledAvgCharge, float scaledCharge) const;
GPUd() void GetClusterErrors2(char sector, int row, float z, float sinPhi, float DzDs, float time, float avgCharge, float charge, float& ErrY2, float& ErrZ2) const;
GPUd() float GetClusterError2(int yz, int type, float zDiff, float angle2, float scaledMult, float scaledAvgInvCharge, float scaledInvCharge) const;
GPUd() void GetClusterErrors2(char sector, int row, float z, float sinPhi, float DzDs, float time, float avgInvCharge, float invCharge, float& ErrY2, float& ErrZ2) const;
GPUd() void UpdateClusterError2ByState(short clusterState, float& ErrY2, float& ErrZ2) const;
GPUd() float GetScaledMult(float time) const;

Expand Down
30 changes: 15 additions & 15 deletions GPU/GPUTracking/Base/GPUParam.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,16 +47,16 @@ GPUdi() void MEM_LG(GPUParam)::Global2Slice(int iSlice, float X, float Y, float
#ifdef GPUCA_TPC_GEOMETRY_O2

MEM_CLASS_PRE()
GPUdi() void MEM_LG(GPUParam)::GetClusterErrorsSeeding2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgCharge, float charge, float& ErrY2, float& ErrZ2) const
GPUdi() void MEM_LG(GPUParam)::GetClusterErrorsSeeding2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgInvCharge, float invCharge, float& ErrY2, float& ErrZ2) const
{
GetClusterErrors2(sector, iRow, z, sinPhi, DzDs, time, avgCharge, charge, ErrY2, ErrZ2);
GetClusterErrors2(sector, iRow, z, sinPhi, DzDs, time, avgInvCharge, invCharge, ErrY2, ErrZ2);
}

MEM_CLASS_PRE()
GPUdi() float MEM_LG(GPUParam)::GetClusterError2(int yz, int type, float z, float angle2, float scaledMult, float scaledAvgCharge, float scaledCharge) const
GPUdi() float MEM_LG(GPUParam)::GetClusterError2(int yz, int type, float zDiff, float angle2, float scaledMult, float scaledInvAvgCharge, float scaledInvCharge) const
{
MakeType(const float*) c = ParamErrors[yz][type];
float v = c[0] + c[1] * angle2 * scaledAvgCharge + c[2] * z * scaledCharge + c[3] * (scaledMult * scaledMult) * (scaledAvgCharge * scaledAvgCharge);
MakeType(const float*) c = ParamErrors[yz][type]; // Note: c[0] = p[0]^2, c[1] = p[1]^2 * padHeight, c[2] = p[2]^2 / tpcLength / padHeight, c[3] = p[3]^2
float v = c[0] + c[1] * angle2 * scaledInvAvgCharge + c[2] * zDiff * scaledInvCharge + c[3] * (scaledMult * scaledMult) * (scaledInvAvgCharge * scaledInvAvgCharge);
v = CAMath::Abs(v);
v *= yz ? rec.tpc.clusterError2CorrectionZ : rec.tpc.clusterError2CorrectionY;
v += yz ? rec.tpc.clusterError2AdditionalZ : rec.tpc.clusterError2AdditionalY;
Expand Down Expand Up @@ -100,16 +100,16 @@ GPUdi() float MEM_LG(GPUParam)::GetSystematicClusterErrorIFC2(float x, float y,
#else // GPUCA_TPC_GEOMETRY_O2

MEM_CLASS_PRE()
GPUdi() float MEM_LG(GPUParam)::GetClusterErrorSeeding(int yz, int type, float z, float angle2) const
GPUdi() float MEM_LG(GPUParam)::GetClusterErrorSeeding(int yz, int type, float zDiff, float angle2) const
{
MakeType(const float*) c = ParamErrorsSeeding0[yz][type];
float v = c[0] + c[1] * z + c[2] * angle2;
float v = c[0] + c[1] * zDiff + c[2] * angle2;
v = CAMath::Abs(v);
return v;
}

MEM_CLASS_PRE()
GPUdi() void MEM_LG(GPUParam)::GetClusterErrorsSeeding2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgCharge, float charge, float& ErrY2, float& ErrZ2) const
GPUdi() void MEM_LG(GPUParam)::GetClusterErrorsSeeding2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgInvCharge, float invCharge, float& ErrY2, float& ErrZ2) const
{
int rowType = tpcGeometry.GetROC(iRow);
z = CAMath::Abs(tpcGeometry.TPCLength() - CAMath::Abs(z));
Expand All @@ -128,10 +128,10 @@ GPUdi() void MEM_LG(GPUParam)::GetClusterErrorsSeeding2(char sector, int iRow, f
}

MEM_CLASS_PRE()
GPUdi() float MEM_LG(GPUParam)::GetClusterError2(int yz, int type, float z, float angle2, float time, float avgCharge, float charge) const
GPUdi() float MEM_LG(GPUParam)::GetClusterError2(int yz, int type, float zDiff, float angle2, float time, float avgInvCharge, float invCharge) const
{
MakeType(const float*) c = ParamS0Par[yz][type];
float v = c[0] + c[1] * z + c[2] * angle2 + c[3] * z * z + c[4] * angle2 * angle2 + c[5] * z * angle2;
float v = c[0] + c[1] * zDiff + c[2] * angle2 + c[3] * zDiff * zDiff + c[4] * angle2 * angle2 + c[5] * zDiff * angle2;
v = CAMath::Abs(v);
if (v < 0.0001f) {
v = 0.0001f;
Expand All @@ -150,7 +150,7 @@ GPUdi() float MEM_LG(GPUParam)::GetSystematicClusterErrorIFC2(float x, float y,
#endif // !GPUCA_TPC_GEOMETRY_O2

MEM_CLASS_PRE()
GPUdi() void MEM_LG(GPUParam)::GetClusterErrors2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgCharge, float charge, float& ErrY2, float& ErrZ2) const
GPUdi() void MEM_LG(GPUParam)::GetClusterErrors2(char sector, int iRow, float z, float sinPhi, float DzDs, float time, float avgInvCharge, float invCharge, float& ErrY2, float& ErrZ2) const
{
// Calibrated cluster error from OCDB for Y and Z
const int rowType = tpcGeometry.GetROC(iRow);
Expand All @@ -164,11 +164,11 @@ GPUdi() void MEM_LG(GPUParam)::GetClusterErrors2(char sector, int iRow, float z,
const float angleZ2 = DzDs * DzDs * sec2; // dz/dx

const float mult = time >= 0.f ? GetScaledMult(time) / tpcGeometry.Row2X(iRow) : 0.f;
const float scaledAvgCharge = avgCharge * rec.tpc.clusterErrorChargeScaler > 0.f ? avgCharge * rec.tpc.clusterErrorChargeScaler : 1.f;
const float scaledCharge = charge * rec.tpc.clusterErrorChargeScaler > 0.f ? charge * rec.tpc.clusterErrorChargeScaler : 1.f;
const float scaledInvAvgCharge = avgInvCharge * rec.tpc.clusterErrorChargeScaler > 0.f ? avgInvCharge * rec.tpc.clusterErrorChargeScaler : 1.f;
const float scaledInvCharge = invCharge * rec.tpc.clusterErrorChargeScaler > 0.f ? invCharge * rec.tpc.clusterErrorChargeScaler : 1.f;

ErrY2 = GetClusterError2(0, rowType, z, angleY2, mult, scaledAvgCharge, scaledCharge);
ErrZ2 = GetClusterError2(1, rowType, z, angleZ2, mult, scaledAvgCharge, scaledCharge);
ErrY2 = GetClusterError2(0, rowType, z, angleY2, mult, scaledInvAvgCharge, scaledInvCharge);
ErrZ2 = GetClusterError2(1, rowType, z, angleZ2, mult, scaledInvAvgCharge, scaledInvCharge);
}

MEM_CLASS_PRE()
Expand Down
4 changes: 2 additions & 2 deletions GPU/GPUTracking/Merger/GPUTPCGMMerger.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,8 +210,8 @@ class GPUTPCGMMerger : public GPUProcessor
std::vector<unsigned int> StreamerOccupancyBin(int iSlice, int iRow, float time) const;
std::vector<float> StreamerUncorrectedZY(int iSlice, int iRow, const GPUTPCGMTrackParam& track, const GPUTPCGMPropagator& prop) const;

void DebugStreamerUpdate(int iTrk, int ihit, float xx, float yy, float zz, const GPUTPCGMMergedTrackHit& cluster, const o2::tpc::ClusterNative& clusterNative, const GPUTPCGMTrackParam& track, const GPUTPCGMPropagator& prop, const gputpcgmmergertypes::InterpolationErrorHit& interpolation, char rejectChi2, bool refit, int retVal, float avgCharge, float charge) const;
static void DebugStreamerReject(float mAlpha, int iRow, float posY, float posZ, short clusterState, char rejectChi2, const gputpcgmmergertypes::InterpolationErrorHit& inter, bool refit, int retVal, float err2Y, float err2Z, const GPUTPCGMTrackParam& track, const GPUParam& param, float time, float avgCharge, float charge);
void DebugStreamerUpdate(int iTrk, int ihit, float xx, float yy, float zz, const GPUTPCGMMergedTrackHit& cluster, const o2::tpc::ClusterNative& clusterNative, const GPUTPCGMTrackParam& track, const GPUTPCGMPropagator& prop, const gputpcgmmergertypes::InterpolationErrorHit& interpolation, char rejectChi2, bool refit, int retVal, float avgInvCharge) const;
static void DebugStreamerReject(float mAlpha, int iRow, float posY, float posZ, short clusterState, char rejectChi2, const gputpcgmmergertypes::InterpolationErrorHit& inter, bool refit, int retVal, float err2Y, float err2Z, const GPUTPCGMTrackParam& track, const GPUParam& param, float time, float avgInvCharge, float invCharge);
#endif

GPUdi() int SliceTrackInfoFirst(int iSlice) const { return mSliceTrackInfoIndex[iSlice]; }
Expand Down
13 changes: 7 additions & 6 deletions GPU/GPUTracking/Merger/GPUTPCGMMergerDump.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -293,12 +293,13 @@ std::vector<float> GPUTPCGMMerger::StreamerUncorrectedZY(int iSlice, int iRow, c
return retVal;
}

void GPUTPCGMMerger::DebugStreamerUpdate(int iTrk, int ihit, float xx, float yy, float zz, const GPUTPCGMMergedTrackHit& cluster, const o2::tpc::ClusterNative& clusterNative, const GPUTPCGMTrackParam& track, const GPUTPCGMPropagator& prop, const gputpcgmmergertypes::InterpolationErrorHit& interpolation, char rejectChi2, bool refit, int retVal, float avgCharge, float charge) const
void GPUTPCGMMerger::DebugStreamerUpdate(int iTrk, int ihit, float xx, float yy, float zz, const GPUTPCGMMergedTrackHit& cluster, const o2::tpc::ClusterNative& clusterNative, const GPUTPCGMTrackParam& track, const GPUTPCGMPropagator& prop, const gputpcgmmergertypes::InterpolationErrorHit& interpolation, char rejectChi2, bool refit, int retVal, float avgInvCharge) const
{
#ifdef DEBUG_STREAMER
float time = clusterNative.getTime();
auto occupancyBins = StreamerOccupancyBin(cluster.slice, cluster.row, time);
auto uncorrectedYZ = StreamerUncorrectedZY(cluster.slice, cluster.row, track, prop);
float invCharge = 1.f / clusterNative.qMax;
o2::utils::DebugStreamer::instance()->getStreamer("debug_update_track", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_update_track").data()
<< "iTrk=" << iTrk
<< "ihit=" << ihit
Expand All @@ -314,13 +315,13 @@ void GPUTPCGMMerger::DebugStreamerUpdate(int iTrk, int ihit, float xx, float yy,
<< "retVal=" << retVal
<< "occupancyBins=" << occupancyBins
<< "trackUncorrectedYZ=" << uncorrectedYZ
<< "avgCharge=" << avgCharge
<< "charge=" << charge
<< "avgInvCharge=" << avgInvCharge
<< "invCharge=" << invCharge
<< "\n";
#endif
}

void GPUTPCGMMerger::DebugStreamerReject(float mAlpha, int iRow, float posY, float posZ, short clusterState, char rejectChi2, const gputpcgmmergertypes::InterpolationErrorHit& inter, bool refit, int retVal, float err2Y, float err2Z, const GPUTPCGMTrackParam& track, const GPUParam& param, float time, float avgCharge, float charge)
void GPUTPCGMMerger::DebugStreamerReject(float mAlpha, int iRow, float posY, float posZ, short clusterState, char rejectChi2, const gputpcgmmergertypes::InterpolationErrorHit& inter, bool refit, int retVal, float err2Y, float err2Z, const GPUTPCGMTrackParam& track, const GPUParam& param, float time, float avgInvCharge, float invCharge)
{
#ifdef DEBUG_STREAMER
float scaledMult = (time >= 0.f ? param.GetScaledMult(time) / param.tpcGeometry.Row2X(iRow) : 0.f);
Expand All @@ -338,8 +339,8 @@ void GPUTPCGMMerger::DebugStreamerReject(float mAlpha, int iRow, float posY, flo
<< "err2Z=" << err2Z
<< "track=" << track
<< "scaledMultiplicity=" << scaledMult
<< "avgCharge=" << avgCharge
<< "charge=" << charge
<< "avgInvCharge=" << avgInvCharge
<< "invCharge=" << invCharge
<< "\n";
#endif
}
6 changes: 3 additions & 3 deletions GPU/GPUTracking/Merger/GPUTPCGMPropagator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -653,18 +653,18 @@ GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, float err2Y
}
}

GPUd() int GPUTPCGMPropagator::Update(float posY, float posZ, int iRow, const GPUParam& GPUrestrict() param, short clusterState, char rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, char sector, float time, float avgCharge, float charge GPUCA_DEBUG_STREAMER_CHECK(, int iTrk))
GPUd() int GPUTPCGMPropagator::Update(float posY, float posZ, int iRow, const GPUParam& GPUrestrict() param, short clusterState, char rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, char sector, float time, float avgInvCharge, float invCharge GPUCA_DEBUG_STREAMER_CHECK(, int iTrk))
{
float err2Y, err2Z;
GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgCharge, charge);
GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgInvCharge, invCharge);

if (rejectChi2 >= rejectInterFill) {
if (rejectChi2 == rejectInterReject && inter->errorY < (GPUCA_MERGER_INTERPOLATION_ERROR_TYPE)0) {
rejectChi2 = rejectDirect;
} else {
int retVal = InterpolateReject(param, posY, posZ, clusterState, rejectChi2, inter, err2Y, err2Z);
GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamRejectCluster, iTrk)) {
GPUTPCGMMerger::DebugStreamerReject(mAlpha, iRow, posY, posZ, clusterState, rejectChi2, *inter, refit, retVal, err2Y, err2Z, *mT, param, time, avgCharge, charge);
GPUTPCGMMerger::DebugStreamerReject(mAlpha, iRow, posY, posZ, clusterState, rejectChi2, *inter, refit, retVal, err2Y, err2Z, *mT, param, time, avgInvCharge, invCharge);
});
if (retVal) {
return retVal;
Expand Down
2 changes: 1 addition & 1 deletion GPU/GPUTracking/Merger/GPUTPCGMPropagator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ class GPUTPCGMPropagator

GPUd() int PropagateToXAlphaBz(float posX, float posAlpha, bool inFlyDirection);

GPUd() int Update(float posY, float posZ, int iRow, const GPUParam& param, short clusterState, char rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, char sideC, float time, float avgCharge, float charge GPUCA_DEBUG_STREAMER_CHECK(, int iTrk = 0));
GPUd() int Update(float posY, float posZ, int iRow, const GPUParam& param, short clusterState, char rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, char sideC, float time, float avgInvCharge, float invCharge GPUCA_DEBUG_STREAMER_CHECK(, int iTrk = 0));
GPUd() int Update(float posY, float posZ, short clusterState, bool rejectChi2, float err2Y, float err2Z, const GPUParam* param = nullptr);
GPUd() int InterpolateReject(const GPUParam& param, float posY, float posZ, short clusterState, char rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z);
GPUd() float PredictChi2(float posY, float posZ, int iRow, const GPUParam& param, short clusterState, char sideC, float time, float avgCharge, float charge) const;
Expand Down
15 changes: 9 additions & 6 deletions GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int iT

for (int iWay = 0; iWay < nWays; iWay++) {
int nMissed = 0, nMissed2 = 0;
float avgCharge = 0.f;
float sumInvSqrtCharge = 0.f;
int nAvgCharge = 0;

if (iWay && storeOuter != 255 && param.rec.tpc.nWaysOuter && outerParam) {
Expand Down Expand Up @@ -322,15 +322,18 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int iT
tup.Fill((float)cluster.row, xx, yy, zz, clAlpha, mX, ImP0, ImP1, mP[2], mP[3], mP[4], ImC0, ImC2, mC[14]);
}
#endif
float time = merger->Param().par.earlyTpcTransform ? -1.f : merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].getTime();
float tmpCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? CAMath::InvSqrt(merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
if (merger->Param().rec.tpc.rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && merger->Param().rejectEdgeClusterByY(uncorrectedY, cluster.row)) {
if (merger->Param().rec.tpc.rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && merger->Param().rejectEdgeClusterByY(uncorrectedY, cluster.row)) { // uncorrectedY > -1e6f implies allowModification
retVal = GPUTPCGMPropagator::updateErrorEdgeCluster;
} else {
retVal = prop.Update(yy, zz, cluster.row, param, clusterState, rejectChi2, &interpolation.hit[ihit], refit, cluster.slice, time, (avgCharge += tmpCharge) / ++nAvgCharge, tmpCharge GPUCA_DEBUG_STREAMER_CHECK(, iTrk)); // TODO: Use avgCharge
const float time = merger->Param().par.earlyTpcTransform ? -1.f : merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].getTime();
const float invSqrtCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? CAMath::InvSqrt(merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
const float invCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? (1.f / merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
float invAvgCharge = (sumInvSqrtCharge += invSqrtCharge) / ++nAvgCharge;
invAvgCharge *= invAvgCharge;
retVal = prop.Update(yy, zz, cluster.row, param, clusterState, rejectChi2, &interpolation.hit[ihit], refit, cluster.slice, time, invAvgCharge, invCharge GPUCA_DEBUG_STREAMER_CHECK(, iTrk));
}
GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamUpdateTrack, iTrk)) {
merger->DebugStreamerUpdate(iTrk, ihit, xx, yy, zz, cluster, merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num], *this, prop, interpolation.hit[ihit], rejectChi2, refit, retVal, avgCharge / nAvgCharge, tmpCharge);
merger->DebugStreamerUpdate(iTrk, ihit, xx, yy, zz, cluster, merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num], *this, prop, interpolation.hit[ihit], rejectChi2, refit, retVal, sumInvSqrtCharge / nAvgCharge * sumInvSqrtCharge / nAvgCharge);
});
}
// clang-format off
Expand Down

0 comments on commit f1c87d2

Please sign in to comment.