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