Updating MUON code with the new class AliMUONTrackExtrap (Philippe Pillot)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 27 Oct 2006 07:59:51 +0000 (07:59 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 27 Oct 2006 07:59:51 +0000 (07:59 +0000)
MUON/AliMUONReconstructor.cxx
MUON/AliMUONTrackK.cxx
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h
MUON/AliMUONTrackReconstructor.cxx
MUON/AliMUONVTrackReconstructor.cxx
MUON/MUONCheck.C
MUON/MUONefficiency.C
MUON/MUONmassPlot_ESD.C
MUON/MUONrecLinkDef.h
MUON/libMUONrec.pkg

index 913d0e3..9be5a6e 100644 (file)
@@ -14,9 +14,9 @@
  **************************************************************************/
 /* $Id$ */
 
-//-----------------------------    
-// Class AliMUONReconstructor                                                                      //
-//-----------------------------                                                                           //
+//-----------------------------
+// Class AliMUONReconstructor  
+//-----------------------------
 // Class for the 
 // MUON track reconstruction
 
@@ -36,6 +36,7 @@
 #include "AliMUONDigitMaker.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliMUONVTrackReconstructor.h"
 #include "AliMUONTrackReconstructor.h"
 #include "AliMUONTrackReconstructorK.h"
@@ -483,7 +484,7 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
    
     if (esdVert->GetNContributors())
-      trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
+      AliMUONTrackExtrap::ExtrapToVertex(trackParam, vertex[0],vertex[1],vertex[2]);
 
     bendingSlope            = trackParam->GetBendingSlope();
     nonBendingSlope         = trackParam->GetNonBendingSlope();
index aa4603a..ce588bc 100644 (file)
@@ -37,6 +37,7 @@
 #include "AliMUONHitForRec.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliMUONEventRecoCombi.h"
 #include "AliMUONDetElement.h"
 #include "AliRun.h"
@@ -696,7 +697,7 @@ void AliMUONTrackK::ParPropagation(Double_t zEnd)
   do {
     step = TMath::Abs(step);
     // Propagate parameters
-    trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+    AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
     //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
     distance = zEnd - vGeant3New[2];
     step *= dZ/(vGeant3New[2]-fPositionNew);
@@ -716,7 +717,7 @@ void AliMUONTrackK::ParPropagation(Double_t zEnd)
     do {
       // binary search
       // Propagate parameters
-      trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+      AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
       //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
       distance = zEnd - vGeant3New[2];
       step /= 2;
@@ -1291,7 +1292,7 @@ void AliMUONTrackK::FillMUONTrack(void)
   for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
     if ((*fChi2Smooth)[i] < 0) {
       // Propagate through last chambers
-      trackParam.ExtrapToZ(((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
+      AliMUONTrackExtrap::ExtrapToZ(&trackParam, ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
     } else {
       // Take saved info
       SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
@@ -1333,8 +1334,7 @@ void AliMUONTrackK::Branson(void)
   */
   SetTrackParam(&trackParam, fTrackPar, fPosition);
 
-  trackParam.ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
-  //trackParam.ExtrapToVertex();
+  AliMUONTrackExtrap::ExtrapToVertex(&trackParam, Double_t(0.), Double_t(0.), Double_t(0.));
 
   (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
   (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
index 8ca824a..d92b72c 100644 (file)
 //
 ///////////////////////////////////////////////////
 
-//#include <Riostream.h>
-#include "AliMUON.h"
-#include "AliMUONTrackParam.h" 
-#include "AliMUONConstants.h"
+#include <Riostream.h>
+
+#include "AliMUONTrackParam.h"
 #include "AliESDMuonTrack.h"
-#include "AliMagF.h" 
-#include "AliLog.h" 
-#include "AliTracker.h"
+#include "AliLog.h"
 #include "AliMUONHitForRec.h"
 
 /// \cond CLASSIMP
@@ -48,13 +45,9 @@ AliMUONTrackParam::AliMUONTrackParam()
     fZ(0.),
     fBendingCoor(0.),
     fNonBendingCoor(0.),
-    fkField(0x0),
     fHitForRecPtr(0x0)
 {
-/// Constructor
-  // get field from outside
-  fkField = AliTracker::GetFieldMap();
-  if (!fkField) AliFatal("No field available");
+  /// Constructor
 }
 
   //_________________________________________________________________________
@@ -66,7 +59,6 @@ AliMUONTrackParam::AliMUONTrackParam(const AliMUONTrackParam& theMUONTrackParam)
     fZ(theMUONTrackParam.fZ),
     fBendingCoor(theMUONTrackParam.fBendingCoor),
     fNonBendingCoor(theMUONTrackParam.fNonBendingCoor),
-    fkField(theMUONTrackParam.fkField),
     fHitForRecPtr(theMUONTrackParam.fHitForRecPtr)
 {
   /// Copy constructor
@@ -88,7 +80,6 @@ AliMUONTrackParam& AliMUONTrackParam::operator=(const AliMUONTrackParam& theMUON
   fZ                      =  theMUONTrackParam.fZ; 
   fBendingCoor            =  theMUONTrackParam.fBendingCoor; 
   fNonBendingCoor         =  theMUONTrackParam.fNonBendingCoor;
-  fkField                 =  theMUONTrackParam.fkField;
   fHitForRecPtr           =  theMUONTrackParam.fHitForRecPtr;
 
   return *this;
@@ -160,417 +151,6 @@ void AliMUONTrackParam::SetParamFor(AliESDMuonTrack& esdMuonTrack)
 }
 
   //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapToZ(Double_t Z)
-{
-  /// Track parameter extrapolation to the plane at "Z".
-  /// On return, the track parameters resulting from the extrapolation
-  /// replace the current track parameters.
-  if (this->fZ == Z) return; // nothing to be done if same Z
-  Double_t forwardBackward; // +1 if forward, -1 if backward
-  if (Z < this->fZ) forwardBackward = 1.0; // spectro. z<0 
-  else forwardBackward = -1.0;
-  Double_t vGeant3[7], vGeant3New[7]; // 7 in parameter ????
-  Int_t iGeant3, stepNumber;
-  Int_t maxStepNumber = 5000; // in parameter ????
-  // For safety: return kTRUE or kFALSE ????
-  // Parameter vector for calling EXTRAP_ONESTEP
-  SetGeant3Parameters(vGeant3, forwardBackward);
-  // sign of charge (sign of fInverseBendingMomentum if forward motion)
-  // must be changed if backward extrapolation
-  Double_t chargeExtrap = forwardBackward *
-    TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
-  Double_t stepLength = 6.0; // in parameter ????
-  // Extrapolation loop
-  stepNumber = 0;
-  while (((-forwardBackward * (vGeant3[2] - Z)) <= 0.0) &&  // spectro. z<0
-        (stepNumber < maxStepNumber)) {
-    stepNumber++;
-    // Option for switching between helix and Runge-Kutta ???? 
-    //ExtrapOneStepRungekutta(chargeExtrap, stepLength, vGeant3, vGeant3New);
-    ExtrapOneStepHelix(chargeExtrap, stepLength, vGeant3, vGeant3New);
-    if ((-forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z spectro. z<0
-    // better use TArray ????
-    for (iGeant3 = 0; iGeant3 < 7; iGeant3++)
-      {vGeant3[iGeant3] = vGeant3New[iGeant3];}
-  }
-  // check maxStepNumber ????
-  // Interpolation back to exact Z (2nd order)
-  // should be in function ???? using TArray ????
-  Double_t dZ12 = vGeant3New[2] - vGeant3[2]; // 1->2
-  if (TMath::Abs(dZ12) > 0) {
-    Double_t dZ1i = Z - vGeant3[2]; // 1-i
-    Double_t dZi2 = vGeant3New[2] - Z; // i->2
-    Double_t xPrime = (vGeant3New[0] - vGeant3[0]) / dZ12;
-    Double_t xSecond =
-      ((vGeant3New[3] / vGeant3New[5]) - (vGeant3[3] / vGeant3[5])) / dZ12;
-    Double_t yPrime = (vGeant3New[1] - vGeant3[1]) / dZ12;
-    Double_t ySecond =
-      ((vGeant3New[4] / vGeant3New[5]) - (vGeant3[4] / vGeant3[5])) / dZ12;
-    vGeant3[0] = vGeant3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
-    vGeant3[1] = vGeant3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
-    vGeant3[2] = Z; // Z
-    Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
-    Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
-    // (PX, PY, PZ)/PTOT assuming forward motion
-    vGeant3[5] =
-      1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
-    vGeant3[3] = xPrimeI * vGeant3[5]; // PX/PTOT
-    vGeant3[4] = yPrimeI * vGeant3[5]; // PY/PTOT
-  } else {
-    AliWarning(Form("Extrap. to Z not reached, Z = %f",Z));    
-  }
-  // Track parameters from Geant3 parameters,
-  // with charge back for forward motion
-  GetFromGeant3Parameters(vGeant3, chargeExtrap * forwardBackward);
-}
-
-  //__________________________________________________________________________
-void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward)
-{
-  /// Set vector of Geant3 parameters pointed to by "VGeant3"
-  /// from track parameters in current AliMUONTrackParam.
-  /// Since AliMUONTrackParam is only geometry, one uses "ForwardBackward"
-  /// to know whether the particle is going forward (+1) or backward (-1).
-  VGeant3[0] = this->fNonBendingCoor; // X
-  VGeant3[1] = this->fBendingCoor; // Y
-  VGeant3[2] = this->fZ; // Z
-  Double_t pYZ = TMath::Abs(1.0 / this->fInverseBendingMomentum);
-  Double_t pZ =
-    pYZ / TMath::Sqrt(1.0 + this->fBendingSlope * this->fBendingSlope);
-  VGeant3[6] =
-    TMath::Sqrt(pYZ * pYZ +
-               pZ * pZ * this->fNonBendingSlope * this->fNonBendingSlope); // PTOT
-  VGeant3[5] = -ForwardBackward * pZ / VGeant3[6]; // PZ/PTOT spectro. z<0
-  VGeant3[3] = this->fNonBendingSlope * VGeant3[5]; // PX/PTOT
-  VGeant3[4] = this->fBendingSlope * VGeant3[5]; // PY/PTOT
-}
-
-  //__________________________________________________________________________
-void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge)
-{
-  /// Get track parameters in current AliMUONTrackParam
-  /// from Geant3 parameters pointed to by "VGeant3",
-  /// assumed to be calculated for forward motion in Z.
-  /// "InverseBendingMomentum" is signed with "Charge".
-  this->fNonBendingCoor = VGeant3[0]; // X
-  this->fBendingCoor = VGeant3[1]; // Y
-  this->fZ = VGeant3[2]; // Z
-  Double_t pYZ = VGeant3[6] * TMath::Sqrt(1.0 - VGeant3[3] * VGeant3[3]);
-  this->fInverseBendingMomentum = Charge / pYZ;
-  this->fBendingSlope = VGeant3[4] / VGeant3[5];
-  this->fNonBendingSlope = VGeant3[3] / VGeant3[5];
-}
-
-  //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam)
-{
-  /// Track parameters extrapolated from current track parameters ("this")
-  /// to both chambers of the station(0..) "Station"
-  /// are returned in the array (dimension 2) of track parameters
-  /// pointed to by "TrackParam" (index 0 and 1 for first and second chambers).
-  Double_t extZ[2], z1, z2;
-  Int_t i1 = -1, i2 = -1; // = -1 to avoid compilation warnings
-  // range of Station to be checked ????
-  z1 = AliMUONConstants::DefaultChamberZ(2 * Station);
-  z2 = AliMUONConstants::DefaultChamberZ(2 * Station + 1);
-  // First and second Z to extrapolate at
-  if ((z1 > this->fZ) && (z2 > this->fZ)) {i1 = 0; i2 = 1;}
-  else if ((z1 < this->fZ) && (z2 < this->fZ)) {i1 = 1; i2 = 0;}
-  else {
-       AliError(Form("Starting Z (%f) in between z1 (%f) and z2 (%f) of station(0..)%d",this->fZ,z1,z2,Station));
-//     cout << "ERROR in AliMUONTrackParam::CreateExtrapSegmentInStation" << endl;
-//     cout << "Starting Z (" << this->fZ << ") in between z1 (" << z1 <<
-//       ") and z2 (" << z2 << ") of station(0..) " << Station << endl;
-  }
-  extZ[i1] = z1;
-  extZ[i2] = z2;
-  // copy of track parameters
-  TrackParam[i1] = *this;
-  // first extrapolation
-  (&(TrackParam[i1]))->ExtrapToZ(extZ[0]);
-  TrackParam[i2] = TrackParam[i1];
-  // second extrapolation
-  (&(TrackParam[i2]))->ExtrapToZ(extZ[1]);
-  return;
-}
-
-  //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapToVertex(Double_t xVtx, Double_t yVtx, Double_t zVtx)
-{
-  /// Extrapolation to the vertex.
-  /// Returns the track parameters resulting from the extrapolation in the current TrackParam.
-  /// Changes parameters according to Branson correction through the absorber 
-  
-  Double_t zAbsorber = -503.0; // to be coherent with the Geant absorber geometry !!!!
-                               // spectro. (z<0) 
-  // Extrapolates track parameters upstream to the "Z" end of the front absorber
-  ExtrapToZ(zAbsorber); // !!!
-  // Makes Branson correction (multiple scattering + energy loss)
-  BransonCorrection(xVtx,yVtx,zVtx);
-  // Makes a simple magnetic field correction through the absorber
-  FieldCorrection(zAbsorber);
-}
-
-
-//  Keep this version for future developments
-  //__________________________________________________________________________
-// void AliMUONTrackParam::BransonCorrection()
-// {
-//   // Branson correction of track parameters
-//   // the entry parameters have to be calculated at the end of the absorber
-//   Double_t zEndAbsorber, zBP, xBP, yBP;
-//   Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
-//   Int_t sign;
-//   // Would it be possible to calculate all that from Geant configuration ????
-//   // and to get the Branson parameters from a function in ABSO module ????
-//   // with an eventual contribution from other detectors like START ????
-//   // Radiation lengths outer part theta > 3 degres
-//   static Double_t x01[9] = { 18.8,    // C (cm)
-//                          10.397,   // Concrete (cm)
-//                          0.56,    // Plomb (cm)
-//                          47.26,   // Polyethylene (cm)
-//                          0.56,   // Plomb (cm)
-//                          47.26,   // Polyethylene (cm)
-//                          0.56,   // Plomb (cm)
-//                          47.26,   // Polyethylene (cm)
-//                          0.56 };   // Plomb (cm)
-//   // inner part theta < 3 degres
-//   static Double_t x02[3] = { 18.8,    // C (cm)
-//                          10.397,   // Concrete (cm)
-//                          0.35 };    // W (cm) 
-//   // z positions of the materials inside the absober outer part theta > 3 degres
-//   static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
-//   // inner part theta < 3 degres
-//   static Double_t z2[4] = { 90, 315, 467, 503 };
-//   static Bool_t first = kTRUE;
-//   static Double_t zBP1, zBP2, rLimit;
-//   // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
-//   if (first) {
-//     first = kFALSE;
-//     Double_t aNBP = 0.0;
-//     Double_t aDBP = 0.0;
-//     Int_t iBound;
-    
-//     for (iBound = 0; iBound < 9; iBound++) {
-//       aNBP = aNBP +
-//     (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
-//      z1[iBound]   * z1[iBound]   * z1[iBound]    ) / x01[iBound];
-//       aDBP = aDBP +
-//     (z1[iBound+1] * z1[iBound+1] - z1[iBound]   * z1[iBound]    ) / x01[iBound];
-//     }
-//     zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
-//     aNBP = 0.0;
-//     aDBP = 0.0;
-//     for (iBound = 0; iBound < 3; iBound++) {
-//       aNBP = aNBP +
-//     (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
-//      z2[iBound]   * z2[iBound ]  * z2[iBound]    ) / x02[iBound];
-//       aDBP = aDBP +
-//     (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
-//     }
-//     zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
-//     rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
-//   }
-
-//   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
-//   sign = 1;      
-//   if (fInverseBendingMomentum < 0) sign = -1;     
-//   pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); 
-//   pX = pZ * fNonBendingSlope; 
-//   pY = pZ * fBendingSlope; 
-//   pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
-//   xEndAbsorber = fNonBendingCoor; 
-//   yEndAbsorber = fBendingCoor; 
-//   radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
-
-//   if (radiusEndAbsorber2 > rLimit*rLimit) {
-//     zEndAbsorber = z1[9];
-//     zBP = zBP1;
-//   } else {
-//     zEndAbsorber = z2[3];
-//     zBP = zBP2;
-//   }
-
-//   xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
-//   yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
-
-//   // new parameters after Branson and energy loss corrections
-//   pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
-//   pX = pZ * xBP / zBP;
-//   pY = pZ * yBP / zBP;
-//   fBendingSlope = pY / pZ;
-//   fNonBendingSlope = pX / pZ;
-  
-//   pT = TMath::Sqrt(pX * pX + pY * pY);      
-//   theta = TMath::ATan2(pT, pZ); 
-//   pTotal =
-//     TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
-
-//   fInverseBendingMomentum = (sign / pTotal) *
-//     TMath::Sqrt(1.0 +
-//             fBendingSlope * fBendingSlope +
-//             fNonBendingSlope * fNonBendingSlope) /
-//     TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
-
-//   // vertex position at (0,0,0)
-//   // should be taken from vertex measurement ???
-//   fBendingCoor = 0.0;
-//   fNonBendingCoor = 0;
-//   fZ= 0;
-// }
-
-void AliMUONTrackParam::BransonCorrection(Double_t xVtx,Double_t yVtx,Double_t zVtx)
-{
-  /// Branson correction of track parameters
-  // the entry parameters have to be calculated at the end of the absorber
-  // simplified version: the z positions of Branson's planes are no longer calculated
-  // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C
-  // to test this correction. 
-  // Would it be possible to calculate all that from Geant configuration ????
-  // and to get the Branson parameters from a function in ABSO module ????
-  // with an eventual contribution from other detectors like START ????
-  // change to take into account the vertex postition (real, reconstruct,....)
-
-  Double_t  zBP, xBP, yBP;
-  Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
-  Int_t sign;
-  static Bool_t first = kTRUE;
-  static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber;
-  // zBP1 for outer part and zBP2 for inner part (only at the first call)
-  if (first) {
-    first = kFALSE;
-  
-    zEndAbsorber = -503;  // spectro (z<0)
-    thetaLimit = 3.0 * (TMath::Pi()) / 180.;
-    rLimit = TMath::Abs(zEndAbsorber) * TMath::Tan(thetaLimit);
-    zBP1 = -450; // values close to those calculated with EvalAbso.C
-    zBP2 = -480;
-  }
-
-  pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
-  sign = 1;      
-  if (fInverseBendingMomentum < 0) sign = -1;  
-  pZ = Pz();
-  pX = Px(); 
-  pY = Py(); 
-  pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
-  xEndAbsorber = fNonBendingCoor; 
-  yEndAbsorber = fBendingCoor; 
-  radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
-
-  if (radiusEndAbsorber2 > rLimit*rLimit) {
-    zBP = zBP1;
-  } else {
-    zBP = zBP2;
-  }
-
-  xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
-  yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
-
-  // new parameters after Branson and energy loss corrections
-//   Float_t zSmear = zBP - gRandom->Gaus(0.,2.);  // !!! possible smearing of Z vertex position
-
-  Float_t zSmear = zBP ;
-  
-   pZ = pTotal * (zSmear-zVtx) / TMath::Sqrt((xBP-xVtx) * (xBP-xVtx) + (yBP-yVtx) * (yBP-yVtx) +( zSmear-zVtx) * (zSmear-zVtx) );
-   pX = pZ * (xBP - xVtx)/ (zSmear-zVtx);
-   pY = pZ * (yBP - yVtx) / (zSmear-zVtx);
-  fBendingSlope = pY / pZ;
-  fNonBendingSlope = pX / pZ;
-
-  
-  pT = TMath::Sqrt(pX * pX + pY * pY);      
-  theta = TMath::ATan2(pT, TMath::Abs(pZ)); 
-  pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta);
-
-  fInverseBendingMomentum = (sign / pTotal) *
-    TMath::Sqrt(1.0 +
-               fBendingSlope * fBendingSlope +
-               fNonBendingSlope * fNonBendingSlope) /
-    TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
-
-  // vertex position at (0,0,0)
-  // should be taken from vertex measurement ???
-
-  fBendingCoor = xVtx;
-  fNonBendingCoor = yVtx;
-  fZ= zVtx;
-
-}
-
-  //__________________________________________________________________________
-Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta)
-{
-  /// Returns the total momentum corrected from energy loss in the front absorber
-  // One can use the macros MUONTestAbso.C and DrawTestAbso.C
-  // to test this correction. 
-  // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002)
-  Double_t deltaP, pTotalCorrected;
-
-   // Parametrization to be redone according to change of absorber material ????
-  // See remark in function BransonCorrection !!!!
-  // The name is not so good, and there are many arguments !!!!
-  if (theta  < thetaLimit ) {
-    if (pTotal < 20) {
-      deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal;
-    } else {
-      deltaP = 3.0714 + 0.011767 *pTotal;
-    }
-    deltaP *= 0.75; // AZ
-  } else {
-    if (pTotal < 20) {
-      deltaP  = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal;
-    } else { 
-      deltaP = 2.6069 + 0.0051705 * pTotal;
-    }
-    deltaP *= 0.9; // AZ
-  }
-  pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
-  return pTotalCorrected;
-}
-
-  //__________________________________________________________________________
-void AliMUONTrackParam::FieldCorrection(Double_t Z)
-{
-  /// Correction of the effect of the magnetic field in the absorber
-  // Assume a constant field along Z axis.
-  Float_t b[3],x[3]; 
-  Double_t bZ;
-  Double_t pYZ,pX,pY,pZ,pT;
-  Double_t pXNew,pYNew;
-  Double_t c;
-
-  pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
-  c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge 
-  pZ = Pz();
-  pX = Px(); 
-  pY = Py();
-  pT = TMath::Sqrt(pX*pX+pY*pY);
-
-  if (TMath::Abs(pZ) <= 0) return;
-  x[2] = Z/2;
-  x[0] = x[2]*fNonBendingSlope;  
-  x[1] = x[2]*fBendingSlope;
-
-  // Take magn. field value at position x.
-  fkField->Field(x, b);
-  bZ =  b[2];
-  // Transverse momentum rotation
-  // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ.
-  Double_t phiShift = c*0.436*0.0003*bZ*Z/pZ; 
- // Rotate momentum around Z axis.
-  pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift);
-  pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift);
-  fBendingSlope = pYNew / pZ;
-  fNonBendingSlope = pXNew / pZ;
-  
-  fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ);
-}
-  //__________________________________________________________________________
 Double_t AliMUONTrackParam::Px() const
 {
   /// return px from track paramaters
@@ -582,6 +162,7 @@ Double_t AliMUONTrackParam::Px() const
   pX = pZ * fNonBendingSlope; 
   return pX;
 }
+
   //__________________________________________________________________________
 Double_t AliMUONTrackParam::Py() const
 {
@@ -594,6 +175,7 @@ Double_t AliMUONTrackParam::Py() const
   pY = pZ * fBendingSlope; 
   return pY;
 }
+
   //__________________________________________________________________________
 Double_t AliMUONTrackParam::Pz() const
 {
@@ -605,6 +187,7 @@ Double_t AliMUONTrackParam::Pz() const
   pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
   return pZ;
 }
+
   //__________________________________________________________________________
 Double_t AliMUONTrackParam::P() const
 {
@@ -619,457 +202,7 @@ Double_t AliMUONTrackParam::P() const
   return p;
   
 }
- //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapOneStepHelix(Double_t charge, Double_t step, 
-                                        Double_t *vect, Double_t *vout) const
-{
-///    ******************************************************************
-///    *                                                                *
-///    *  Performs the tracking of one step in a magnetic field         *
-///    *  The trajectory is assumed to be a helix in a constant field   *
-///    *  taken at the mid point of the step.                           *
-///    *  Parameters:                                                   *
-///    *   input                                                        *
-///    *     STEP =arc length of the step asked                         *
-///    *     VECT =input vector (position,direction cos and momentum)   *
-///    *     CHARGE=  electric charge of the particle                   *
-///    *   output                                                       *
-///    *     VOUT = same as VECT after completion of the step           *
-///    *                                                                *
-///    *    ==>Called by : <USER>, GUSWIM                               *
-///    *       Author    m.hansroul  *********                          *
-///    *       modified  s.egli, s.v.levonian                           *
-///    *       modified  v.perevoztchikov
-///    *                                                                *
-///    ******************************************************************
-
-// modif: everything in double precision
-
-    Double_t xyz[3], h[4], hxp[3];
-    Double_t h2xy, hp, rho, tet;
-    Double_t sint, sintt, tsint, cos1t;
-    Double_t f1, f2, f3, f4, f5, f6;
-
-    const Int_t kix  = 0;
-    const Int_t kiy  = 1;
-    const Int_t kiz  = 2;
-    const Int_t kipx = 3;
-    const Int_t kipy = 4;
-    const Int_t kipz = 5;
-    const Int_t kipp = 6;
-
-    const Double_t kec = 2.9979251e-4;
-    //
-    //    ------------------------------------------------------------------
-    //
-    //       units are kgauss,centimeters,gev/c
-    //
-    vout[kipp] = vect[kipp];
-    if (TMath::Abs(charge) < 0.00001) {
-      for (Int_t i = 0; i < 3; i++) {
-       vout[i] = vect[i] + step * vect[i+3];
-       vout[i+3] = vect[i+3];
-      }
-      return;
-    }
-    xyz[0]    = vect[kix] + 0.5 * step * vect[kipx];
-    xyz[1]    = vect[kiy] + 0.5 * step * vect[kipy];
-    xyz[2]    = vect[kiz] + 0.5 * step * vect[kipz];
-
-    //cmodif: call gufld (xyz, h) changed into:
-    GetField (xyz, h);
-    h2xy = h[0]*h[0] + h[1]*h[1];
-    h[3] = h[2]*h[2]+ h2xy;
-    if (h[3] < 1.e-12) {
-      for (Int_t i = 0; i < 3; i++) {
-       vout[i] = vect[i] + step * vect[i+3];
-       vout[i+3] = vect[i+3];
-      }
-      return;
-    }
-    if (h2xy < 1.e-12*h[3]) {
-      ExtrapOneStepHelix3(charge*h[2], step, vect, vout);
-      return;
-    }
-    h[3] = TMath::Sqrt(h[3]);
-    h[0] /= h[3];
-    h[1] /= h[3];
-    h[2] /= h[3];
-    h[3] *= kec;
-
-    hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy];
-    hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz];
-    hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx];
-    hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz];
-
-    rho = -charge*h[3]/vect[kipp];
-    tet = rho * step;
-
-    if (TMath::Abs(tet) > 0.15) {
-      sint = TMath::Sin(tet);
-      sintt = (sint/tet);
-      tsint = (tet-sint)/tet;
-      cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet;
-    } else {
-      tsint = tet*tet/36.;
-      sintt = (1. - tsint);
-      sint = tet*sintt;
-      cos1t = 0.5*tet;
-    }
-
-    f1 = step * sintt;
-    f2 = step * cos1t;
-    f3 = step * tsint * hp;
-    f4 = -tet*cos1t;
-    f5 = sint;
-    f6 = tet * cos1t * hp;
-    vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0];
-    vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1];
-    vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2];
-    vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0];
-    vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1];
-    vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2];
-    return;
-}
-
- //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapOneStepHelix3(Double_t field, Double_t step, 
-                                              Double_t *vect, Double_t *vout) const
-{
-///    ******************************************************************
-///    *                                                                *
-///    *       Tracking routine in a constant field oriented            *
-///    *       along axis 3                                             *
-///    *       Tracking is performed with a conventional                *
-///    *       helix step method                                        *
-///    *                                                                *
-///    *    ==>Called by : <USER>, GUSWIM                               *
-///    *       Authors    R.Brun, M.Hansroul  *********                 *
-///    *       Rewritten  V.Perevoztchikov
-///    *                                                                *
-///    ******************************************************************
-
-    Double_t hxp[3];
-    Double_t h4, hp, rho, tet;
-    Double_t sint, sintt, tsint, cos1t;
-    Double_t f1, f2, f3, f4, f5, f6;
-
-    const Int_t kix  = 0;
-    const Int_t kiy  = 1;
-    const Int_t kiz  = 2;
-    const Int_t kipx = 3;
-    const Int_t kipy = 4;
-    const Int_t kipz = 5;
-    const Int_t kipp = 6;
-
-    const Double_t kec = 2.9979251e-4;
-
-// 
-//     ------------------------------------------------------------------
-// 
-//       units are kgauss,centimeters,gev/c
-// 
-    vout[kipp] = vect[kipp];
-    h4 = field * kec;
-
-    hxp[0] = - vect[kipy];
-    hxp[1] = + vect[kipx];
-    hp = vect[kipz];
-
-    rho = -h4/vect[kipp];
-    tet = rho * step;
-    if (TMath::Abs(tet) > 0.15) {
-      sint = TMath::Sin(tet);
-      sintt = (sint/tet);
-      tsint = (tet-sint)/tet;
-      cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
-    } else {
-      tsint = tet*tet/36.;
-      sintt = (1. - tsint);
-      sint = tet*sintt;
-      cos1t = 0.5*tet;
-    }
-
-    f1 = step * sintt;
-    f2 = step * cos1t;
-    f3 = step * tsint * hp;
-    f4 = -tet*cos1t;
-    f5 = sint;
-    f6 = tet * cos1t * hp;
-    vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
-    vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
-    vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
-    vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
-    vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
-    vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
-
-    return;
-}
- //__________________________________________________________________________
-void AliMUONTrackParam::ExtrapOneStepRungekutta(Double_t charge, Double_t step, 
-                                                    Double_t* vect, Double_t* vout) const
-{
-///    ******************************************************************
-///    *                                                                *
-///    *  Runge-Kutta method for tracking a particle through a magnetic *
-///    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
-///    *  Standards, procedure 25.5.20)                                 *
-///    *                                                                *
-///    *  Input parameters                                              *
-///    *       CHARGE    Particle charge                                *
-///    *       STEP      Step size                                      *
-///    *       VECT      Initial co-ords,direction cosines,momentum     *
-///    *  Output parameters                                             *
-///    *       VOUT      Output co-ords,direction cosines,momentum      *
-///    *  User routine called                                           *
-///    *       CALL GUFLD(X,F)                                          *
-///    *                                                                *
-///    *    ==>Called by : <USER>, GUSWIM                               *
-///    *       Authors    R.Brun, M.Hansroul  *********                 *
-///    *                  V.Perevoztchikov (CUT STEP implementation)    *
-///    *                                                                *
-///    *                                                                *
-///    ******************************************************************
-
-    Double_t h2, h4, f[4];
-    Double_t xyzt[3], a, b, c, ph,ph2;
-    Double_t secxs[4],secys[4],seczs[4],hxp[3];
-    Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
-    Double_t est, at, bt, ct, cba;
-    Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
-    
-    Double_t x;
-    Double_t y;
-    Double_t z;
-    
-    Double_t xt;
-    Double_t yt;
-    Double_t zt;
-
-    Double_t maxit = 1992;
-    Double_t maxcut = 11;
-
-    const Double_t kdlt   = 1e-4;
-    const Double_t kdlt32 = kdlt/32.;
-    const Double_t kthird = 1./3.;
-    const Double_t khalf  = 0.5;
-    const Double_t kec = 2.9979251e-4;
-
-    const Double_t kpisqua = 9.86960440109;
-    const Int_t kix  = 0;
-    const Int_t kiy  = 1;
-    const Int_t kiz  = 2;
-    const Int_t kipx = 3;
-    const Int_t kipy = 4;
-    const Int_t kipz = 5;
-  
-    // *.
-    // *.    ------------------------------------------------------------------
-    // *.
-    // *             this constant is for units cm,gev/c and kgauss
-    // *
-    Int_t iter = 0;
-    Int_t ncut = 0;
-    for(Int_t j = 0; j < 7; j++)
-      vout[j] = vect[j];
 
-    Double_t  pinv   = kec * charge / vect[6];
-    Double_t tl = 0.;
-    Double_t h = step;
-    Double_t rest;
-
-    do {
-      rest  = step - tl;
-      if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
-      //cmodif: call gufld(vout,f) changed into:
-
-      GetField(vout,f);
-
-      // *
-      // *             start of integration
-      // *
-      x      = vout[0];
-      y      = vout[1];
-      z      = vout[2];
-      a      = vout[3];
-      b      = vout[4];
-      c      = vout[5];
-
-      h2     = khalf * h;
-      h4     = khalf * h2;
-      ph     = pinv * h;
-      ph2    = khalf * ph;
-      secxs[0] = (b * f[2] - c * f[1]) * ph2;
-      secys[0] = (c * f[0] - a * f[2]) * ph2;
-      seczs[0] = (a * f[1] - b * f[0]) * ph2;
-      ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
-      if (ang2 > kpisqua) break;
-
-      dxt    = h2 * a + h4 * secxs[0];
-      dyt    = h2 * b + h4 * secys[0];
-      dzt    = h2 * c + h4 * seczs[0];
-      xt     = x + dxt;
-      yt     = y + dyt;
-      zt     = z + dzt;
-      // *
-      // *              second intermediate point
-      // *
-
-      est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
-      if (est > h) {
-       if (ncut++ > maxcut) break;
-       h *= khalf;
-       continue;
-      }
-      xyzt[0] = xt;
-      xyzt[1] = yt;
-      xyzt[2] = zt;
-
-      //cmodif: call gufld(xyzt,f) changed into:
-      GetField(xyzt,f);
-
-      at     = a + secxs[0];
-      bt     = b + secys[0];
-      ct     = c + seczs[0];
-
-      secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
-      secys[1] = (ct * f[0] - at * f[2]) * ph2;
-      seczs[1] = (at * f[1] - bt * f[0]) * ph2;
-      at     = a + secxs[1];
-      bt     = b + secys[1];
-      ct     = c + seczs[1];
-      secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
-      secys[2] = (ct * f[0] - at * f[2]) * ph2;
-      seczs[2] = (at * f[1] - bt * f[0]) * ph2;
-      dxt    = h * (a + secxs[2]);
-      dyt    = h * (b + secys[2]);
-      dzt    = h * (c + seczs[2]);
-      xt     = x + dxt;
-      yt     = y + dyt;
-      zt     = z + dzt;
-      at     = a + 2.*secxs[2];
-      bt     = b + 2.*secys[2];
-      ct     = c + 2.*seczs[2];
-
-      est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
-      if (est > 2.*TMath::Abs(h)) {
-       if (ncut++ > maxcut) break;
-       h *= khalf;
-       continue;
-      }
-      xyzt[0] = xt;
-      xyzt[1] = yt;
-      xyzt[2] = zt;
-
-      //cmodif: call gufld(xyzt,f) changed into:
-      GetField(xyzt,f);
-
-      z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
-      y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
-      x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
-
-      secxs[3] = (bt*f[2] - ct*f[1])* ph2;
-      secys[3] = (ct*f[0] - at*f[2])* ph2;
-      seczs[3] = (at*f[1] - bt*f[0])* ph2;
-      a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
-      b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
-      c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
-
-      est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
-       + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
-       + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
-
-      if (est > kdlt && TMath::Abs(h) > 1.e-4) {
-       if (ncut++ > maxcut) break;
-       h *= khalf;
-       continue;
-      }
-
-      ncut = 0;
-      // *               if too many iterations, go to helix
-      if (iter++ > maxit) break;
-
-      tl += h;
-      if (est < kdlt32) 
-       h *= 2.;
-      cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
-      vout[0] = x;
-      vout[1] = y;
-      vout[2] = z;
-      vout[3] = cba*a;
-      vout[4] = cba*b;
-      vout[5] = cba*c;
-      rest = step - tl;
-      if (step < 0.) rest = -rest;
-      if (rest < 1.e-5*TMath::Abs(step)) return;
-
-    } while(1);
-
-    // angle too big, use helix
-
-    f1  = f[0];
-    f2  = f[1];
-    f3  = f[2];
-    f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
-    rho = -f4*pinv;
-    tet = rho * step;
-    hnorm = 1./f4;
-    f1 = f1*hnorm;
-    f2 = f2*hnorm;
-    f3 = f3*hnorm;
-
-    hxp[0] = f2*vect[kipz] - f3*vect[kipy];
-    hxp[1] = f3*vect[kipx] - f1*vect[kipz];
-    hxp[2] = f1*vect[kipy] - f2*vect[kipx];
-    hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
-
-    rho1 = 1./rho;
-    sint = TMath::Sin(tet);
-    cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
-
-    g1 = sint*rho1;
-    g2 = cost*rho1;
-    g3 = (tet-sint) * hp*rho1;
-    g4 = -cost;
-    g5 = sint;
-    g6 = cost * hp;
-    vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
-    vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
-    vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
-    vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
-    vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
-    vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
-
-    return;
-}
-//___________________________________________________________
- void  AliMUONTrackParam::GetField(Double_t *Position, Double_t *Field) const
-{
-    /// interface for arguments in double precision (Why ? ChF)
-    Float_t x[3], b[3];
-
-    x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2];
-
-    fkField->Field(x, b);
-    Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2];
-
-    return;
-  }
 //_____________________________________________-
 void AliMUONTrackParam::Print(Option_t* opt) const
 {
index 6601aed..fb71200 100644 (file)
@@ -18,7 +18,6 @@
 #include "AliMUONHitForRec.h"
 
 class AliESDMuonTrack;
-class AliMagF;
 
 class AliMUONTrackParam : public TObject 
 {
@@ -32,7 +31,6 @@ class AliMUONTrackParam : public TObject
   void GetParamFrom(const AliESDMuonTrack& esdMuonTrack);
   void SetParamFor(AliESDMuonTrack& esdMuonTrack);
 
-
   // Get and Set methods for data
        /// return inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
   Double_t GetInverseBendingMomentum(void) const {return fInverseBendingMomentum;}
@@ -73,29 +71,9 @@ class AliMUONTrackParam : public TObject
        /// "Compare" function for sorting
   Int_t Compare(const TObject* TrackParam) const;
 
-  void ExtrapToZ(Double_t Z);
-  void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
-  void ExtrapToVertex(Double_t xVtx, Double_t yVtx, Double_t zVtx);  // extrapolation to vertex through the absorber (with true vertex) 
-   void BransonCorrection(Double_t xVtx, Double_t yVtx, Double_t zVtx); // makes Branson correction with true vertex  
-   // returns total momentum after energy loss correction in the absorber
-  Double_t TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta);
-  void FieldCorrection(Double_t Z); // makes simple magnetic field correction through the absorber 
-
-  void ExtrapOneStepHelix(Double_t charge, Double_t step, 
-                         Double_t *vect, Double_t *vout) const;
-  void ExtrapOneStepHelix3(Double_t field, Double_t step, 
-                          Double_t *vect, Double_t *vout) const;
-
-  void ExtrapOneStepRungekutta(Double_t charge, Double_t step, 
-                              Double_t* vect, Double_t* vout) const;
-  
   virtual void Print(Option_t* opt="") const;
  
-       /// set field map
-  void SetField(const AliMagF* magField) {fkField = magField;}
 
-
- protected:
  private:
   Double_t fInverseBendingMomentum; ///< Inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
   Double_t fBendingSlope; ///< Bending slope (cm ** -1)
@@ -104,16 +82,9 @@ class AliMUONTrackParam : public TObject
   Double_t fBendingCoor; ///< bending coordinate (cm)
   Double_t fNonBendingCoor; ///< non bending coordinate (cm)
 
-  const AliMagF* fkField;     //!< field map
-
-  void SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward);
-  void GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge);
-
-  void GetField(Double_t *Position, Double_t *Field) const;
-  
   AliMUONHitForRec *fHitForRecPtr; //!< Pointer to associated HitForRec if any
   
   ClassDef(AliMUONTrackParam, 2) // Track parameters in ALICE dimuon spectrometer
