Changes to EventReconstructor...:
[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
16/*
17$Log$
04b5ea16 18Revision 1.2 2000/06/15 07:58:49 morsch
19Code from MUON-dev joined
20
a9e2aefa 21Revision 1.1.2.3 2000/06/09 21:03:09 morsch
22Make includes consistent with new file structure.
23
24Revision 1.1.2.2 2000/06/09 12:58:05 gosset
25Removed comment beginnings in Log sections of .cxx files
26Suppressed most violations of coding rules
27
28Revision 1.1.2.1 2000/06/07 14:44:53 gosset
29Addition of files for track reconstruction in C++
30*/
31
32//__________________________________________________________________________
33//
34// Track parameters in ALICE dimuon spectrometer
35//__________________________________________________________________________
36
37#include <iostream.h>
38
39#include "AliCallf77.h"
40#include "AliMUON.h"
41#include "AliMUONHitForRec.h"
42#include "AliMUONSegment.h"
43#include "AliMUONTrackParam.h"
44#include "AliMUONChamber.h"
45#include "AliRun.h"
46
47ClassImp(AliMUONTrackParam) // Class implementation in ROOT context
48
49#ifndef WIN32
50# define reco_ghelix reco_ghelix_
51#else
52# define reco_ghelix RECO_GHELIX
53#endif
54
55extern "C"
56{
57void type_of_call reco_ghelix(Double_t &Charge, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
58}
59
60// Inline functions for Get and Set: inline removed because it does not work !!!!
61Double_t AliMUONTrackParam::GetInverseBendingMomentum(void) {
62 // Get fInverseBendingMomentum
63 return fInverseBendingMomentum;}
64void AliMUONTrackParam::SetInverseBendingMomentum(Double_t InverseBendingMomentum) {
65 // Set fInverseBendingMomentum
66 fInverseBendingMomentum = InverseBendingMomentum;}
67Double_t AliMUONTrackParam::GetBendingSlope(void) {
68 // Get fBendingSlope
69 return fBendingSlope;}
70void AliMUONTrackParam::SetBendingSlope(Double_t BendingSlope) {
71 // Set fBendingSlope
72 fBendingSlope = BendingSlope;}
73Double_t AliMUONTrackParam::GetNonBendingSlope(void) {
74 // Get fNonBendingSlope
75 return fNonBendingSlope;}
76void AliMUONTrackParam::SetNonBendingSlope(Double_t NonBendingSlope) {
77 // Set fNonBendingSlope
78 fNonBendingSlope = NonBendingSlope;}
79Double_t AliMUONTrackParam::GetZ(void) {
80 // Get fZ
81 return fZ;}
82void AliMUONTrackParam::SetZ(Double_t Z) {
83 // Set fZ
84 fZ = Z;}
85Double_t AliMUONTrackParam::GetBendingCoor(void) {
86 // Get fBendingCoor
87 return fBendingCoor;}
88void AliMUONTrackParam::SetBendingCoor(Double_t BendingCoor) {
89 // Set fBendingCoor
90 fBendingCoor = BendingCoor;}
91Double_t AliMUONTrackParam::GetNonBendingCoor(void) {
92 // Get fNonBendingCoor
93 return fNonBendingCoor;}
94void AliMUONTrackParam::SetNonBendingCoor(Double_t NonBendingCoor) {
95 // Set fNonBendingCoor
96 fNonBendingCoor = NonBendingCoor;}
97
98 //__________________________________________________________________________
99void AliMUONTrackParam::ExtrapToZ(Double_t Z)
100{
101 // Track parameter extrapolation to the plane at "Z".
102 // On return, the track parameters resulting from the extrapolation
103 // replace the current track parameters.
104 // Use "reco_ghelix" which should be replaced by something else !!!!
105 if (this->fZ == Z) return; // nothing to be done if same Z
106 Double_t forwardBackward; // +1 if forward, -1 if backward
107 if (Z > this->fZ) forwardBackward = 1.0;
108 else forwardBackward = -1.0;
109 Double_t temp, vGeant3[7], vGeant3New[7]; // 7 in parameter ????
110 Int_t iGeant3, stepNumber;
111 Int_t maxStepNumber = 5000; // in parameter ????
112 // For safety: return kTRUE or kFALSE ????
113 // Parameter vector for calling GHELIX in Geant3
114 SetGeant3Parameters(vGeant3, forwardBackward);
115 // For use of reco_ghelix...: invert X and Y, PX/PTOT and PY/PTOT !!!!
116 temp = vGeant3[0]; vGeant3[0] = vGeant3[1]; vGeant3[1] = temp;
117 temp = vGeant3[3]; vGeant3[3] = vGeant3[4]; vGeant3[4] = temp;
118 // charge must be changed with momentum for backward motion
119 Double_t charge =
120 forwardBackward * TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
121 Double_t stepLength = 6.0; // in parameter ????
122 // Extrapolation loop
123 stepNumber = 0;
124 while (((forwardBackward * (vGeant3[2] - Z)) <= 0.0) &&
125 (stepNumber < maxStepNumber)) {
126 stepNumber++;
127 // call Geant3 "ghelix" subroutine through a copy in "reco_muon.F":
128 // the true function should be called, but how ???? and remove prototyping ...
129 reco_ghelix(charge, stepLength, vGeant3, vGeant3New);
130 if ((forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z
131 // better use TArray ????
132 for (iGeant3 = 0; iGeant3 < 7; iGeant3++)
133 {vGeant3[iGeant3] = vGeant3New[iGeant3];}
134 }
135 // check maxStepNumber ????
136 // For use of reco_ghelix...:
137 // invert back X and Y, PX/PTOT and PY/PTOT, both for vGeant3 and vGeant3New !!!!
138 temp = vGeant3[0]; vGeant3[0] = vGeant3[1]; vGeant3[1] = temp;
139 temp = vGeant3New[0]; vGeant3New[0] = vGeant3New[1]; vGeant3New[1] = temp;
140 temp = vGeant3[3]; vGeant3[3] = vGeant3[4]; vGeant3[4] = temp;
141 temp = vGeant3New[3]; vGeant3New[3] = vGeant3New[4]; vGeant3New[4] = temp;
142 // Interpolation back to exact Z (2nd order)
143 // should be in function ???? using TArray ????
144 Double_t dZ12 = vGeant3New[2] - vGeant3[2]; // 1->2
145 Double_t dZ1i = Z - vGeant3[2]; // 1-i
146 Double_t dZi2 = vGeant3New[2] - Z; // i->2
147 Double_t xPrime = (vGeant3New[0] - vGeant3[0]) / dZ12;
148 Double_t xSecond =
149 ((vGeant3New[3] / vGeant3New[5]) - (vGeant3[3] / vGeant3[5])) / dZ12;
150 Double_t yPrime = (vGeant3New[1] - vGeant3[1]) / dZ12;
151 Double_t ySecond =
152 ((vGeant3New[4] / vGeant3New[5]) - (vGeant3[4] / vGeant3[5])) / dZ12;
153 vGeant3[0] = vGeant3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
154 vGeant3[1] = vGeant3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
155 vGeant3[2] = Z; // Z
156 Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
157 Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
158 vGeant3[5] =
159 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
160 vGeant3[3] = xPrimeI * vGeant3[5]; // PX/PTOT
161 vGeant3[4] = yPrimeI * vGeant3[5]; // PY/PTOT
162 // Track parameters from Geant3 parameters
163 GetFromGeant3Parameters(vGeant3, charge);
164}
165
166 //__________________________________________________________________________
167void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward)
168{
169 // Set vector of Geant3 parameters pointed to by "VGeant3"
170 // from track parameters in current AliMUONTrackParam.
171 // Since AliMUONTrackParam is only geometry, one uses "ForwardBackward"
172 // to know whether the particle is going forward (+1) or backward (-1).
173 VGeant3[0] = this->fNonBendingCoor; // X
174 VGeant3[1] = this->fBendingCoor; // Y
175 VGeant3[2] = this->fZ; // Z
176 Double_t pYZ = TMath::Abs(1.0 / this->fInverseBendingMomentum);
177 Double_t pZ =
178 pYZ / TMath::Sqrt(1.0 + this->fBendingSlope * this->fBendingSlope);
179 VGeant3[6] =
180 TMath::Sqrt(pYZ * pYZ +
181 pZ * pZ * this->fNonBendingSlope * this->fNonBendingSlope); // PTOT
182 VGeant3[5] = ForwardBackward * pZ / VGeant3[6]; // PZ/PTOT
183 VGeant3[3] = this->fNonBendingSlope * VGeant3[5]; // PX/PTOT
184 VGeant3[4] = this->fBendingSlope * VGeant3[5]; // PY/PTOT
185}
186
187 //__________________________________________________________________________
188void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge)
189{
190 // Get track parameters in current AliMUONTrackParam
191 // from Geant3 parameters pointed to by "VGeant3".
192 // "InverseBendingMomentum" is signed with "Charge".
193 this->fNonBendingCoor = VGeant3[0]; // X
194 this->fBendingCoor = VGeant3[1]; // Y
195 this->fZ = VGeant3[2]; // Z
196 Double_t pYZ = VGeant3[6] * TMath::Sqrt(1.0 - VGeant3[3] * VGeant3[3]);
197 this->fInverseBendingMomentum = Charge / pYZ;
198 this->fBendingSlope = VGeant3[4] / VGeant3[5];
199 this->fNonBendingSlope = VGeant3[3] / VGeant3[5];
200}
201
202 //__________________________________________________________________________
203void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam)
204{
205 // Track parameters extrapolated from current track parameters ("this")
206 // to both chambers of the station(0..) "Station"
207 // are returned in the array (dimension 2) of track parameters
208 // pointed to by "TrackParam" (index 0 and 1 for first and second chambers).
209 Double_t extZ[2], z1, z2;
210 Int_t i1, i2;
211 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
212 // range of Station to be checked ????
213 z1 = (&(pMUON->Chamber(2 * Station)))->Z(); // Z of first chamber
214 z2 = (&(pMUON->Chamber(2 * Station + 1)))->Z(); // Z of second chamber
215 // First and second Z to extrapolate at
216 if ((z1 > this->fZ) && (z2 > this->fZ)) {i1 = 0; i2 = 1;}
217 else if ((z1 < this->fZ) && (z2 < this->fZ)) {i1 = 1; i2 = 0;}
218 else {
219 cout << "ERROR in AliMUONTrackParam::CreateExtrapSegmentInStation" << endl;
220 cout << "Starting Z (" << this->fZ << ") in between z1 (" << z1 <<
221 ") and z2 (" << z2 << ") of station(0..) " << Station << endl;
222 }
223 extZ[i1] = z1;
224 extZ[i2] = z2;
225 // copy of track parameters
226 TrackParam[i1] = *this;
227 // first extrapolation
228 (&(TrackParam[i1]))->ExtrapToZ(extZ[0]);
229 TrackParam[i2] = TrackParam[i1];
230 // second extrapolation
231 (&(TrackParam[i2]))->ExtrapToZ(extZ[1]);
232 return;
233}
234
04b5ea16 235 //__________________________________________________________________________
236void AliMUONTrackParam::ExtrapToVertex()
237{
238 // Extrapolation to the vertex.
239 // Returns the track parameters resulting from the extrapolation,
240 // in the current TrackParam.
241 // Changes parameters according to branson correction through the absorber
242
243 Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
244 // Extrapolates track parameters upstream to the "Z" end of the front absorber
245 ExtrapToZ(zAbsorber);
246 // Makes Branson correction (multiple scattering + energy loss)
247 BransonCorrection();
248}
249
250 //__________________________________________________________________________
251void AliMUONTrackParam::BransonCorrection()
252{
253 // Branson correction of track parameters
254 // the entry parameters have to be calculated at the end of the absorber
255 Double_t zEndAbsorber, zBP, xBP, yBP;
256 Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
257 Int_t sign;
258 // Would it be possible to calculate all that from Geant configuration ????
259 // Radiation lengths outer part theta > 3 degres
260 static Double_t x01[9] = { 18.8, // C (cm)
261 10.397, // Concrete (cm)
262 0.56, // Plomb (cm)
263 47.26, // Polyethylene (cm)
264 0.56, // Plomb (cm)
265 47.26, // Polyethylene (cm)
266 0.56, // Plomb (cm)
267 47.26, // Polyethylene (cm)
268 0.56 }; // Plomb (cm)
269 // inner part theta < 3 degres
270 static Double_t x02[3] = { 18.8, // C (cm)
271 10.397, // Concrete (cm)
272 0.35 }; // W (cm)
273 // z positions of the materials inside the absober outer part theta > 3 degres
274 static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
275 // inner part theta < 3 degres
276 static Double_t z2[4] = { 90, 315, 467, 503 };
277 static Bool_t first = kTRUE;
278 static Double_t zBP1, zBP2, rLimit;
279 // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
280 if (first) {
281 first = kFALSE;
282 Double_t aNBP = 0.0;
283 Double_t aDBP = 0.0;
284 for (Int_t iBound = 0; iBound < 9; iBound++) {
285 aNBP = aNBP +
286 (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
287 z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound];
288 aDBP = aDBP +
289 (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound];
290 }
291 zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
292 aNBP = 0.0;
293 aDBP = 0.0;
294 for (Int_t iBound = 0; iBound < 3; iBound++) {
295 aNBP = aNBP +
296 (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
297 z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound];
298 aDBP = aDBP +
299 (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
300 }
301 zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
302 rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
303 }
304
305 pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
306 sign = 1;
307 if (fInverseBendingMomentum < 0) sign = -1;
308 pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
309 pX = pZ * fNonBendingSlope;
310 pY = pZ * fBendingSlope;
311 pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
312 xEndAbsorber = fNonBendingCoor;
313 yEndAbsorber = fBendingCoor;
314 radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
315
316 if (radiusEndAbsorber2 > rLimit*rLimit) {
317 zEndAbsorber = z1[9];
318 zBP = zBP1;
319 } else {
320 zEndAbsorber = z2[2];
321 zBP = zBP2;
322 }
323
324 xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
325 yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
326
327 // new parameters after Branson and energy loss corrections
328 pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
329 pX = pZ * xBP / zBP;
330 pY = pZ * yBP / zBP;
331 fBendingSlope = pY / pZ;
332 fNonBendingSlope = pX / pZ;
333
334 pT = TMath::Sqrt(pX * pX + pY * pY);
335 theta = TMath::ATan2(pT, pZ);
336 pTotal =
337 TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
338
339 fInverseBendingMomentum = (sign / pTotal) *
340 TMath::Sqrt(1.0 +
341 fBendingSlope * fBendingSlope +
342 fNonBendingSlope * fNonBendingSlope) /
343 TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
344
345 // vertex position at (0,0,0)
346 // should be taken from vertex measurement ???
347 fBendingCoor = 0.0;
348 fNonBendingCoor = 0;
349 fZ= 0;
350}
351
352 //__________________________________________________________________________
353Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
354{
355 // Returns the total momentum corrected from energy loss in the front absorber
356 Double_t deltaP, pTotalCorrected;
357
358 Double_t radiusEndAbsorber2 =
359 xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
360 // Parametrization to be redone according to change of absorber material ????
361 // The name is not so good, and there are many arguments !!!!
362 if (radiusEndAbsorber2 < rLimit * rLimit) {
363 if (pTotal < 15) {
364 deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
365 } else {
366 deltaP = 3.0643 + 0.01346 *pTotal;
367 }
368 } else {
369 if (pTotal < 15) {
370 deltaP = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
371 } else {
372 deltaP = 2.407 + 0.00702 * pTotal;
373 }
374 }
375 pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
376 return pTotalCorrected;
377}
378