@@ -39,6 +39,21 @@ using namespace o2::tpc;
3939using GTrackID = o2::dataformats::GlobalTrackID;
4040using DetID = o2::detectors::DetID;
4141
42+ float UnbinnedResid::getAlpha () const
43+ {
44+ return row < 180 ? o2::math_utils::sector2Angle (sec) : getITSAlpha ();
45+ }
46+
47+ float UnbinnedResid::getTOFX () const
48+ {
49+ int det[5 ];
50+ o2::tof::Geo::getVolumeIndices (channel, det);
51+ float pos[3 ];
52+ o2::tof::Geo::getPos (det, pos);
53+ o2::tof::Geo::rotateToSector (pos, c.getSector ());
54+ return pos[2 ]; // coordinates in sector frame: note that the rotation above puts z in pos[1], the radial coordinate in pos[2], and the tangent coordinate in pos[0] (this is to match the TOF residual system, where we don't use the radial component), so we swap their positions.
55+ }
56+
4257void TrackInterpolation::init (o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap)
4358{
4459 // perform initialization
@@ -386,6 +401,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
386401 trackData.gid = mGIDs [iSeed];
387402 trackData.par = mSeeds [iSeed];
388403 auto & trkWork = mSeeds [iSeed];
404+ auto trkInner = trkWork;
389405 // reset the cache array (sufficient to set cluster available to zero)
390406 for (auto & elem : mCache ) {
391407 elem.clAvailable = 0 ;
@@ -468,10 +484,10 @@ void TrackInterpolation::interpolateTrack(int iSeed)
468484 std::array<float , 2 > trkltTRDYZ{};
469485 std::array<float , 3 > trkltTRDCov{};
470486 int res = processTRDLayer (trkTRD, iLayer, trkWork, &trkltTRDYZ, &trkltTRDCov);
471- if (res == 0 ) { // no TRD tracklet in this layer
487+ if (res == - 1 ) { // no TRD tracklet in this layer
472488 continue ;
473489 }
474- if (res < 0 ) { // failed to reach this layer
490+ if (res < - 1 ) { // failed to reach this layer
475491 return ;
476492 }
477493 if (!trkWork.update (trkltTRDYZ, trkltTRDCov)) {
@@ -484,6 +500,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
484500 if (mDumpTrackPoints ) {
485501 (*trackDataExtended).trkOuter = trkWork;
486502 }
503+ auto trkOuter = trkWork; // outer param
487504
488505 // go back through the TPC and store updated track positions
489506 bool outerParamStored = false ;
@@ -576,13 +593,94 @@ void TrackInterpolation::interpolateTrack(int iSeed)
576593 }
577594 }
578595 trackData.clIdx .setEntries (nClValidated);
596+
597+ // do we have TRD residuals to add?
598+ bool stopPropagation = false ;
599+ trkWork = trkOuter;
600+ if (gidTable[GTrackID::TRD].isIndexSet ()) {
601+ const auto & trkTRD = mRecoCont ->getITSTPCTRDTrack <o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]);
602+ for (int iLayer = 0 ; iLayer < o2::trd::constants::NLAYER; iLayer++) {
603+ std::array<float , 2 > trkltTRDYZ{};
604+ int res = processTRDLayer (trkTRD, iLayer, trkWork, &trkltTRDYZ);
605+ if (res == -1 ) { // no traklet on this layer
606+ continue ;
607+ }
608+ if (res < -1 ) { // failed to reach this layer
609+ stopPropagation = true ;
610+ break ;
611+ }
612+
613+ float tgPhi = trkWork.getSnp () / std::sqrt ((1 .f - trkWork.getSnp ()) * (1 .f + trkWork.getSnp ()));
614+ auto dy = trkltTRDYZ[0 ] - trkWork.getY ();
615+ auto dz = trkltTRDYZ[1 ] - trkWork.getZ ();
616+ const auto sec = clusterResiduals[iCl].sec ;
617+ if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWork.getY ()) < param::MaxY) && (std::abs (trkWork.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
618+ mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 160 + iLayer, o2::math_utils::angle2Sector (trkWork.getAlpha ()), (short )res);
619+ trackData.nTrkltsTRD ++;
620+ trackData.nExtDetResid ++;
621+ }
622+ }
623+ }
624+
625+ // do we have TOF residual to add?
626+ while (gidTable[GTrackID::TOF].isIndexSet () && !stopPropagation) {
627+ const auto & clTOF = mRecoCont ->getTOFClusters ()[gidTable[GTrackID::TOF]];
628+ const float clTOFAlpha = o2::math_utils::sector2Angle (clTOF.getCount ());
629+ if (trkWork.getAlpha () != clTOFAlpha && !trkWork.rotate (clTOFAlpha)) {
630+ LOG (debug) << " Failed to rotate into TOF cluster sector frame" ;
631+ stopPropagation = true ;
632+ break ;
633+ }
634+ if (!propagator->PropagateToXBxByBz (trkWork, clTOF.getX (), mParams ->maxSnp , mParams ->maxStep , mMatCorr )) {
635+ LOG (debug) << " Failed final propagation to TOF radius" ;
636+ break ;
637+ }
638+
639+ float tgPhi = trkWork.getSnp () / std::sqrt ((1 .f - trkWork.getSnp ()) * (1 .f + trkWork.getSnp ()));
640+ auto dy = clTOF.getY () - trkWork.getY ();
641+ auto dz = clTOF.getZ () - trkWork.getZ ();
642+ if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWork.getY ()) < param::MaxY) && (std::abs (trkWork.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
643+ mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 170 , clTOF.getCount (), clTOF.getPadInSector ());
644+ trackData.nExtDetResid ++;
645+ }
646+ break ;
647+ }
648+
649+ // add ITS residuals
650+ while (!stopPropagation) {
651+ trkWork = trkInner; // this is ITS outer param
652+ auto nClITS = trkITS.getNumberOfClusters ();
653+ auto clEntryITS = trkITS.getFirstClusterEntry ();
654+ o2::its::GeometryTGeo* geom = GeometryTGeo::Instance ();
655+ for (int iCl = 0 ; iCl < nCl; iCl++) { // clusters are stored from outer to inner layers
656+ const auto & clsITS = mITSClustersArray [mITSTrackClusIdx [clEntry + iCl]];
657+ int chip = clsITS.getSensorID ();
658+ float chipX, chipAlpha;
659+ geom->getSensorXAlphaRefPlane (clsITS.getSensorID (), chipX, chipAlpha);
660+ if (!trkWorkITS.rotate (chipAlpha) || !propagator->PropagateToXBxByBz (trkWorkITS, chipX, mParams ->maxSnp , mParams ->maxStep , mMatCorr )) {
661+ LOGP (debug, " Failed final propagation to ITS X={} alpha={}" , chipX, chipAlpha);
662+ stopPropagation = true ;
663+ break ;
664+ }
665+ float tgPhi = trkWorkITS.getSnp () / std::sqrt ((1 .f - trkWorkITS.getSnp ()) * (1 .f + trkWorkITS.getSnp ()));
666+ auto dy = clsITS.getY () - trkWorkITS.getY ();
667+ auto dz = clsITS.getZ () - trkWorkITS.getZ ();
668+ if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWorkITS.getY ()) < param::MaxY) && (std::abs (trkWorkITS.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
669+ mClRes .emplace_back (dy, dz, tgPhi, trkWorkITS.getY (), trkWorkITS.getZ (), 180 + geom->getLayer (clsITS.getSensorID ()), -1 , clsITS.getSensorID ());
670+ trackData.nExtDetResid ++;
671+ }
672+ }
673+ break ;
674+ }
675+
579676 mTrackData .push_back (std::move (trackData));
580677 if (mDumpTrackPoints ) {
581678 (*trackDataExtended).clIdx .setEntries (nClValidated);
679+ (*trackDataExtended).nExtDetResid = trackData.nExtDetResid ;
582680 mTrackDataExtended .push_back (std::move (*trackDataExtended));
583681 }
584682 mGIDsSuccess .push_back (mGIDs [iSeed]);
585- mTrackDataCompact .emplace_back (mClRes .size () - nClValidated, nClValidated, mGIDs [iSeed].getSource ());
683+ mTrackDataCompact .emplace_back (mClRes .size () - nClValidated, nClValidated, mGIDs [iSeed].getSource (), trackData. nExtDetResid );
586684 }
587685 if (mParams ->writeUnfiltered ) {
588686 TrackData trkDataTmp = trackData;
@@ -596,10 +694,10 @@ void TrackInterpolation::interpolateTrack(int iSeed)
596694int TrackInterpolation::processTRDLayer (const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TracParCov& trkWork,
597695 std::array<float , 2 >* trkltTRDYZ, std::array<float , 3 >* trkltTRDCov)
598696{
599- // return 1 in case of successful processing, 0 if there is no TRD tracklet at given layer, -1 if processing failed
697+ // return chamber ID (0:539) in case of successful processing, -1 if there is no TRD tracklet at given layer, -2 if processing failed
600698 int trkltIdx = trkTRD.getTrackletIndex (iLayer);
601699 if (trkltIdx < 0 ) {
602- return 0 ; // no TRD tracklet in this layer
700+ return - 1 ; // no TRD tracklet in this layer
603701 }
604702 const auto & trdSP = mRecoCont ->getTRDCalibratedTracklets ()[trkltIdx];
605703 const auto & trdTrklt = mRecoCont ->getTRDTracklets ()[trkltIdx];
@@ -608,12 +706,12 @@ int TrackInterpolation::processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLa
608706 if (trkltSec != o2::math_utils::angle2Sector (trkWork.getAlpha ())) {
609707 if (!trkWork.rotate (o2::math_utils::sector2Angle (trkltSec))) {
610708 LOG (debug) << " Track could not be rotated in TRD tracklet coordinate system in layer " << iLayer;
611- return -1 ;
709+ return -2 ;
612710 }
613711 }
614712 if (!propagator->PropagateToXBxByBz (trkWork, trdSP.getX (), mParams ->maxSnp , mParams ->maxStep , mMatCorr )) {
615713 LOG (debug) << " Failed propagation to TRD layer " << iLayer;
616- return -1 ;
714+ return -2 ;
617715 }
618716 if (trkltTRDYZ) {
619717 const auto * pad = mGeoTRD ->getPadPlane (trkltDet);
@@ -630,7 +728,7 @@ int TrackInterpolation::processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLa
630728 mRecoParam .recalcTrkltCov (tilt, trkWork.getSnp (), pad->getRowSize (trdTrklt.getPadRow ()), *trkltTRDCov);
631729 }
632730 }
633- return 1 ;
731+ return trkltDet ;
634732}
635733
636734void TrackInterpolation::extrapolateTrack (int iSeed)
@@ -755,10 +853,10 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
755853 for (int iLayer = 0 ; iLayer < o2::trd::constants::NLAYER; iLayer++) {
756854 std::array<float , 2 > trkltTRDYZ{};
757855 int res = processTRDLayer (trkTRD, iLayer, trkWork, &trkltTRDYZ);
758- if (res == 0 ) { // no traklet on this layer
856+ if (res == - 1 ) { // no traklet on this layer
759857 continue ;
760858 }
761- if (res < 0 ) { // failed to reach this layer
859+ if (res < - 1 ) { // failed to reach this layer
762860 stopPropagation = true ;
763861 break ;
764862 }
@@ -768,10 +866,9 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
768866 auto dz = trkltTRDYZ[1 ] - trkWork.getZ ();
769867 const auto sec = clusterResiduals[iCl].sec ;
770868 if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWork.getY ()) < param::MaxY) && (std::abs (trkWork.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
771- mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 160 + iLayer, o2::math_utils::angle2Sector (trkWork.getAlpha ()));
869+ mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 160 + iLayer, o2::math_utils::angle2Sector (trkWork.getAlpha ()), ( short )res );
772870 trackData.nTrkltsTRD ++;
773871 trackData.nExtDetResid ++;
774- // NEED TO STORE: TRD X
775872 }
776873 }
777874 }
@@ -794,10 +891,9 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
794891 auto dy = clTOF.getY () - trkWork.getY ();
795892 auto dz = clTOF.getZ () - trkWork.getZ ();
796893 if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWork.getY ()) < param::MaxY) && (std::abs (trkWork.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
797- mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 170 , clTOF.getCount ());
894+ mClRes .emplace_back (dy, dz, tgPhi, trkWork.getY (), trkWork.getZ (), 170 , clTOF.getCount (), clTOF. getPadInSector () );
798895 trackData.clAvailTOF = 1 ;
799896 trackData.nExtDetResid ++;
800- // NEED TO STORE: TOF X
801897 }
802898 break ;
803899 }
@@ -822,18 +918,19 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
822918 auto dy = clsITS.getY () - trkWorkITS.getY ();
823919 auto dz = clsITS.getZ () - trkWorkITS.getZ ();
824920 if ((std::abs (dy) < param::MaxResid) && (std::abs (dz) < param::MaxResid) && (std::abs (trkWorkITS.getY ()) < param::MaxY) && (std::abs (trkWorkITS.getZ ()) < param::MaxZ) && (std::abs (tgPhi) < param::MaxTgSlp)) {
825- mClRes .emplace_back (dy, dz, tgPhi, trkWorkITS.getY (), trkWorkITS.getZ (), 180 + geom->getLayer (clsITS.getSensorID ()), 0 );
921+ mClRes .emplace_back (dy, dz, tgPhi, trkWorkITS.getY (), trkWorkITS.getZ (), 180 + geom->getLayer (clsITS.getSensorID ()), - 1 , clsITS. getSensorID () );
826922 trackData.nExtDetResid ++;
827- // NEED TO STORE: ITS X, alpha
828923 }
829924 }
830925 break ;
831926 }
927+
832928 mTrackData .push_back (std::move (trackData));
833929 mGIDsSuccess .push_back (mGIDs [iSeed]);
834- mTrackDataCompact .emplace_back (mClRes .size () - nClValidated, nClValidated, mGIDs [iSeed].getSource ()); // RS Check
930+ mTrackDataCompact .emplace_back (mClRes .size () - nClValidated, nClValidated, mGIDs [iSeed].getSource (), trackData. nExtDetResid );
835931 if (mDumpTrackPoints ) {
836932 (*trackDataExtended).clIdx .setEntries (nClValidated);
933+ (*trackDataExtended).nExtDetResid = trackData.nExtDetResid ;
837934 mTrackDataExtended .push_back (std::move (*trackDataExtended));
838935 }
839936 }
0 commit comments