-    };
+};
        
 #endif
index 1cb6010..2c365ff 100644 (file)
@@ -37,6 +37,8 @@
 #include "AliMUONHitForRec.h"
 #include "AliMUONSegment.h"
 #include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliLog.h"
 #include <TVirtualFitter.h>
 
@@ -408,7 +410,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
       // Track parameters at first track hit (smallest Z)
       trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
       // extrapolation to station
-      trackParam1->ExtrapToStation(station, trackParam);
+      AliMUONTrackExtrap::ExtrapToStation(trackParam1, station, trackParam);
       extrapSegment = new AliMUONSegment(); //  empty segment
       // multiple scattering factor corresponding to one chamber
       // and momentum in bending plane (not total)
@@ -443,8 +445,8 @@ void AliMUONTrackReconstructor::FollowTracks(void)
        segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
        // correction of corrected segment (fBendingCoor and fNonBendingCoor)
        // according to real Z value of "segment" and slopes of "extrapSegment"
-       trackParam[0].ExtrapToZ(segment->GetZ());
-       trackParam[1].ExtrapToZ(segment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
+       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), segment->GetZ());
+       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), segment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
        extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
        extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
        extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
@@ -458,9 +460,9 @@ void AliMUONTrackReconstructor::FollowTracks(void)
       }
       if (bestSegment) {
        // best segment found: add it to track candidate
-       trackParam[0].ExtrapToZ(bestSegment->GetZ());
+       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), bestSegment->GetZ());
        track->AddTrackParamAtHit(&(trackParam[0]),bestSegment->GetHitForRec1());
-       trackParam[1].ExtrapToZ(bestSegment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
+       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), bestSegment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
        track->AddTrackParamAtHit(&(trackParam[1]),bestSegment->GetHitForRec2());
        AliDebug(3, Form("FollowTracks: track candidate(0..): %d  Added segment in station(0..): %d", trackIndex, station));
        if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
@@ -481,7 +483,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
          for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; iHit < fIndexOfFirstHitForRecPerChamber[ch]+fNHitsForRecPerChamber[ch]; iHit++) {
            hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
            // coordinates of extrapolated hit
-           trackParam[chInStation].ExtrapToZ(hit->GetZ());
+           AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chInStation]), hit->GetZ());
            extrapHit->SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
            extrapHit->SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
            // resolutions from "extrapSegment"
@@ -500,7 +502,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
        }
        if (bestHit) {
          // best hit found: add it to track candidate
-         trackParam[chBestHit].ExtrapToZ(bestHit->GetZ());
+         AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chBestHit]), bestHit->GetZ());
          track->AddTrackParamAtHit(&(trackParam[chBestHit]),bestHit);
          if (AliLog::GetGlobalDebugLevel() > 2) {
            cout << "FollowTracks: track candidate(0..): " << trackIndex
@@ -562,7 +564,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
       // after going through the first station
       if (station == 0) {
        trackParamVertex = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()));
-       (&trackParamVertex)->ExtrapToVertex(0.,0.,0.);
+       AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, 0., 0., 0.);
        track->SetTrackParamAtVertex(&trackParamVertex);
        if (AliLog::GetGlobalDebugLevel() > 0) {
          cout << "FollowTracks: track candidate(0..): " << trackIndex
@@ -699,7 +701,7 @@ void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t
   while (trackParamAtHit) {
     hitForRec = trackParamAtHit->GetHitForRecPtr();
     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
-    param1.ExtrapToZ(hitForRec->GetZ());
+    AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
     // update track parameters of the current hit
     trackParamAtHit->SetTrackParam(param1);
     // Increment Chi2
@@ -757,7 +759,7 @@ void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double
     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
     hitForRec = trackParamAtHit->GetHitForRecPtr();
     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
-    param1.ExtrapToZ(hitForRec->GetZ());
+    AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
     // update track parameters of the current hit
     trackParamAtHit->SetTrackParam(param1);
     // square of multiple scattering angle at current hit, with one chamber
index 02c956a..90bb6ab 100644 (file)
@@ -44,6 +44,8 @@
 #include "AliMUONGlobalTrigger.h"
 #include "AliMUONSegment.h"
 #include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliMagF.h"
 #include "AliLog.h"
 #include "AliTracker.h"
@@ -124,6 +126,9 @@ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data)
   fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]);
   // See how to get fSimple(BValue, BLength, BPosition)
   // automatically calculated from the actual magnetic field ????
+  
+  // set the magnetic field for track extrapolations
+  AliMUONTrackExtrap::SetField(kField);
 }
 
   //__________________________________________________________________________
@@ -417,7 +422,7 @@ void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void)
     chi2MatchTrigger = 0.;
     
     trackParam = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last()));
-    trackParam.ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
+    AliMUONTrackExtrap::ExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
     
     xTrack = trackParam.GetNonBendingCoor();
     yTrack = trackParam.GetBendingCoor();
index e73fb35..8aefd6a 100644 (file)
@@ -50,6 +50,7 @@
 #include "AliMUONLocalTrigger.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliMUONTriggerCircuit.h"
 #include "AliMUONTriggerCrateStore.h"
 
@@ -630,11 +631,13 @@ void MUONrectracks (Int_t event2Check=0, char * filename="galice.root"){
     AliMUONTrack* recTrack = 0;
     AliMUONTrackParam* trackParam = 0;
 
+    // set the magnetic field for track extrapolations
+    AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
     for (Int_t iRecTracks = 0; iRecTracks <  nrectracks;  iRecTracks++) {
    //  // reading info from tracks
        recTrack = (AliMUONTrack*) RecTracks->At(iRecTracks);
        trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
-       trackParam->ExtrapToZ(0.0);
+       AliMUONTrackExtrap::ExtrapToZ(trackParam,0.);
       recTrack->Print("full");
     } // end loop tracks
 
index 79a5c39..e6adc8f 100644 (file)
@@ -63,6 +63,7 @@
 
 // MUON includes
 #include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliESDMuonTrack.h"
 #endif
 
@@ -314,6 +315,8 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
       cout << " number of tracks: " << nTracks  <<endl;
     }
 
+    // set the magnetic field for track extrapolations
+    AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
     // loop over all reconstructed tracks (also first track of combination)
     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 
@@ -322,7 +325,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
       if (!Vertex->GetNContributors()) {
        //re-extrapolate to vertex, if not kown before.
        trackParam.GetParamFrom(*muonTrack);
-       trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex);
+       AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
        trackParam.SetParamFor(*muonTrack);
       }
 
@@ -388,7 +391,7 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven
           
          if (!Vertex->GetNContributors()) {
            trackParam.GetParamFrom(*muonTrack);
-           trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex);
+           AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
            trackParam.SetParamFor(*muonTrack);
          }
 
index f1755bd..de460fd 100644 (file)
@@ -24,6 +24,7 @@
 
 // MUON includes
 #include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
 #include "AliESDMuonTrack.h"
 #endif
 //
@@ -196,6 +197,8 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t
     //    printf("\n Nb of events analysed: %d\r",iEvent);
     //      cout << " number of tracks: " << nTracks  <<endl;
   
+    // set the magnetic field for track extrapolations
+    AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
     // loop over all reconstructed tracks (also first track of combination)
     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
 
@@ -204,7 +207,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t
       if (!Vertex->GetNContributors()) {
        //re-extrapolate to vertex, if not kown before.
        trackParam.GetParamFrom(*muonTrack);
-       trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex);
+       AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
        trackParam.SetParamFor(*muonTrack);
       }
       thetaX = muonTrack->GetThetaX();
@@ -259,7 +262,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t
 
          if (!Vertex->GetNContributors()) {
            trackParam.GetParamFrom(*muonTrack);
-           trackParam.ExtrapToVertex(fXVertex, fYVertex, fZVertex);
+           AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
            trackParam.SetParamFor(*muonTrack);
          }
 
index 1e5360b..425a14f 100644 (file)
@@ -17,6 +17,7 @@
 #pragma link C++ class AliMUONTrack+; 
 #pragma link C++ class AliMUONTrackK+; 
 #pragma link C++ class AliMUONTrackParam+; 
+#pragma link C++ class AliMUONTrackExtrap+; 
 #pragma link C++ class AliMUONTriggerTrack+; 
 #pragma link C++ class AliMUONRecoTrack+; 
 #pragma link C++ class AliMUONHitForRec+; 
index 65ed2f0..972faae 100644 (file)
@@ -13,6 +13,7 @@ SRCS:= AliMUONClusterReconstructor.cxx \
        AliMUONTrack.cxx \
        AliMUONTrackK.cxx \
        AliMUONTrackParam.cxx \
+       AliMUONTrackExtrap.cxx \
        AliMUONTriggerTrack.cxx \
        AliMUONRecoTrack.cxx \
        AliMUONHitForRec.cxx \