X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrackExtrap.cxx;h=c2ec6aa2deba05c201cafc4d71772f358900ba4a;hb=690d2205d260c6401c46797c529839b4304c506b;hp=42fe6ffb4db9f55e755c34f4e88ae4ce059bafd1;hpb=4284483e6e9366546ff7c5fd437912accdf263f5;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrackExtrap.cxx b/MUON/AliMUONTrackExtrap.cxx index 42fe6ffb4db..c2ec6aa2deb 100644 --- a/MUON/AliMUONTrackExtrap.cxx +++ b/MUON/AliMUONTrackExtrap.cxx @@ -15,31 +15,27 @@ /* $Id$ */ -/////////////////////////////////////////////////// -// -// Tools -// for -// track -// extrapolation -// in -// ALICE -// dimuon -// spectrometer -// -/////////////////////////////////////////////////// - -#include -#include -#include +//----------------------------------------------------------------------------- +// Class AliMUONTrackExtrap +// ------------------------ +// Tools for track extrapolation in ALICE dimuon spectrometer +// Author: Philippe Pillot +//----------------------------------------------------------------------------- #include "AliMUONTrackExtrap.h" #include "AliMUONTrackParam.h" #include "AliMUONConstants.h" + #include "AliMagF.h" -#include "AliLog.h" -#include "AliTracker.h" +#include +#include + +#include + +/// \cond CLASSIMP ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context +/// \endcond const AliMagF* AliMUONTrackExtrap::fgkField = 0x0; const Bool_t AliMUONTrackExtrap::fgkUseHelix = kFALSE; @@ -47,7 +43,7 @@ const Int_t AliMUONTrackExtrap::fgkMaxStepNumber = 5000; const Double_t AliMUONTrackExtrap::fgkHelixStepLength = 6.; const Double_t AliMUONTrackExtrap::fgkRungeKuttaMaxResidue = 0.002; - //__________________________________________________________________________ +//__________________________________________________________________________ Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum) { /// Returns impact parameter at vertex in bending plane (cm), @@ -70,7 +66,7 @@ Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingM return (-0.0003 * simpleBValue * simpleBLength * simpleBPosition / bendingMomentum); } - //__________________________________________________________________________ +//__________________________________________________________________________ Double_t AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(Double_t impactParam) { /// Returns signed bending momentum in bending plane (GeV/c), @@ -93,7 +89,35 @@ Double_t AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(Double_t impactPa return (-0.0003 * simpleBValue * simpleBLength * simpleBPosition / impactParam); } - //__________________________________________________________________________ +//__________________________________________________________________________ +void AliMUONTrackExtrap::LinearExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd) +{ + /// Track parameters (and their covariances if any) linearly extrapolated to the plane at "zEnd". + /// On return, results from the extrapolation are updated in trackParam. + + if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z + + // Compute track parameters + Double_t dZ = zEnd - trackParam->GetZ(); + trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ); + trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ); + trackParam->SetZ(zEnd); + + // Update track parameters covariances if any + if (trackParam->CovariancesExist()) { + TMatrixD paramCov(trackParam->GetCovariances()); + paramCov(0,0) += dZ * dZ * paramCov(1,1) + 2. * dZ * paramCov(0,1); + paramCov(0,1) += dZ * paramCov(1,1); + paramCov(1,0) = paramCov(0,1); + paramCov(2,2) += dZ * dZ * paramCov(3,3) + 2. * dZ * paramCov(2,3); + paramCov(2,3) += dZ * paramCov(3,3); + paramCov(3,2) = paramCov(2,3); + trackParam->SetCovariances(paramCov); + } + +} + +//__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd) { /// Interface to track parameter extrapolation to the plane at "Z" using Helix or Rungekutta algorithm. @@ -102,7 +126,7 @@ void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd) else AliMUONTrackExtrap::ExtrapToZRungekutta(trackParam,zEnd); } - //__________________________________________________________________________ +//__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapToZHelix(AliMUONTrackParam* trackParam, Double_t zEnd) { /// Track parameter extrapolation to the plane at "Z" using Helix algorithm. @@ -125,7 +149,7 @@ void AliMUONTrackExtrap::ExtrapToZHelix(AliMUONTrackParam* trackParam, Double_t stepNumber++; ExtrapOneStepHelix(chargeExtrap, fgkHelixStepLength, v3, v3New); if ((-forwardBackward * (v3New[2] - zEnd)) > 0.0) break; // one is beyond Z spectro. z<0 - // better use TArray ???? + // better use TArray ???? for (i3 = 0; i3 < 7; i3++) {v3[i3] = v3New[i3];} } // check fgkMaxStepNumber ???? @@ -155,7 +179,7 @@ void AliMUONTrackExtrap::ExtrapToZHelix(AliMUONTrackParam* trackParam, Double_t RecoverTrackParam(v3, chargeExtrap * forwardBackward, trackParam); } - //__________________________________________________________________________ +//__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Double_t zEnd) { /// Track parameter extrapolation to the plane at "Z" using Rungekutta algorithm. @@ -177,7 +201,7 @@ void AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Doub dZ = zEnd - trackParam->GetZ(); // step lenght assuming linear trajectory step = dZ * TMath::Sqrt(1.0 + trackParam->GetBendingSlope()*trackParam->GetBendingSlope() + - trackParam->GetNonBendingSlope()*trackParam->GetNonBendingSlope()); + trackParam->GetNonBendingSlope()*trackParam->GetNonBendingSlope()); ConvertTrackParamForExtrap(trackParam, forwardBackward, v3); do { // reduce step lenght while zEnd oversteped if (stepNumber > fgkMaxStepNumber) { @@ -199,7 +223,7 @@ void AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Doub trackParam->SetZ(zEnd); } - //__________________________________________________________________________ +//__________________________________________________________________________ void AliMUONTrackExtrap::ConvertTrackParamForExtrap(AliMUONTrackParam* trackParam, Double_t forwardBackward, Double_t *v3) { /// Set vector of Geant3 parameters pointed to by "v3" from track parameters in trackParam. @@ -216,7 +240,7 @@ void AliMUONTrackExtrap::ConvertTrackParamForExtrap(AliMUONTrackParam* trackPara v3[4] = trackParam->GetBendingSlope() * v3[5]; // PY/PTOT } - //__________________________________________________________________________ +//__________________________________________________________________________ void AliMUONTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMUONTrackParam* trackParam) { /// Set track parameters in trackParam from Geant3 parameters pointed to by "v3", @@ -231,169 +255,342 @@ void AliMUONTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMUO trackParam->SetNonBendingSlope(v3[3]/v3[5]); } - //__________________________________________________________________________ -void AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd) +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator) { /// Track parameters and their covariances extrapolated to the plane at "zEnd". /// On return, results from the extrapolation are updated in trackParam. if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z + // No need to propagate the covariance matrix if it does not exist + if (!trackParam->CovariancesExist()) { + cout<<"W-AliMUONTrackExtrap::ExtrapToZCov: Covariance matrix does not exist"<GetCovariances(); // Extrapolate track parameters to "zEnd" ExtrapToZ(trackParam,zEnd); - Double_t extrapNonBendingCoor = trackParam->GetNonBendingCoor(); - Double_t extrapNonBendingSlope = trackParam->GetNonBendingSlope(); - Double_t extrapBendingCoor = trackParam->GetBendingCoor(); - Double_t extrapBendingSlope = trackParam->GetBendingSlope(); - Double_t extrapInverseBendingMomentum = trackParam->GetInverseBendingMomentum(); - // Get the pointer to the parameter covariance matrix - if (!trackParam->CovariancesExist()) { - //cout<<"W-AliMUONTrackExtrap::ExtrapToZCov: track parameter covariance matrix does not exist"< nothing to extrapolate !!"<GetCovariances(); + // Get reference to the extrapolated parameters + const TMatrixD& extrapParam = trackParam->GetParameters(); // Calculate the jacobian related to the track parameters extrapolation to "zEnd" TMatrixD jacob(5,5); - jacob = 0.; - Double_t dParam[5]; + jacob.Zero(); + TMatrixD dParam(5,1); for (Int_t i=0; i<5; i++) { // Skip jacobian calculation for parameters with no associated error - if ((*paramCov)(i,i) == 0.) continue; + if (kParamCov(i,i) == 0.) continue; + // Small variation of parameter i only for (Int_t j=0; j<5; j++) { if (j==i) { - dParam[j] = TMath::Sqrt((*paramCov)(i,i)); - if (j == 4) dParam[j] *= TMath::Sign(1.,-inverseBendingMomentum); // variation always in the same direction - } else dParam[j] = 0.; + dParam(j,0) = TMath::Sqrt(kParamCov(i,i)); + if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramSave(4,0)); // variation always in the same direction + } else dParam(j,0) = 0.; } + // Set new parameters - trackParamSave.SetNonBendingCoor (nonBendingCoor + dParam[0]); - trackParamSave.SetNonBendingSlope (nonBendingSlope + dParam[1]); - trackParamSave.SetBendingCoor (bendingCoor + dParam[2]); - trackParamSave.SetBendingSlope (bendingSlope + dParam[3]); - trackParamSave.SetInverseBendingMomentum(inverseBendingMomentum + dParam[4]); - trackParamSave.SetZ (zBegin); + trackParamSave.SetParameters(paramSave); + trackParamSave.AddParameters(dParam); + trackParamSave.SetZ(zBegin); + // Extrapolate new track parameters to "zEnd" ExtrapToZ(&trackParamSave,zEnd); + // Calculate the jacobian - jacob(0,i) = (trackParamSave.GetNonBendingCoor() - extrapNonBendingCoor ) / dParam[i]; - jacob(1,i) = (trackParamSave.GetNonBendingSlope() - extrapNonBendingSlope ) / dParam[i]; - jacob(2,i) = (trackParamSave.GetBendingCoor() - extrapBendingCoor ) / dParam[i]; - jacob(3,i) = (trackParamSave.GetBendingSlope() - extrapBendingSlope ) / dParam[i]; - jacob(4,i) = (trackParamSave.GetInverseBendingMomentum() - extrapInverseBendingMomentum) / dParam[i]; + TMatrixD jacobji(trackParamSave.GetParameters(),TMatrixD::kMinus,extrapParam); + jacobji *= 1. / dParam(i,0); + jacob.SetSub(0,i,jacobji); } // Extrapolate track parameter covariances to "zEnd" - TMatrixD tmp((*paramCov),TMatrixD::kMultTranspose,jacob); - (*paramCov) = TMatrixD(jacob,TMatrixD::kMult,tmp); + TMatrixD tmp(kParamCov,TMatrixD::kMultTranspose,jacob); + TMatrixD tmp2(jacob,TMatrixD::kMult,tmp); + trackParam->SetCovariances(tmp2); + // Update the propagator if required + if (updatePropagator) trackParam->UpdatePropagator(jacob); } - //__________________________________________________________________________ -void AliMUONTrackExtrap::ExtrapToStation(AliMUONTrackParam* trackParamIn, Int_t station, AliMUONTrackParam *trackParamOut) +//__________________________________________________________________________ +void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2) { - /// Track parameters extrapolated from "trackParamIn" to both chambers of the station(0..) "station" - /// are returned in the array (dimension 2) of track parameters pointed to by "TrackParamOut" - /// (index 0 and 1 for first and second chambers). - Double_t extZ[2], z1, z2; - Int_t i1 = -1, i2 = -1; // = -1 to avoid compilation warnings - // range of station to be checked ???? - z1 = AliMUONConstants::DefaultChamberZ(2 * station); - z2 = AliMUONConstants::DefaultChamberZ(2 * station + 1); - // First and second Z to extrapolate at - if ((z1 > trackParamIn->GetZ()) && (z2 > trackParamIn->GetZ())) {i1 = 0; i2 = 1;} - else if ((z1 < trackParamIn->GetZ()) && (z2 < trackParamIn->GetZ())) {i1 = 1; i2 = 0;} - else { - cout<<"E-AliMUONTrackExtrap::ExtrapToStation: Starting Z ("<GetZ() - <<") in between z1 ("<GetBendingSlope(); + Double_t nonBendingSlope = param->GetNonBendingSlope(); + Double_t inverseBendingMomentum = param->GetInverseBendingMomentum(); + Double_t alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) / + (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1 + Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2); + Double_t covCorrSlope = alpha2 * (pathLength * f0 - f1); + Double_t varSlop = alpha2 * f0; + + // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeX + Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope); + Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope / + (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope); + + // Set MCS covariance matrix + TMatrixD newParamCov(param->GetCovariances()); + // Non bending plane + newParamCov(0,0) += varCoor; newParamCov(0,1) += covCorrSlope; + newParamCov(1,0) += covCorrSlope; newParamCov(1,1) += varSlop; + // Bending plane + newParamCov(2,2) += varCoor; newParamCov(2,3) += covCorrSlope; + newParamCov(3,2) += covCorrSlope; newParamCov(3,3) += varSlop; + // Inverse bending momentum (due to dependences with bending and non bending slopes) + newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope; + newParamCov(4,1) += dqPxydSlopeX * varSlop; newParamCov(1,4) += dqPxydSlopeX * varSlop; + newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope; + newParamCov(4,3) += dqPxydSlopeY * varSlop; newParamCov(3,4) += dqPxydSlopeY * varSlop; + newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop; + + // Set new covariances + param->SetCovariances(newParamCov); } - //__________________________________________________________________________ -void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx) +//__________________________________________________________________________ +void AliMUONTrackExtrap::CorrectMCSEffectInAbsorber(AliMUONTrackParam* param, + Double_t xVtx, Double_t yVtx, Double_t zVtx, + Double_t errXVtx, Double_t errYVtx, + Double_t absZBeg, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2) { - /// Extrapolation to the vertex (at the z position "zVtx") without Branson and Field correction. - /// Returns the track parameters resulting from the extrapolation in the current TrackParam. - /// Include multiple Coulomb scattering effects in trackParam covariances. + /// Correct parameters and corresponding covariances using Branson correction + /// - input param are parameters and covariances at the end of absorber + /// - output param are parameters and covariances at vertex + /// Absorber correction parameters are supposed to be calculated at the current track z-position - if (trackParam->GetZ() > zVtx) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<GetZ() - <<") upstream the vertex (zVtx = "<0.) ? absZBeg - f2/f1 : 0.; - if (zVtx < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Ending Z ("<GetNonBendingCoor() - xVtx) / (zB - zVtx); + newParam(2,0) = yVtx; + newParam(3,0) = (param->GetBendingCoor() - yVtx) / (zB - zVtx); + newParam(4,0) = param->GetCharge() / param->P() * + TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) / + TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0)); + + // Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system + TMatrixD paramCovP(param->GetCovariances()); + Cov2CovP(param->GetParameters(),paramCovP); + + // Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system + TMatrixD paramCovVtx(5,5); + paramCovVtx.Zero(); + paramCovVtx(0,0) = errXVtx * errXVtx; + paramCovVtx(1,1) = paramCovP(0,0); + paramCovVtx(2,2) = errYVtx * errYVtx; + paramCovVtx(3,3) = paramCovP(2,2); + paramCovVtx(4,4) = paramCovP(4,4); + paramCovVtx(1,3) = paramCovP(0,2); + paramCovVtx(3,1) = paramCovP(2,0); + paramCovVtx(1,4) = paramCovP(0,4); + paramCovVtx(4,1) = paramCovP(4,0); + paramCovVtx(3,4) = paramCovP(2,4); + paramCovVtx(4,3) = paramCovP(4,2); + + // Jacobian of the transformation (XVtx, X, YVtx, Y, q*PTot) -> (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx) + TMatrixD jacob(5,5); + jacob.UnitMatrix(); + jacob(1,0) = - 1. / (zB - zVtx); + jacob(1,1) = 1. / (zB - zVtx); + jacob(3,2) = - 1. / (zB - zVtx); + jacob(3,3) = 1. / (zB - zVtx); + + // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx) coordinate system + TMatrixD tmp(paramCovVtx,TMatrixD::kMultTranspose,jacob); + TMatrixD newParamCov(jacob,TMatrixD::kMult,tmp); + + // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q/PyzVtx) coordinate system + CovP2Cov(newParam,newParamCov); + + // Set parameters and covariances at vertex + param->SetParameters(newParam); + param->SetZ(zVtx); + param->SetCovariances(newParamCov); +} + +//__________________________________________________________________________ +void AliMUONTrackExtrap::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss, Double_t sigmaELoss2) +{ + /// Correct parameters for energy loss and add energy loss fluctuation effect to covariances + + // Get parameter covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system + TMatrixD newParamCov(param->GetCovariances()); + Cov2CovP(param->GetParameters(),newParamCov); + + // Add effects of energy loss fluctuation to covariances + newParamCov(4,4) += sigmaELoss2; + + // Compute new parameters corrected for energy loss + Double_t nonBendingSlope = param->GetNonBendingSlope(); + Double_t bendingSlope = param->GetBendingSlope(); + param->SetInverseBendingMomentum(param->GetCharge() / (param->P() + eLoss) * + TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) / + TMath::Sqrt(1.0 + bendingSlope*bendingSlope)); + + // Get new parameter covariances in (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system + CovP2Cov(param->GetParameters(),newParamCov); + + // Set new parameter covariances + param->SetCovariances(newParamCov); +} + +//__________________________________________________________________________ +void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t pTotal, + Double_t &pathLength, Double_t &f0, Double_t &f1, Double_t &f2, + Double_t &meanRho, Double_t &totalELoss, Double_t &sigmaELoss2) +{ + /// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber + /// Calculated assuming a linear propagation from trackXYZIn to trackXYZOut (order is important) + // pathLength: path length between trackXYZIn and trackXYZOut (cm) + // f0: 0th moment of z calculated with the inverse radiation-length distribution + // f1: 1st moment of z calculated with the inverse radiation-length distribution + // f2: 2nd moment of z calculated with the inverse radiation-length distribution + // meanRho: average density of crossed material (g/cm3) + // totalELoss: total energy loss in absorber + + // Reset absorber's parameters + pathLength = 0.; + f0 = 0.; + f1 = 0.; + f2 = 0.; + meanRho = 0.; + totalELoss = 0.; + sigmaELoss2 = 0.; + + // Check whether the geometry is available + if (!gGeoManager) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<GetZ() < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0) - ExtrapToZCov(trackParam,AliMUONConstants::ZAbsorberEnd()); - } else { - cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<GetZ() - <<") upstream or inside the front absorber (zAbsorberEnd = "<InitTrack(trackXYZIn, b); + if (!currentnode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<=0; iElement--) { - zElement = AliMUONConstants::ZAbsorberElement(iElement); - z = trackParam->GetZ(); - if (z > zElement) continue; // spectro. (z<0) - nonBendingCoor = trackParam->GetNonBendingCoor(); - bendingCoor = trackParam->GetBendingCoor(); - r0Norm = nonBendingCoor * nonBendingCoor + bendingCoor * bendingCoor; - r0Norm = TMath::Sqrt(r0Norm) / TMath::Abs(trackParam->GetZ()) / tan3; - if (r0Norm > 1.) x0 = AliMUONConstants::X0AbsorberOut(iElement); // outer part of the absorber - else x0 = AliMUONConstants::X0AbsorberIn(iElement); // inner part of the absorber + // loop over absorber slices and calculate absorber's parameters + Double_t rho = 0.; // material density (g/cm3) + Double_t x0 = 0.; // radiation-length (cm-1) + Double_t atomicA = 0.; // A of material + Double_t atomicZ = 0.; // Z of material + Double_t localPathLength = 0; + Double_t remainingPathLength = pathLength; + Double_t zB = trackXYZIn[2]; + Double_t zE, dzB, dzE; + do { + // Get material properties + TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial(); + rho = material->GetDensity(); + x0 = material->GetRadLen(); + if (!material->IsMixture()) x0 /= rho; // different normalization in the modeler for mixture + atomicA = material->GetA(); + atomicZ = material->GetZ(); - if (zVtx > zElement) { - ExtrapToZCov(trackParam,zElement); // extrapolate to absorber element "iElement" - dZ = zElement - z; - AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances - } else { - ExtrapToZCov(trackParam,zVtx); // extrapolate to zVtx - dZ = zVtx - z; - AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances - break; + // Get path length within this material + gGeoManager->FindNextBoundary(remainingPathLength); + localPathLength = gGeoManager->GetStep() + 1.e-6; + // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step + if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength; + else { + currentnode = gGeoManager->Step(); + if (!currentnode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<IsEntering()) { + // make another small step to try to enter in new absorber slice + gGeoManager->SetStep(0.001); + currentnode = gGeoManager->Step(); + if (!gGeoManager->IsEntering() || !currentnode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"< TGeoShape::Tolerance()); - // finally go to the vertex - ExtrapToZCov(trackParam,zVtx); + meanRho /= pathLength; +} + +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::GetMCSAngle2(const AliMUONTrackParam& param, Double_t dZ, Double_t x0) +{ + /// Return the angular dispersion square due to multiple Coulomb scattering + /// through a material of thickness "dZ" and of radiation length "x0" + /// assuming linear propagation and using the small angle approximation. + + Double_t bendingSlope = param.GetBendingSlope(); + Double_t nonBendingSlope = param.GetNonBendingSlope(); + Double_t inverseTotalMomentum2 = param.GetInverseBendingMomentum() * param.GetInverseBendingMomentum() * + (1.0 + bendingSlope * bendingSlope) / + (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); + // Path length in the material + Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope); + // relativistic velocity + Double_t velo = 1.; + // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory + Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0)); + return theta02 * theta02 * inverseTotalMomentum2 * pathLength / x0; } - //__________________________________________________________________________ -void AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, Double_t dZ, Double_t x0) +//__________________________________________________________________________ +void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0) { /// Add to the track parameter covariances the effects of multiple Coulomb scattering /// through a material of thickness "dZ" and of radiation length "x0" @@ -401,9 +598,10 @@ void AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, D Double_t bendingSlope = param->GetBendingSlope(); Double_t nonBendingSlope = param->GetNonBendingSlope(); - Double_t inverseTotalMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum() * - (1.0 + bendingSlope * bendingSlope) / - (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); + Double_t inverseBendingMomentum = param->GetInverseBendingMomentum(); + Double_t inverseTotalMomentum2 = inverseBendingMomentum * inverseBendingMomentum * + (1.0 + bendingSlope * bendingSlope) / + (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // Path length in the material Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope); Double_t pathLength2 = pathLength * pathLength; @@ -413,300 +611,311 @@ void AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, D Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0)); theta02 *= theta02 * inverseTotalMomentum2 * pathLength / x0; - // Add effects of multiple Coulomb scattering in track parameter covariances - TMatrixD* paramCov = param->GetCovariances(); Double_t varCoor = pathLength2 * theta02 / 3.; Double_t varSlop = theta02; Double_t covCorrSlope = pathLength * theta02 / 2.; + + // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeX + Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope); + Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope / + (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope); + + // Set MCS covariance matrix + TMatrixD newParamCov(param->GetCovariances()); // Non bending plane - (*paramCov)(0,0) += varCoor; (*paramCov)(0,1) += covCorrSlope; - (*paramCov)(1,0) += covCorrSlope; (*paramCov)(1,1) += varSlop; + newParamCov(0,0) += varCoor; newParamCov(0,1) += covCorrSlope; + newParamCov(1,0) += covCorrSlope; newParamCov(1,1) += varSlop; // Bending plane - (*paramCov)(2,2) += varCoor; (*paramCov)(2,3) += covCorrSlope; - (*paramCov)(3,2) += covCorrSlope; (*paramCov)(3,3) += varSlop; + newParamCov(2,2) += varCoor; newParamCov(2,3) += covCorrSlope; + newParamCov(3,2) += covCorrSlope; newParamCov(3,3) += varSlop; + // Inverse bending momentum (due to dependences with bending and non bending slopes) + newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope; + newParamCov(4,1) += dqPxydSlopeX * varSlop; newParamCov(1,4) += dqPxydSlopeX * varSlop; + newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope; + newParamCov(4,3) += dqPxydSlopeY * varSlop; newParamCov(3,4) += dqPxydSlopeY * varSlop; + newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop; + // Set new covariances + param->SetCovariances(newParamCov); } - //__________________________________________________________________________ -void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx) -{ - /// Extrapolation to the vertex. - /// Returns the track parameters resulting from the extrapolation in the current TrackParam. - /// Changes parameters according to Branson correction through the absorber - - // Extrapolates track parameters upstream to the "Z" end of the front absorber - ExtrapToZ(trackParam,AliMUONConstants::ZAbsorberEnd()); // !!! - // Makes Branson correction (multiple scattering + energy loss) - BransonCorrection(trackParam,xVtx,yVtx,zVtx); - // Makes a simple magnetic field correction through the absorber - FieldCorrection(trackParam,AliMUONConstants::ZAbsorberEnd()); -} - - -// Keep this version for future developments - //__________________________________________________________________________ -// void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam) -// { -// /// Branson correction of track parameters -// // the entry parameters have to be calculated at the end of the absorber -// Double_t zEndAbsorber, zBP, xBP, yBP; -// Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta; -// Int_t sign; -// // Would it be possible to calculate all that from Geant configuration ???? -// // and to get the Branson parameters from a function in ABSO module ???? -// // with an eventual contribution from other detectors like START ???? -// // Radiation lengths outer part theta > 3 degres -// static Double_t x01[9] = { 18.8, // C (cm) -// 10.397, // Concrete (cm) -// 0.56, // Plomb (cm) -// 47.26, // Polyethylene (cm) -// 0.56, // Plomb (cm) -// 47.26, // Polyethylene (cm) -// 0.56, // Plomb (cm) -// 47.26, // Polyethylene (cm) -// 0.56 }; // Plomb (cm) -// // inner part theta < 3 degres -// static Double_t x02[3] = { 18.8, // C (cm) -// 10.397, // Concrete (cm) -// 0.35 }; // W (cm) -// // z positions of the materials inside the absober outer part theta > 3 degres -// static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 }; -// // inner part theta < 3 degres -// static Double_t z2[4] = { 90, 315, 467, 503 }; -// static Bool_t first = kTRUE; -// static Double_t zBP1, zBP2, rLimit; -// // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call) -// if (first) { -// first = kFALSE; -// Double_t aNBP = 0.0; -// Double_t aDBP = 0.0; -// Int_t iBound; -// -// for (iBound = 0; iBound < 9; iBound++) { -// aNBP = aNBP + -// (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] - -// z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound]; -// aDBP = aDBP + -// (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound]; -// } -// zBP1 = (2.0 * aNBP) / (3.0 * aDBP); -// aNBP = 0.0; -// aDBP = 0.0; -// for (iBound = 0; iBound < 3; iBound++) { -// aNBP = aNBP + -// (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] - -// z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound]; -// aDBP = aDBP + -// (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound]; -// } -// zBP2 = (2.0 * aNBP) / (3.0 * aDBP); -// rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.); -// } -// -// pYZ = TMath::Abs(1.0 / trackParam->GetInverseBendingMomentum()); -// sign = 1; -// if (trackParam->GetInverseBendingMomentum() < 0) sign = -1; -// pZ = pYZ / (TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope())); -// pX = pZ * trackParam->GetNonBendingSlope(); -// pY = pZ * trackParam->GetBendingSlope(); -// pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX); -// xEndAbsorber = trackParam->GetNonBendingCoor(); -// yEndAbsorber = trackParam->GetBendingCoor(); -// radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber; -// -// if (radiusEndAbsorber2 > rLimit*rLimit) { -// zEndAbsorber = z1[9]; -// zBP = zBP1; -// } else { -// zEndAbsorber = z2[3]; -// zBP = zBP2; -// } -// -// xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP); -// yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP); -// -// // new parameters after Branson and energy loss corrections -// pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP); -// pX = pZ * xBP / zBP; -// pY = pZ * yBP / zBP; -// trackParam->SetBendingSlope(pY/pZ); -// trackParam->SetNonBendingSlope(pX/pZ); -// -// pT = TMath::Sqrt(pX * pX + pY * pY); -// theta = TMath::ATan2(pT, pZ); -// pTotal = TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber); -// -// trackParam->SetInverseBendingMomentum((sign / pTotal) * -// TMath::Sqrt(1.0 + -// trackParam->GetBendingSlope() * trackParam->GetBendingSlope() + -// trackParam->GetNonBendingSlope() * trackParam->GetNonBendingSlope()) / -// TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope())); -// -// // vertex position at (0,0,0) -// // should be taken from vertex measurement ??? -// trackParam->SetBendingCoor(0.); -// trackParam->SetNonBendingCoor(0.); -// trackParam->SetZ(0.); -// } - -void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx) +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, + Double_t xVtx, Double_t yVtx, Double_t zVtx, + Double_t errXVtx, Double_t errYVtx, + Bool_t correctForMCS, Bool_t correctForEnergyLoss) { - /// Branson correction of track parameters - // the entry parameters have to be calculated at the end of the absorber - // simplified version: the z positions of Branson's planes are no longer calculated - // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C - // to test this correction. - // Would it be possible to calculate all that from Geant configuration ???? - // and to get the Branson parameters from a function in ABSO module ???? - // with an eventual contribution from other detectors like START ???? - // change to take into account the vertex postition (real, reconstruct,....) - - Double_t zBP, xBP, yBP; - Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta; - Int_t sign; - static Bool_t first = kTRUE; - static Double_t zBP1, zBP2, rLimit, thetaLimit; - // zBP1 for outer part and zBP2 for inner part (only at the first call) - if (first) { - first = kFALSE; - - thetaLimit = 3.0 * (TMath::Pi()) / 180.; - rLimit = TMath::Abs(AliMUONConstants::ZAbsorberEnd()) * TMath::Tan(thetaLimit); - zBP1 = -450; // values close to those calculated with EvalAbso.C - zBP2 = -480; + /// Main method for extrapolation to the vertex: + /// Returns the track parameters and covariances resulting from the extrapolation of the current trackParam + /// Changes parameters and covariances according to multiple scattering and energy loss corrections: + /// if correctForMCS=kTRUE: compute parameters using Branson correction and add correction resolution to covariances + /// if correctForMCS=kFALSE: add parameter dispersion due to MCS in parameter covariances + /// if correctForEnergyLoss=kTRUE: correct parameters for energy loss and add energy loss fluctuation to covariances + /// if correctForEnergyLoss=kFALSE: do nothing about energy loss + + if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex + + if (trackParam->GetZ() > zVtx) { // spectro. (z<0) + cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") upstream the vertex (zVtx = "<GetInverseBendingMomentum()); - sign = 1; - if (trackParam->GetInverseBendingMomentum() < 0) sign = -1; - pZ = trackParam->Pz(); - pX = trackParam->Px(); - pY = trackParam->Py(); - pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX); - xEndAbsorber = trackParam->GetNonBendingCoor(); - yEndAbsorber = trackParam->GetBendingCoor(); - radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber; - - if (radiusEndAbsorber2 > rLimit*rLimit) { - zBP = zBP1; + + // Check the vertex position relatively to the absorber + if (zVtx < AliMUONConstants::AbsZBeg() && zVtx > AliMUONConstants::AbsZEnd()) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<CovariancesExist()) ExtrapToZCov(trackParam,zVtx); + else ExtrapToZ(trackParam,zVtx); + return; + } + + // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed + if (trackParam->GetZ() > AliMUONConstants::AbsZBeg()) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") upstream the front absorber (zAbsorberBegin = "<CovariancesExist()) ExtrapToZCov(trackParam,zVtx); + else ExtrapToZ(trackParam,zVtx); + return; + } else if (trackParam->GetZ() > AliMUONConstants::AbsZEnd()) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") inside the front absorber ("<CovariancesExist()) ExtrapToZCov(trackParam,AliMUONConstants::AbsZEnd()); + else ExtrapToZ(trackParam,AliMUONConstants::AbsZEnd()); } - - xBP = xEndAbsorber - (pX / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP); - yBP = yEndAbsorber - (pY / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP); - - // new parameters after Branson and energy loss corrections -// Float_t zSmear = zBP - gRandom->Gaus(0.,2.); // !!! possible smearing of Z vertex position - - Float_t zSmear = zBP ; - pZ = pTotal * (zSmear-zVtx) / TMath::Sqrt((xBP-xVtx) * (xBP-xVtx) + (yBP-yVtx) * (yBP-yVtx) +( zSmear-zVtx) * (zSmear-zVtx) ); - pX = pZ * (xBP - xVtx)/ (zSmear-zVtx); - pY = pZ * (yBP - yVtx) / (zSmear-zVtx); - trackParam->SetBendingSlope(pY/pZ); - trackParam->SetNonBendingSlope(pX/pZ); - + // Get absorber correction parameters assuming linear propagation in absorber + Double_t trackXYZOut[3]; + trackXYZOut[0] = trackParam->GetNonBendingCoor(); + trackXYZOut[1] = trackParam->GetBendingCoor(); + trackXYZOut[2] = trackParam->GetZ(); + Double_t trackXYZIn[3]; + if (correctForMCS) { // assume linear propagation until the vertex + trackXYZIn[2] = TMath::Min(zVtx, AliMUONConstants::AbsZBeg()); // spectro. (z<0) + trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]); + trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]); + } else { + AliMUONTrackParam trackParamIn(*trackParam); + ExtrapToZ(&trackParamIn, TMath::Min(zVtx, AliMUONConstants::AbsZBeg())); + trackXYZIn[0] = trackParamIn.GetNonBendingCoor(); + trackXYZIn[1] = trackParamIn.GetBendingCoor(); + trackXYZIn[2] = trackParamIn.GetZ(); + } + Double_t pTot = trackParam->P(); + Double_t pathLength = 0.; + Double_t f0 = 0.; + Double_t f1 = 0.; + Double_t f2 = 0.; + Double_t meanRho = 0.; + Double_t deltaP = 0.; + Double_t sigmaDeltaP2 = 0.; + GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,deltaP,sigmaDeltaP2); - pT = TMath::Sqrt(pX * pX + pY * pY); - theta = TMath::ATan2(pT, TMath::Abs(pZ)); - pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta); + // Compute track parameters and covariances at vertex according to correctForMCS and correctForEnergyLoss flags + if (correctForMCS) { + + if (correctForEnergyLoss) { + + // Correct for multiple scattering and energy loss + CorrectELossEffectInAbsorber(trackParam, 0.5*deltaP, 0.5*sigmaDeltaP2); + CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, + trackXYZIn[2], pathLength, f0, f1, f2); + CorrectELossEffectInAbsorber(trackParam, 0.5*deltaP, 0.5*sigmaDeltaP2); + + } else { + + // Correct for multiple scattering + CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, + trackXYZIn[2], pathLength, f0, f1, f2); + } + + } else { + + if (correctForEnergyLoss) { + + // Correct for energy loss + CorrectELossEffectInAbsorber(trackParam, 0.5*deltaP, 0.5*sigmaDeltaP2); + AddMCSEffectInAbsorber(trackParam, pathLength, f0, f1, f2); + ExtrapToZCov(trackParam, trackXYZIn[2]); + CorrectELossEffectInAbsorber(trackParam, 0.5*deltaP, 0.5*sigmaDeltaP2); + ExtrapToZCov(trackParam, zVtx); + + } else { + + // Correct for multiple scattering and energy loss + AddMCSEffectInAbsorber(trackParam, pathLength, f0, f1, f2); + ExtrapToZCov(trackParam, zVtx); + + } + + } + +} - trackParam->SetInverseBendingMomentum((sign / pTotal) * - TMath::Sqrt(1.0 + - trackParam->GetBendingSlope() * trackParam->GetBendingSlope() + - trackParam->GetNonBendingSlope() * trackParam->GetNonBendingSlope()) / - TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope())); +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, + Double_t xVtx, Double_t yVtx, Double_t zVtx, + Double_t errXVtx, Double_t errYVtx) +{ + /// Extrapolate track parameters to vertex, corrected for multiple scattering and energy loss effects + /// Add branson correction resolution and energy loss fluctuation to parameter covariances + ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kTRUE); +} - // vertex position at (0,0,0) - // should be taken from vertex measurement ??? +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToVertexWithoutELoss(AliMUONTrackParam* trackParam, + Double_t xVtx, Double_t yVtx, Double_t zVtx, + Double_t errXVtx, Double_t errYVtx) +{ + /// Extrapolate track parameters to vertex, corrected for multiple scattering effects only + /// Add branson correction resolution to parameter covariances + ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kFALSE); +} - trackParam->SetBendingCoor(xVtx); - trackParam->SetNonBendingCoor(yVtx); - trackParam->SetZ(zVtx); +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(AliMUONTrackParam* trackParam, Double_t zVtx) +{ + /// Extrapolate track parameters to vertex, corrected for energy loss effects only + /// Add dispersion due to multiple scattering and energy loss fluctuation to parameter covariances + ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kTRUE); +} +//__________________________________________________________________________ +void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx) +{ + /// Extrapolate track parameters to vertex without multiple scattering and energy loss corrections + /// Add dispersion due to multiple scattering to parameter covariances + ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kFALSE); } - //__________________________________________________________________________ -Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta) +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx) { - /// Returns the total momentum corrected from energy loss in the front absorber - // One can use the macros MUONTestAbso.C and DrawTestAbso.C - // to test this correction. - // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002) - Double_t deltaP, pTotalCorrected; - - // Parametrization to be redone according to change of absorber material ???? - // See remark in function BransonCorrection !!!! - // The name is not so good, and there are many arguments !!!! - if (theta < thetaLimit ) { - if (pTotal < 20) { - deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal; - } else { - deltaP = 3.0714 + 0.011767 *pTotal; - } - deltaP *= 0.75; // AZ - } else { - if (pTotal < 20) { - deltaP = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal; - } else { - deltaP = 2.6069 + 0.0051705 * pTotal; - } - deltaP *= 0.9; // AZ + /// Calculate the total momentum energy loss in-between the track position and the vertex assuming a linear propagation + + if (trackParam->GetZ() == zVtx) return 0.; // nothing to be done if already at vertex + + // Check whether the geometry is available + if (!gGeoManager) { + cout<<"E-AliMUONTrackExtrap::TotalMomentumEnergyLoss: no TGeo"<GetNonBendingCoor(); + trackXYZOut[1] = trackParam->GetBendingCoor(); + trackXYZOut[2] = trackParam->GetZ(); + Double_t trackXYZIn[3]; + trackXYZIn[0] = xVtx; + trackXYZIn[1] = yVtx; + trackXYZIn[2] = zVtx; + Double_t pTot = trackParam->P(); + Double_t pathLength = 0.; + Double_t f0 = 0.; + Double_t f1 = 0.; + Double_t f2 = 0.; + Double_t meanRho = 0.; + Double_t totalELoss = 0.; + Double_t sigmaELoss2 = 0.; + GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2); + + return totalELoss; } - //__________________________________________________________________________ -void AliMUONTrackExtrap::FieldCorrection(AliMUONTrackParam *trackParam, Double_t zEnd) +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ) { - /// Correction of the effect of the magnetic field in the absorber - // Assume a constant field along Z axis. - Float_t b[3],x[3]; - Double_t bZ; - Double_t pYZ,pX,pY,pZ,pT; - Double_t pXNew,pYNew; - Double_t c; - - pYZ = TMath::Abs(1.0 / trackParam->GetInverseBendingMomentum()); - c = TMath::Sign(1.0,trackParam->GetInverseBendingMomentum()); // particle charge - - pZ = trackParam->Pz(); - pX = trackParam->Px(); - pY = trackParam->Py(); - pT = TMath::Sqrt(pX*pX+pY*pY); + /// Returns the mean total momentum energy loss of muon with total momentum='pTotal' + /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ' + Double_t muMass = 0.105658369; // GeV + Double_t eMass = 0.510998918e-3; // GeV + Double_t k = 0.307075e-3; // GeV.g^-1.cm^2 + Double_t i = 9.5e-9; // mean exitation energy per atomic Z (GeV) + Double_t p2=pTotal*pTotal; + Double_t beta2=p2/(p2 + muMass*muMass); + + Double_t w = k * rho * pathLength * atomicZ / atomicA / beta2; + + if (beta2/(1-beta2)>3.5*3.5) + return w * (log(2.*eMass*3.5/(i*atomicZ)) + 0.5*log(beta2/(1-beta2)) - beta2); + + return w * (log(2.*eMass*beta2/(1-beta2)/(i*atomicZ)) - beta2); +} - if (TMath::Abs(pZ) <= 0) return; - x[2] = zEnd/2; - x[0] = x[2]*trackParam->GetNonBendingSlope(); - x[1] = x[2]*trackParam->GetBendingSlope(); +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ) +{ + /// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal' + /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ' + Double_t muMass = 0.105658369; // GeV + //Double_t eMass = 0.510998918e-3; // GeV + Double_t k = 0.307075e-3; // GeV.g^-1.cm^2 + Double_t p2=pTotal*pTotal; + Double_t beta2=p2/(p2 + muMass*muMass); + + Double_t fwhm = 2. * k * rho * pathLength * atomicZ / atomicA / beta2; // FWHM of the energy loss Landau distribution + Double_t sigma2 = fwhm * fwhm / (8.*log(2.)); // gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma) + + //sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; // sigma2 of the energy loss gaussian distribution + + return sigma2; +} - // Take magn. field value at position x. - if (fgkField) fgkField->Field(x,b); - else { - cout<<"F-AliMUONTrackExtrap::FieldCorrection: fgkField = 0x0"<SetBendingSlope(pYNew/pZ); - trackParam->SetNonBendingSlope(pXNew/pZ); +//__________________________________________________________________________ +void AliMUONTrackExtrap::Cov2CovP(const TMatrixD ¶m, TMatrixD &cov) +{ + /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot) + /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system - trackParam->SetInverseBendingMomentum(c/TMath::Sqrt(pYNew*pYNew+pZ*pZ)); - + // charge * total momentum + Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) / + TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0); + + // Jacobian of the opposite transformation + TMatrixD jacob(5,5); + jacob.UnitMatrix(); + jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)); + jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) / + (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)); + jacob(4,4) = - qPTot / param(4,0); + + // compute covariances in new coordinate system + TMatrixD tmp(cov,TMatrixD::kMultTranspose,jacob); + cov.Mult(jacob,tmp); +} + +//__________________________________________________________________________ +void AliMUONTrackExtrap::CovP2Cov(const TMatrixD ¶m, TMatrixD &covP) +{ + /// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz) + /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system + + // charge * total momentum + Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) / + TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0); + + // Jacobian of the transformation + TMatrixD jacob(5,5); + jacob.UnitMatrix(); + jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)); + jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) / + (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)); + jacob(4,4) = - param(4,0) / qPTot; + + // compute covariances in new coordinate system + TMatrixD tmp(covP,TMatrixD::kMultTranspose,jacob); + covP.Mult(jacob,tmp); } //__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t *vect, Double_t *vout) { +/// 
 ///    ******************************************************************
 ///    *                                                                *
 ///    *  Performs the tracking of one step in a magnetic field         *
@@ -720,12 +929,13 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Doub
 ///    *   output                                                       *
 ///    *     VOUT = same as VECT after completion of the step           *
 ///    *                                                                *
-///    *    ==>Called by : , GUSWIM                               *
+///    *    ==>Called by : USER, GUSWIM                               *
 ///    *       Author    m.hansroul  *********                          *
 ///    *       modified  s.egli, s.v.levonian                           *
 ///    *       modified  v.perevoztchikov
 ///    *                                                                *
 ///    ******************************************************************
+///
// modif: everything in double precision @@ -824,6 +1034,7 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Doub //__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout) { +/// 
 ///	******************************************************************
 ///	*								 *
 ///	*	Tracking routine in a constant field oriented		 *
@@ -831,11 +1042,12 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Doub
 ///	*	Tracking is performed with a conventional		 *
 ///	*	helix step method					 *
 ///	*								 *
-///	*    ==>Called by : , GUSWIM				 *
+///	*    ==>Called by : USER, GUSWIM				 *
 ///	*	Authors    R.Brun, M.Hansroul  *********		 *
 ///	*	Rewritten  V.Perevoztchikov
 ///	*								 *
 ///	******************************************************************
+///
Double_t hxp[3]; Double_t h4, hp, rho, tet; @@ -896,9 +1108,11 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Doub return; } + //__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout) { +/// 
 ///	******************************************************************
 ///	*								 *
 ///	*  Runge-Kutta method for tracking a particle through a magnetic *
@@ -914,12 +1128,13 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step,
 ///	*  User routine called  					 *
 ///	*	CALL GUFLD(X,F) 					 *
 ///	*								 *
-///	*    ==>Called by : , GUSWIM				 *
+///	*    ==>Called by : USER, GUSWIM				 *
 ///	*	Authors    R.Brun, M.Hansroul  *********		 *
 ///	*		   V.Perevoztchikov (CUT STEP implementation)	 *
 ///	*								 *
 ///	*								 *
 ///	******************************************************************
+///
Double_t h2, h4, f[4]; Double_t xyzt[3], a, b, c, ph,ph2; @@ -1139,14 +1354,15 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, return; } + //___________________________________________________________ - void AliMUONTrackExtrap::GetField(Double_t *Position, Double_t *Field) +void AliMUONTrackExtrap::GetField(Double_t *Position, Double_t *Field) { /// interface for arguments in double precision (Why ? ChF) Float_t x[3], b[3]; - + x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2]; - + if (fgkField) fgkField->Field(x,b); else { cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<