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