diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index edb44394317eb..02a3eb08612a4 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -153,9 +153,9 @@ void TRDGlobalTracking::updateTimeDependentParams(ProcessingContext& pc) mBase->init(pc); mBase->setLocalGainFactors(pc.inputs().get("localgainfactors").get()); } - + pc.inputs().get*>("chamberstatus"); // called to trigger finaliseCCDB - pc.inputs().get("padstatus"); // called to trigger finaliseCCDB + pc.inputs().get("padstatus"); // called to trigger finaliseCCDB } const auto& trackTune = TrackTuneParams::Instance(); @@ -207,9 +207,8 @@ void TRDGlobalTracking::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { if ((*chamberStatus)[iDet] == 3) { mTracker->SetChamberStatus(iDet, false); // chamber is good - } - else { - mTracker->SetChamberStatus(iDet, true); // chamber is bad + } else { + mTracker->SetChamberStatus(iDet, true); // chamber is bad } } return; @@ -219,14 +218,13 @@ void TRDGlobalTracking::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) const o2::trd::PadStatus* padStatus = (const o2::trd::PadStatus*)obj; for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { for (int iCol = 0; iCol < constants::NCOLUMN; iCol++) { - for (int iRow = 0; iRow < ((iDet%30)/6 == 2 ? constants::NROWC0 : constants::NROWC1); iRow++) { + for (int iRow = 0; iRow < ((iDet % 30) / 6 == 2 ? constants::NROWC0 : constants::NROWC1); iRow++) { if (padStatus->isMasked(iDet, iCol, iRow) || padStatus->isNotConnected(iDet, iCol, iRow)) { mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, true); // pad is masked - } - else { + } else { mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, false); // pad is not masked } - } + } } } return; @@ -533,7 +531,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc) if (trdTrack.getChi2() / trdTrack.getNtracklets() > mTracker->Param().rec.trd.maxChi2Red) { continue; } - + // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets float maxProb = 0.f; // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability @@ -806,7 +804,7 @@ bool TRDGlobalTracking::refitTPCTRDTrack(TrackTRD& trk, float timeTRD, o2::globa return false; } if (pileUpOn) { // account pileup time uncertainty in Z errors - //timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS(); + // timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS(); timeZErr = mTPCVdrift * mTPCVdrift * mTErrPileUp2; outerParam.updateCov(timeZErr, o2::track::CovLabels::kSigZ2); } @@ -930,7 +928,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards, int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet); float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber)); - + // Correction of y position based on angular pull if (mRec->GetParam().rec.trd.useAngularPull == 3 || mRec->GetParam().rec.trd.useAngularPull == 4) { float corrPull = -angularPull * mRecoParam.getCorrYDy(trkParam->getSnp()); @@ -1046,7 +1044,7 @@ DataProcessorSpec getTRDGlobalTrackingSpec(bool useMC, GTrackID::mask_t src, boo // request calibration data inputs.emplace_back("localgainfactors", "TRD", "LOCALGAINFACTORS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/LocalGainFactor")); } - + // request list of bad chambers and masked pads to estimate better the number of findable tracklets inputs.emplace_back("chamberstatus", "TRD", "CHAMBERSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/DCSDPsFedChamberStatus")); inputs.emplace_back("padstatus", "TRD", "PADSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/PadStatus")); diff --git a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h index e75c1e77a2b67..61ed436c5c665 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h @@ -136,32 +136,31 @@ GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy); } -GPUdi() double GPUTRDRecoParam::getDyLikelihood(float snp, float slope, int occupancy) const -{ +GPUdi() double GPUTRDRecoParam::getDyLikelihood(float snp, float slope, int occupancy) const +{ // Gaussian + left exponential + right exponential - double likelihood = CAMath::Exp(- 0.5f * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) / GPUTRDRecoParam::getDyRes(snp, occupancy)); - + double likelihood = CAMath::Exp(-0.5f * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) / GPUTRDRecoParam::getDyRes(snp, occupancy)); + // Normalization for the exponential is parametrized with respect to the Gaussian, it is smaller at lorentz angle double expNorm = CAMath::Sqrt(mDyExpNormA + mDyExpNormC * (snp - mLorentzAngle) * (snp - mLorentzAngle)); if (slope < convertAngleToDy(snp)) { // left tail in this case, larger tail for large positive snp likelihood += expNorm * CAMath::Exp((mDyExpA + mDyExpC * (snp - mLorentzAngle)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp))); - } - else { + } else { // right tail, larger tail for large negative snp likelihood += expNorm * CAMath::Exp(-(mDyExpA - mDyExpC * (snp - mLorentzAngle)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp))); } - + // Normalized such that likelihood is 1 when track angle and tracklet slope agree, and always lower than 1 otherwise (such that -log(likelihood) is always positive) likelihood /= (1.f + expNorm); return likelihood; } -GPUdi() double GPUTRDRecoParam::getZLikelihood(float deltaZ, float padLength, float sigmaZtrk) const +GPUdi() double GPUTRDRecoParam::getZLikelihood(float deltaZ, float padLength, float sigmaZtrk) const { // logistic function as approximation of convolution between uniform tracklet Z and gaussian track, depends on sqrt(3)/pi times the track resolution // normalized so that maximum is 1 - double lmax = 1.f/(1.f + CAMath::Exp(- 0.5f * padLength / 0.5513f / sigmaZtrk)) - 1.f/(1.f + CAMath::Exp(0.5f * padLength / 0.5513f / sigmaZtrk)); + double lmax = 1.f / (1.f + CAMath::Exp(-0.5f * padLength / 0.5513f / sigmaZtrk)) - 1.f / (1.f + CAMath::Exp(0.5f * padLength / 0.5513f / sigmaZtrk)); return 1.f / lmax / (1.f + CAMath::Exp((deltaZ - 0.5f * padLength) / 0.5513f / sigmaZtrk)) - 1.f / lmax / (1.f + CAMath::Exp((deltaZ + 0.5f * padLength) / 0.5513f / sigmaZtrk)); } diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 749667712226e..b48a2368319a0 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -635,14 +635,13 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (probBC > maxProb) { maxProb = probBC; yCorrPileUp = -slopeFactor * deltaBC; - zShiftTrkPileUp = - deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide; + zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide; } } if (sumProb > 1e-6f) { yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; } } - // number of tracklets within the chamber is the current TRD occupancy estimator int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; @@ -652,13 +651,13 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK float zPosCorr = spacePoints[trkltIdx].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl(); float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp; zPosCorr -= zShiftTrk + zShiftTrkPileUp; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision - + // Correction of y position based on angular pull if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) { float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp()); yPosCorr += corrPull; } - + float deltaY = yPosCorr - projY; float deltaZ = zPosCorr - projZ; @@ -677,18 +676,18 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK trkltCovTmpWithDy[0] += trkWork->getSigmaY2(); trkltCovTmpWithDy[1] += trkWork->getSigmaZY(); trkltCovTmpWithDy[2] += trkWork->getSigmaZ2(); - + if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) { // In this case the correlation between tracklet y and tracklet dy is already corrected for so we can cancel it trkltCovTmpWithDy[3] = 0.; trkltCovTmpWithDy[4] = 0.; } - + // We add the correlation between track y and track dy trkltCovTmpWithDy[3] += trkWork->getSigmaSnpY() * mGeo->GetCdrHght(); trkltCovTmpWithDy[4] += trkWork->getSigmaSnpZ() * mGeo->GetCdrHght(); - - // keep these parameters before matrix inversion + + // keep these parameters before matrix inversion float sigmaZ2 = trkltCovTmpWithDy[2]; float sigmaDy2 = trkltCovTmpWithDy[5]; @@ -698,13 +697,15 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (Param().rec.trd.addDeflectionInChi2 == 2 || Param().rec.trd.addDeflectionInChi2 == 3) { // In this case we take into account the full likelihood, so we replace (deltaDy/sigmaDy)^2 by -2*ln(likelihood), which is the same in the default Gaussian case double likelihood = mRecoParam->getDyLikelihood(trkWork->getSnp(), spacePoints[trkltIdx].getDy() + dyTiltCorr, nTrackletsChamber); - if (likelihood < 1e-6f) continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation + if (likelihood < 1e-6f) + continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation deltaDy = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaDy2) * (deltaDy > 0.f ? 1.f : -1.f); } if (Param().rec.trd.addDeflectionInChi2 == 3) { // We do the same for deltaZ double likelihood = mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2())); - if (likelihood < 1e-6f) continue; + if (likelihood < 1e-6f) + continue; deltaZ = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaZ2) * (deltaZ > 0.f ? 1.f : -1.f); } chi2 = deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy; @@ -724,8 +725,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK Hypothesis hypoNoUpdate(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2() + Param().rec.trd.penaltyChi2); InsertHypothesis(hypoNoUpdate, nCurrHypothesis, hypothesisIdxOffset); isOK = true; - } // end candidate loop - + } // end candidate loop mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->getChi2(), iLayer); // only meaningful for ONE candidate!!! mDebug->SetRoad(roadY, roadZ, iLayer); // only meaningful for ONE candidate @@ -791,9 +791,9 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK pad = mGeo->GetPadPlane(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()); float tiltCorrUp = tilt * (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ()); float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght(); - + float yPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp; - + float zPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl(); zPosCorrUp -= zShiftTrk; float padLength = pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()); @@ -825,21 +825,21 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (probBC > maxProb) { maxProb = probBC; yCorrPileUp = -slopeFactor * deltaBC; - zShiftTrkPileUp = - deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide; + zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide; } } if (sumProb > 1e-6f) { yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp; } } - + zPosCorrUp -= zShiftTrkPileUp; yPosCorrUp += yCorrPileUp; const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector(); int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet]; float angularPull = GetAngularPull(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber); - + // Correction of y position based on angular pull if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) { float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp()); @@ -938,8 +938,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK return false; } } // end layer loop - - // -------------------------------------------------------------------------------- // add some debug information (compare labels of attached tracklets to track label) @@ -1288,7 +1286,7 @@ GPUd() bool GPUTRDTracker_t::IsGeoFindable(const TRDTRK* t, const if (!mGeo->ChamberInGeometry(det)) { return false; } - + // reject tracks in bad chambers if (mChamberStatus[det]) { return false; @@ -1310,12 +1308,12 @@ GPUd() bool GPUTRDTracker_t::IsGeoFindable(const TRDTRK* t, const if (!((zTrk > zMin + epsZ) && (zTrk < zMax - epsZ))) { return false; } - + // reject tracks if the corresponding pad and neighboring pads (to take into account track uncertainty and charge sharing with neighbours) are masked int32_t padrow = pp->GetPadRowNumber(zTrk); int32_t padcol = pp->GetPadColNumber(t->getY()); int32_t idxPad = det * kNPadColumns * kNPadRows + padcol * kNPadRows + padrow; - if (padrow != -1 && padcol != -1 && mPadStatus[idxPad] && (padrow == 0 || mPadStatus[idxPad - kNPadRows]) && (padrow == kNPadRows-1 || mPadStatus[idxPad + kNPadRows])) { + if (padrow != -1 && padcol != -1 && mPadStatus[idxPad] && (padrow == 0 || mPadStatus[idxPad - kNPadRows]) && (padrow == kNPadRows - 1 || mPadStatus[idxPad + kNPadRows])) { return false; } diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index 70b44ee21130a..246cd3ab3d0c1 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -72,7 +72,7 @@ class GPUTRDTracker_t : public GPUProcessor kNChambers = 540, kNPadRows = 16, kNPadColumns = 144, - kNPads = 1244160}; + kNPads = 1244160 }; struct HelperTrackAttributes { // additional TRD track attributes which are transient