X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrackExtrap.cxx;h=e7762acbff43e0f3e0e75dd3297a308ecaf9475b;hb=0d106d88a2aa0a0f6605b8c7f4c46aaa9fa05e2c;hp=2565ba14fa262594ec26c8ffb869f1df13d93e84;hpb=7dc68cf49f56629ebaeb98000a28382de2edcb14;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrackExtrap.cxx b/MUON/AliMUONTrackExtrap.cxx index 2565ba14fa2..e7762acbff4 100644 --- a/MUON/AliMUONTrackExtrap.cxx +++ b/MUON/AliMUONTrackExtrap.cxx @@ -15,41 +15,65 @@ /* $Id$ */ -/////////////////////////////////////////////////// -// -// Tools -// for -// track -// extrapolation -// in -// ALICE -// dimuon -// spectrometer -// -/////////////////////////////////////////////////// +//----------------------------------------------------------------------------- +// Class AliMUONTrackExtrap +// ------------------------ +// Tools for track extrapolation in ALICE dimuon spectrometer +// Author: Philippe Pillot +//----------------------------------------------------------------------------- #include "AliMUONTrackExtrap.h" #include "AliMUONTrackParam.h" #include "AliMUONConstants.h" +#include "AliMUONReconstructor.h" +#include "AliMUONRecoParam.h" #include "AliMagF.h" -#include #include -#include #include +#include + /// \cond CLASSIMP ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context /// \endcond const AliMagF* AliMUONTrackExtrap::fgkField = 0x0; +const Double_t AliMUONTrackExtrap::fgkSimpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ()); +const Double_t AliMUONTrackExtrap::fgkSimpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL()); + Double_t AliMUONTrackExtrap::fgSimpleBValue = 0.; + Bool_t AliMUONTrackExtrap::fgFieldON = kFALSE; const Bool_t AliMUONTrackExtrap::fgkUseHelix = kFALSE; const Int_t AliMUONTrackExtrap::fgkMaxStepNumber = 5000; const Double_t AliMUONTrackExtrap::fgkHelixStepLength = 6.; const Double_t AliMUONTrackExtrap::fgkRungeKuttaMaxResidue = 0.002; - //__________________________________________________________________________ +//__________________________________________________________________________ +void AliMUONTrackExtrap::SetField(const AliMagF* magField) +{ + /// set magnetic field + + // set field map + fgkField = magField; + if (!fgkField) { + cout<<"E-AliMUONTrackExtrap::SetField: fgkField = 0x0"<Factor() == 0.) ? kFALSE : kTRUE; + + // set field at the centre of the dipole + if (fgFieldON) { + Float_t b[3] = {0.,0.,0.}, x[3] = {50.,50.,(Float_t) fgkSimpleBPosition}; + fgkField->Field(x,b); + fgSimpleBValue = (Double_t) b[0]; + } else fgSimpleBValue = 0.; + +} + +//__________________________________________________________________________ Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum) { /// Returns impact parameter at vertex in bending plane (cm), @@ -59,21 +83,19 @@ Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingM if (bendingMomentum == 0.) return 1.e10; - Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ()); - Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL()); - Float_t b[3], x[3] = {0.,0.,(Float_t) simpleBPosition}; - if (fgkField) fgkField->Field(x,b); - else { + if (!fgkField) { cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<Field(x,b); - else { + if (!fgkField) { cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<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); + + // Update the propagator if required + if (updatePropagator) { + TMatrixD jacob(5,5); + jacob.UnitMatrix(); + jacob(0,1) = dZ; + jacob(2,3) = dZ; + trackParam->UpdatePropagator(jacob); + } + + } + +} + +//__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd) { /// Interface to track parameter extrapolation to the plane at "Z" using Helix or Rungekutta algorithm. /// On return, the track parameters resulting from the extrapolation are updated in trackParam. - if (fgkUseHelix) AliMUONTrackExtrap::ExtrapToZHelix(trackParam,zEnd); + if (!fgFieldON) AliMUONTrackExtrap::LinearExtrapToZ(trackParam,zEnd); + else if (fgkUseHelix) AliMUONTrackExtrap::ExtrapToZHelix(trackParam,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. @@ -127,7 +192,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 ???? @@ -157,7 +222,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. @@ -179,7 +244,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) { @@ -201,7 +266,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. @@ -218,7 +283,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", @@ -233,227 +298,230 @@ 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 + if (!fgFieldON) { // linear extrapolation if no magnetic field + AliMUONTrackExtrap::LinearExtrapToZ(trackParam,zEnd,updatePropagator); + return; + } + + // 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 energy loss corrections. - /// 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) return; // nothing to be done if already at vertex + // Position of the Branson plane (spectro. (z<0)) + Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.; - if (trackParam->GetZ() > zVtx) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<GetZ() - <<") upstream the vertex (zVtx = "<GetVolume("ALIC")->GetNode("ABSM_1"); - if (!absNode) { - cout<<"E-AliMUONTrackExtrap::ExtrapToVertexUncorrected: failed to get absorber node"<GetVolume()->GetShape()->GetAxisRange(3,kZAbsBeg,kZAbsEnd); - const Double_t *absPos = absNode->GetMatrix()->GetTranslation(); - kZAbsBeg = absPos[2] - kZAbsBeg; // spectro. (z<0) - kZAbsEnd = absPos[2] - kZAbsEnd; // spectro. (z<0) -*/ - static const Double_t kZAbsBeg = -90.; - static const Double_t kZAbsEnd = -505.; - - // Check the vertex position relatively to the absorber - if (zVtx < kZAbsBeg && zVtx > kZAbsEnd) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<GetZ() > kZAbsBeg) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() - <<") upstream the front absorber (zAbsorberBegin = "<GetZ() > kZAbsEnd) { // spectro. (z<0) - cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() - <<") inside the front absorber ("<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)); - // Then add MCS effect in absorber to the parameters covariances - AliMUONTrackParam trackParamIn(*trackParam); - ExtrapToZ(&trackParamIn, TMath::Min(zVtx, kZAbsBeg)); - Double_t trackXYZIn[3]; - trackXYZIn[0] = trackParamIn.GetNonBendingCoor(); - trackXYZIn[1] = trackParamIn.GetBendingCoor(); - trackXYZIn[2] = trackParamIn.GetZ(); - Double_t trackXYZOut[3]; - trackXYZOut[0] = trackParam->GetNonBendingCoor(); - trackXYZOut[1] = trackParam->GetBendingCoor(); - trackXYZOut[2] = trackParam->GetZ(); - Double_t pathLength = 0.; - Double_t f0 = 0.; - Double_t f1 = 0.; - Double_t f2 = 0.; - Double_t meanRho = 0.; - GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho); - AddMCSEffectInAbsorber(trackParam,pathLength,f0,f1,f2); + // 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); - // finally go to the vertex - ExtrapToZCov(trackParam,zVtx); + // 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::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2) +//__________________________________________________________________________ +void AliMUONTrackExtrap::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss, Double_t sigmaELoss2) { - /// Add to the track parameter covariances the effects of multiple Coulomb scattering - /// at the end of the front absorber using the absorber correction parameters + /// Correct parameters for energy loss and add energy loss fluctuation effect to covariances - // absorber related covariance parameters - Double_t bendingSlope = param->GetBendingSlope(); + // 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 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; + 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)); - TMatrixD* paramCov = param->GetCovariances(); - // Non bending plane - (*paramCov)(0,0) += varCoor; (*paramCov)(0,1) += covCorrSlope; - (*paramCov)(1,0) += covCorrSlope; (*paramCov)(1,1) += varSlop; - // Bending plane - (*paramCov)(2,2) += varCoor; (*paramCov)(2,3) += covCorrSlope; - (*paramCov)(3,2) += covCorrSlope; (*paramCov)(3,3) += varSlop; + // 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 &pathLength, - Double_t &f0, Double_t &f1, Double_t &f2, Double_t &meanRho) +//__________________________________________________________________________ +Bool_t 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 between track positions trackXYZIn and trackXYZOut + /// 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.; @@ -461,18 +529,20 @@ void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Doub f1 = 0.; f2 = 0.; meanRho = 0.; + totalELoss = 0.; + sigmaELoss2 = 0.; // Check whether the geometry is available if (!gGeoManager) { cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<InitTrack(trackXYZIn, b); if (!currentnode) { cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<GetDensity(); x0 = material->GetRadLen(); if (!material->IsMixture()) x0 /= rho; // different normalization in the modeler for mixture + atomicA = material->GetA(); + atomicZ = material->GetZ(); // Get path length within this material gGeoManager->FindNextBoundary(remainingPathLength); @@ -506,8 +580,8 @@ void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Doub currentnode = gGeoManager->Step(); if (!currentnode) { cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<IsEntering()) { // make another small step to try to enter in new absorber slice @@ -515,8 +589,8 @@ void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Doub currentnode = gGeoManager->Step(); if (!gGeoManager->IsEntering() || !currentnode) { cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"< TGeoShape::Tolerance()); meanRho /= pathLength; + + return kTRUE; } - //__________________________________________________________________________ +//__________________________________________________________________________ +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::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0) { /// Add to the track parameter covariances the effects of multiple Coulomb scattering @@ -548,9 +648,10 @@ void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Dou 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; @@ -560,154 +661,187 @@ void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Dou 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, - Bool_t CorrectForMCS, Bool_t CorrectForEnergyLoss) +//__________________________________________________________________________ +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) { - /// Extrapolation to the vertex. - /// Returns the track parameters resulting from the extrapolation of the current TrackParam. - /// Changes parameters according to Branson correction through the absorber and energy loss + /// 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<<"F-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() - <<") upstream the vertex (zVtx = "<GetVolume("ALIC")->GetNode("ABSM_1"); - if (!absNode) { - cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: failed to get absorber node"<GetZ() + <<") upstream the vertex (zVtx = "<GetVolume()->GetShape()->GetAxisRange(3,kZAbsBeg,kZAbsEnd); - const Double_t *absPos = absNode->GetMatrix()->GetTranslation(); - kZAbsBeg = absPos[2] - kZAbsBeg; // spectro. (z<0) - kZAbsEnd = absPos[2] - kZAbsEnd; // spectro. (z<0) -*/ - static const Double_t kZAbsBeg = -90.; - static const Double_t kZAbsEnd = -505.; // Check the vertex position relatively to the absorber - if (zVtx < kZAbsBeg && zVtx > kZAbsEnd) { // spectro. (z<0) + 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() > kZAbsBeg) { // spectro. (z<0) + 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() > kZAbsEnd) { // spectro. (z<0) + } 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()); } - // Get absorber correction parameters assuming linear propagation from vertex to the track position + // 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]; - trackXYZIn[2] = TMath::Min(zVtx, kZAbsBeg); // 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]); - Double_t pathLength = 0.; - Double_t f0 = 0.; - Double_t f1 = 0.; - Double_t f2 = 0.; - Double_t meanRho = 0.; - GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho); - - // Calculate energy loss + 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 charge = TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum()); - Double_t deltaP = TotalMomentumEnergyLoss(pTot,pathLength,meanRho); - - // Correct for half of energy loss - Double_t nonBendingSlope, bendingSlope; - if (CorrectForEnergyLoss) { - pTot += 0.5 * deltaP; - nonBendingSlope = trackParam->GetNonBendingSlope(); - bendingSlope = trackParam->GetBendingSlope(); - trackParam->SetInverseBendingMomentum(charge / pTot * - TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) / - TMath::Sqrt(1.0 + bendingSlope*bendingSlope)); + Double_t pathLength, f0, f1, f2, meanRho, deltaP, sigmaDeltaP2; + if (!GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,deltaP,sigmaDeltaP2)) { + cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: Unable to take into account the absorber effects"<CovariancesExist()) ExtrapToZCov(trackParam,zVtx); + else ExtrapToZ(trackParam,zVtx); + return; } - if (CorrectForMCS) { - // Position of the Branson plane (spectro. (z<0)) - Double_t zB = (f1>0.) ? trackXYZIn[2] - f2/f1 : 0.; - - // Get track position in the Branson plane corrected for magnetic field effect - ExtrapToZ(trackParam,zVtx); - Double_t xB = trackParam->GetNonBendingCoor() + (zB - zVtx) * trackParam->GetNonBendingSlope(); - Double_t yB = trackParam->GetBendingCoor() + (zB - zVtx) * trackParam->GetBendingSlope(); + // Compute track parameters and covariances at vertex according to correctForMCS and correctForEnergyLoss flags + if (correctForMCS) { - // Get track slopes corrected for multiple scattering (spectro. (z<0)) - nonBendingSlope = (zB<0.) ? (xB - xVtx) / (zB - zVtx) : trackParam->GetNonBendingSlope(); - bendingSlope = (zB<0.) ? (yB - yVtx) / (zB - zVtx) : trackParam->GetBendingSlope(); + 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); + } - // Set track parameters at vertex - trackParam->SetNonBendingCoor(xVtx); - trackParam->SetBendingCoor(yVtx); - trackParam->SetZ(zVtx); - trackParam->SetNonBendingSlope(nonBendingSlope); - trackParam->SetBendingSlope(bendingSlope); } else { - ExtrapToZ(trackParam,zVtx); - nonBendingSlope = trackParam->GetNonBendingSlope(); - bendingSlope = trackParam->GetBendingSlope(); + + if (correctForEnergyLoss) { + + // Correct for energy loss add multiple scattering dispersion in covariance matrix + 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 { + + // add multiple scattering dispersion in covariance matrix + AddMCSEffectInAbsorber(trackParam, pathLength, f0, f1, f2); + ExtrapToZCov(trackParam, zVtx); + + } + } - // Correct for second half of energy loss - if (CorrectForEnergyLoss) pTot += 0.5 * deltaP; - - // Set track parameters at vertex - trackParam->SetInverseBendingMomentum(charge / pTot * - TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) / - TMath::Sqrt(1.0 + bendingSlope*bendingSlope)); - } - //__________________________________________________________________________ +//__________________________________________________________________________ +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); +} + +//__________________________________________________________________________ +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); +} + +//__________________________________________________________________________ +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(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx) { /// Calculate the total momentum energy loss in-between the track position and the vertex assuming a linear propagation @@ -729,41 +863,96 @@ Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(AliMUONTrackParam* trackPar trackXYZIn[0] = xVtx; trackXYZIn[1] = yVtx; trackXYZIn[2] = zVtx; - Double_t pathLength = 0.; - Double_t f0 = 0.; - Double_t f1 = 0.; - Double_t f2 = 0.; - Double_t meanRho = 0.; - GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho); - - // Calculate energy loss Double_t pTot = trackParam->P(); - return TotalMomentumEnergyLoss(pTot,pathLength,meanRho); + Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2; + GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2); + + return totalELoss; } - //__________________________________________________________________________ -Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(Double_t pTotal, Double_t pathLength, Double_t rho) +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ) { - /// Returns the total momentum energy loss in the front absorber - Double_t muMass = 0.105658369; + /// 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 dE=ApproximateBetheBloch(beta2)*pathLength*rho; - return dE; + 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); } - //__________________________________________________________________________ -Double_t AliMUONTrackExtrap::ApproximateBetheBloch(Double_t beta2) +//__________________________________________________________________________ +Double_t AliMUONTrackExtrap::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicA, Double_t atomicZ) { -/// This is an approximation of the Bethe-Bloch formula with -/// the density effect taken into account at beta*gamma > 3.5 -/// (the approximation is reasonable only for solid materials) + /// 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; +} - if (beta2/(1-beta2)>3.5*3.5) - return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2); +//__________________________________________________________________________ +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 + + // 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); +} - return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2); +//__________________________________________________________________________ +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); } //__________________________________________________________________________ @@ -783,7 +972,7 @@ 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 @@ -896,7 +1085,7 @@ 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 /// * * @@ -982,7 +1171,7 @@ 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) * /// * * @@ -1210,13 +1399,13 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, } //___________________________________________________________ - 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"<