MUON/AliMUONTrackParam.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////
  19 //
  20 // Track parameters
  21 // in
  22 // ALICE
  23 // dimuon
  24 // spectrometer
  25 //
  26 ///////////////////////////////////////////////////
  27
  28 //#include <Riostream.h>
  29 #include "AliMUON.h"
  30 #include "AliMUONTrackParam.h"
  31 #include "AliMUONConstants.h"
  32 #include "AliESDMuonTrack.h"
  33 #include "AliMagF.h"
  34 #include "AliLog.h"
  35 #include "AliTracker.h"
  36 #include "AliMUONHitForRec.h"
  37
  38 ClassImp(AliMUONTrackParam) // Class implementation in ROOT context
  39
  40   //_________________________________________________________________________
  41 AliMUONTrackParam::AliMUONTrackParam()
  42   : TObject(),
  43     fInverseBendingMomentum(0.),
  44     fBendingSlope(0.),
  45     fNonBendingSlope(0.),
  46     fZ(0.),
  47     fBendingCoor(0.),
  48     fNonBendingCoor(0.),
  49     fkField(0x0),
  50     fHitForRecPtr(0x0)
  51 {
  52 /// Constructor
  53   // get field from outside
  54   fkField = AliTracker::GetFieldMap();
  55   if (!fkField) AliFatal("No field available");
  56 }
  57
  58   //_________________________________________________________________________
  59 AliMUONTrackParam::AliMUONTrackParam(const AliMUONTrackParam& theMUONTrackParam)
  60   : TObject(theMUONTrackParam),
  61     fInverseBendingMomentum(theMUONTrackParam.fInverseBendingMomentum),
  62     fBendingSlope(theMUONTrackParam.fBendingSlope),
  63     fNonBendingSlope(theMUONTrackParam.fNonBendingSlope),
  64     fZ(theMUONTrackParam.fZ),
  65     fBendingCoor(theMUONTrackParam.fBendingCoor),
  66     fNonBendingCoor(theMUONTrackParam.fNonBendingCoor),
  67     fkField(theMUONTrackParam.fkField),
  68     fHitForRecPtr(theMUONTrackParam.fHitForRecPtr)
  69 {
  70   /// Copy constructor
  71 }
  72
  73   //_________________________________________________________________________
  74 AliMUONTrackParam& AliMUONTrackParam::operator=(const AliMUONTrackParam& theMUONTrackParam)
  75 {
  76   /// Asignment operator
  77   if (this == &theMUONTrackParam)
  78     return *this;
  79
  80   // base class assignement
  81   TObject::operator=(theMUONTrackParam);
  82
  83   fInverseBendingMomentum =  theMUONTrackParam.fInverseBendingMomentum;
  84   fBendingSlope           =  theMUONTrackParam.fBendingSlope;
  85   fNonBendingSlope        =  theMUONTrackParam.fNonBendingSlope;
  86   fZ                      =  theMUONTrackParam.fZ;
  87   fBendingCoor            =  theMUONTrackParam.fBendingCoor;
  88   fNonBendingCoor         =  theMUONTrackParam.fNonBendingCoor;
  89   fkField                 =  theMUONTrackParam.fkField;
  90   fHitForRecPtr           =  theMUONTrackParam.fHitForRecPtr;
  91
  92   return *this;
  93 }
  94
  95   //__________________________________________________________________________
  96 AliMUONTrackParam::~AliMUONTrackParam()
  97 {
  98 /// Destructor
  99 /// Update the number of TrackHit's connected to the attached HitForRec if any
 100   if (fHitForRecPtr) fHitForRecPtr->SetNTrackHits(fHitForRecPtr->GetNTrackHits() - 1); // decrement NTrackHits of hit
 101 }
 102
 103   //__________________________________________________________________________
 104 void AliMUONTrackParam::SetTrackParam(AliMUONTrackParam& theMUONTrackParam)
 105 {
 106   /// Set track parameters from "TrackParam" leaving pointer to fHitForRecPtr unchanged
 107   fInverseBendingMomentum =  theMUONTrackParam.fInverseBendingMomentum;
 108   fBendingSlope           =  theMUONTrackParam.fBendingSlope;
 109   fNonBendingSlope        =  theMUONTrackParam.fNonBendingSlope;
 110   fZ                      =  theMUONTrackParam.fZ;
 111   fBendingCoor            =  theMUONTrackParam.fBendingCoor;
 112   fNonBendingCoor         =  theMUONTrackParam.fNonBendingCoor;
 113
 114 }
 115
 116   //__________________________________________________________________________
 117 AliMUONHitForRec* AliMUONTrackParam::GetHitForRecPtr(void) const
 118 {
 119 /// return pointer to HitForRec attached to the current TrackParam
 120 /// this method should not be called when fHitForRecPtr == NULL
 121   if (!fHitForRecPtr) AliWarning("AliMUONTrackParam::GetHitForRecPtr: fHitForRecPtr == NULL");
 122   return fHitForRecPtr;
 123 }
 124
 125   //__________________________________________________________________________
 126 Int_t AliMUONTrackParam::Compare(const TObject* TrackParam) const
 127 {
 128 /// "Compare" function to sort with decreasing Z (spectro. muon Z <0).
 129 /// Returns 1 (0, -1) if Z of current TrackHit
 130 /// is smaller than (equal to, larger than) Z of TrackHit
 131   if (fHitForRecPtr->GetZ() < ((AliMUONTrackParam*)TrackParam)->fHitForRecPtr->GetZ()) return(1);
 132   else if (fHitForRecPtr->GetZ() == ((AliMUONTrackParam*)TrackParam)->fHitForRecPtr->GetZ()) return(0);
 133   else return(-1);
 134 }
 135
 136   //_________________________________________________________________________
 137 void AliMUONTrackParam::GetParamFrom(const AliESDMuonTrack& esdMuonTrack)
 138 {
 139   /// assigned value form ESD track.
 140   fInverseBendingMomentum =  esdMuonTrack.GetInverseBendingMomentum();
 141   fBendingSlope           =  TMath::Tan(esdMuonTrack.GetThetaY());
 142   fNonBendingSlope        =  TMath::Tan(esdMuonTrack.GetThetaX());
 143   fZ                      =  esdMuonTrack.GetZ();
 144   fBendingCoor            =  esdMuonTrack.GetBendingCoor();
 145   fNonBendingCoor         =  esdMuonTrack.GetNonBendingCoor();
 146 }
 147
 148   //_________________________________________________________________________
 149 void AliMUONTrackParam::SetParamFor(AliESDMuonTrack& esdMuonTrack)
 150 {
 151   /// assigned value form ESD track.
 152   esdMuonTrack.SetInverseBendingMomentum(fInverseBendingMomentum);
 153   esdMuonTrack.SetThetaX(TMath::ATan(fNonBendingSlope));
 154   esdMuonTrack.SetThetaY(TMath::ATan(fBendingSlope));
 155   esdMuonTrack.SetZ(fZ);
 156   esdMuonTrack.SetBendingCoor(fBendingCoor);
 157   esdMuonTrack.SetNonBendingCoor(fNonBendingCoor);
 158 }
 159
 160   //__________________________________________________________________________
 161 void AliMUONTrackParam::ExtrapToZ(Double_t Z)
 162 {
 163   /// Track parameter extrapolation to the plane at "Z".
 164   /// On return, the track parameters resulting from the extrapolation
 165   /// replace the current track parameters.
 166   if (this->fZ == Z) return; // nothing to be done if same Z
 167   Double_t forwardBackward; // +1 if forward, -1 if backward
 168   if (Z < this->fZ) forwardBackward = 1.0; // spectro. z<0
 169   else forwardBackward = -1.0;
 170   Double_t vGeant3[7], vGeant3New[7]; // 7 in parameter ????
 171   Int_t iGeant3, stepNumber;
 172   Int_t maxStepNumber = 5000; // in parameter ????
 173   // For safety: return kTRUE or kFALSE ????
 174   // Parameter vector for calling EXTRAP_ONESTEP
 175   SetGeant3Parameters(vGeant3, forwardBackward);
 176   // sign of charge (sign of fInverseBendingMomentum if forward motion)
 177   // must be changed if backward extrapolation
 178   Double_t chargeExtrap = forwardBackward *
 179     TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
 180   Double_t stepLength = 6.0; // in parameter ????
 181   // Extrapolation loop
 182   stepNumber = 0;
 183   while (((-forwardBackward * (vGeant3[2] - Z)) <= 0.0) &&  // spectro. z<0
 184          (stepNumber < maxStepNumber)) {
 185     stepNumber++;
 186     // Option for switching between helix and Runge-Kutta ????
 187     //ExtrapOneStepRungekutta(chargeExtrap, stepLength, vGeant3, vGeant3New);
 188     ExtrapOneStepHelix(chargeExtrap, stepLength, vGeant3, vGeant3New);
 189     if ((-forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z spectro. z<0
 190     // better use TArray ????
 191     for (iGeant3 = 0; iGeant3 < 7; iGeant3++)
 192       {vGeant3[iGeant3] = vGeant3New[iGeant3];}
 193   }
 194   // check maxStepNumber ????
 195   // Interpolation back to exact Z (2nd order)
 196   // should be in function ???? using TArray ????
 197   Double_t dZ12 = vGeant3New[2] - vGeant3[2]; // 1->2
 198   if (TMath::Abs(dZ12) > 0) {
 199     Double_t dZ1i = Z - vGeant3[2]; // 1-i
 200     Double_t dZi2 = vGeant3New[2] - Z; // i->2
 201     Double_t xPrime = (vGeant3New[0] - vGeant3[0]) / dZ12;
 202     Double_t xSecond =
 203       ((vGeant3New[3] / vGeant3New[5]) - (vGeant3[3] / vGeant3[5])) / dZ12;
 204     Double_t yPrime = (vGeant3New[1] - vGeant3[1]) / dZ12;
 205     Double_t ySecond =
 206       ((vGeant3New[4] / vGeant3New[5]) - (vGeant3[4] / vGeant3[5])) / dZ12;
 207     vGeant3[0] = vGeant3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
 208     vGeant3[1] = vGeant3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
 209     vGeant3[2] = Z; // Z
 210     Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
 211     Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
 212     // (PX, PY, PZ)/PTOT assuming forward motion
 213     vGeant3[5] =
 214       1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
 215     vGeant3[3] = xPrimeI * vGeant3[5]; // PX/PTOT
 216     vGeant3[4] = yPrimeI * vGeant3[5]; // PY/PTOT
 217   } else {
 218     AliWarning(Form("Extrap. to Z not reached, Z = %f",Z));
 219   }
 220   // Track parameters from Geant3 parameters,
 221   // with charge back for forward motion
 222   GetFromGeant3Parameters(vGeant3, chargeExtrap * forwardBackward);
 223 }
 224
 225   //__________________________________________________________________________
 226 void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward)
 227 {
 228   /// Set vector of Geant3 parameters pointed to by "VGeant3"
 229   /// from track parameters in current AliMUONTrackParam.
 230   /// Since AliMUONTrackParam is only geometry, one uses "ForwardBackward"
 231   /// to know whether the particle is going forward (+1) or backward (-1).
 232   VGeant3[0] = this->fNonBendingCoor; // X
 233   VGeant3[1] = this->fBendingCoor; // Y
 234   VGeant3[2] = this->fZ; // Z
 235   Double_t pYZ = TMath::Abs(1.0 / this->fInverseBendingMomentum);
 236   Double_t pZ =
 237     pYZ / TMath::Sqrt(1.0 + this->fBendingSlope * this->fBendingSlope);
 238   VGeant3[6] =
 239     TMath::Sqrt(pYZ * pYZ +
 240                 pZ * pZ * this->fNonBendingSlope * this->fNonBendingSlope); // PTOT
 241   VGeant3[5] = -ForwardBackward * pZ / VGeant3[6]; // PZ/PTOT spectro. z<0
 242   VGeant3[3] = this->fNonBendingSlope * VGeant3[5]; // PX/PTOT
 243   VGeant3[4] = this->fBendingSlope * VGeant3[5]; // PY/PTOT
 244 }
 245
 246   //__________________________________________________________________________
 247 void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge)
 248 {
 249   /// Get track parameters in current AliMUONTrackParam
 250   /// from Geant3 parameters pointed to by "VGeant3",
 251   /// assumed to be calculated for forward motion in Z.
 252   /// "InverseBendingMomentum" is signed with "Charge".
 253   this->fNonBendingCoor = VGeant3[0]; // X
 254   this->fBendingCoor = VGeant3[1]; // Y
 255   this->fZ = VGeant3[2]; // Z
 256   Double_t pYZ = VGeant3[6] * TMath::Sqrt(1.0 - VGeant3[3] * VGeant3[3]);
 257   this->fInverseBendingMomentum = Charge / pYZ;
 258   this->fBendingSlope = VGeant3[4] / VGeant3[5];
 259   this->fNonBendingSlope = VGeant3[3] / VGeant3[5];
 260 }
 261
 262   //__________________________________________________________________________
 263 void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam)
 264 {
 265   /// Track parameters extrapolated from current track parameters ("this")
 266   /// to both chambers of the station(0..) "Station"
 267   /// are returned in the array (dimension 2) of track parameters
 268   /// pointed to by "TrackParam" (index 0 and 1 for first and second chambers).
 269   Double_t extZ[2], z1, z2;
 270   Int_t i1 = -1, i2 = -1; // = -1 to avoid compilation warnings
 271   // range of Station to be checked ????
 272   z1 = AliMUONConstants::DefaultChamberZ(2 * Station);
 273   z2 = AliMUONConstants::DefaultChamberZ(2 * Station + 1);
 274   // First and second Z to extrapolate at
 275   if ((z1 > this->fZ) && (z2 > this->fZ)) {i1 = 0; i2 = 1;}
 276   else if ((z1 < this->fZ) && (z2 < this->fZ)) {i1 = 1; i2 = 0;}
 277   else {
 278         AliError(Form("Starting Z (%f) in between z1 (%f) and z2 (%f) of station(0..)%d",this->fZ,z1,z2,Station));
 279 //     cout << "ERROR in AliMUONTrackParam::CreateExtrapSegmentInStation" << endl;
 280 //     cout << "Starting Z (" << this->fZ << ") in between z1 (" << z1 <<
 281 //       ") and z2 (" << z2 << ") of station(0..) " << Station << endl;
 282   }
 283   extZ[i1] = z1;
 284   extZ[i2] = z2;
 285   // copy of track parameters
 286   TrackParam[i1] = *this;
 287   // first extrapolation
 288   (&(TrackParam[i1]))->ExtrapToZ(extZ[0]);
 289   TrackParam[i2] = TrackParam[i1];
 290   // second extrapolation
 291   (&(TrackParam[i2]))->ExtrapToZ(extZ[1]);
 292   return;
 293 }
 294
 295   //__________________________________________________________________________
 296 void AliMUONTrackParam::ExtrapToVertex(Double_t xVtx, Double_t yVtx, Double_t zVtx)
 297 {
 298   /// Extrapolation to the vertex.
 299   /// Returns the track parameters resulting from the extrapolation in the current TrackParam.
 300   /// Changes parameters according to Branson correction through the absorber
 301
 302   Double_t zAbsorber = -503.0; // to be coherent with the Geant absorber geometry !!!!
 303                                // spectro. (z<0)
 304   // Extrapolates track parameters upstream to the "Z" end of the front absorber
 305   ExtrapToZ(zAbsorber); // !!!
 306   // Makes Branson correction (multiple scattering + energy loss)
 307   BransonCorrection(xVtx,yVtx,zVtx);
 308   // Makes a simple magnetic field correction through the absorber
 309   FieldCorrection(zAbsorber);
 310 }
 311
 312
 313 //  Keep this version for future developments
 314   //__________________________________________________________________________
 315 // void AliMUONTrackParam::BransonCorrection()
 316 // {
 317 //   // Branson correction of track parameters
 318 //   // the entry parameters have to be calculated at the end of the absorber
 319 //   Double_t zEndAbsorber, zBP, xBP, yBP;
 320 //   Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
 321 //   Int_t sign;
 322 //   // Would it be possible to calculate all that from Geant configuration ????
 323 //   // and to get the Branson parameters from a function in ABSO module ????
 324 //   // with an eventual contribution from other detectors like START ????
 325 //   // Radiation lengths outer part theta > 3 degres
 326 //   static Double_t x01[9] = { 18.8,    // C (cm)
 327 //                           10.397,   // Concrete (cm)
 328 //                           0.56,    // Plomb (cm)
 329 //                           47.26,   // Polyethylene (cm)
 330 //                           0.56,   // Plomb (cm)
 331 //                           47.26,   // Polyethylene (cm)
 332 //                           0.56,   // Plomb (cm)
 333 //                           47.26,   // Polyethylene (cm)
 334 //                           0.56 };   // Plomb (cm)
 335 //   // inner part theta < 3 degres
 336 //   static Double_t x02[3] = { 18.8,    // C (cm)
 337 //                           10.397,   // Concrete (cm)
 338 //                           0.35 };    // W (cm)
 339 //   // z positions of the materials inside the absober outer part theta > 3 degres
 340 //   static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
 341 //   // inner part theta < 3 degres
 342 //   static Double_t z2[4] = { 90, 315, 467, 503 };
 343 //   static Bool_t first = kTRUE;
 344 //   static Double_t zBP1, zBP2, rLimit;
 345 //   // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
 346 //   if (first) {
 347 //     first = kFALSE;
 348 //     Double_t aNBP = 0.0;
 349 //     Double_t aDBP = 0.0;
 350 //     Int_t iBound;
 351
 352 //     for (iBound = 0; iBound < 9; iBound++) {
 353 //       aNBP = aNBP +
 354 //      (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
 355 //       z1[iBound]   * z1[iBound]   * z1[iBound]    ) / x01[iBound];
 356 //       aDBP = aDBP +
 357 //      (z1[iBound+1] * z1[iBound+1] - z1[iBound]   * z1[iBound]    ) / x01[iBound];
 358 //     }
 359 //     zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
 360 //     aNBP = 0.0;
 361 //     aDBP = 0.0;
 362 //     for (iBound = 0; iBound < 3; iBound++) {
 363 //       aNBP = aNBP +
 364 //      (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
 365 //       z2[iBound]   * z2[iBound ]  * z2[iBound]    ) / x02[iBound];
 366 //       aDBP = aDBP +
 367 //      (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
 368 //     }
 369 //     zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
 370 //     rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
 371 //   }
 372
 373 //   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 374 //   sign = 1;
 375 //   if (fInverseBendingMomentum < 0) sign = -1;
 376 //   pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
 377 //   pX = pZ * fNonBendingSlope;
 378 //   pY = pZ * fBendingSlope;
 379 //   pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
 380 //   xEndAbsorber = fNonBendingCoor;
 381 //   yEndAbsorber = fBendingCoor;
 382 //   radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
 383
 384 //   if (radiusEndAbsorber2 > rLimit*rLimit) {
 385 //     zEndAbsorber = z1[9];
 386 //     zBP = zBP1;
 387 //   } else {
 388 //     zEndAbsorber = z2[3];
 389 //     zBP = zBP2;
 390 //   }
 391
 392 //   xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
 393 //   yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
 394
 395 //   // new parameters after Branson and energy loss corrections
 396 //   pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
 397 //   pX = pZ * xBP / zBP;
 398 //   pY = pZ * yBP / zBP;
 399 //   fBendingSlope = pY / pZ;
 400 //   fNonBendingSlope = pX / pZ;
 401
 402 //   pT = TMath::Sqrt(pX * pX + pY * pY);
 403 //   theta = TMath::ATan2(pT, pZ);
 404 //   pTotal =
 405 //     TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
 406
 407 //   fInverseBendingMomentum = (sign / pTotal) *
 408 //     TMath::Sqrt(1.0 +
 409 //              fBendingSlope * fBendingSlope +
 410 //              fNonBendingSlope * fNonBendingSlope) /
 411 //     TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
 412
 413 //   // vertex position at (0,0,0)
 414 //   // should be taken from vertex measurement ???
 415 //   fBendingCoor = 0.0;
 416 //   fNonBendingCoor = 0;
 417 //   fZ= 0;
 418 // }
 419
 420 void AliMUONTrackParam::BransonCorrection(Double_t xVtx,Double_t yVtx,Double_t zVtx)
 421 {
 422   /// Branson correction of track parameters
 423   // the entry parameters have to be calculated at the end of the absorber
 424   // simplified version: the z positions of Branson's planes are no longer calculated
 425   // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C
 426   // to test this correction.
 427   // Would it be possible to calculate all that from Geant configuration ????
 428   // and to get the Branson parameters from a function in ABSO module ????
 429   // with an eventual contribution from other detectors like START ????
 430   // change to take into account the vertex postition (real, reconstruct,....)
 431
 432   Double_t  zBP, xBP, yBP;
 433   Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
 434   Int_t sign;
 435   static Bool_t first = kTRUE;
 436   static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber;
 437   // zBP1 for outer part and zBP2 for inner part (only at the first call)
 438   if (first) {
 439     first = kFALSE;
 440
 441     zEndAbsorber = -503;  // spectro (z<0)
 442     thetaLimit = 3.0 * (TMath::Pi()) / 180.;
 443     rLimit = TMath::Abs(zEndAbsorber) * TMath::Tan(thetaLimit);
 444     zBP1 = -450; // values close to those calculated with EvalAbso.C
 445     zBP2 = -480;
 446   }
 447
 448   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 449   sign = 1;
 450   if (fInverseBendingMomentum < 0) sign = -1;
 451   pZ = Pz();
 452   pX = Px();
 453   pY = Py();
 454   pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
 455   xEndAbsorber = fNonBendingCoor;
 456   yEndAbsorber = fBendingCoor;
 457   radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
 458
 459   if (radiusEndAbsorber2 > rLimit*rLimit) {
 460     zBP = zBP1;
 461   } else {
 462     zBP = zBP2;
 463   }
 464
 465   xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
 466   yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
 467
 468   // new parameters after Branson and energy loss corrections
 469 //   Float_t zSmear = zBP - gRandom->Gaus(0.,2.);  // !!! possible smearing of Z vertex position
 470
 471   Float_t zSmear = zBP ;
 472
 473    pZ = pTotal * (zSmear-zVtx) / TMath::Sqrt((xBP-xVtx) * (xBP-xVtx) + (yBP-yVtx) * (yBP-yVtx) +( zSmear-zVtx) * (zSmear-zVtx) );
 474    pX = pZ * (xBP - xVtx)/ (zSmear-zVtx);
 475    pY = pZ * (yBP - yVtx) / (zSmear-zVtx);
 476   fBendingSlope = pY / pZ;
 477   fNonBendingSlope = pX / pZ;
 478
 479
 480   pT = TMath::Sqrt(pX * pX + pY * pY);
 481   theta = TMath::ATan2(pT, TMath::Abs(pZ));
 482   pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta);
 483
 484   fInverseBendingMomentum = (sign / pTotal) *
 485     TMath::Sqrt(1.0 +
 486                 fBendingSlope * fBendingSlope +
 487                 fNonBendingSlope * fNonBendingSlope) /
 488     TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
 489
 490   // vertex position at (0,0,0)
 491   // should be taken from vertex measurement ???
 492
 493   fBendingCoor = xVtx;
 494   fNonBendingCoor = yVtx;
 495   fZ= zVtx;
 496
 497 }
 498
 499   //__________________________________________________________________________
 500 Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta)
 501 {
 502   /// Returns the total momentum corrected from energy loss in the front absorber
 503   // One can use the macros MUONTestAbso.C and DrawTestAbso.C
 504   // to test this correction.
 505   // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002)
 506   Double_t deltaP, pTotalCorrected;
 507
 508    // Parametrization to be redone according to change of absorber material ????
 509   // See remark in function BransonCorrection !!!!
 510   // The name is not so good, and there are many arguments !!!!
 511   if (theta  < thetaLimit ) {
 512     if (pTotal < 20) {
 513       deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal;
 514     } else {
 515       deltaP = 3.0714 + 0.011767 *pTotal;
 516     }
 517     deltaP *= 0.75; // AZ
 518   } else {
 519     if (pTotal < 20) {
 520       deltaP  = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal;
 521     } else {
 522       deltaP = 2.6069 + 0.0051705 * pTotal;
 523     }
 524     deltaP *= 0.9; // AZ
 525   }
 526   pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
 527   return pTotalCorrected;
 528 }
 529
 530   //__________________________________________________________________________
 531 void AliMUONTrackParam::FieldCorrection(Double_t Z)
 532 {
 533   /// Correction of the effect of the magnetic field in the absorber
 534   // Assume a constant field along Z axis.
 535   Float_t b[3],x[3];
 536   Double_t bZ;
 537   Double_t pYZ,pX,pY,pZ,pT;
 538   Double_t pXNew,pYNew;
 539   Double_t c;
 540
 541   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 542   c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge
 543
 544   pZ = Pz();
 545   pX = Px();
 546   pY = Py();
 547   pT = TMath::Sqrt(pX*pX+pY*pY);
 548
 549   if (TMath::Abs(pZ) <= 0) return;
 550   x[2] = Z/2;
 551   x[0] = x[2]*fNonBendingSlope;
 552   x[1] = x[2]*fBendingSlope;
 553
 554   // Take magn. field value at position x.
 555   fkField->Field(x, b);
 556   bZ =  b[2];
 557
 558   // Transverse momentum rotation
 559   // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ.
 560   Double_t phiShift = c*0.436*0.0003*bZ*Z/pZ;
 561  // Rotate momentum around Z axis.
 562   pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift);
 563   pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift);
 564
 565   fBendingSlope = pYNew / pZ;
 566   fNonBendingSlope = pXNew / pZ;
 567
 568   fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ);
 569
 570 }
 571   //__________________________________________________________________________
 572 Double_t AliMUONTrackParam::Px() const
 573 {
 574   /// return px from track paramaters
 575   Double_t pYZ, pZ, pX;
 576   pYZ = 0;
 577   if (  TMath::Abs(fInverseBendingMomentum) > 0 )
 578     pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 579   pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
 580   pX = pZ * fNonBendingSlope;
 581   return pX;
 582 }
 583   //__________________________________________________________________________
 584 Double_t AliMUONTrackParam::Py() const
 585 {
 586   /// return px from track paramaters
 587   Double_t pYZ, pZ, pY;
 588   pYZ = 0;
 589   if (  TMath::Abs(fInverseBendingMomentum) > 0 )
 590     pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 591   pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
 592   pY = pZ * fBendingSlope;
 593   return pY;
 594 }
 595   //__________________________________________________________________________
 596 Double_t AliMUONTrackParam::Pz() const
 597 {
 598   /// return px from track paramaters
 599   Double_t pYZ, pZ;
 600   pYZ = 0;
 601   if (  TMath::Abs(fInverseBendingMomentum) > 0 )
 602     pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 603   pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
 604   return pZ;
 605 }
 606   //__________________________________________________________________________
 607 Double_t AliMUONTrackParam::P() const
 608 {
 609   /// return p from track paramaters
 610   Double_t  pYZ, pZ, p;
 611   pYZ = 0;
 612   if (  TMath::Abs(fInverseBendingMomentum) > 0 )
 613     pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
 614   pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
 615   p = TMath::Abs(pZ) *
 616     TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope + fNonBendingSlope * fNonBendingSlope);
 617   return p;
 618
 619 }
 620  //__________________________________________________________________________
 621 void AliMUONTrackParam::ExtrapOneStepHelix(Double_t charge, Double_t step,
 622                                          Double_t *vect, Double_t *vout) const
 623 {
 624 ///    ******************************************************************
 625 ///    *                                                                *
 626 ///    *  Performs the tracking of one step in a magnetic field         *
 627 ///    *  The trajectory is assumed to be a helix in a constant field   *
 628 ///    *  taken at the mid point of the step.                           *
 629 ///    *  Parameters:                                                   *
 630 ///    *   input                                                        *
 631 ///    *     STEP =arc length of the step asked                         *
 632 ///    *     VECT =input vector (position,direction cos and momentum)   *
 633 ///    *     CHARGE=  electric charge of the particle                   *
 634 ///    *   output                                                       *
 635 ///    *     VOUT = same as VECT after completion of the step           *
 636 ///    *                                                                *
 637 ///    *    ==>Called by : <USER>, GUSWIM                               *
 638 ///    *       Author    m.hansroul  *********                          *
 639 ///    *       modified  s.egli, s.v.levonian                           *
 640 ///    *       modified  v.perevoztchikov
 641 ///    *                                                                *
 642 ///    ******************************************************************
 643
 644 // modif: everything in double precision
 645
 646     Double_t xyz[3], h[4], hxp[3];
 647     Double_t h2xy, hp, rho, tet;
 648     Double_t sint, sintt, tsint, cos1t;
 649     Double_t f1, f2, f3, f4, f5, f6;
 650
 651     const Int_t kix  = 0;
 652     const Int_t kiy  = 1;
 653     const Int_t kiz  = 2;
 654     const Int_t kipx = 3;
 655     const Int_t kipy = 4;
 656     const Int_t kipz = 5;
 657     const Int_t kipp = 6;
 658
 659     const Double_t kec = 2.9979251e-4;
 660     //
 661     //    ------------------------------------------------------------------
 662     //
 663     //       units are kgauss,centimeters,gev/c
 664     //
 665     vout[kipp] = vect[kipp];
 666     if (TMath::Abs(charge) < 0.00001) {
 667       for (Int_t i = 0; i < 3; i++) {
 668         vout[i] = vect[i] + step * vect[i+3];
 669         vout[i+3] = vect[i+3];
 670       }
 671       return;
 672     }
 673     xyz[0]    = vect[kix] + 0.5 * step * vect[kipx];
 674     xyz[1]    = vect[kiy] + 0.5 * step * vect[kipy];
 675     xyz[2]    = vect[kiz] + 0.5 * step * vect[kipz];
 676
 677     //cmodif: call gufld (xyz, h) changed into:
 678     GetField (xyz, h);
 679
 680     h2xy = h[0]*h[0] + h[1]*h[1];
 681     h[3] = h[2]*h[2]+ h2xy;
 682     if (h[3] < 1.e-12) {
 683       for (Int_t i = 0; i < 3; i++) {
 684         vout[i] = vect[i] + step * vect[i+3];
 685         vout[i+3] = vect[i+3];
 686       }
 687       return;
 688     }
 689     if (h2xy < 1.e-12*h[3]) {
 690       ExtrapOneStepHelix3(charge*h[2], step, vect, vout);
 691       return;
 692     }
 693     h[3] = TMath::Sqrt(h[3]);
 694     h[0] /= h[3];
 695     h[1] /= h[3];
 696     h[2] /= h[3];
 697     h[3] *= kec;
 698
 699     hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy];
 700     hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz];
 701     hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx];
 702
 703     hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz];
 704
 705     rho = -charge*h[3]/vect[kipp];
 706     tet = rho * step;
 707
 708     if (TMath::Abs(tet) > 0.15) {
 709       sint = TMath::Sin(tet);
 710       sintt = (sint/tet);
 711       tsint = (tet-sint)/tet;
 712       cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet;
 713     } else {
 714       tsint = tet*tet/36.;
 715       sintt = (1. - tsint);
 716       sint = tet*sintt;
 717       cos1t = 0.5*tet;
 718     }
 719
 720     f1 = step * sintt;
 721     f2 = step * cos1t;
 722     f3 = step * tsint * hp;
 723     f4 = -tet*cos1t;
 724     f5 = sint;
 725     f6 = tet * cos1t * hp;
 726
 727     vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0];
 728     vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1];
 729     vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2];
 730
 731     vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0];
 732     vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1];
 733     vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2];
 734
 735     return;
 736 }
 737
 738  //__________________________________________________________________________
 739 void AliMUONTrackParam::ExtrapOneStepHelix3(Double_t field, Double_t step,
 740                                                Double_t *vect, Double_t *vout) const
 741 {
 742 ///     ******************************************************************
 743 ///     *                                                                *
 744 ///     *       Tracking routine in a constant field oriented            *
 745 ///     *       along axis 3                                             *
 746 ///     *       Tracking is performed with a conventional                *
 747 ///     *       helix step method                                        *
 748 ///     *                                                                *
 749 ///     *    ==>Called by : <USER>, GUSWIM                               *
 750 ///     *       Authors    R.Brun, M.Hansroul  *********                 *
 751 ///     *       Rewritten  V.Perevoztchikov
 752 ///     *                                                                *
 753 ///     ******************************************************************
 754
 755     Double_t hxp[3];
 756     Double_t h4, hp, rho, tet;
 757     Double_t sint, sintt, tsint, cos1t;
 758     Double_t f1, f2, f3, f4, f5, f6;
 759
 760     const Int_t kix  = 0;
 761     const Int_t kiy  = 1;
 762     const Int_t kiz  = 2;
 763     const Int_t kipx = 3;
 764     const Int_t kipy = 4;
 765     const Int_t kipz = 5;
 766     const Int_t kipp = 6;
 767
 768     const Double_t kec = 2.9979251e-4;
 769
 770 //
 771 //     ------------------------------------------------------------------
 772 //
 773 //       units are kgauss,centimeters,gev/c
 774 //
 775     vout[kipp] = vect[kipp];
 776     h4 = field * kec;
 777
 778     hxp[0] = - vect[kipy];
 779     hxp[1] = + vect[kipx];
 780
 781     hp = vect[kipz];
 782
 783     rho = -h4/vect[kipp];
 784     tet = rho * step;
 785     if (TMath::Abs(tet) > 0.15) {
 786       sint = TMath::Sin(tet);
 787       sintt = (sint/tet);
 788       tsint = (tet-sint)/tet;
 789       cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
 790     } else {
 791       tsint = tet*tet/36.;
 792       sintt = (1. - tsint);
 793       sint = tet*sintt;
 794       cos1t = 0.5*tet;
 795     }
 796
 797     f1 = step * sintt;
 798     f2 = step * cos1t;
 799     f3 = step * tsint * hp;
 800     f4 = -tet*cos1t;
 801     f5 = sint;
 802     f6 = tet * cos1t * hp;
 803
 804     vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
 805     vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
 806     vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
 807
 808     vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
 809     vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
 810     vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
 811
 812     return;
 813 }
 814  //__________________________________________________________________________
 815 void AliMUONTrackParam::ExtrapOneStepRungekutta(Double_t charge, Double_t step,
 816                                                      Double_t* vect, Double_t* vout) const
 817 {
 818 ///     ******************************************************************
 819 ///     *                                                                *
 820 ///     *  Runge-Kutta method for tracking a particle through a magnetic *
 821 ///     *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
 822 ///     *  Standards, procedure 25.5.20)                                 *
 823 ///     *                                                                *
 824 ///     *  Input parameters                                              *
 825 ///     *       CHARGE    Particle charge                                *
 826 ///     *       STEP      Step size                                      *
 827 ///     *       VECT      Initial co-ords,direction cosines,momentum     *
 828 ///     *  Output parameters                                             *
 829 ///     *       VOUT      Output co-ords,direction cosines,momentum      *
 830 ///     *  User routine called                                           *
 831 ///     *       CALL GUFLD(X,F)                                          *
 832 ///     *                                                                *
 833 ///     *    ==>Called by : <USER>, GUSWIM                               *
 834 ///     *       Authors    R.Brun, M.Hansroul  *********                 *
 835 ///     *                  V.Perevoztchikov (CUT STEP implementation)    *
 836 ///     *                                                                *
 837 ///     *                                                                *
 838 ///     ******************************************************************
 839
 840     Double_t h2, h4, f[4];
 841     Double_t xyzt[3], a, b, c, ph,ph2;
 842     Double_t secxs[4],secys[4],seczs[4],hxp[3];
 843     Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
 844     Double_t est, at, bt, ct, cba;
 845     Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
 846
 847     Double_t x;
 848     Double_t y;
 849     Double_t z;
 850
 851     Double_t xt;
 852     Double_t yt;
 853     Double_t zt;
 854
 855     Double_t maxit = 1992;
 856     Double_t maxcut = 11;
 857
 858     const Double_t kdlt   = 1e-4;
 859     const Double_t kdlt32 = kdlt/32.;
 860     const Double_t kthird = 1./3.;
 861     const Double_t khalf  = 0.5;
 862     const Double_t kec = 2.9979251e-4;
 863
 864     const Double_t kpisqua = 9.86960440109;
 865     const Int_t kix  = 0;
 866     const Int_t kiy  = 1;
 867     const Int_t kiz  = 2;
 868     const Int_t kipx = 3;
 869     const Int_t kipy = 4;
 870     const Int_t kipz = 5;
 871
 872     // *.
 873     // *.    ------------------------------------------------------------------
 874     // *.
 875     // *             this constant is for units cm,gev/c and kgauss
 876     // *
 877     Int_t iter = 0;
 878     Int_t ncut = 0;
 879     for(Int_t j = 0; j < 7; j++)
 880       vout[j] = vect[j];
 881
 882     Double_t  pinv   = kec * charge / vect[6];
 883     Double_t tl = 0.;
 884     Double_t h = step;
 885     Double_t rest;
 886
 887
 888     do {
 889       rest  = step - tl;
 890       if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
 891       //cmodif: call gufld(vout,f) changed into:
 892
 893       GetField(vout,f);
 894
 895       // *
 896       // *             start of integration
 897       // *
 898       x      = vout[0];
 899       y      = vout[1];
 900       z      = vout[2];
 901       a      = vout[3];
 902       b      = vout[4];
 903       c      = vout[5];
 904
 905       h2     = khalf * h;
 906       h4     = khalf * h2;
 907       ph     = pinv * h;
 908       ph2    = khalf * ph;
 909       secxs[0] = (b * f[2] - c * f[1]) * ph2;
 910       secys[0] = (c * f[0] - a * f[2]) * ph2;
 911       seczs[0] = (a * f[1] - b * f[0]) * ph2;
 912       ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
 913       if (ang2 > kpisqua) break;
 914
 915       dxt    = h2 * a + h4 * secxs[0];
 916       dyt    = h2 * b + h4 * secys[0];
 917       dzt    = h2 * c + h4 * seczs[0];
 918       xt     = x + dxt;
 919       yt     = y + dyt;
 920       zt     = z + dzt;
 921       // *
 922       // *              second intermediate point
 923       // *
 924
 925       est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
 926       if (est > h) {
 927         if (ncut++ > maxcut) break;
 928         h *= khalf;
 929         continue;
 930       }
 931
 932       xyzt[0] = xt;
 933       xyzt[1] = yt;
 934       xyzt[2] = zt;
 935
 936       //cmodif: call gufld(xyzt,f) changed into:
 937       GetField(xyzt,f);
 938
 939       at     = a + secxs[0];
 940       bt     = b + secys[0];
 941       ct     = c + seczs[0];
 942
 943       secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
 944       secys[1] = (ct * f[0] - at * f[2]) * ph2;
 945       seczs[1] = (at * f[1] - bt * f[0]) * ph2;
 946       at     = a + secxs[1];
 947       bt     = b + secys[1];
 948       ct     = c + seczs[1];
 949       secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
 950       secys[2] = (ct * f[0] - at * f[2]) * ph2;
 951       seczs[2] = (at * f[1] - bt * f[0]) * ph2;
 952       dxt    = h * (a + secxs[2]);
 953       dyt    = h * (b + secys[2]);
 954       dzt    = h * (c + seczs[2]);
 955       xt     = x + dxt;
 956       yt     = y + dyt;
 957       zt     = z + dzt;
 958       at     = a + 2.*secxs[2];
 959       bt     = b + 2.*secys[2];
 960       ct     = c + 2.*seczs[2];
 961
 962       est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
 963       if (est > 2.*TMath::Abs(h)) {
 964         if (ncut++ > maxcut) break;
 965         h *= khalf;
 966         continue;
 967       }
 968
 969       xyzt[0] = xt;
 970       xyzt[1] = yt;
 971       xyzt[2] = zt;
 972
 973       //cmodif: call gufld(xyzt,f) changed into:
 974       GetField(xyzt,f);
 975
 976       z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
 977       y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
 978       x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
 979
 980       secxs[3] = (bt*f[2] - ct*f[1])* ph2;
 981       secys[3] = (ct*f[0] - at*f[2])* ph2;
 982       seczs[3] = (at*f[1] - bt*f[0])* ph2;
 983       a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
 984       b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
 985       c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
 986
 987       est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
 988         + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
 989         + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
 990
 991       if (est > kdlt && TMath::Abs(h) > 1.e-4) {
 992         if (ncut++ > maxcut) break;
 993         h *= khalf;
 994         continue;
 995       }
 996
 997       ncut = 0;
 998       // *               if too many iterations, go to helix
 999       if (iter++ > maxit) break;
1000
1001       tl += h;
1002       if (est < kdlt32)
1003         h *= 2.;
1004       cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
1005       vout[0] = x;
1006       vout[1] = y;
1007       vout[2] = z;
1008       vout[3] = cba*a;
1009       vout[4] = cba*b;
1010       vout[5] = cba*c;
1011       rest = step - tl;
1012       if (step < 0.) rest = -rest;
1013       if (rest < 1.e-5*TMath::Abs(step)) return;
1014
1015     } while(1);
1016
1017     // angle too big, use helix
1018
1019     f1  = f[0];
1020     f2  = f[1];
1021     f3  = f[2];
1022     f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1023     rho = -f4*pinv;
1024     tet = rho * step;
1025
1026     hnorm = 1./f4;
1027     f1 = f1*hnorm;
1028     f2 = f2*hnorm;
1029     f3 = f3*hnorm;
1030
1031     hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1032     hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1033     hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1034
1035     hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1036
1037     rho1 = 1./rho;
1038     sint = TMath::Sin(tet);
1039     cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1040
1041     g1 = sint*rho1;
1042     g2 = cost*rho1;
1043     g3 = (tet-sint) * hp*rho1;
1044     g4 = -cost;
1045     g5 = sint;
1046     g6 = cost * hp;
1047
1048     vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1049     vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1050     vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1051
1052     vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1053     vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1054     vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1055
1056     return;
1057 }
1058 //___________________________________________________________
1059  void  AliMUONTrackParam::GetField(Double_t *Position, Double_t *Field) const
1060 {
1061     /// interface for arguments in double precision (Why ? ChF)
1062     Float_t x[3], b[3];
1063
1064     x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2];
1065
1066     fkField->Field(x, b);
1067     Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2];
1068
1069     return;
1070   }
1071 //_____________________________________________-
1072 void AliMUONTrackParam::Print(Option_t* opt) const
1073 {
1074   /// Printing TrackParam information
1075   /// "full" option for printing all the information about the TrackParam
1076   TString sopt(opt);
1077   sopt.ToUpper();
1078
1079   if ( sopt.Contains("FULL") ) {
1080     cout << "<AliMUONTrackParam> Bending P=" << setw(5) << setprecision(3)  << 1./GetInverseBendingMomentum() <<
1081       ", NonBendSlope=" << setw(5) << setprecision(3)  << GetNonBendingSlope()*180./TMath::Pi() <<
1082       ", BendSlope=" << setw(5) << setprecision(3)     << GetBendingSlope()*180./TMath::Pi()  <<
1083       ", (x,y,z)_IP=(" <<  setw(5) << setprecision(3) << GetNonBendingCoor() <<
1084       "," <<  setw(5) << setprecision(3) << GetBendingCoor() <<
1085       "," <<  setw(5) << setprecision(3) << GetZ() <<
1086       ") cm, (px,py,pz)=(" << setw(5) << setprecision(3) << Px() <<
1087       "," << setw(5) << setprecision(3) << Py() <<
1088       "," << setw(5) << setprecision(3) << Pz() << ") GeV/c" << endl;
1089   }
1090   else {
1091     cout << "<AliMUONTrackParam>"  << endl;
1092   }
1093
1094 }