#include <Riostream.h>
+using std::endl;
+using std::cout;
/// \cond CLASSIMP
ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
/// \endcond
Double_t b[3] = {0.,0.,0.};
TGeoGlobalMagField::Instance()->Field(x,b);
fgSimpleBValue = b[0];
- fgFieldON = fgSimpleBValue ? kTRUE : kFALSE;
+ fgFieldON = (TMath::Abs(fgSimpleBValue) > 1.e-10) ? kTRUE : kFALSE;
}
// Extrapolation loop (until within tolerance or the track turn around)
Double_t residue = zEnd - trackParam->GetZ();
Bool_t uturn = kFALSE;
+ Bool_t trackingFailed = kFALSE;
Bool_t tooManyStep = kFALSE;
while (TMath::Abs(residue) > fgkRungeKuttaMaxResidue && stepNumber <= fgkMaxStepNumber) {
}
stepNumber ++;
step = TMath::Abs(step);
- AliMUONTrackExtrap::ExtrapOneStepRungekutta(chargeExtrap,step,v3,v3New);
+ if (!AliMUONTrackExtrap::ExtrapOneStepRungekutta(chargeExtrap,step,v3,v3New)) {
+ trackingFailed = kTRUE;
+ break;
+ }
residue = zEnd - v3New[2];
step *= dZ/(v3New[2]-trackParam->GetZ());
} while (residue*dZ < 0 && TMath::Abs(residue) > fgkRungeKuttaMaxResidue);
- if (v3New[5]*v3[5] < 0) { // the track turned around
+ if (trackingFailed) break;
+ else if (v3New[5]*v3[5] < 0) { // the track turned around
cout<<"W-AliMUONTrackExtrap::ExtrapToZRungekutta: The track turned around"<<endl;
uturn = kTRUE;
break;
}
// terminate the extropolation with a straight line up to the exact "zEnd" value
- if (uturn) {
+ if (trackingFailed || uturn) {
// track ends +-100 meters away in the bending direction
dZ = zEnd - v3[2];
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2)
+void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t signedPathLength, Double_t f0, Double_t f1, Double_t f2)
{
/// Add to the track parameter covariances the effects of multiple Coulomb scattering
- /// The absorber correction parameters are supposed to be calculated at the current track z-position
+ /// signedPathLength must have the sign of (zOut - zIn) where all other parameters are assumed to be given at zOut.
// absorber related covariance parameters
Double_t bendingSlope = param->GetBendingSlope();
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 pathLength = TMath::Abs(signedPathLength);
Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
- Double_t covCorrSlope = alpha2 * (pathLength * f0 - f1);
+ Double_t covCorrSlope = TMath::Sign(1.,signedPathLength) * alpha2 * (pathLength * f0 - f1);
Double_t varSlop = alpha2 * f0;
// Set MCS covariance matrix
// Position of the Branson plane (spectro. (z<0))
Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
- // Add MCS effects to current parameter covariances
- AddMCSEffectInAbsorber(param, pathLength, f0, f1, f2);
+ // Add MCS effects to current parameter covariances (spectro. (z<0))
+ AddMCSEffectInAbsorber(param, -pathLength, f0, f1, f2);
// Get track parameters and covariances in the Branson plane corrected for magnetic field effect
ExtrapToZCov(param,zVtx);
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"
+ /// through a material of thickness "Abs(dZ)" and of radiation length "x0"
/// assuming linear propagation and using the small angle approximation.
+ /// dZ = zOut - zIn (sign is important) and "param" is assumed to be given zOut.
+ /// If x0 <= 0., assume dZ = pathLength/x0 and consider the material thickness as negligible.
Double_t bendingSlope = param->GetBendingSlope();
Double_t nonBendingSlope = param->GetNonBendingSlope();
(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;
+ Double_t signedPathLength = dZ * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope);
+ Double_t pathLengthOverX0 = (x0 > 0.) ? TMath::Abs(signedPathLength) / x0 : TMath::Abs(signedPathLength);
// 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));
- theta02 *= theta02 * inverseTotalMomentum2 * pathLength / x0;
+ Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
+ theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
- Double_t varCoor = pathLength2 * theta02 / 3.;
+ Double_t varCoor = (x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
Double_t varSlop = theta02;
- Double_t covCorrSlope = pathLength * theta02 / 2.;
+ Double_t covCorrSlope = (x0 > 0.) ? signedPathLength * theta02 / 2. : 0.;
// Set MCS covariance matrix
TMatrixD newParamCov(param->GetCovariances());
// Correct for energy loss add multiple scattering dispersion in covariance matrix
CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
- AddMCSEffectInAbsorber(trackParam, pathLength, f0, f1, f2);
+ AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
ExtrapToZCov(trackParam, trackXYZIn[2]);
CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
ExtrapToZCov(trackParam, zVtx);
} else {
// add multiple scattering dispersion in covariance matrix
- AddMCSEffectInAbsorber(trackParam, pathLength, f0, f1, f2);
+ AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
ExtrapToZCov(trackParam, zVtx);
}
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t *vect, Double_t *vout)
+void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
{
/// <pre>
/// ******************************************************************
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout)
+void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
{
/// <pre>
/// ******************************************************************
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout)
+Bool_t AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t* vect, Double_t* vout)
{
/// <pre>
/// ******************************************************************
/// </pre>
Double_t h2, h4, f[4];
- Double_t xyzt[3], a, b, c, ph,ph2;
+ Double_t xyzt[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
+ Double_t a, b, c, ph,ph2;
Double_t secxs[4],secys[4],seczs[4],hxp[3];
Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
Double_t est, at, bt, ct, cba;
vout[5] = cba*c;
rest = step - tl;
if (step < 0.) rest = -rest;
- if (rest < 1.e-5*TMath::Abs(step)) return;
+ if (rest < 1.e-5*TMath::Abs(step)) return kTRUE;
} while(1);
// angle too big, use helix
+ cout<<"W-AliMUONTrackExtrap::ExtrapOneStepRungekutta: Ruge-Kutta failed: switch to helix"<<endl;
f1 = f[0];
f2 = f[1];
f3 = f[2];
f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
+ if (f4 < 1.e-10) {
+ cout<<"E-AliMUONTrackExtrap::ExtrapOneStepRungekutta: magnetic field at (";
+ cout<<xyzt[0]<<", "<<xyzt[1]<<", "<<xyzt[2]<<") = "<<f4<<": giving up"<<endl;
+ return kFALSE;
+ }
rho = -f4*pinv;
tet = rho * step;
vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
- return;
+ return kTRUE;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff