Lightweight cascade checks -> code implementation

[u/mrichter/AliRoot.git] / MUON / AliMUONTrackExtrap.cxx

diff --git a/MUON/AliMUONTrackExtrap.cxx b/MUON/AliMUONTrackExtrap.cxx
index a9cc2a2f50040157782c04cc0dec0cd5ead692b8..45dcc35c0c02e311e4bd3aba9bec84d1b19baf5f 100644 (file)
--- a/MUON/AliMUONTrackExtrap.cxx
+++ b/MUON/AliMUONTrackExtrap.cxx
@@ -37,6 +37,8 @@
 
 #include <Riostream.h>
 
+using std::endl;
+using std::cout;
 /// \cond CLASSIMP
 ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
 /// \endcond
@@ -59,7 +61,7 @@ void AliMUONTrackExtrap::SetField()
   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;
   
 }
 
@@ -239,6 +241,7 @@ Bool_t AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Do
   // 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) {
     
@@ -256,12 +259,16 @@ Bool_t AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Do
       }
       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;
@@ -270,7 +277,7 @@ Bool_t AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Do
   }
   
   // 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];
@@ -609,6 +616,7 @@ Bool_t AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Do
   Double_t atomicZoverA = 0.; // Z/A of material
   Double_t localPathLength = 0;
   Double_t remainingPathLength = pathLength;
+  Double_t sigmaELoss = 0.;
   Double_t zB = trackXYZIn[2];
   Double_t zE, dzB, dzE;
   do {
@@ -664,7 +672,7 @@ Bool_t AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Do
     f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
     meanRho += localPathLength * rho;
     totalELoss += BetheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
-    sigmaELoss2 += EnergyLossFluctuation2(pTotal, localPathLength, rho, atomicZoverA);
+    sigmaELoss += EnergyLossFluctuation(pTotal, localPathLength, rho, atomicZoverA);
     
     // prepare next step
     zB = zE;
@@ -672,6 +680,7 @@ Bool_t AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Do
   } while (remainingPathLength > TGeoShape::Tolerance());
   
   meanRho /= pathLength;
+  sigmaELoss2 = sigmaELoss*sigmaELoss;
   
   return kTRUE;
 }
@@ -956,7 +965,7 @@ Double_t AliMUONTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Do
 }
 
 //__________________________________________________________________________
-Double_t AliMUONTrackExtrap::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZoverA)
+Double_t AliMUONTrackExtrap::EnergyLossFluctuation(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZoverA)
 {
   /// 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'
@@ -967,11 +976,11 @@ Double_t AliMUONTrackExtrap::EnergyLossFluctuation2(Double_t pTotal, Double_t pa
   Double_t beta2=p2/(p2 + muMass*muMass);
   
   Double_t fwhm = 2. * k * rho * pathLength * atomicZoverA / 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)
+  Double_t sigma = fwhm / TMath::Sqrt(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;
+  return sigma;
 }
 
 //__________________________________________________________________________
@@ -1021,7 +1030,7 @@ void AliMUONTrackExtrap::CovP2Cov(const TMatrixD &param, TMatrixD &covP)
 }
 
  //__________________________________________________________________________
-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>
 ///    ******************************************************************
@@ -1140,7 +1149,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)
+void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
 {
 /// <pre>
 ///    ******************************************************************
@@ -1218,7 +1227,7 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Doub
 }
 
  //__________________________________________________________________________
-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>
 ///    ******************************************************************
@@ -1245,7 +1254,8 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step,
 /// </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;
@@ -1416,16 +1426,22 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step,
       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;
  
@@ -1459,6 +1475,6 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t 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;
 }