anchor.diff:  !++  to simple pencil beam
   !
             IF (ipencil .eq. 0) bool_create=.true.
  -!
  -!! pencil beam stuffffffffffffffffffffffffffffffffff  (~80 lines)
  -!c ----------------------------pencil beam generatio at the crystal--------------------------------------
  -!          IF ( (ICOLL.EQ.IPENCIL
  -!     1           .AND. ITURN.EQ.1)
  -!     2           ) THEN
  -!            X    = pencil_dx(ICOLL)
  -!            XP = cry_tilt0 !valentina (if the beam is generated @ crystal, I want it to have the natural beam divergence)
  -!            bool_create=.true.
  -!c
  -!            AMPLZ = RAN_GAUSS(3.)!*pencil_spread(ICOLL)
  -!            Z     = AMPLZ
  -!            ZP   = 0.
  -!            dpop = 0.
  -!c
  -!            if (write_c_out) then
  -!c++             I want to write the original coordinates in a file... I have to transform back all coordinates....Include collimator tilt
  -!c
  -!              IF (tiltangle.GT.0.) THEN
  -!                x_print  = X  + tiltangle*C_LENGTH
  -!                XP_print = XP + tiltangle
  -!              ELSEIF (tiltangle.LT.0.) THEN
  -!                x_print  = X + tiltangle*C_LENGTH
  -!                XP_print = XP + tiltangle
  -!                x_print  = X - SIN(tiltangle) * C_LENGTH
  -!              ELSE
  -!                x_print = x
  -!                xp_print = xp
  -!              ENDIF
  -!
  -!c++  Transform back to particle coordinates with opening and offset
  -!c
  -!              x_print = x_print + C_APERTURE/2 + MIRROR*C_OFFSET
  -!c++  Now mirror at the horizontal axis for negative X offset
  -!c
  -!              x_print    = MIRROR * x_print
  -!              XP_print   = MIRROR * xp_print
  -!
  -!c++  Last do rotation into collimator frame
  -!c
  -!              Y_print  = Z  *COS(-1.*C_ROTATION) -
  -!     1                     X_print  *SIN(-1.*C_ROTATION)
  -!              x_print  = x_print  *COS(-1.*C_ROTATION) +
  -!     1                     Z  *SIN(-1.*C_ROTATION)
  -!              YP_print = ZP *COS(-1.*C_ROTATION) -
  -!     1                     XP_print *SIN(-1.*C_ROTATION)
  -!              XP_print = XP_print *COS(-1.*C_ROTATION) +
  -!     1                     ZP *SIN(-1.*C_ROTATION)
  -!c              WRITE(*,'(i4,2x,i4,2x,a,2x,5(f15.8,2x))')
  -!c     2          ITURN,ICOLL,C_MATERIAL,x_print,XP_print,Y_print,YP_print
  -!c     3          ,P
  -!C
  -!c++ Then drift forward of half lenght
  -!c
  -!              x_print= x_print+xp_print*C_length/2        !valentina drift removed
  -!              y_print= y_print+yp_print*C_length/2
  -!              WRITE(881,'(i4,2x,i4,2x,a,2x,5(f15.8,2x))')
  -!     2        ITURN,ICOLL,C_MATERIAL,x_print,XP_print,Y_print,YP_print
  -!     3        ,P
  -!
  -!              X_NORM=x_print/ SQRT(bx)/sqrt(EMITX0)
  -!              XP_NORM=(x_print*ax+xp_print*bx)/SQRT(bx)/sqrt(EMITX0)
  -!              Y_NORM=y_print/ SQRT(by)/sqrt(EMITY0)
  -!              YP_NORM=(y_print*AY+yp_print*by)/SQRT(by)/sqrt(EMITY0)
  -!c
  -!              WRITE(883,'(i4,2x,i4,2x,a,2x,7(f15.8,2x))')
  -!     2          ITURN,ICOLL,C_MATERIAL,X_NORM,XP_NORM,Y_NORM,YP_NORM,
  -!     3          SQRT(X_NORM**2+XP_NORM**2),
  -!     4          SQRT(Y_NORM**2+YP_NORM**2),P
  -!c
  -!c
  -!            endif
  -!          ENDIF
  -!c.--------------------end of pencil beam stuff-------------------------------
  -
   
   !-valentina for the first impact file
   !        s_in0(j)   = s_in(j)                         !daniele SEE COMMENTS ABOVE
  
positive.code:       use mod_ranlux
        use mod_funlux
        use floatPrecision
  
        IMPLICIT none
        integer IS
  
negative.code:                                             !(in case of dechanneling)
        real(kind=fPrec) Am_length            !Amorphous length
        real(kind=fPrec) Length_xs, Length_ys !Amorphous length
        real(kind=fPrec) L_chan, tchan
        real(kind=fPrec) xp_rel               !xp-miscut angle in mrad
  !      REAL RNDM                           !random numbers
  !      REAL      rndm4
  !      REAL      RAN_GAUSS
        real(kind=fPrec) eUm(4)                !maximum potential
  
        real(kind=fPrec) rho(4),z(4),Ime(4) !Daniele: material parameters for dE/dX calculation
        real(kind=fPrec) k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
  