@@ -406,6 +406,25 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
406406 }
407407 LOGF (info, " %i tracks are loaded into the TRD tracker. Out of those %i ITS-TPC tracks and %i TPC tracks" , nTracksLoadedITSTPC + nTracksLoadedTPC, nTracksLoadedITSTPC, nTracksLoadedTPC);
408408
409+ // Load the FT0 triggered BCs if this is requested
410+
411+ if (mTrkMask [GTrackID::FT0 ]) { // pile-up tagging was requested
412+ auto ft0recPoints = inputTracks.getFT0RecPoints ();
413+ uint32_t firstOrbit = 0 ;
414+ for (size_t ft0id = 0 ; ft0id < ft0recPoints.size (); ft0id++) {
415+ const auto & f0rec = ft0recPoints[ft0id];
416+ if (ft0id == 0 ) {
417+ firstOrbit = f0rec.getInteractionRecord ().orbit ;
418+ }
419+ if (o2::ft0::InteractionTag::Instance ().isSelected (f0rec)) {
420+ uint32_t currentOrbit = f0rec.getInteractionRecord ().orbit ;
421+ mTriggeredBCFT0 .push_back (f0rec.getInteractionRecord ().bc + (currentOrbit - firstOrbit) * o2::constants::lhc::LHCMaxBunches);
422+ }
423+ }
424+ }
425+
426+ mTracker ->SetFT0TriggeredBC (mTriggeredBCFT0 .data (), mTriggeredBCFT0 .size ());
427+
409428 // start the tracking
410429 // mTracker->DumpTracks();
411430 mChainTracking ->DoTRDGPUTracking <GPUTRDTrackerKernels::o2Version>(mTracker );
@@ -788,6 +807,46 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
788807 }
789808 }
790809
810+ // Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets
811+ float tCorrPileUp = 0 .;
812+ float tErrPileUp2 = 0 ;
813+ float maxProb = 0 .f ;
814+ // The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
815+ float sumCorr = 0 .f ;
816+ float sumCorr2 = 0 .f ;
817+ float sumProb = 0 .f ;
818+ for (int iBC = 0 ; iBC < mTriggeredBCFT0 .size (); iBC++) {
819+ int deltaBC = roundf (mTriggeredBCFT0 [iBC] - mChainTracking ->mIOPtrs .trdTriggerTimes [trk.getCollisionId ()] / o2::constants::lhc::LHCBunchSpacingMUS);
820+ if (deltaBC <= mRecoParam .getPileUpRangeBefore () || deltaBC >= mRecoParam .getPileUpRangeAfter ()) {
821+ continue ;
822+ }
823+ // collect the charges
824+ std::array<int , 6 > q0;
825+ std::array<int , 6 > q1;
826+ for (int iLy = 0 ; iLy < NLAYER ; iLy++) {
827+ int trkltId = trk.getTrackletIndex (iLy);
828+ if (trkltId < 0 ) {
829+ q0[iLy] = -1 ;
830+ q1[iLy] = -1 ;
831+ } else {
832+ q0[iLy] = mTrackletsRaw [trkltId].getQ0 ();
833+ q1[iLy] = mTrackletsRaw [trkltId].getQ1 ();
834+ }
835+ }
836+ // get pile-up probability
837+ float probBC = mRecoParam .getPileUpProbTrack (deltaBC, q0, q1);
838+ sumCorr += probBC * deltaBC;
839+ sumCorr2 += probBC * deltaBC * deltaBC;
840+ sumProb += probBC;
841+ if (probBC > maxProb) {
842+ maxProb = probBC;
843+ tCorrPileUp = -deltaBC;
844+ }
845+ }
846+ if (sumProb > 1e-6 ) {
847+ tErrPileUp2 = sumCorr2 / sumProb - 2 * tCorrPileUp * sumCorr / sumProb + tCorrPileUp * tCorrPileUp;
848+ }
849+
791850 if (inwards) {
792851 // reset covariance to something big for inwards refit
793852 trkParam->resetCovariance (100 );
@@ -811,16 +870,26 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
811870 }
812871 const PadPlane* pad = Geometry::instance ()->getPadPlane (trkltDet);
813872 float tilt = tan (TMath::DegToRad () * pad->getTiltingAngle ()); // tilt is signed! and returned in degrees
873+ float dyTiltCorr = tilt * trkParam->getTgl () * Geometry::instance ()->cdrHght ();
814874 float tiltCorrUp = tilt * (mTrackletsCalib [trkltId].getZ () - trkParam->getZ ());
815875 float zPosCorrUp = mTrackletsCalib [trkltId].getZ () + mRecoParam .getZCorrCoeffNRC () * trkParam->getTgl ();
816876 float padLength = pad->getRowSize (mTrackletsRaw [trkltId].getPadRow ());
817877 if (!((trkParam->getSigmaZ2 () < (padLength * padLength / 12 .f )) && (std::fabs (mTrackletsCalib [trkltId].getZ () - trkParam->getZ ()) < padLength))) {
818878 tiltCorrUp = 0 .f ;
819879 }
820880
821- std::array<float , 2 > trkltPosUp{mTrackletsCalib [trkltId].getY () - tiltCorrUp, zPosCorrUp};
881+ // conversion from slope in pad per time bin to slope in cm per BC = tracklets[trkltIdx].getSlopeFloat() * padWidth / BCperTimeBin
882+ float slopeFactor = mTrackletsRaw [trkltId].getSlopeFloat () * pad->getWidthIPad () / 4 .f ;
883+ float yCorrPileUp = tCorrPileUp * slopeFactor;
884+ float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
885+
886+ int nTrackletsChamber = mTracker ->GetNtrackletsChamber (trk.getCollisionId (), trkltDet);
887+ float angularPull = (mTrackletsCalib [trkltId].getDy () + dyTiltCorr - mRecoParam .convertAngleToDy (trkParam->getSnp ())) / std::sqrt (mRecoParam .getDyRes (trkParam->getSnp (), nTrackletsChamber));
888+
889+ std::array<float , 2 > trkltPosUp{mTrackletsCalib [trkltId].getY () - tiltCorrUp + yCorrPileUp, zPosCorrUp};
822890 std::array<float , 3 > trkltCovUp;
823- mRecoParam .recalcTrkltCov (tilt, trkParam->getSnp (), pad->getRowSize (mTrackletsRaw [trkltId].getPadRow ()), trkltCovUp);
891+ mRecoParam .recalcTrkltCov (tilt, trkParam->getSnp (), pad->getRowSize (mTrackletsRaw [trkltId].getPadRow ()), trkltCovUp, (mRec ->GetParam ().rec .trd .useAngularPull != 0 ? angularPull : 0 .), nTrackletsChamber);
892+ trkltCovUp[0 ] += yAddErrPileUp2;
824893
825894 chi2 += trkParam->getPredictedChi2 (trkltPosUp, trkltCovUp);
826895 if (!trkParam->update (trkltPosUp, trkltCovUp)) {
0 commit comments