1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.2 2000/06/15 07:58:49 morsch
19 Code from MUON-dev joined
21 Revision 1.1.2.3 2000/06/09 21:03:09 morsch
22 Make includes consistent with new file structure.
24 Revision 1.1.2.2 2000/06/09 12:58:05 gosset
25 Removed comment beginnings in Log sections of .cxx files
26 Suppressed most violations of coding rules
28 Revision 1.1.2.1 2000/06/07 14:44:53 gosset
29 Addition of files for track reconstruction in C++
32 //__________________________________________________________________________
34 // Track parameters in ALICE dimuon spectrometer
35 //__________________________________________________________________________
39 #include "AliCallf77.h"
41 #include "AliMUONHitForRec.h"
42 #include "AliMUONSegment.h"
43 #include "AliMUONTrackParam.h"
44 #include "AliMUONChamber.h"
47 ClassImp(AliMUONTrackParam) // Class implementation in ROOT context
50 # define reco_ghelix reco_ghelix_
52 # define reco_ghelix RECO_GHELIX
57 void type_of_call reco_ghelix(Double_t &Charge, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
60 // Inline functions for Get and Set: inline removed because it does not work !!!!
61 Double_t AliMUONTrackParam::GetInverseBendingMomentum(void) {
62 // Get fInverseBendingMomentum
63 return fInverseBendingMomentum;}
64 void AliMUONTrackParam::SetInverseBendingMomentum(Double_t InverseBendingMomentum) {
65 // Set fInverseBendingMomentum
66 fInverseBendingMomentum = InverseBendingMomentum;}
67 Double_t AliMUONTrackParam::GetBendingSlope(void) {
69 return fBendingSlope;}
70 void AliMUONTrackParam::SetBendingSlope(Double_t BendingSlope) {
72 fBendingSlope = BendingSlope;}
73 Double_t AliMUONTrackParam::GetNonBendingSlope(void) {
74 // Get fNonBendingSlope
75 return fNonBendingSlope;}
76 void AliMUONTrackParam::SetNonBendingSlope(Double_t NonBendingSlope) {
77 // Set fNonBendingSlope
78 fNonBendingSlope = NonBendingSlope;}
79 Double_t AliMUONTrackParam::GetZ(void) {
82 void AliMUONTrackParam::SetZ(Double_t Z) {
85 Double_t AliMUONTrackParam::GetBendingCoor(void) {
88 void AliMUONTrackParam::SetBendingCoor(Double_t BendingCoor) {
90 fBendingCoor = BendingCoor;}
91 Double_t AliMUONTrackParam::GetNonBendingCoor(void) {
92 // Get fNonBendingCoor
93 return fNonBendingCoor;}
94 void AliMUONTrackParam::SetNonBendingCoor(Double_t NonBendingCoor) {
95 // Set fNonBendingCoor
96 fNonBendingCoor = NonBendingCoor;}
98 //__________________________________________________________________________
99 void AliMUONTrackParam::ExtrapToZ(Double_t Z)
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
120 forwardBackward * TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
121 Double_t stepLength = 6.0; // in parameter ????
122 // Extrapolation loop
124 while (((forwardBackward * (vGeant3[2] - Z)) <= 0.0) &&
125 (stepNumber < maxStepNumber)) {
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];}
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;
149 ((vGeant3New[3] / vGeant3New[5]) - (vGeant3[3] / vGeant3[5])) / dZ12;
150 Double_t yPrime = (vGeant3New[1] - vGeant3[1]) / dZ12;
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
156 Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
157 Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
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);
166 //__________________________________________________________________________
167 void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward)
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);
178 pYZ / TMath::Sqrt(1.0 + this->fBendingSlope * this->fBendingSlope);
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
187 //__________________________________________________________________________
188 void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge)
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];
202 //__________________________________________________________________________
203 void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam)
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;
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;}
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;
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]);
235 //__________________________________________________________________________
236 void AliMUONTrackParam::ExtrapToVertex()
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
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)
250 //__________________________________________________________________________
251 void AliMUONTrackParam::BransonCorrection()
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;
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)
263 47.26, // Polyethylene (cm)
265 47.26, // Polyethylene (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)
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)
284 for (Int_t iBound = 0; iBound < 9; iBound++) {
286 (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
287 z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound];
289 (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound];
291 zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
294 for (Int_t iBound = 0; iBound < 3; iBound++) {
296 (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
297 z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound];
299 (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
301 zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
302 rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
305 pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
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;
316 if (radiusEndAbsorber2 > rLimit*rLimit) {
317 zEndAbsorber = z1[9];
320 zEndAbsorber = z2[2];
324 xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
325 yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
327 // new parameters after Branson and energy loss corrections
328 pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
331 fBendingSlope = pY / pZ;
332 fNonBendingSlope = pX / pZ;
334 pT = TMath::Sqrt(pX * pX + pY * pY);
335 theta = TMath::ATan2(pT, pZ);
337 TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
339 fInverseBendingMomentum = (sign / pTotal) *
341 fBendingSlope * fBendingSlope +
342 fNonBendingSlope * fNonBendingSlope) /
343 TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
345 // vertex position at (0,0,0)
346 // should be taken from vertex measurement ???
352 //__________________________________________________________________________
353 Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
355 // Returns the total momentum corrected from energy loss in the front absorber
356 Double_t deltaP, pTotalCorrected;
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) {
364 deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
366 deltaP = 3.0643 + 0.01346 *pTotal;
370 deltaP = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
372 deltaP = 2.407 + 0.00702 * pTotal;
375 pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
376 return pTotalCorrected;