| 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 | |
| c04e3238 |
31 | #include "AliMUONTrackExtrap.h" |
| 32 | #include "AliMUONTrackParam.h" |
| 33 | #include "AliMUONConstants.h" |
| 8cde4af5 |
34 | |
| c04e3238 |
35 | #include "AliMagF.h" |
| 8cde4af5 |
36 | |
| 37 | #include <Riostream.h> |
| 38 | #include <TMath.h> |
| 39 | #include <TMatrixD.h> |
| 40 | #include <TGeoManager.h> |
| c04e3238 |
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 | { |
| 8cde4af5 |
338 | /// Extrapolation to the vertex (at the z position "zVtx") without Branson and energy loss corrections. |
| 208f139e |
339 | /// Returns the track parameters resulting from the extrapolation in the current TrackParam. |
| 340 | /// Include multiple Coulomb scattering effects in trackParam covariances. |
| 341 | |
| 8cde4af5 |
342 | if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex |
| 343 | |
| 208f139e |
344 | if (trackParam->GetZ() > zVtx) { // spectro. (z<0) |
| 345 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<<trackParam->GetZ() |
| 346 | <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl; |
| 347 | exit(-1); |
| 348 | } |
| 349 | |
| 8cde4af5 |
350 | // Check whether the geometry is available and get absorber boundaries |
| 351 | if (!gGeoManager) { |
| 352 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl; |
| 353 | return; |
| 354 | } |
| 355 | TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1"); |
| 356 | if (!absNode) { |
| 357 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<<endl; |
| 358 | return; |
| 359 | } |
| 360 | Double_t zAbsBeg, zAbsEnd; |
| 361 | absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd); |
| 362 | const Double_t *absPos = absNode->GetMatrix()->GetTranslation(); |
| 363 | zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0) |
| 364 | zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0) |
| 365 | |
| 366 | // Check the vertex position relatively to the absorber |
| 367 | if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0) |
| 368 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx |
| 369 | <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl; |
| 370 | } else if (zVtx < zAbsEnd ) { // spectro. (z<0) |
| 371 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx |
| 372 | <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl; |
| 208f139e |
373 | ExtrapToZCov(trackParam,zVtx); |
| 374 | return; |
| 375 | } |
| 376 | |
| 8cde4af5 |
377 | // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed |
| 378 | if (trackParam->GetZ() > zAbsBeg) { // spectro. (z<0) |
| 379 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ() |
| 380 | <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl; |
| 381 | ExtrapToZCov(trackParam,zVtx); |
| 382 | return; |
| 383 | } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0) |
| 384 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ() |
| 385 | <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl; |
| 208f139e |
386 | } else { |
| 8cde4af5 |
387 | ExtrapToZCov(trackParam,zAbsEnd); |
| 208f139e |
388 | } |
| 389 | |
| 8cde4af5 |
390 | // Then add MCS effect in absorber to the parameters covariances |
| 391 | AliMUONTrackParam trackParamIn(*trackParam); |
| 392 | ExtrapToZ(&trackParamIn, TMath::Min(zVtx, zAbsBeg)); |
| 393 | Double_t trackXYZIn[3]; |
| 394 | trackXYZIn[0] = trackParamIn.GetNonBendingCoor(); |
| 395 | trackXYZIn[1] = trackParamIn.GetBendingCoor(); |
| 396 | trackXYZIn[2] = trackParamIn.GetZ(); |
| 397 | Double_t trackXYZOut[3]; |
| 398 | trackXYZOut[0] = trackParam->GetNonBendingCoor(); |
| 399 | trackXYZOut[1] = trackParam->GetBendingCoor(); |
| 400 | trackXYZOut[2] = trackParam->GetZ(); |
| 401 | Double_t pathLength = 0.; |
| 402 | Double_t f0 = 0.; |
| 403 | Double_t f1 = 0.; |
| 404 | Double_t f2 = 0.; |
| 405 | Double_t meanRho = 0.; |
| 406 | GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho); |
| 407 | AddMCSEffectInAbsorber(trackParam,pathLength,f0,f1,f2); |
| 208f139e |
408 | |
| 409 | // finally go to the vertex |
| 410 | ExtrapToZCov(trackParam,zVtx); |
| 411 | |
| 412 | } |
| 413 | |
| 414 | //__________________________________________________________________________ |
| 8cde4af5 |
415 | void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2) |
| 416 | { |
| 417 | /// Add to the track parameter covariances the effects of multiple Coulomb scattering |
| 418 | /// at the end of the front absorber using the absorber correction parameters |
| 419 | |
| 420 | // absorber related covariance parameters |
| 421 | Double_t bendingSlope = param->GetBendingSlope(); |
| 422 | Double_t nonBendingSlope = param->GetNonBendingSlope(); |
| 423 | Double_t inverseBendingMomentum = param->GetInverseBendingMomentum(); |
| 424 | Double_t alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) / |
| 425 | (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1 |
| 426 | Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2); |
| 427 | Double_t covCorrSlope = alpha2 * (pathLength * f0 - f1); |
| 428 | Double_t varSlop = alpha2 * f0; |
| 429 | |
| 430 | TMatrixD* paramCov = param->GetCovariances(); |
| 431 | // Non bending plane |
| 432 | (*paramCov)(0,0) += varCoor; (*paramCov)(0,1) += covCorrSlope; |
| 433 | (*paramCov)(1,0) += covCorrSlope; (*paramCov)(1,1) += varSlop; |
| 434 | // Bending plane |
| 435 | (*paramCov)(2,2) += varCoor; (*paramCov)(2,3) += covCorrSlope; |
| 436 | (*paramCov)(3,2) += covCorrSlope; (*paramCov)(3,3) += varSlop; |
| 437 | |
| 438 | } |
| 439 | |
| 440 | //__________________________________________________________________________ |
| 441 | void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t &pathLength, |
| 442 | Double_t &f0, Double_t &f1, Double_t &f2, Double_t &meanRho) |
| 443 | { |
| 444 | /// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber |
| 445 | /// Calculated assuming a linear propagation between track positions trackXYZIn and trackXYZOut |
| 446 | // pathLength: path length between trackXYZIn and trackXYZOut (cm) |
| 447 | // f0: 0th moment of z calculated with the inverse radiation-length distribution |
| 448 | // f1: 1st moment of z calculated with the inverse radiation-length distribution |
| 449 | // f2: 2nd moment of z calculated with the inverse radiation-length distribution |
| 450 | // meanRho: average density of crossed material (g/cm3) |
| 451 | |
| 452 | // Reset absorber's parameters |
| 453 | pathLength = 0.; |
| 454 | f0 = 0.; |
| 455 | f1 = 0.; |
| 456 | f2 = 0.; |
| 457 | meanRho = 0.; |
| 458 | |
| 459 | // Check whether the geometry is available |
| 460 | if (!gGeoManager) { |
| 461 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl; |
| 462 | return; |
| 463 | } |
| 464 | |
| 465 | // Initialize starting point and direction |
| 466 | pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+ |
| 467 | (trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+ |
| 468 | (trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2])); |
| 469 | if (pathLength < TGeoShape::Tolerance()) return; |
| 470 | Double_t b[3]; |
| 471 | b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength; |
| 472 | b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength; |
| 473 | b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength; |
| 474 | TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b); |
| 475 | if (!currentnode) { |
| 476 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<<endl; |
| 477 | return; |
| 478 | } |
| 479 | |
| 480 | // loop over absorber slices and calculate absorber's parameters |
| 481 | Double_t rho = 0.; // material density (g/cm3) |
| 482 | Double_t x0 = 0.; // radiation-length (cm-1) |
| 483 | Double_t localPathLength = 0; |
| 484 | Double_t remainingPathLength = pathLength; |
| 485 | Double_t zB = trackXYZIn[2]; |
| 486 | Double_t zE, dzB, dzE; |
| 487 | do { |
| 488 | // Get material properties |
| 489 | TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial(); |
| 490 | rho = material->GetDensity(); |
| 491 | x0 = material->GetRadLen(); |
| 492 | if (!material->IsMixture()) x0 /= rho; // different normalization in the modeler for mixture |
| 493 | |
| 494 | // Get path length within this material |
| 495 | gGeoManager->FindNextBoundary(remainingPathLength); |
| 496 | localPathLength = gGeoManager->GetStep() + 1.e-6; |
| 497 | // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step |
| 498 | if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength; |
| 499 | else { |
| 500 | currentnode = gGeoManager->Step(); |
| 501 | if (!currentnode) { |
| 502 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl; |
| 503 | f0 = f1 = f2 = meanRho = 0.; |
| 504 | return; |
| 505 | } |
| 506 | if (!gGeoManager->IsEntering()) { |
| 507 | // make another small step to try to enter in new absorber slice |
| 508 | gGeoManager->SetStep(0.001); |
| 509 | currentnode = gGeoManager->Step(); |
| 510 | if (!gGeoManager->IsEntering() || !currentnode) { |
| 511 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl; |
| 512 | f0 = f1 = f2 = meanRho = 0.; |
| 513 | return; |
| 514 | } |
| 515 | localPathLength += 0.001; |
| 516 | } |
| 517 | } |
| 518 | |
| 519 | // calculate absorber's parameters |
| 520 | zE = b[2] * localPathLength + zB; |
| 521 | dzB = zB - trackXYZIn[2]; |
| 522 | dzE = zE - trackXYZIn[2]; |
| 523 | f0 += localPathLength / x0; |
| 524 | f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.; |
| 525 | f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.; |
| 526 | meanRho += localPathLength * rho; |
| 527 | |
| 528 | // prepare next step |
| 529 | zB = zE; |
| 530 | remainingPathLength -= localPathLength; |
| 531 | } while (remainingPathLength > TGeoShape::Tolerance()); |
| 532 | |
| 533 | meanRho /= pathLength; |
| 534 | } |
| 535 | |
| 536 | //__________________________________________________________________________ |
| 537 | void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0) |
| 208f139e |
538 | { |
| 539 | /// Add to the track parameter covariances the effects of multiple Coulomb scattering |
| 540 | /// through a material of thickness "dZ" and of radiation length "x0" |
| 541 | /// assuming linear propagation and using the small angle approximation. |
| 542 | |
| 543 | Double_t bendingSlope = param->GetBendingSlope(); |
| 544 | Double_t nonBendingSlope = param->GetNonBendingSlope(); |
| 545 | Double_t inverseTotalMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum() * |
| 546 | (1.0 + bendingSlope * bendingSlope) / |
| 547 | (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); |
| 548 | // Path length in the material |
| 549 | Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope); |
| 550 | Double_t pathLength2 = pathLength * pathLength; |
| 551 | // relativistic velocity |
| 552 | Double_t velo = 1.; |
| 553 | // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory |
| 554 | Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0)); |
| 555 | theta02 *= theta02 * inverseTotalMomentum2 * pathLength / x0; |
| 556 | |
| 557 | // Add effects of multiple Coulomb scattering in track parameter covariances |
| 558 | TMatrixD* paramCov = param->GetCovariances(); |
| 559 | Double_t varCoor = pathLength2 * theta02 / 3.; |
| 560 | Double_t varSlop = theta02; |
| 561 | Double_t covCorrSlope = pathLength * theta02 / 2.; |
| 562 | // Non bending plane |
| 563 | (*paramCov)(0,0) += varCoor; (*paramCov)(0,1) += covCorrSlope; |
| 564 | (*paramCov)(1,0) += covCorrSlope; (*paramCov)(1,1) += varSlop; |
| 565 | // Bending plane |
| 566 | (*paramCov)(2,2) += varCoor; (*paramCov)(2,3) += covCorrSlope; |
| 567 | (*paramCov)(3,2) += covCorrSlope; (*paramCov)(3,3) += varSlop; |
| 568 | |
| c04e3238 |
569 | } |
| 570 | |
| 571 | //__________________________________________________________________________ |
| 572 | void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx) |
| 573 | { |
| 574 | /// Extrapolation to the vertex. |
| 575 | /// Returns the track parameters resulting from the extrapolation in the current TrackParam. |
| 8cde4af5 |
576 | /// Changes parameters according to Branson correction through the absorber and energy loss |
| c04e3238 |
577 | |
| 8cde4af5 |
578 | if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex |
| c04e3238 |
579 | |
| 8cde4af5 |
580 | if (trackParam->GetZ() > zVtx) { // spectro. (z<0) |
| 581 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ() |
| 582 | <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl; |
| 583 | exit(-1); |
| c04e3238 |
584 | } |
| 8cde4af5 |
585 | |
| 586 | // Check whether the geometry is available and get absorber boundaries |
| 587 | if (!gGeoManager) { |
| 588 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl; |
| 589 | return; |
| 590 | } |
| 591 | TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1"); |
| 592 | if (!absNode) { |
| 593 | cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<<endl; |
| 594 | return; |
| 595 | } |
| 596 | Double_t zAbsBeg, zAbsEnd; |
| 597 | absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd); |
| 598 | const Double_t *absPos = absNode->GetMatrix()->GetTranslation(); |
| 599 | zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0) |
| 600 | zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0) |
| 601 | |
| 602 | // Check the vertex position relatively to the absorber |
| 603 | if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0) |
| 604 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx |
| 605 | <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl; |
| 606 | } else if (zVtx < zAbsEnd ) { // spectro. (z<0) |
| 607 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx |
| 608 | <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl; |
| 609 | ExtrapToZ(trackParam,zVtx); |
| 610 | return; |
| 611 | } |
| 612 | |
| 613 | // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed |
| 614 | if (trackParam->GetZ() > zAbsBeg) { // spectro. (z<0) |
| 615 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ() |
| 616 | <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl; |
| 617 | ExtrapToZ(trackParam,zVtx); |
| 618 | return; |
| 619 | } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0) |
| 620 | cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ() |
| 621 | <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl; |
| c04e3238 |
622 | } else { |
| 8cde4af5 |
623 | ExtrapToZ(trackParam,zAbsEnd); |
| c04e3238 |
624 | } |
| c04e3238 |
625 | |
| 8cde4af5 |
626 | // Get absorber correction parameters assuming linear propagation from vertex to the track position |
| 627 | Double_t trackXYZOut[3]; |
| 628 | trackXYZOut[0] = trackParam->GetNonBendingCoor(); |
| 629 | trackXYZOut[1] = trackParam->GetBendingCoor(); |
| 630 | trackXYZOut[2] = trackParam->GetZ(); |
| 631 | Double_t trackXYZIn[3]; |
| 632 | trackXYZIn[2] = TMath::Min(zVtx, zAbsBeg); // spectro. (z<0) |
| 633 | trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]); |
| 634 | trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]); |
| 635 | Double_t pathLength = 0.; |
| 636 | Double_t f0 = 0.; |
| 637 | Double_t f1 = 0.; |
| 638 | Double_t f2 = 0.; |
| 639 | Double_t meanRho = 0.; |
| 640 | GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho); |
| c04e3238 |
641 | |
| 8cde4af5 |
642 | // Calculate energy loss |
| 643 | Double_t pTot = trackParam->P(); |
| 644 | Double_t charge = TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum()); |
| 645 | Double_t deltaP = TotalMomentumEnergyLoss(pTot,pathLength,meanRho); |
| 646 | |
| 647 | // Correct for half of energy loss |
| 648 | pTot += 0.5 * deltaP; |
| 649 | |
| 650 | // Position of the Branson plane (spectro. (z<0)) |
| 651 | Double_t zB = (f1>0.) ? trackXYZIn[2] - f2/f1 : 0.; |
| 652 | |
| 653 | // Get track position in the Branson plane corrected for magnetic field effect |
| 654 | ExtrapToZ(trackParam,zVtx); |
| 655 | Double_t xB = trackParam->GetNonBendingCoor() + (zB - zVtx) * trackParam->GetNonBendingSlope(); |
| 656 | Double_t yB = trackParam->GetBendingCoor() + (zB - zVtx) * trackParam->GetBendingSlope(); |
| 657 | |
| 658 | // Get track slopes corrected for multiple scattering (spectro. (z<0)) |
| 659 | Double_t nonBendingSlope = (zB<0.) ? (xB - xVtx) / (zB - zVtx) : trackParam->GetNonBendingSlope(); |
| 660 | Double_t bendingSlope = (zB<0.) ? (yB - yVtx) / (zB - zVtx) : trackParam->GetBendingSlope(); |
| 661 | |
| 662 | // Correct for second half of energy loss |
| 663 | pTot += 0.5 * deltaP; |
| 664 | |
| 665 | // Set track parameters at vertex |
| 666 | trackParam->SetNonBendingCoor(xVtx); |
| 667 | trackParam->SetBendingCoor(yVtx); |
| c04e3238 |
668 | trackParam->SetZ(zVtx); |
| 8cde4af5 |
669 | trackParam->SetNonBendingSlope(nonBendingSlope); |
| 670 | trackParam->SetBendingSlope(bendingSlope); |
| 671 | trackParam->SetInverseBendingMomentum(charge / pTot * |
| 672 | TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) / |
| 673 | TMath::Sqrt(1.0 + bendingSlope*bendingSlope)); |
| 674 | |
| c04e3238 |
675 | } |
| 676 | |
| 677 | //__________________________________________________________________________ |
| 8cde4af5 |
678 | Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(Double_t pTotal, Double_t pathLength, Double_t rho) |
| c04e3238 |
679 | { |
| 8cde4af5 |
680 | /// Returns the total momentum energy loss in the front absorber |
| 681 | Double_t muMass = 0.10566; |
| 682 | Double_t p2=pTotal*pTotal; |
| 683 | Double_t beta2=p2/(p2 + muMass*muMass); |
| 684 | Double_t dE=ApproximateBetheBloch(beta2)*pathLength*rho; |
| 685 | |
| 686 | return dE; |
| c04e3238 |
687 | } |
| 688 | |
| 689 | //__________________________________________________________________________ |
| 8cde4af5 |
690 | Double_t AliMUONTrackExtrap::ApproximateBetheBloch(Double_t beta2) { |
| 691 | //------------------------------------------------------------------ |
| 692 | // This is an approximation of the Bethe-Bloch formula with |
| 693 | // the density effect taken into account at beta*gamma > 3.5 |
| 694 | // (the approximation is reasonable only for solid materials) |
| 695 | //------------------------------------------------------------------ |
| 696 | if (beta2/(1-beta2)>3.5*3.5) |
| 697 | return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2); |
| 698 | |
| 699 | return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2); |
| c04e3238 |
700 | } |
| 701 | |
| 702 | //__________________________________________________________________________ |
| 703 | void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t *vect, Double_t *vout) |
| 704 | { |
| 705 | /// ****************************************************************** |
| 706 | /// * * |
| 707 | /// * Performs the tracking of one step in a magnetic field * |
| 708 | /// * The trajectory is assumed to be a helix in a constant field * |
| 709 | /// * taken at the mid point of the step. * |
| 710 | /// * Parameters: * |
| 711 | /// * input * |
| 712 | /// * STEP =arc length of the step asked * |
| 713 | /// * VECT =input vector (position,direction cos and momentum) * |
| 714 | /// * CHARGE= electric charge of the particle * |
| 715 | /// * output * |
| 716 | /// * VOUT = same as VECT after completion of the step * |
| 717 | /// * * |
| 718 | /// * ==>Called by : <USER>, GUSWIM * |
| 719 | /// * Author m.hansroul ********* * |
| 720 | /// * modified s.egli, s.v.levonian * |
| 721 | /// * modified v.perevoztchikov |
| 722 | /// * * |
| 723 | /// ****************************************************************** |
| 724 | |
| 725 | // modif: everything in double precision |
| 726 | |
| 727 | Double_t xyz[3], h[4], hxp[3]; |
| 728 | Double_t h2xy, hp, rho, tet; |
| 729 | Double_t sint, sintt, tsint, cos1t; |
| 730 | Double_t f1, f2, f3, f4, f5, f6; |
| 731 | |
| 732 | const Int_t kix = 0; |
| 733 | const Int_t kiy = 1; |
| 734 | const Int_t kiz = 2; |
| 735 | const Int_t kipx = 3; |
| 736 | const Int_t kipy = 4; |
| 737 | const Int_t kipz = 5; |
| 738 | const Int_t kipp = 6; |
| 739 | |
| 740 | const Double_t kec = 2.9979251e-4; |
| 741 | // |
| 742 | // ------------------------------------------------------------------ |
| 743 | // |
| 744 | // units are kgauss,centimeters,gev/c |
| 745 | // |
| 746 | vout[kipp] = vect[kipp]; |
| 747 | if (TMath::Abs(charge) < 0.00001) { |
| 748 | for (Int_t i = 0; i < 3; i++) { |
| 749 | vout[i] = vect[i] + step * vect[i+3]; |
| 750 | vout[i+3] = vect[i+3]; |
| 751 | } |
| 752 | return; |
| 753 | } |
| 754 | xyz[0] = vect[kix] + 0.5 * step * vect[kipx]; |
| 755 | xyz[1] = vect[kiy] + 0.5 * step * vect[kipy]; |
| 756 | xyz[2] = vect[kiz] + 0.5 * step * vect[kipz]; |
| 757 | |
| 758 | //cmodif: call gufld (xyz, h) changed into: |
| 759 | GetField (xyz, h); |
| 760 | |
| 761 | h2xy = h[0]*h[0] + h[1]*h[1]; |
| 762 | h[3] = h[2]*h[2]+ h2xy; |
| 763 | if (h[3] < 1.e-12) { |
| 764 | for (Int_t i = 0; i < 3; i++) { |
| 765 | vout[i] = vect[i] + step * vect[i+3]; |
| 766 | vout[i+3] = vect[i+3]; |
| 767 | } |
| 768 | return; |
| 769 | } |
| 770 | if (h2xy < 1.e-12*h[3]) { |
| 771 | ExtrapOneStepHelix3(charge*h[2], step, vect, vout); |
| 772 | return; |
| 773 | } |
| 774 | h[3] = TMath::Sqrt(h[3]); |
| 775 | h[0] /= h[3]; |
| 776 | h[1] /= h[3]; |
| 777 | h[2] /= h[3]; |
| 778 | h[3] *= kec; |
| 779 | |
| 780 | hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy]; |
| 781 | hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz]; |
| 782 | hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx]; |
| 783 | |
| 784 | hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz]; |
| 785 | |
| 786 | rho = -charge*h[3]/vect[kipp]; |
| 787 | tet = rho * step; |
| 788 | |
| 789 | if (TMath::Abs(tet) > 0.15) { |
| 790 | sint = TMath::Sin(tet); |
| 791 | sintt = (sint/tet); |
| 792 | tsint = (tet-sint)/tet; |
| 793 | cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet; |
| 794 | } else { |
| 795 | tsint = tet*tet/36.; |
| 796 | sintt = (1. - tsint); |
| 797 | sint = tet*sintt; |
| 798 | cos1t = 0.5*tet; |
| 799 | } |
| 800 | |
| 801 | f1 = step * sintt; |
| 802 | f2 = step * cos1t; |
| 803 | f3 = step * tsint * hp; |
| 804 | f4 = -tet*cos1t; |
| 805 | f5 = sint; |
| 806 | f6 = tet * cos1t * hp; |
| 807 | |
| 808 | vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0]; |
| 809 | vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1]; |
| 810 | vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2]; |
| 811 | |
| 812 | vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0]; |
| 813 | vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1]; |
| 814 | vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2]; |
| 815 | |
| 816 | return; |
| 817 | } |
| 818 | |
| 819 | //__________________________________________________________________________ |
| 820 | void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout) |
| 821 | { |
| 822 | /// ****************************************************************** |
| 823 | /// * * |
| 824 | /// * Tracking routine in a constant field oriented * |
| 825 | /// * along axis 3 * |
| 826 | /// * Tracking is performed with a conventional * |
| 827 | /// * helix step method * |
| 828 | /// * * |
| 829 | /// * ==>Called by : <USER>, GUSWIM * |
| 830 | /// * Authors R.Brun, M.Hansroul ********* * |
| 831 | /// * Rewritten V.Perevoztchikov |
| 832 | /// * * |
| 833 | /// ****************************************************************** |
| 834 | |
| 835 | Double_t hxp[3]; |
| 836 | Double_t h4, hp, rho, tet; |
| 837 | Double_t sint, sintt, tsint, cos1t; |
| 838 | Double_t f1, f2, f3, f4, f5, f6; |
| 839 | |
| 840 | const Int_t kix = 0; |
| 841 | const Int_t kiy = 1; |
| 842 | const Int_t kiz = 2; |
| 843 | const Int_t kipx = 3; |
| 844 | const Int_t kipy = 4; |
| 845 | const Int_t kipz = 5; |
| 846 | const Int_t kipp = 6; |
| 847 | |
| 848 | const Double_t kec = 2.9979251e-4; |
| 849 | |
| 850 | // |
| 851 | // ------------------------------------------------------------------ |
| 852 | // |
| 853 | // units are kgauss,centimeters,gev/c |
| 854 | // |
| 855 | vout[kipp] = vect[kipp]; |
| 856 | h4 = field * kec; |
| 857 | |
| 858 | hxp[0] = - vect[kipy]; |
| 859 | hxp[1] = + vect[kipx]; |
| 860 | |
| 861 | hp = vect[kipz]; |
| 862 | |
| 863 | rho = -h4/vect[kipp]; |
| 864 | tet = rho * step; |
| 865 | if (TMath::Abs(tet) > 0.15) { |
| 866 | sint = TMath::Sin(tet); |
| 867 | sintt = (sint/tet); |
| 868 | tsint = (tet-sint)/tet; |
| 869 | cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet; |
| 870 | } else { |
| 871 | tsint = tet*tet/36.; |
| 872 | sintt = (1. - tsint); |
| 873 | sint = tet*sintt; |
| 874 | cos1t = 0.5*tet; |
| 875 | } |
| 876 | |
| 877 | f1 = step * sintt; |
| 878 | f2 = step * cos1t; |
| 879 | f3 = step * tsint * hp; |
| 880 | f4 = -tet*cos1t; |
| 881 | f5 = sint; |
| 882 | f6 = tet * cos1t * hp; |
| 883 | |
| 884 | vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0]; |
| 885 | vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1]; |
| 886 | vout[kiz] = vect[kiz] + f1*vect[kipz] + f3; |
| 887 | |
| 888 | vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0]; |
| 889 | vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1]; |
| 890 | vout[kipz] = vect[kipz] + f4*vect[kipz] + f6; |
| 891 | |
| 892 | return; |
| 893 | } |
| 8cde4af5 |
894 | |
| c04e3238 |
895 | //__________________________________________________________________________ |
| 896 | void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout) |
| 897 | { |
| 898 | /// ****************************************************************** |
| 899 | /// * * |
| 900 | /// * Runge-Kutta method for tracking a particle through a magnetic * |
| 901 | /// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * |
| 902 | /// * Standards, procedure 25.5.20) * |
| 903 | /// * * |
| 904 | /// * Input parameters * |
| 905 | /// * CHARGE Particle charge * |
| 906 | /// * STEP Step size * |
| 907 | /// * VECT Initial co-ords,direction cosines,momentum * |
| 908 | /// * Output parameters * |
| 909 | /// * VOUT Output co-ords,direction cosines,momentum * |
| 910 | /// * User routine called * |
| 911 | /// * CALL GUFLD(X,F) * |
| 912 | /// * * |
| 913 | /// * ==>Called by : <USER>, GUSWIM * |
| 914 | /// * Authors R.Brun, M.Hansroul ********* * |
| 915 | /// * V.Perevoztchikov (CUT STEP implementation) * |
| 916 | /// * * |
| 917 | /// * * |
| 918 | /// ****************************************************************** |
| 919 | |
| 920 | Double_t h2, h4, f[4]; |
| 921 | Double_t xyzt[3], a, b, c, ph,ph2; |
| 922 | Double_t secxs[4],secys[4],seczs[4],hxp[3]; |
| 923 | Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt; |
| 924 | Double_t est, at, bt, ct, cba; |
| 925 | Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost; |
| 926 | |
| 927 | Double_t x; |
| 928 | Double_t y; |
| 929 | Double_t z; |
| 930 | |
| 931 | Double_t xt; |
| 932 | Double_t yt; |
| 933 | Double_t zt; |
| 934 | |
| 935 | Double_t maxit = 1992; |
| 936 | Double_t maxcut = 11; |
| 937 | |
| 938 | const Double_t kdlt = 1e-4; |
| 939 | const Double_t kdlt32 = kdlt/32.; |
| 940 | const Double_t kthird = 1./3.; |
| 941 | const Double_t khalf = 0.5; |
| 942 | const Double_t kec = 2.9979251e-4; |
| 943 | |
| 944 | const Double_t kpisqua = 9.86960440109; |
| 945 | const Int_t kix = 0; |
| 946 | const Int_t kiy = 1; |
| 947 | const Int_t kiz = 2; |
| 948 | const Int_t kipx = 3; |
| 949 | const Int_t kipy = 4; |
| 950 | const Int_t kipz = 5; |
| 951 | |
| 952 | // *. |
| 953 | // *. ------------------------------------------------------------------ |
| 954 | // *. |
| 955 | // * this constant is for units cm,gev/c and kgauss |
| 956 | // * |
| 957 | Int_t iter = 0; |
| 958 | Int_t ncut = 0; |
| 959 | for(Int_t j = 0; j < 7; j++) |
| 960 | vout[j] = vect[j]; |
| 961 | |
| 962 | Double_t pinv = kec * charge / vect[6]; |
| 963 | Double_t tl = 0.; |
| 964 | Double_t h = step; |
| 965 | Double_t rest; |
| 966 | |
| 967 | |
| 968 | do { |
| 969 | rest = step - tl; |
| 970 | if (TMath::Abs(h) > TMath::Abs(rest)) h = rest; |
| 971 | //cmodif: call gufld(vout,f) changed into: |
| 972 | |
| 973 | GetField(vout,f); |
| 974 | |
| 975 | // * |
| 976 | // * start of integration |
| 977 | // * |
| 978 | x = vout[0]; |
| 979 | y = vout[1]; |
| 980 | z = vout[2]; |
| 981 | a = vout[3]; |
| 982 | b = vout[4]; |
| 983 | c = vout[5]; |
| 984 | |
| 985 | h2 = khalf * h; |
| 986 | h4 = khalf * h2; |
| 987 | ph = pinv * h; |
| 988 | ph2 = khalf * ph; |
| 989 | secxs[0] = (b * f[2] - c * f[1]) * ph2; |
| 990 | secys[0] = (c * f[0] - a * f[2]) * ph2; |
| 991 | seczs[0] = (a * f[1] - b * f[0]) * ph2; |
| 992 | ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]); |
| 993 | if (ang2 > kpisqua) break; |
| 994 | |
| 995 | dxt = h2 * a + h4 * secxs[0]; |
| 996 | dyt = h2 * b + h4 * secys[0]; |
| 997 | dzt = h2 * c + h4 * seczs[0]; |
| 998 | xt = x + dxt; |
| 999 | yt = y + dyt; |
| 1000 | zt = z + dzt; |
| 1001 | // * |
| 1002 | // * second intermediate point |
| 1003 | // * |
| 1004 | |
| 1005 | est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt); |
| 1006 | if (est > h) { |
| 1007 | if (ncut++ > maxcut) break; |
| 1008 | h *= khalf; |
| 1009 | continue; |
| 1010 | } |
| 1011 | |
| 1012 | xyzt[0] = xt; |
| 1013 | xyzt[1] = yt; |
| 1014 | xyzt[2] = zt; |
| 1015 | |
| 1016 | //cmodif: call gufld(xyzt,f) changed into: |
| 1017 | GetField(xyzt,f); |
| 1018 | |
| 1019 | at = a + secxs[0]; |
| 1020 | bt = b + secys[0]; |
| 1021 | ct = c + seczs[0]; |
| 1022 | |
| 1023 | secxs[1] = (bt * f[2] - ct * f[1]) * ph2; |
| 1024 | secys[1] = (ct * f[0] - at * f[2]) * ph2; |
| 1025 | seczs[1] = (at * f[1] - bt * f[0]) * ph2; |
| 1026 | at = a + secxs[1]; |
| 1027 | bt = b + secys[1]; |
| 1028 | ct = c + seczs[1]; |
| 1029 | secxs[2] = (bt * f[2] - ct * f[1]) * ph2; |
| 1030 | secys[2] = (ct * f[0] - at * f[2]) * ph2; |
| 1031 | seczs[2] = (at * f[1] - bt * f[0]) * ph2; |
| 1032 | dxt = h * (a + secxs[2]); |
| 1033 | dyt = h * (b + secys[2]); |
| 1034 | dzt = h * (c + seczs[2]); |
| 1035 | xt = x + dxt; |
| 1036 | yt = y + dyt; |
| 1037 | zt = z + dzt; |
| 1038 | at = a + 2.*secxs[2]; |
| 1039 | bt = b + 2.*secys[2]; |
| 1040 | ct = c + 2.*seczs[2]; |
| 1041 | |
| 1042 | est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt); |
| 1043 | if (est > 2.*TMath::Abs(h)) { |
| 1044 | if (ncut++ > maxcut) break; |
| 1045 | h *= khalf; |
| 1046 | continue; |
| 1047 | } |
| 1048 | |
| 1049 | xyzt[0] = xt; |
| 1050 | xyzt[1] = yt; |
| 1051 | xyzt[2] = zt; |
| 1052 | |
| 1053 | //cmodif: call gufld(xyzt,f) changed into: |
| 1054 | GetField(xyzt,f); |
| 1055 | |
| 1056 | z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h; |
| 1057 | y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h; |
| 1058 | x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h; |
| 1059 | |
| 1060 | secxs[3] = (bt*f[2] - ct*f[1])* ph2; |
| 1061 | secys[3] = (ct*f[0] - at*f[2])* ph2; |
| 1062 | seczs[3] = (at*f[1] - bt*f[0])* ph2; |
| 1063 | a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird; |
| 1064 | b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird; |
| 1065 | c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird; |
| 1066 | |
| 1067 | est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2])) |
| 1068 | + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2])) |
| 1069 | + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2])); |
| 1070 | |
| 1071 | if (est > kdlt && TMath::Abs(h) > 1.e-4) { |
| 1072 | if (ncut++ > maxcut) break; |
| 1073 | h *= khalf; |
| 1074 | continue; |
| 1075 | } |
| 1076 | |
| 1077 | ncut = 0; |
| 1078 | // * if too many iterations, go to helix |
| 1079 | if (iter++ > maxit) break; |
| 1080 | |
| 1081 | tl += h; |
| 1082 | if (est < kdlt32) |
| 1083 | h *= 2.; |
| 1084 | cba = 1./ TMath::Sqrt(a*a + b*b + c*c); |
| 1085 | vout[0] = x; |
| 1086 | vout[1] = y; |
| 1087 | vout[2] = z; |
| 1088 | vout[3] = cba*a; |
| 1089 | vout[4] = cba*b; |
| 1090 | vout[5] = cba*c; |
| 1091 | rest = step - tl; |
| 1092 | if (step < 0.) rest = -rest; |
| 1093 | if (rest < 1.e-5*TMath::Abs(step)) return; |
| 1094 | |
| 1095 | } while(1); |
| 1096 | |
| 1097 | // angle too big, use helix |
| 1098 | |
| 1099 | f1 = f[0]; |
| 1100 | f2 = f[1]; |
| 1101 | f3 = f[2]; |
| 1102 | f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3); |
| 1103 | rho = -f4*pinv; |
| 1104 | tet = rho * step; |
| 1105 | |
| 1106 | hnorm = 1./f4; |
| 1107 | f1 = f1*hnorm; |
| 1108 | f2 = f2*hnorm; |
| 1109 | f3 = f3*hnorm; |
| 1110 | |
| 1111 | hxp[0] = f2*vect[kipz] - f3*vect[kipy]; |
| 1112 | hxp[1] = f3*vect[kipx] - f1*vect[kipz]; |
| 1113 | hxp[2] = f1*vect[kipy] - f2*vect[kipx]; |
| 1114 | |
| 1115 | hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz]; |
| 1116 | |
| 1117 | rho1 = 1./rho; |
| 1118 | sint = TMath::Sin(tet); |
| 1119 | cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet); |
| 1120 | |
| 1121 | g1 = sint*rho1; |
| 1122 | g2 = cost*rho1; |
| 1123 | g3 = (tet-sint) * hp*rho1; |
| 1124 | g4 = -cost; |
| 1125 | g5 = sint; |
| 1126 | g6 = cost * hp; |
| 1127 | |
| 1128 | vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1; |
| 1129 | vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2; |
| 1130 | vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3; |
| 1131 | |
| 1132 | vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1; |
| 1133 | vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2; |
| 1134 | vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3; |
| 1135 | |
| 1136 | return; |
| 1137 | } |
| 8cde4af5 |
1138 | |
| c04e3238 |
1139 | //___________________________________________________________ |
| 1140 | void AliMUONTrackExtrap::GetField(Double_t *Position, Double_t *Field) |
| 1141 | { |
| 1142 | /// interface for arguments in double precision (Why ? ChF) |
| 1143 | Float_t x[3], b[3]; |
| 1144 | |
| 1145 | x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2]; |
| 1146 | |
| 1147 | if (fgkField) fgkField->Field(x,b); |
| 1148 | else { |
| 1149 | cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<<endl; |
| 1150 | exit(-1); |
| 1151 | } |
| 1152 | |
| 1153 | Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2]; |
| 1154 | |
| 1155 | return; |
| 1156 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame