[PWGDQ] add code to emulate a shift of the dipole

PWGDQ/Tasks/muonGlobalAlignment.cxx

@@ -150,6 +150,7 @@
   Configurable<uint32_t> fMftTracksMultiplicityMax{"cfgMftTracksMultiplicityMax", 0, "Maximum number of MFT tracks to be processed per event (zero means no limit)"};
 
   Configurable<float> fVertexZshift{"cfgVertexZshift", 0.0f, "Correction to the vertex z position"};
+  Configurable<float> fDipoleZshift{"cfgDipoleZshift", 0.0f, "Correction to the dipole z position"};
 
   ////   Variables for MFT alignment corrections
   struct : ConfigurableGroup {
@@ -272,7 +273,7 @@
                       std::map<uint64_t, CollisionInfo>& collisionInfos)
   {
     // fill collision information for global muon tracks (MFT-MCH-MID matches)
     for (auto muonTrack : muonTracks) {
       if (!muonTrack.has_collision())
         continue;
 
@@ -320,8 +321,8 @@
       return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT());
     };
 
     for (auto& [collisionIndex, collisionInfo] : collisionInfos) {
       for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) {
         std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2);
       }
     }
@@ -336,7 +337,7 @@
     InitCollisions(collisions, bcs, muonTracks, collisionInfos);
 
     // fill collision information for MFT standalone tracks
     for (auto mftTrack : mftTracks) {
       if (!mftTrack.has_collision())
         continue;
 
@@ -1094,8 +1095,8 @@
     track.setBendingCoor(y + yCorrection);
     track.setBendingSlope(ySlope + ySlopeCorrection);
     /*
     std::cout << std::format("[TOTO] MFT position:    pos={:0.3f},{:0.3f}", x, y) << std::endl;
     std::cout << std::format("[TOTO] MFT corrections: pos={:0.3f},{:0.3f}  slope={:0.12f},{:0.12f}  angle={:0.12f},{:0.12f}",
         xCorrection, yCorrection, xSlopeCorrection, ySlopeCorrection,
         std::atan2(xSlopeCorrection, 1), std::atan2(ySlopeCorrection, 1)) << std::endl;
     */
@@ -1133,8 +1134,8 @@
   template <typename T>
   T UpdateTrackMomentum(const T& track, const double p, int sign)
   {
     double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
     double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
     double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
 
     SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
@@ -1154,8 +1155,8 @@
   template <typename T>
   T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom)
   {
     double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
     double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
     double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
     double sign = track4mom.getCharge();
 
@@ -1200,6 +1201,15 @@
     if (mchTrack.getZ() < absBack && z > absFront) {
       // extrapolation through the absorber in the upstream direction
       o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, z);
+    } else if (z < absBack) {
+      // extrapolation downstream of the absorber, correct for dipole longitudinal shift if needed
+      if (fDipoleZshift.value != 0) {
+        mchTrack.setZ(mchTrack.getZ() + fDipoleZshift.value);
+        o2::mch::TrackExtrap::extrapToZCov(mchTrack, z + fDipoleZshift.value);
+        mchTrack.setZ(mchTrack.getZ() - fDipoleZshift.value);
+      } else {
+        o2::mch::TrackExtrap::extrapToZCov(mchTrack, z);
+      }
     } else {
       // all other cases
       o2::mch::TrackExtrap::extrapToZCov(mchTrack, z);
@@ -1389,7 +1399,24 @@
       o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -466.f);
       UpdateTrackMomentum(mftTrackProp, mchTrackPar);
     }
-    o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z);
+
+    if (fDipoleZshift.value != 0) {
+      // extrapolate to the back of the absorber, taking into account the dipole shift,
+      // to avoid that the correction bring the track starting point back into the absorber
+      if (fDipoleZshift.value < 0) {
+        o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -505.f);
+      } else if (fDipoleZshift.value > 0) {
+        o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -505.f - fDipoleZshift.value);
+      }
+      // shift the track starting point
+      mftTrackProp.setZ(mftTrackProp.getZ() + fDipoleZshift.value);
+      // extrapolate to the final z, corrected for the dipole shift
+      o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z + fDipoleZshift.value);
+      // remove the shift from the extrapolated track
+      mftTrackProp.setZ(mftTrackProp.getZ() - fDipoleZshift.value);
+    } else {
+      o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z);
+    }
 
     return MCHtoFwd(mftTrackProp);
   }
@@ -1595,6 +1622,11 @@
         transformNew[cluster.deId()].LocalToMaster(local, master);
       }
 
+      // shift the clusters to correct the longitudinal shift of the dipole
+      if (fDipoleZshift.value != 0) {
+        master.SetZ(master.z() + fDipoleZshift.value);
+      }
+
       if (applyCorrections) {
         auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId());
         if (correctionsIt != mMchAlignmentCorrections.end()) {
@@ -1628,6 +1660,12 @@
       LOGF(fatal, "Muon track %d has no associated clusters.", mchTrack.globalIndex());
     }
 
+    // subtract the longitudinal shift of the dipole from the track z
+    if (fDipoleZshift.value != 0) {
+      auto& trackParam = *(convertedTrack.begin());
+      trackParam.setZ(trackParam.getZ() - fDipoleZshift.value);
+    }
+
     return !removable;
   }