/* $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 "AliMagF.h"
-#include <Riostream.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TGeoManager.h>
+#include <Riostream.h>
+
/// \cond CLASSIMP
ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
/// \endcond
Double_t simpleBValue = (Double_t) b[0];
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::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"<<endl;
+ // Extrapolate track parameters to "zEnd"
+ ExtrapToZ(trackParam,zEnd);
+ return;
+ }
+
// Save the actual track parameters
AliMUONTrackParam trackParamSave(*trackParam);
- Double_t nonBendingCoor = trackParamSave.GetNonBendingCoor();
- Double_t nonBendingSlope = trackParamSave.GetNonBendingSlope();
- Double_t bendingCoor = trackParamSave.GetBendingCoor();
- Double_t bendingSlope = trackParamSave.GetBendingSlope();
- Double_t inverseBendingMomentum = trackParamSave.GetInverseBendingMomentum();
- Double_t zBegin = trackParamSave.GetZ();
+ TMatrixD paramSave(trackParamSave.GetParameters());
+ Double_t zBegin = trackParamSave.GetZ();
+
+ // Get reference to the parameter covariance matrix
+ const TMatrixD& kParamCov = trackParam->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"<<endl;
- //cout<<" -> nothing to extrapolate !!"<<endl;
- return;
- }
- TMatrixD* paramCov = trackParam->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);
}
exit(-1);
}
- // Check whether the geometry is available and get absorber boundaries
- if (!gGeoManager) {
- cout<<"E-AliMUONTrackExtrap::ExtrapToVertexUncorrected: no TGeo"<<endl;
- return;
- }
- TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1");
- if (!absNode) {
- cout<<"E-AliMUONTrackExtrap::ExtrapToVertexUncorrected: failed to get absorber node"<<endl;
- return;
- }
- Double_t zAbsBeg, zAbsEnd;
- absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd);
- const Double_t *absPos = absNode->GetMatrix()->GetTranslation();
- zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0)
- zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0)
-
// Check the vertex position relatively to the absorber
- if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0)
+ if (zVtx < AliMUONConstants::AbsZBeg() && zVtx > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
- <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
- } else if (zVtx < zAbsEnd ) { // spectro. (z<0)
+ <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
+ } else if (zVtx < AliMUONConstants::AbsZEnd() ) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
- <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl;
- ExtrapToZCov(trackParam,zVtx);
+ <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::AbsZEnd()<<")"<<endl;
+ if (trackParam->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() > zAbsBeg) { // spectro. (z<0)
+ if (trackParam->GetZ() > AliMUONConstants::AbsZBeg()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
- <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl;
- ExtrapToZCov(trackParam,zVtx);
+ <<") upstream the front absorber (zAbsorberBegin = "<<AliMUONConstants::AbsZBeg()<<")"<<endl;
+ if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
+ else ExtrapToZ(trackParam,zVtx);
return;
- } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0)
+ } else if (trackParam->GetZ() > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
- <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
+ <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
} else {
- ExtrapToZCov(trackParam,zAbsEnd);
+ if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,AliMUONConstants::AbsZEnd());
+ else ExtrapToZ(trackParam,AliMUONConstants::AbsZEnd());
}
// Then add MCS effect in absorber to the parameters covariances
AliMUONTrackParam trackParamIn(*trackParam);
- ExtrapToZ(&trackParamIn, TMath::Min(zVtx, zAbsBeg));
+ ExtrapToZ(&trackParamIn, TMath::Min(zVtx, AliMUONConstants::AbsZBeg()));
Double_t trackXYZIn[3];
trackXYZIn[0] = trackParamIn.GetNonBendingCoor();
trackXYZIn[1] = trackParamIn.GetBendingCoor();
Double_t f1 = 0.;
Double_t f2 = 0.;
Double_t meanRho = 0.;
- GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho);
+ Double_t totalELoss = 0.;
+ GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,trackParam->P(),pathLength,f0,f1,f2,meanRho,totalELoss);
AddMCSEffectInAbsorber(trackParam,pathLength,f0,f1,f2);
// finally go to the vertex
Double_t covCorrSlope = alpha2 * (pathLength * f0 - f1);
Double_t varSlop = alpha2 * f0;
- TMatrixD* paramCov = param->GetCovariances();
+ 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;
+
+ // Set new 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)
+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)
{
/// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber
/// Calculated assuming a linear propagation between track positions trackXYZIn and trackXYZOut
// 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.;
f1 = 0.;
f2 = 0.;
meanRho = 0.;
+ totalELoss = 0.;
// Check whether the geometry is available
if (!gGeoManager) {
// 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];
rho = material->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);
f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.;
f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
meanRho += localPathLength * rho;
+ totalELoss += BetheBloch(pTotal, localPathLength, rho, atomicA, atomicZ);
// prepare next step
zB = zE;
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::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0)
{
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.;
+
+ // Add effects of multiple Coulomb scattering in track parameter covariances
+ 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;
+ // Set new covariances
+ param->SetCovariances(newParamCov);
}
//__________________________________________________________________________
return;
}
- // Check whether the geometry is available and get absorber boundaries
- if (!gGeoManager) {
- cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: no TGeo"<<endl;
- return;
- }
- TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1");
- if (!absNode) {
- cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: failed to get absorber node"<<endl;
- return;
- }
- Double_t zAbsBeg, zAbsEnd;
- absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd);
- const Double_t *absPos = absNode->GetMatrix()->GetTranslation();
- zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0)
- zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0)
-
// Check the vertex position relatively to the absorber
- if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0)
+ if (zVtx < AliMUONConstants::AbsZBeg() && zVtx > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
- <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
- } else if (zVtx < zAbsEnd ) { // spectro. (z<0)
+ <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
+ } else if (zVtx < AliMUONConstants::AbsZEnd() ) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
- <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl;
+ <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::AbsZEnd()<<")"<<endl;
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() > zAbsBeg) { // spectro. (z<0)
+ if (trackParam->GetZ() > AliMUONConstants::AbsZBeg()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
- <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl;
+ <<") upstream the front absorber (zAbsorberBegin = "<<AliMUONConstants::AbsZBeg()<<")"<<endl;
ExtrapToZ(trackParam,zVtx);
return;
- } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0)
+ } else if (trackParam->GetZ() > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
- <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
+ <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
} else {
- ExtrapToZ(trackParam,zAbsEnd);
+ ExtrapToZ(trackParam,AliMUONConstants::AbsZEnd());
}
// Get absorber correction parameters assuming linear propagation from vertex to the track position
trackXYZOut[1] = trackParam->GetBendingCoor();
trackXYZOut[2] = trackParam->GetZ();
Double_t trackXYZIn[3];
- trackXYZIn[2] = TMath::Min(zVtx, zAbsBeg); // spectro. (z<0)
+ 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]);
+ 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.;
- GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho);
-
- // Calculate energy loss
- Double_t pTot = trackParam->P();
- Double_t charge = TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
- Double_t deltaP = TotalMomentumEnergyLoss(pTot,pathLength,meanRho);
+ Double_t deltaP = 0.;
+ GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,deltaP);
// Correct for half of energy loss
Double_t nonBendingSlope, bendingSlope;
+ Double_t charge = TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
if (CorrectForEnergyLoss) {
pTot += 0.5 * deltaP;
nonBendingSlope = trackParam->GetNonBendingSlope();
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.;
- GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho);
+ Double_t totalELoss = 0.;
+ GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss);
- // Calculate energy loss
- Double_t pTot = trackParam->P();
- return TotalMomentumEnergyLoss(pTot,pathLength,meanRho);
+ 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 AliMUONTrackExtrap::ApproximateBetheBloch(Double_t beta2) {
- //------------------------------------------------------------------
- // 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)
- //------------------------------------------------------------------
+ Double_t w = k * rho * pathLength * atomicZ / atomicA / beta2;
+
if (beta2/(1-beta2)>3.5*3.5)
- return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
+ return w * (log(2.*eMass*3.5/(i*atomicZ)) + 0.5*log(beta2/(1-beta2)) - beta2);
- return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
+ return w * (log(2.*eMass*beta2/(1-beta2)/(i*atomicZ)) - beta2);
}
//__________________________________________________________________________
void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t *vect, Double_t *vout)
{
+/// <pre>
/// ******************************************************************
/// * *
/// * Performs the tracking of one step in a magnetic field *
/// * output *
/// * VOUT = same as VECT after completion of the step *
/// * *
-/// * ==>Called by : <USER>, GUSWIM *
+/// * ==>Called by : USER, GUSWIM *
/// * Author m.hansroul ********* *
/// * modified s.egli, s.v.levonian *
/// * modified v.perevoztchikov
/// * *
/// ******************************************************************
+/// </pre>
// modif: everything in double precision
//__________________________________________________________________________
void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout)
{
+/// <pre>
/// ******************************************************************
/// * *
/// * Tracking routine in a constant field oriented *
/// * Tracking is performed with a conventional *
/// * helix step method *
/// * *
-/// * ==>Called by : <USER>, GUSWIM *
+/// * ==>Called by : USER, GUSWIM *
/// * Authors R.Brun, M.Hansroul ********* *
/// * Rewritten V.Perevoztchikov
/// * *
/// ******************************************************************
+/// </pre>
Double_t hxp[3];
Double_t h4, hp, rho, tet;
//__________________________________________________________________________
void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout)
{
+/// <pre>
/// ******************************************************************
/// * *
/// * Runge-Kutta method for tracking a particle through a magnetic *
/// * User routine called *
/// * CALL GUFLD(X,F) *
/// * *
-/// * ==>Called by : <USER>, GUSWIM *
+/// * ==>Called by : USER, GUSWIM *
/// * Authors R.Brun, M.Hansroul ********* *
/// * V.Perevoztchikov (CUT STEP implementation) *
/// * *
/// * *
/// ******************************************************************
+/// </pre>
Double_t h2, h4, f[4];
Double_t xyzt[3], a, b, c, ph,ph2;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff