From [email protected] and [email protected]

in PANDA we prepared some additional bug fixes for GEANE

There are two corrections:
1) the propagation to length used to contain a bug so that if the
propagation to length was chosen but the given length was too small
(or 0, as limit case) it did not work.
We fixed it in such a way that if the track length is too small but not 0
the step size is set to PREC, if it is <= 0 it is set to 0 and the start
point is coincident with the end point.
The flag IERR has been provided of the additional value 5: in this case,
when start and end position coincide, IERR is set to 5 and this means that
the error matrix does not have to be updated. In this way if track length
is 0 the particle just does not move.
2) it has been added the possibility to track along the z axis, though
with an approximation. Up to now, due to the matrix elements with
1/cos(lambda), the propagation exactly along z axis (or below a certain
limit) was forbidden. The approximation we introduced is in TGeant3gu.gxx,
in eustep() function: if the theta angle drops below 1e-4 deg it is forced
to this limit value in order to be able to go on tracking. Obviously this
means that if you try to track exactly along z (theta = 0) you will be
forced to track at 1e-4 deg instead. 

We checked the pulls and residual values for a plane @12M and @26m and 
they both with and without magnetic field and they are acceptable.

TGeant3/TGeant3gu.cxx

@@ -622,6 +622,51 @@ void eustep(){
    if(geant3->Gctrak()->sleng == 0) icc = 0;
    cflag=geant3->Gcmore()->iclose;
    if (geant3->Gcflag()->idebug * geant3->Gcflag()->iswit[2] != 0)geant3->Erxyzc();
+
+   // ---------------------------------------
+   // A. Rotondi and L. Lavezzi (sept 2010)
+   // changes to be able to approximate a solution when propagating
+   // along z axis. In that case lambda is 90 deg and 1/cos(lambda)
+   // becomes inf. Here when the theta angle drops below 1e-4 deg
+   // we set it to 1e-4 and this allows to propagate (with a small
+   // approximation) in the forward direction.
+
+   // direction vector set unitary
+   TVector3 dir(geant3->Gctrak()->vect[3], 
+		geant3->Gctrak()->vect[4],
+		geant3->Gctrak()->vect[5]);
+   dir.SetMag(1.);
+
+   if(dir.Mag() != 0.) {
+
+     // if the direction is REALLY along z, (0., 0., 1), change it
+     // a little
+     if(dir.Z() == 1) dir.SetMagThetaPhi(1., 1e-9, gRandom->Uniform(0., 2 * TMath::Pi())); 
+
+     Double_t cosLam0 = TMath::Sin(dir.Theta());
+
+     // cos(lambda) < 2.e-6 means theta < 1e-4 deg
+     if(fabs(cosLam0) < 2.e-6) {
+
+       // cos(lambda) @ limit angle
+       Double_t cosLam = TMath::Sign(1., cosLam0) * 2.e-6;
+       Double_t sinLam = TMath::Sqrt(1 - cosLam * cosLam);
+
+       // px, py, pz reset to limit angle
+       dir.SetX(dir.X() * (cosLam/cosLam0));
+       dir.SetY(dir.Y() * (cosLam/cosLam0));
+       dir.SetZ(sinLam);
+       dir.SetMag(1.);
+       geant3->Gctrak()->vect[3] = dir.X();
+       geant3->Gctrak()->vect[4] = dir.Y();
+       geant3->Gctrak()->vect[5] = dir.Z();
+
+     }
+   }
+   // ---------------------------------------
+
+
+
    if(cflag==1){
 // distance between the track point and the point
       if(icc==0) prdist2=geant3->Gconst()->big;

erdecks/ertrch.F

@@ -208,7 +208,7 @@ SUBROUTINE ERTRCH
       ASCL1  = BIG
       DO 20 IPR = 1,NEPRED
          STEPE  = BIG
-         IF (LELENG) STEPE = ERLENG(IPR) - SLENG
+c         IF (LELENG) STEPE = ERLENG(IPR) - SLENG
          IF (LEPLAN) THEN
             SCAL1 = 0.
             SCAL2 = 0.
@@ -219,6 +219,24 @@ SUBROUTINE ERTRCH
             STEPE = SCAL1/SCAL2
          ENDIF
          IF (STEPE.LE.PREC) STEPE = BIG
+c *******************************
+c     A. Rotondi and L. Lavezzi (sept 2010)
+c     when leleng (track length propagation) 
+c     if track length (ERLENG) <= 0 then the
+c     particle must not move; else if the step
+c     is <= PREC it is set = to it
+         IF (LELENG) THEN
+            STEPE = ERLENG(IPR) - SLENG
+            IF(ERLENG(IPR).LE.0.) THEN
+               STEPE = 0
+            ELSE IF(STEPE.LE.PREC) THEN
+               STEPE = PREC
+            ENDIF
+         ENDIF
+c *******************************
+
+
+
          IF (STEPE.LT.STEPER) THEN
             STEPER = STEPE
             INLIST = IPR

erdecks/ertrnt.F

@@ -37,7 +37,7 @@ SUBROUTINE ERTRNT
       STEPER = BIG
       DO 20 IPR = 1,NEPRED
          STEPE  = BIG
-         IF (LELENG) STEPE = ERLENG(IPR) - SLENG
+c         IF (LELENG) STEPE = ERLENG(IPR) - SLENG
          IF (LEPLAN) THEN
             SCAL1 = 0.
             SCAL2 = 0.
@@ -48,6 +48,22 @@ SUBROUTINE ERTRNT
             STEPE = SCAL1/SCAL2
          ENDIF
          IF (STEPE.LE.PREC) STEPE = BIG
+c *********************
+c     A. Rotondi and L. Lavezzi (sept 2010)
+c     when leleng (track length propagation) 
+c     if track length (ERLENG) <= 0 then the
+c     particle must not move; else if the step
+c     is <= PREC it is set = to it
+         IF (LELENG) THEN
+            STEPE = ERLENG(IPR) - SLENG
+            IF(ERLENG(IPR).LE.0.) THEN
+               STEPE = 0
+            ELSE IF(STEPE.LE.PREC) THEN
+               STEPE = PREC
+            ENDIF
+         ENDIF
+c *********************
+
          IF (STEPE.LT.STEPER) THEN
             STEPER = STEPE
             INLIST = IPR

erpremc/trprfn.F

@@ -65,6 +65,7 @@ SUBROUTINE TRPRFN(X1,P1,H1,X2,P2,H2,CH,XL,R,MVAR,IFLAG,ITRAN,IERR)
 *               2    MOMENTUM IS ZERO
 *               3    H*ALFA/P AT X1 AND X2 DIFFER TOO MUCH
 *               4    PARTICLE MOVES IN Z - DIRECTION
+*               5    INITIAL AND FINAL POINTS ARE THE SAME
 *
 ************************************************************************
 *
@@ -120,10 +121,21 @@ SUBROUTINE TRPRFN(X1,P1,H1,X2,P2,H2,CH,XL,R,MVAR,IFLAG,ITRAN,IERR)
    11 IF(MVAR.NE.0) GO TO 901
    19 GO TO 900
 *
-* *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
-
+* *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES 
+*
+*
 *
-   20 PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
+*     A. Rotondi and L. Lavezzi (sept 2010)
+*     if the first and last point coincide
+*     set IERR to 5 and exit
+ 20   IF((X1(1).EQ.X2(1)).AND.(X1(2).EQ.X2(2))
+     +     .AND.(X1(3).EQ.X2(3))) THEN
+         IERR = 5
+         GOTO 999
+      ENDIF
+*     ***********************************
+
+      PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
       PA2=SQRT(P2(1)**2+P2(2)**2+P2(3)**2)
       IF(PA1*PA2.EQ.0.) GO TO 902
 C      DPA = PA2 - PA1