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