source: ZHANGProjects/ICOSIM/CPP/trunk/source/crystal_dan_FINAL_VERSION_CHvsVR.f @ 17

Last change on this file since 17 was 17, checked in by zhangj, 10 years ago

Added comments...
File size: 63.2 KB
Line 
1C.**************************************************************************
2C     SUBROUTINE FOR THE MOVEMENTS OF THE PARTICLES IN THE CRYSTAL
3C.**************************************************************************
4      SUBROUTINE CRYST(Length,layerflag,PROC,
5     +          L_chan,tchan,Alayer,ymin,ymax,Rcurv,s_length,IS,NAM, 
6     +          xpcrit0,xpcrit,Rcrit,ratio,C_orient,miscut,Cry_length, 
7     +          ZN,C_xmax,C_ymax,x,xp,y,yp,PC,s,Chann,xp_rel,Ang_rms,
8     +          Ang_avr,Vcapt,Dechan,DES,DLYI,DLRI,eUm,AI)
9             
10
11
12c      SUBROUTINE CRYST(Length,layerflag,crypara,partpara,crysparaDes,
13c     +          crysparaDlri,
14c     +          crysparaDlyi,crysparaeUmIS,crysparaAIIS,PROC)
15C      SUBROUTINE CRYST(IS,x,xp,y,yp,PC,Length,p1,p2,p3,p4,p5,p6,p7,PROC,LAYERFLAG,Chann,Dechan,L_chan)
16
17c
18C     Simple tranport protons in crystal 2
19C-----------------------------------------------------------C
20C      J -   number of element                              C
21C      S - longitudinal coordinate     
22C      IS -   number of substance 1-4: Si,W,C,Ge(110)       C
23C      x,xp,y,yp - coordinates at input of crystal          C
24C      PC -   momentum of particle*c [GeV]                  C
25C      W  -   weigth of particle                            C
26C-----------------------------------------------------------C
27C
28C
29      IMPLICIT none
30c     
31C
32C
33      integer LAYERFLAG
34      double precision p1, p2,p3,p4 !common block: Cry_length, Rcurv,C_xmax,C_ymax
35      double precision p5 !common block: Alayer
36      integer p6 !common block: C_orient
37C
38      double precision p7 !common block: miscut
39      CHARACTER*50 p8 !common block: PROC
40
41     
42     
43     
44      double precision Rcurv,Length,C_xmax,C_ymax !crystal geometrical parameters
45                                                  ! [m],[m],[m],[m],[rad]
46      double precision ymax,ymin       !crystal geometrical parameters
47      double precision s_length             !element length along s
48      double precision Alayer               !amorphous layer [m]
49      integer C_orient                      !crystal orientation
50      integer IS                            !index of the material
51c      integer counter
52      double precision  DLRI(4),DLYI(4),AI(4),DES(4)!cry parameters:see line~270
53      double precision DESt                  ! Daniele: changed energy loss by ionization now calculated and not tabulated
54      integer NAM,ZN                        !switch on/off the nuclear
55                                            !interaction (NAM) and the MCS (ZN)
56      double precision x,xp,y,yp,PC         !coordinates of the particle
57                                            ![m],[rad],[m],[rad],[GeV]
58      double precision x0,y0                !coordinates of the particle [m]
59      double precision s                    !long coordinates of the particle [m]
60      double precision a_eq,b_eq,c_eq,Delta !second order equation param.
61      double precision Ang_rms, Ang_avr     !Volume reflection mean angle [rad]
62      double precision c_v1                 !fitting coefficient
63      double precision c_v2                 !fitting coefficient
64      double precision Dechan               !probability for dechanneling
65      double precision Lrefl, Srefl         !distance of the reflection point [m]
66      double precision Vcapt                !volume capture probability
67      double precision Chann                !channeling probability
68      double precision N_atom               !probability for entering
69                                            !channeling near atomic planes
70      double precision Dxp                  !variation in angle
71      double precision xpcrit               !critical angle for curved crystal[rad]
72      double precision xpcrit0              !critical angle for str. crystal [rad]
73      double precision Rcrit                !critical curvature radius [m]
74      double precision ratio                !x=Rcurv/Rcrit
75      double precision Cry_length
76      double precision TLdech2              !tipical dechanneling length(1) [m]
77      double precision TLdech1              !tipical dechanneling length(2) [m]
78      double precision tdech, Ldech,Sdech   !angle, lenght, and S coordinate
79                                            !of dechanneling point 
80      double precision Rlength, Red_S       !reduced length/s coordinate
81                                            !(in case of dechanneling)
82      double precision Am_length            !Amorphous length
83      double precision Length_xs, Length_ys !Amorphous length
84      double precision miscut               !miscut angle in rad
85      double precision L_chan, tchan
86      double precision xp_rel               !xp-miscut angle in mrad
87c      REAL*4 RNDM                           !random numbers
88c      REAL*4      rndm4
89c      REAL*4      RAN_GAUSS
90      double precision RNDM                           !random numbers
91      double precision  rndm4
92      double precision  RAN_GAUSS
93      double precision eUm(4)                !maximum potential
94      CHARACTER*50 PROC        !string that contains the physical process 
95
96      double precision rho(4),z(4),Ime(4) !Daniele: material parameters for dE/dX calculation 
97      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
98      double precision enr,mom,betar,gammar,bgr !Daniele: energy,momentum,beta relativistic, gamma relativistic
99      double precision Tmax,plen !Daniele: maximum energy tranfer in single collision, plasma energy (see pdg)
100      double precision anuc_cry2(4),aTF,dP,const_dech,u1,xpin,ypin
101
102
103      real emr_cry(4)
104
105C     global variables shared by this file and the main file
106C     "crystaltrack_dan.f".
107C      common /Par_Cry1/ Cry_length, Rcurv,C_xmax,C_ymax,Alayer,C_orient
108C      common /miscut/ miscut
109C      common /Proc2/PROC
110
111C     global variables shared by the subroutines in this file.
112C      COMMON/NPC/     NAM,ZN                 
113C      COMMON/CRYS/    DLRI,DLYI,AI,DES       
114C      COMMON/eUc/     eUm  !
115      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
116      common/ion2/enr,mom,gammar,betar,bgr,Tmax,plen
117      common/dech/aTF,dP,u1
118
119c     global variables shared by this file and the main file
120c     "crystaltrack_dan.f"
121c      Cry_length=p1
122c      Rcurv=p2
123c      C_xmax=p3
124c      C_ymax=p4
125c     Alayer=p5
126c      C_orient=p6
127c      miscut=p7
128
129      IS = IS + 1; 
130   
131c      common/utils/ counter
132C
133c
134      NAM=1 !switch on/off the nuclear interaction (NAM) and the MCS (ZN)
135      ZN=1
136
137!-------------Daniele: dE/dX and dechanneling length calculation--------------------
138
139
140      enr=PC*1.0d3 !GeV -> MeV
141      mom=(enr*enr-mp*mp)**0.5 !p momentum [MeV/c]
142      gammar=enr/mp
143      betar=mom/enr
144      bgr=betar*gammar
145
146      Tmax=(2.0d0*me*bgr**2)/(1.0d0+2*gammar*me/mp+(me/mp)**2) ![MeV]
147
148      plen=((rho(IS)*z(IS)/anuc_cry2(IS))**0.5)*28.816d-6 ![MeV]
149
150c      DESt=((k*z(IS))/(anuc_cry2(IS)*betar**2))*
151c     + (0.5*log((2.0d0*me*bgr*bgr*Tmax)/(Ime(IS)*Ime(IS)))
152c     + -betar**2.0-log(plen/Ime(IS))-log(bgr)+0.5);
153
154c      DESt=DESt*rho(IS)*0.1 ![GeV/m]
155
156      const_dech=(256.0/(9.0*(4.D0*DATAN(1.D0))**2))*
157     + (1.0/(log(2.0*me*gammar/Ime(IS))-1.0))*((aTF*dP)/(re*me))  ![m/MeV]
158
159      const_dech=const_dech*1.0d3    ![m/GeV]
160
161!      write(*,*)DESt, const_dech
162
163!----------------------------------------------------------
164
165
166c      miscut=0.001000
167c
168c      write(*,*)"last miscut angle =",miscut
169c      write(*,*) 'enter crystal subroutine'
170c      write(*,*) 'particle energy Gev :', PC
171c      write(*,*) 'x_initial :', x
172c      write(*,*) 'Length [m]:', Length
173c      write(*,*) 'Random:', rndm4()
174c      write(*,*)'xp',xp,'x',x , 's', s
175      s=0
176      s_length=Rcurv*(sin(length/Rcurv)) !
177      L_chan=length
178
179C******************************************
180C Start to call the crystal routine
181C  the same as: SimCrys::bases() in
182C  simcrys.cc
183C
184C   Commented by Jianfeng Zhang @ LAL, 29/04/2014
185C
186C******************************************
187C     
188C-----------SimCrys::bases() in simcrys.cc
189C
190C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
191
192      if ( miscut .lt. 0 
193     &     .and. x .gt. 0 !should be useless
194     &     .and. x .lt. -length*tan(miscut)) then
195            L_chan=-x/sin(miscut)
196      endif
197      tchan=L_chan/Rcurv
198      xp_rel=xp-miscut
199C----------------------------
200
201
202C     
203C-----------Crystal::Crystal() in crystal.cc
204C
205C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
206
207c  FIRST CASE: p don't interact with crystal
208      ymin = - C_ymax/2
209      ymax =  C_ymax/2
210C--------------------------------------------
211
212C     
213C-----------SimCrys::cross() in simcrys.cc
214C
215C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
216
217      IF (y.LT.ymin .or. y.GT.ymax .or. x.gt.C_xmax) THEN
218     
219        write(*,*) "y = ", y
220        write(*,*) "ymin = ", ymin
221        write(*,*) "ymax = ", ymax
222        write(*,*) "x = ", x
223        write(*,*) "C_xmax = ", C_xmax
224       
225        x = x+xp*s_length
226        y = y+yp*s_length
227        PROC='out'
228        GOTO 111
229
230C---------------------------
231
232
233C--------------------------------------------
234
235C     
236C-----------SimCrys::layer() in simcrys.cc
237C
238C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
239
240c SECOND CASE: p hits the amorphous layer       
241      ELSEIF ( (x.LT.Alayer) .or.  ((y-ymin).LT.Alayer) .or.
242     1  ((ymax-y).lt.Alayer)  ) THEN
243        x0=x
244        y0=y
245        a_eq=(1+(xp)**2)
246        b_eq=(2*(x)*(xp)-2*(xp)*Rcurv)
247        c_eq=(x)**2-2*(x)*Rcurv
248        Delta=b_eq**2-4*a_eq*c_eq
249        s=((-b_eq+sqrt(Delta))/(2*a_eq))
250
251        if (s .ge. s_length) s=s_length
252
253        x=(xp)*s+x0
254        Length_xs=sqrt((x-x0)**2+s**2)
255
256        if ( (yp .ge.0 .and. (y+yp*s).le.ymax)) then
257          Length_ys = yp*Length_xs
258        elseif (yp.lt.0 .and. (y+yp*s).ge. ymin) then
259          Length_ys = yp*Length_xs
260        else
261          s=(ymax-y)/yp
262          Length_ys = sqrt((ymax-y)**2+s**2)
263          x=x0+xp*s
264          Length_xs=sqrt((x-x0)**2+s**2)
265        endif
266        Am_length   = sqrt(Length_xs**2+Length_ys**2)       
267        s=s/2
268        x=x0+xp*s
269        y=y0+yp*s
270        PROC='AM'
271!        CALL MOVE_AM_(IS,NAM,Am_Length,DES(IS),DLYi(IS),DLRi(IS),xp,yp
272!     + ,PC)
273        CALL CALC_ION_LOSS(IS,PC,AM_Length,DESt)
274        CALL MOVE_AM_(IS,NAM,Am_Length,DESt,DLYi(IS),DLRi(IS),xp,yp
275     + ,PC)
276        x=x+xp*(s_length-s)
277        y=y+yp*(s_length-s)
278
279        LAYERFLAG = 0;
280        GOTO 111
281      ELSEIF ((x.GT.(C_xmax-Alayer)) .and. x.LT.(C_xmax)  ) THEN
282        PROC='AM'
283!        CALL MOVE_AM_(IS,NAM,s_length,DES(IS),DLYi(IS),DLRi(IS), xp,yp
284!     + ,PC)
285        CALL CALC_ION_LOSS(IS,PC,s_length,DESt)
286        CALL MOVE_AM_(IS,NAM,s_length,DESt,DLYi(IS),DLRi(IS), xp,yp
287     + ,PC)
288        WRITE(*,*)'Fix here!'
289        LAYERFLAG = 0;
290        GOTO 111
291      END IF
292
293
294C--------------------------------------------
295
296C     
297C-----------SimCrys::parameters() in simcrys.cc
298C
299C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
300C
301c
302c THIRD CASE: the p interacts with the crystal.     
303C. Define typical angles/probabilities for orientation 110
304c
305      xpcrit0 = (2.e-9*eUm(IS)/PC)**0.5       ! critical angle (rad) for
306                                              ! straight crystals
307      Rcrit  = PC/(2.e-6*eUm(IS))*AI(IS)      ! critical curvature radius [m]
308                                              ! if R>Rcritical=>no channeling is
309                                              ! possible (ratio<1)
310      ratio = Rcurv/Rcrit                     ! parameter Rcry/Rcritical
311c      write(*,*) "Critical Radius: ",Rcrit
312c       write(*,*) "eUm(IS): ", eUm(IS)
313c       write(*,*)  "PC: ", PC
314c       write(*,*) "xpcrit0: ", xpcrit0
315      xpcrit = xpcrit0*(Rcurv-Rcrit)/Rcurv    ! critical angle for curved crystal
316c----------------valentina approx-----------   
317c      xpcrit = xpcrit0*(1-(Rcrit/Rcurv))**0.5
318c      write(*,*)(Rcurv-Rcrit)/Rcurv,(1-(Rcrit/Rcurv))**0.5 
319
320                                              ! NB: if ratio<1 => xpcrit<0
321      c_v1 = 1.7                              ! fitting coefficient ??!
322      c_v2 = -1.5                             ! fitting coefficient ???
323      if (ratio .le. 1.) then                 ! case 1:no possibile channeling
324        Ang_rms = c_v1*0.42*xpcrit0*sin(1.4*ratio)  ! rms scattering
325        Ang_avr = c_v2*xpcrit0*0.05*ratio           ! average angle reflection
326        Vcapt = 0.0                                 ! probability of VC
327        elseif (ratio .le. 3) then              ! case 2: strongly bent xstal
328          Ang_rms = c_v1*0.42*xpcrit0*sin(1.571*0.3*ratio+0.85)! rms scattering
329          Ang_avr = c_v2*xpcrit0*(0.1972*ratio-0.1472)  ! avg angle reflection
330c          Vcapt   = 0.01*(ratio-0.7)/(PC**2)
331          Vcapt   = 0.0007*(ratio-0.7)/PC**0.2 !correction by sasha drozdin/armen   
332          !K=0.00070 is taken based on simulations using CATCH.f (V.Biryukov)   
333        else                                       ! case 3: Rcry >> Rcrit
334          Ang_rms = c_v1*xpcrit0*(1./ratio)        !   
335          Ang_avr = c_v2*xpcrit0*(1.-1.6667/ratio) ! average angle for VR
336c          Vcapt = 0.01*(ratio-0.7)/(PC**2)        ! probability for VC 
337          Vcapt = 0.0007*(ratio-0.7)/PC**0.2  !correction by sasha drozdin/armen
338          ! K=0.0007 is taken based on simulations using CATCH.f (V.Biryukov)
339      endif
340cc----------------valentina approx-----------   
341c      Ang_avr=-(xpcrit+xpcrit0) 
342c      Ang_rms=(xpcrit0-xpcrit)/2 
343cc-----------end valentina approx--------------
344
345c      write(*,*) "Rcrit" , Rcrit,"Rcurv",Rcurv,
346c     c "Ratio: ",ratio,"average VR angle:", ang_avr*1e6,"+-",           
347c     c ang_rms*1e6, "ang crit:", xpcrit0*1e6,xpcrit*1e6
348c     
349      if(C_orient .eq. 2) then
350        Ang_avr = Ang_avr * 0.93                     ! for (111)
351        Ang_rms = Ang_rms * 1.05
352        xpcrit  = xpcrit * 0.98
353      endif
354
355C---------------------------
356
357
358
359
360C--------------------------------------------
361
362C     
363C-----------SimCrys::channel() in simcrys.cc
364C
365C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
366c
367C. case 3-1: channeling
368      IF (abs(xp_rel) .lt. xpcrit) THEN              ! if R' < R'c (ok CH) (1)
369        Chann  = (xpcrit**2-xp_rel**2)**0.5/xpcrit0  ! probability of CH/VC
370        N_atom = 0.1                                 ! probability of entering
371                                                     ! close to atomic planes
372        IF (rndm4() .le. Chann) then      ! if they can channel: 2 options
373                                          ! option 1:channeling
374c          TLdech1= 0.00054*PC*(1.-1./ratio)**2 ! calculate dechanneling length
375!          TLdech1= 0.0005*PC*(1.-1./ratio)**2 !calculate tipical dech. length(m)
376          TLdech1= const_dech*PC*(1.-1./ratio)**2 !Daniele: updated calculate tipical dech. length(m)
377          IF (rndm4() .le. N_atom) then
378c            TLdech1= 0.000004*PC*(1.-1./ratio)**2! calculate tipical dechanneling
379                                                  !length near atomic planes(m)
380           !next line new from sasha
381!            TLdech1= 2.0e-6*PC*(1.-1./ratio)**2  ! dechanneling length (m)
382            TLdech1= (const_dech/200.d0)*PC*(1.-1./ratio)**2  ! Daniele: updated dechanneling length (m)
383
384                               !for short crystal for high amplitude particles
385          ENDIF 
386
387c          TLdech1=TLdech1/100 !!!!CHECK
388         
389          Dechan = -log(rndm4())                 ! probability of dechanneling
390          Ldech  = TLdech1*Dechan                ! actual dechan. length
391                     ! careful: the dechanneling lentgh is along the trajectory
392                     ! of the particle -not along the longitudinal coordinate...   
393          if(Ldech .LT. L_chan) THEN   
394           PROC='DC'                  !                               
395            Dxp= Ldech/Rcurv             ! change angle from channeling [mrad]
396            Sdech=Ldech*cos(miscut+0.5*Dxp)
397
398            x  = x+ Ldech*(sin(0.5*Dxp+miscut))   ! trajectory at channeling exit
399            xp = xp + Dxp + 2.0*(rndm4()-0.5)*xpcrit
400            y= y + yp * Sdech
401            CALL CALC_ION_LOSS(IS,PC,Ldech,DESt)
402            PC = PC - 0.5*DESt*Ldech          ! energy loss to ionization while in CH [GeV]
403
404            x = x + 0.5*(s_length-Sdech)*xp
405            y = y + 0.5*(s_length-Sdech)*yp
406
407            CALL CALC_ION_LOSS(IS,PC,s_length-Sdech,DESt)
408            CALL 
409!     +MOVE_AM_(IS,NAM,s_length-Sdech,DES(IS),DLYi(IS),DLRi(IS),xp,yp,PC)
410     +MOVE_AM_(IS,NAM,s_length-Sdech,DESt,DLYi(IS),DLRi(IS),xp,yp,PC)
411           !next line new from sasha
412!            PC = PC - 0.5*DES(IS)*y          ! energy loss to ionization [GeV]
413            x = x + 0.5*(s_length-Sdech)*xp
414            y = y + 0.5*(s_length-Sdech)*yp
415          else   !channeling                 
416            PROC='CH'
417            xpin=XP
418            ypin=YP
419
420c            write(*,*) PROC
421
422            CALL MOVE_CH_(IS,NAM,L_chan,X,XP,YP,PC,Rcurv,Rcrit)  !daniele:check if a nuclear interaction happen while in CH
423
424c            write(*,*) PROC
425
426            if(PROC(1:2).ne.'CH')then             !daniele: if an nuclear interaction happened, move until the middle with initial xp,yp
427            x = x + 0.5 * L_chan * xpin           !then propagate until the "crystal exit" with the new xp,yp
428            y = y + 0.5 * L_chan * ypin           !accordingly with the rest of the code in "thin lens approx"
429            x = x + 0.5 * L_chan * XP
430            y = y + 0.5 * L_chan * YP
431            CALL CALC_ION_LOSS(IS,PC,Length,DESt)
432            PC = PC - DESt*Length       ! energy loss to ionization [GeV]
433            else
434            Dxp= L_chan/Rcurv + 0.5*RAN_GAUSS(1.)*xpcrit ! change angle[rad]
435            xp = Dxp
436            !next line new from sasha
437            x  = x+ L_chan*(sin(0.5*Dxp+miscut)) ! trajectory at channeling exit
438c            xp = xp + Dxp + 2.0*(rndm4()-0.5)*xpcrit 
439            y = y + s_length * yp
440c            !next line new from sasha
441!            PC = PC - 0.5*DES(IS)*Length       ! energy loss to ionization [GeV]
442            CALL CALC_ION_LOSS(IS,PC,Length,DESt)
443            PC = PC - 0.5*DESt*Length       ! energy loss to ionization [GeV]
444            endif
445          endif                     
446        ELSE                                   !option 2: VR
447                                               ! good for channeling
448                                               ! but don't channel         (1-2)
449          PROC='VR'                            !volume reflection at the surface
450c          Dxp=0.5*(xp_rel/xpcrit+1)*Ang_avr
451c          xp=xp+Dxp+Ang_rms*RAN_GAUSS(1.)
452            !next line new from sasha
453         
454c          write(*,*) "y = ", y
455c           write(*,*) "before enter VR: x = ", x, "  xp = ",xp
456c           write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
457           
458          xp=xp+0.45*(xp/xpcrit+1)*Ang_avr
459          x = x + 0.5*s_length * xp
460          y = y + 0.5*s_length * yp
461!          CALL MOVE_AM_(IS,NAM,s_length,DES(IS),DLYi(IS),DLRi(IS),
462!     +      xp ,yp,PC)
463          CALL CALC_ION_LOSS(IS,PC,s_length,DESt)
464          CALL MOVE_AM_(IS,NAM,s_length,DESt,DLYi(IS),DLRi(IS), 
465     +      xp ,yp,PC)
466     
467c          write(*,*) "after enter VR: x = ", x, "  xp = ",xp
468c          write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
469         
470          x = x + 0.5*s_length * xp
471          y = y + 0.5*s_length * yp
472        ENDIF       
473
474C------------------------------                             !           
475C
476C
477C--------------------------------------------
478C     
479C-----------SimCrys::reflection() in simcrys.cc
480C
481C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
482
483c case 3-2: no good for channeling. check if the  can VR
484      ELSE                                       
485        Lrefl =  (xp_rel)*Rcurv                  ! distance of refl. point [m]
486c        Srefl = sin(xp) * Lrefl
487        Srefl = sin(xp_rel/2+miscut) * Lrefl
488        if(Lrefl .gt. 0. .and. Lrefl .lt. Length) then 
489                ! VR point inside
490                !2 options: volume capture and volume reflection
491          IF (rndm4() .gt. Vcapt .or. ZN. eq. 0.) THEN   !opt. 1: VR
492           PROC='VR'
493           
494            x = x + xp * Srefl
495            y = y + yp * Srefl
496            Dxp= Ang_avr
497            xp = xp + Dxp + Ang_rms*RAN_GAUSS(1.)
498            x = x + 0.5* xp * (s_length - Srefl)
499            y = y + 0.5* yp * (s_length - Srefl)     !       
500!            CALL MOVE_AM_(IS,NAM,s_length-Srefl,DES(IS),DLYi(IS),
501!     +          DLRi(IS),xp ,yp,PC)
502
503c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
504c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
505           
506            CALL CALC_ION_LOSS(IS,PC,s_length-Srefl,DESt)
507           
508c            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
509c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
510c            write(*,*) "IS = ", IS, " NAM = ", NAM
511c            write(*,*)  "s_length-Srefl = ",s_length-Srefl
512c            write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
513c            write(*,*) "DLRi(IS) = ", DLRi(IS)
514           
515            CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS),
516     +          DLRi(IS),xp ,yp,PC)
517     
518c            write(*,*) "after enter VR: x = ", x, "  xp = ",xp
519c            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
520           
521            x = x + 0.5 * xp * (s_length - Srefl)
522            y = y + 0.5 * yp * (s_length - Srefl)   
523
524        ELSE                                      !opt 2: VC
525            x = x + xp * Srefl
526            y = y + yp * Srefl
527c            TLdech2= 0.00011*PC**0.25*(1.-1./ratio)**2 ! dechanneling length(m)
528c            Dechan = log(1.-rndm4())
529c            Ldech  = -TLdech2*Dechan
530           !next 2 lines new from sasha - different dechanneling
531           !probability
532 !           TLdech2= 0.01*PC*(1.-1./ratio)**2   ! typical dechanneling length(m)
533 !           Ldech  = 0.005*TLdech2*(sqrt(0.01-log(rndm4())) -0.1)**2 ! DC length
534            TLdech2= (const_dech/10.0d0)*PC*(1.-1./ratio)**2   ! Daniele: updated typical dechanneling length(m)
535            Ldech  = TLdech2*(sqrt(0.01-log(rndm4())) -0.1)**2 ! daniele: updated DC length
536            tdech=Ldech/Rcurv
537            Sdech=Ldech*cos(xp+0.5*tdech)
538            IF(Ldech .LT. (Length-Lrefl)) then   
539              PROC='DC'
540              Dxp= Ldech/Rcurv + 0.5*ran_gauss(1)*xpcrit
541              x  = x+ Ldech*(sin(0.5*Dxp+xp))   ! trajectory at channeling exit
542              y = y + Sdech * yp
543              xp =  Dxp
544              Red_S = s_length-Srefl -Sdech       
545              x = x + 0.5 * xp * Red_S
546              y = y + 0.5 * yp * Red_S
547              CALL CALC_ION_LOSS(IS,PC,Srefl,DESt)                     
548              PC=PC - DESt * Srefl !Daniele: "added" energy loss before capture
549              CALL CALC_ION_LOSS(IS,PC,Sdech,DESt)
550              PC=PC - 0.5 * DESt * Sdech !Daniele: "added" energy loss while captured
551!              CALL MOVE_AM_(IS,NAM,Red_S,DES(IS),DLYi(IS),DLRi(IS),
552!     +          xp,yp,PC)
553              CALL CALC_ION_LOSS(IS,PC,Red_S,DESt)
554              CALL MOVE_AM_(IS,NAM,Red_S,DESt,DLYi(IS),DLRi(IS),
555     +          xp,yp,PC)
556              x = x + 0.5 * xp * Red_S
557              y = y + 0.5 * yp * Red_S                   
558            else                               
559              PROC='VC'
560              Rlength = Length-Lrefl
561              tchan = Rlength / Rcurv
562              Red_S=Rlength*cos(xp+0.5*tchan)
563              CALL CALC_ION_LOSS(IS,PC,Lrefl,DESt)
564              PC=PC - DESt*Lrefl  !Daniele: "added" energy loss before capture   
565              xpin=XP
566              ypin=YP
567              CALL MOVE_CH_(IS,NAM,Rlength,X,XP,YP,PC,Rcurv,Rcrit)  !daniele:check if a nuclear interaction happen while in CH
568
569              if(PROC(1:2).ne.'VC')then             !daniele: if an nuclear interaction happened, move until the middle with initial xp,yp
570              x = x + 0.5 * Rlength * xpin           !then propagate until the "crystal exit" with the new xp,yp
571              y = y + 0.5 * Rlength * ypin           !accordingly with the rest of the code in "thin lens approx"
572              x = x + 0.5 * Rlength * XP
573              y = y + 0.5 * Rlength * YP
574              CALL CALC_ION_LOSS(IS,PC,Rlength,DESt)
575              PC=PC - DESt*Rlength
576              else
577              Dxp = (Length-Lrefl)/Rcurv
578              x  = x+ sin(0.5*Dxp+xp)*Rlength     ! trajectory at channeling exit
579              y = y + red_S * yp
580              xp =  Length/Rcurv + 0.5*ran_gauss(1)*xpcrit ! [mrad]
581              CALL CALC_ION_LOSS(IS,PC,Rlength,DESt)
582              PC=PC - 0.5*DESt*Rlength  !Daniele: "added" energy loss once captured   
583              endif
584            endif                       
585          ENDIF                       
586C.  case 3-3: move in amorphous substance (big input angles)---------------
587        else                         
588           PROC='AM'
589           x = x + 0.5 * s_length * xp
590           y = y + 0.5 * s_length * yp
591          if(ZN .gt. 0) then                       
592!           CALL MOVE_AM_(IS,NAM,s_length,DES(IS),DLYi(IS),DLRi(IS),
593!     +          xp,yp,PC)
594           CALL CALC_ION_LOSS(IS,PC,s_length,DESt)
595           CALL MOVE_AM_(IS,NAM,s_length,DESt,DLYi(IS),DLRi(IS), 
596     +          xp,yp,PC)
597          endif
598          x = x + 0.5 * s_length * xp
599          y = y + 0.5 * s_length * yp           
600        endif                                 
601       ENDIF                                 
602C---------------------------
603C
604C
605C     print the crystal parameters to an external file
606       open(unit=833,file='crys_para_check.dat')
607     
608c      if (counter .eq. 0) then
609111   write(833,*)'crystal parameters:\n Length:',Length, '\n Rcurv:'
610     + , Rcurv ,'\n Critical Radius:', Rcrit, 'ratio',ratio             +
611     +, '\n Critical angle for straight:',                              +
612     + xpcrit0,'\n critical angle for curved crystal:', xpcrit,' \n Leng+
613     +th:', Length, '\n xmax:', C_xmax, ' ymax:', ymax, '  C_orient: '  +
614     +, C_orient                                                        +
615     +, '\n Avg angle reflection:', Ang_avr, '\n full channeling angle: +
616     +',(Length/Rcurv) 
617c      counter=1
618c      endif
619c     
620c      WRITE(*,*)'Crystal process: ',PROC
621c      write(*,*)'xp_final :', xp
622c      WRITE(*,*)'Crystal process: ',PROC,'Chann Angle',Ch_angle/1000,   +
623c     1'Critical angle: ', xpcrit/1000
624c     2 DLRI(IS),DLYI(IS),AI(IS),DES(IS),eUm(IS),IS,ZN,NAM,C_orient     !,W
625
626      close(833)
627      END
628
629C.**************************************************************************
630C     subroutine for the calculazion of the energy loss by ionization
631C.**************************************************************************
632      SUBROUTINE CALC_ION_LOSS(IS,PC,DZ,EnLo)
633
634      IMPLICIT none
635      integer IS
636      double precision PC,DZ,EnLo
637      double precision rho(4),z(4),Ime(4) !Daniele: material parameters for dE/dX calculation 
638      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
639      double precision enr,mom,betar,gammar,bgr !Daniele: energy,momentum,beta relativistic, gamma relativistic
640      double precision Tmax,plen !Daniele: maximum energy tranfer in single collision, plasma energy (see pdg)
641      double precision anuc_cry2(4)
642      real emr_cry(4)
643      double precision thl,Tt,cs_tail,prob_tail
644      double precision ranc
645c      REAL*4 RNDM4
646      double precision RNDM4
647
648      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
649      common/ion2/enr,mom,gammar,betar,bgr,Tmax,plen
650
651
652       thl= 4.0d0*k*z(IS)*DZ*100.0d0*rho(IS)/(anuc_cry2(IS)*betar**2) ![MeV]
653
654       EnLo=((k*z(IS))/(anuc_cry2(IS)*betar**2))*
655     + (0.5*log((2.0d0*me*bgr*bgr*Tmax)/(Ime(IS)*Ime(IS)))
656     + -betar**2.0-log(plen/Ime(IS))-log(bgr)+0.5);
657
658       EnLo=EnLo*rho(IS)*0.1*DZ ![GeV]
659
660       Tt=EnLo*1000.0d0+thl  ![MeV]
661
662       cs_tail=((k*z(IS))/(anuc_cry2(IS)*betar**2))*
663     + ((0.5*((1.0d0/Tt)-(1.0d0/Tmax)))-
664     + (log(Tmax/Tt)*(betar**2)/(2.0d0*Tmax))+
665     + ((Tmax-Tt)/(4.0d0*(gammar**2)*(mp**2))))
666
667       prob_tail=cs_tail*rho(IS)*DZ*100.0d0;
668
669       ranc=dble(rndm4())
670
671       if(ranc.lt.prob_tail)then
672         EnLo=((k*z(IS))/(anuc_cry2(IS)*betar**2))*
673     +   (0.5*log((2.0d0*me*bgr*bgr*Tmax)/(Ime(IS)*Ime(IS)))
674     +   -betar**2.0-log(plen/Ime(IS))-log(bgr)+0.5+
675     +   (TMax**2)/(8.0d0*(gammar**2)*(mp**2)));
676
677         EnLo=EnLo*rho(IS)*0.1 ![GeV/m]
678
679       else
680         EnLo=EnLo/DZ  ![GeV/m]
681       endif
682     
683c      write(*,*)cs_tail,prob_tail,ranc,EnLo*DZ
684
685      RETURN
686      END
687
688
689
690C--------------------------------------------
691
692C     
693C-----------SimCrys::move_am() in simcrys.cc
694C
695C The two routines are NOT the same!!!!!!!
696C
697C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
698
699C.**************************************************************************
700C     subroutine for the movement in the amorphous
701C.**************************************************************************
702      SUBROUTINE MOVE_AM_(IS,NAM,DZ,DEI,DLY,DLr, XP,YP,PC)
703C. Moving in amorphous substance...........................
704      IMPLICIT none
705      integer IS,NAM
706      double precision DZ,DEI,DLY,DLr, XP,YP,PC
707      double precision DLAI(4),SAI(4),DES(4)
708      double precision DLRI(4),DLYI(4),AI(4)
709      double precision AM(30),QP(30),NPAI
710      double precision Dc(4),eUm(4)
711      double precision DYA,W_p
712c      REAL*4 RNDM4
713c         REAL*4      RAN_GAUSS
714      double precision RNDM4
715      double precision      RAN_GAUSS
716      CHARACTER*50 PROC              !string that contains the physical process
717      COMMON /ALAST/DLAI,SAI
718      COMMON/CRYS/ DLRI,DLYI,AI,DES
719      common /Proc2/PROC
720
721
722!------- adding variables --------
723
724
725      double precision cprob_cry(0:5,1:4)
726      double precision cs(0:5,1:4)
727      double precision csref_cry(0:5,1:4)
728      double precision freep(4)
729      double precision anuc_cry(4)
730      double precision collnt_cry(4)
731      double precision bn(4)
732      double precision bnref_cry(4)
733
734      double precision freeco_cry,ppel,ppsd,pptot,pptref_cry,pptco_cry
735      double precision pperef_cry,sdcoe_cry,pref_cry,ppeco_cry
736
737      double precision ecmsq,xln15s,bpp,xm2,bsd,t,teta,va,vb,va2,vb2
738      double precision tx,tz,r2,zlm,f
739
740      integer i,ichoix
741      double precision aran
742
743      common/scat_cry/cprob_cry,csref_cry,collnt_cry,bnref_cry,anuc_cry
744      common/scat_cry/pptref_cry,pperef_cry,sdcoe_cry,pref_cry
745      common/scat_cry/pptco_cry,ppeco_cry,freeco_cry
746
747      double precision tlcut_cry,hcut_cry(4)
748      real cgen_cry(200,4),tlow,thigh,ruth_cry
749      external ruth_cry
750
751      common/ruth_scat_cry/tlcut_cry,hcut_cry
752      common/ruth_scat_cry/cgen_cry,mcurr_cry
753
754      integer length_cry,mcurr_cry
755      real xran_cry(1)
756
757!-------daniele--------------
758! adding new variables to study the energy loss when the routine is called
759!---------------------------
760
761      double precision PC_in_dan     !daniele     
762      integer PROC_dan    !daniele
763      double precision xp_in,yp_in,kxmcs,kymcs
764       double precision  ran_x, ran_y
765
766      PC_in_dan=PC                   !daniele
767
768      xp_in=XP
769      yp_in=YP
770
771
772!---------------daniele------------
773! NEW TREATMENT OF SCATTERING ROUTINE BASED ON STANDARD SIXTRACK ROUTINE
774!----------------------------------
775
776
777!------- useful calculations for cross-section and event topology calculation --------------
778
779      ecmsq = 2 * 0.93828d0 * PC
780      xln15s=log(0.15*ecmsq)
781! pp(pn) data
782      pptot = pptref_cry *(PC / pref_cry)** pptco_cry
783      ppel = pperef_cry *(PC / pref_cry)** ppeco_cry
784      ppsd = sdcoe_cry * log(0.15d0 * ecmsq)
785      bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq))
786
787!------ distribution for Ruth. scatt.---------
788
789      tlow=tlcut_cry
790      mcurr_cry=IS
791      thigh=hcut_cry(IS)
792!      write(*,*)tlow,thigh
793      call funlxp(ruth_cry,cgen_cry(1,IS),tlow,thigh)
794
795
796
797!---------- cross-section calculation -----------
798
799!
800! freep: number of nucleons involved in single scattering
801        freep(IS) = freeco_cry * anuc_cry(IS)**(1d0/3d0)
802! compute pp and pn el+single diff contributions to cross-section
803! (both added : quasi-elastic or qel later)
804        cs(3,IS) = freep(IS) * ppel
805        cs(4,IS) = freep(IS) * ppsd
806!
807! correct TOT-CSec for energy dependence of qel
808! TOT CS is here without a Coulomb contribution
809        cs(0,IS) = csref_cry(0,IS) + freep(IS) * (pptot - pptref_cry)
810        bn(IS) = bnref_cry(IS) * cs(0,IS) / csref_cry(0,IS)               
811! also correct inel-CS
812        cs(1,IS) = csref_cry(1,IS) * cs(0,IS) / csref_cry(0,IS)
813!
814! Nuclear Elastic is TOT-inel-qel ( see definition in RPP)
815        cs(2,IS) = cs(0,IS) - cs(1,IS) - cs(3,IS) - cs(4,IS)
816        cs(5,IS) = csref_cry(5,IS)
817! Now add Coulomb
818        cs(0,IS) = cs(0,IS) + cs(5,IS)
819! Interaction length in meter
820!        xintl(ma) = 0.01d0*anuc(ma)/(fnavo * rho(ma)*cs(0,ma)*1d-24)  !don't need at the moment, take it from pdg
821
822
823! Calculate cumulative probability
824
825        do i=1,4
826          cprob_cry(i,IS)=cprob_cry(i-1,IS)+cs(i,IS)/cs(0,IS)
827        end do
828C------------------------------
829
830C--------------------------------------------
831
832C     
833C-----------SimCrys::cross() in simcrys.cc
834C   ----The physics are the same.
835C  Commented by Jianfeng Zhang @ LAL, 29/04/2014
836
837!--------- Multiple Coulomb Scattering ---------
838
839     
840c      write(*,*) "In Move_AM(): "
841
842
843      xp=xp*1000
844      yp=yp*1000
845c      write(*,*)'xp initial:', xp, 'yp initial', yp
846C. DEI - dE/dx stoping energy
847      PC    = PC - DEI*DZ    ! energy lost because of ionization process[GeV]
848
849C. DYA - rms of coloumb scattering
850      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
851c      write(*,*)'dya= ',dya
852     
853      ran_x = RAN_GAUSS(1.)
854c      write(*,*) "ran_x = ", ran_x
855     
856      ran_y = RAN_GAUSS(1.)
857c      write(*,*) "ran_y = ", ran_y
858     
859      kxmcs = DYA*ran_x
860      kymcs = DYA*ran_y
861     
862c      kxmcs = DYA*RAN_GAUSS(1.)
863c     kymcs = DYA*RAN_GAUSS(1.)
864
865      XP = xp+kxmcs
866      yp = yp+kymcs
867     
868c       write(*,*) "ran_x = ", ran_x, "ran_y = ", ran_y
869c      write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs
870c      write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP
871     
872cc     XP    = XP+DYA*RAN_GAUSS(1.)
873cc      YP    = YP+DYA*RAN_GAUSS(1.)
874
875      if(NAM .eq. 0) return                    !turn on/off nuclear interactions
876
877
878!--------- Can nuclear interaction happen? -----
879
880
881      zlm=-collnt_cry(IS)*log(dble(rndm4()))
882
883c      zlm=0.0 
884
885      if(zlm.lt.DZ) then
886
887!--------- Choose nuclear interaction --------
888
889
890      aran=dble(rndm4())
891      i=1
892  10  if ( aran.gt.cprob_cry(i,IS) ) then
893          i=i+1
894          goto 10
895      endif
896      ichoix=i
897
898
899!-------- Do the interaction ----------
900
901c      ichoix=2
902
903      if (ichoix.eq.1) then
904
905        PROC = 'absorbed'            !deep inelastic, impinging p disappeared
906
907      endif
908
909      if (ichoix.eq.2) then          ! p-n elastic
910        PROC = 'pne'
911        t = -log(dble(rndm4()))/bn(IS)
912      endif
913
914      if ( ichoix .eq. 3 ) then   !p-p elastic
915        PROC='ppe'
916        t = -log(dble(rndm4()))/bpp
917      endif
918
919      if ( ichoix .eq. 4 ) then   !single diffractive
920        PROC='diff'
921        xm2 = exp( dble(rndm4()) * xln15s )
922        PC = PC  *(1.d0 - xm2/ecmsq)
923          if ( xm2 .lt. 2.d0 ) then
924             bsd = 2.d0 * bpp
925           elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then
926                bsd = (106.d0-17.d0*xm2) *  bpp / 26.d0
927           elseif ( xm2 .gt. 5.d0 ) then
928                bsd = 7.d0 * bpp / 12.d0
929          endif
930        t = -log(dble(rndm4()))/bsd
931      endif
932
933      if ( ichoix .eq. 5 ) then
934        PROC='ruth'
935           length_cry=1
936           call funlux( cgen_cry(1,IS) ,xran_cry,length_cry)
937           t=xran_cry(1)
938      endif
939
940c      write(*,*) "ichoix = ", ichoix
941c       write(*,*)'xp final:', xp, 'yp final', yp
942       
943!---------- calculate the related kick -----------
944
945
946      if ( ichoix .eq. 4) then     
947        teta = sqrt(t)/PC_in_dan                !DIFF has changed PC!!!
948      else
949        teta = sqrt(t)/PC
950      endif
951
952
953c 100  va  =2d0*rndm4()-1d0 
954c      vb = dble(rndm4())
955c      va2 = va*va
956c      vb2 = vb*vb
957c      r2 = va2 + vb2
958c      if ( r2.gt.1.d0) go to 100
959c      tx = teta * (2.d0*va*vb) / r2
960c      tz = teta * (va2 - vb2) / r2
961
962cc 100  va  =2d0*rndm4()-1d0 
963cc      vb  =2d0*rndm4()-1d0 
964cc      va2 = va*va
965cc      vb2 = vb*vb
966cc      r2 = va2 + vb2
967cc      if ( r2.gt.1.d0) go to 100
968cc      f = SQRT(-2d0*LOG(r2)/r2)
969cc      tx = teta*f*va
970cc      tz = teta*f*vb
971
972
973      tx= teta * RAN_GAUSS(1.)
974      tz= teta * RAN_GAUSS(1.)
975
976      tx = tx * 1000
977      tz = tz * 1000
978
979!---------- change p angle --------
980
981
982c      write(*,*) "tx = ", tx, " tz = ", tz
983     
984      XP = XP + tx
985      YP = YP + tz
986
987      end if                !close the if(zlm.lt.DZ)
988
989      xp=xp/1000
990      yp=yp/1000
991
992c      write(*,*)'xp final:', xp, 'yp final', yp
993!-----------------------------------------------
994! print out the energy loss and process experienced
995!-------------------
996
997               if (PROC(1:2).eq.'AM')then
998                 PROC_dan=1
999               elseif (PROC(1:2).eq.'VR') then
1000                 PROC_dan=2
1001               elseif (PROC(1:2).eq.'CH')then
1002                 PROC_dan=3
1003               elseif (PROC(1:2).eq.'VC') then
1004                 PROC_dan=4
1005               elseif (PROC(1:3).eq.'out')then
1006                 PROC_dan=-1
1007               elseif (PROC(1:8).eq.'absorbed') then
1008                 PROC_dan=5
1009               elseif (PROC(1:2).eq.'DC')then
1010                 PROC_dan=6
1011               elseif (PROC(1:3).eq.'pne')then               
1012                 PROC_dan=7
1013               elseif (PROC(1:3).eq.'ppe')then               
1014                 PROC_dan=8
1015               elseif (PROC(1:4).eq.'diff')then 
1016                 PROC_dan=9
1017               elseif (PROC(1:4).eq.'ruth')then 
1018                 PROC_dan=10
1019
1020               endif
1021
1022c      write(889,*) PC_in_dan, PC, PROC_dan, tx,tz,
1023c     & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) 
1024
1025
1026      RETURN
1027      END
1028
1029
1030!--------- Multiple Coulomb Scattering ---------
1031
1032cc      xp=xp*1000
1033cc      yp=yp*1000
1034c      write(*,*)'xp initial:', xp, 'yp initial', yp
1035C. DEI - dE/dx stoping energy
1036cc      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1037
1038C. DYA - rms of coloumb scattering
1039cc      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1040c      write(*,*)'dya=',dya
1041
1042cc      kxmcs = DYA*RAN_GAUSS(1.)
1043cc      kymcs = DYA*RAN_GAUSS(1.)
1044
1045cc      XP = xp+kxmcs
1046cc      yp = yp+kymcs
1047 
1048c      XP    = XP+DYA*RAN_GAUSS(1.)
1049c      YP    = YP+DYA*RAN_GAUSS(1.)
1050
1051cc      xp = xp/1000
1052cc      yp = yp/1000 
1053
1054!--------END DANIELE------------------------------------
1055
1056
1057C.**************************************************************************
1058C     DANIELE subroutine for check if a nuclear interaction happen while in channeling
1059C.**************************************************************************
1060      SUBROUTINE MOVE_CH_(IS,NAM,DZ,X,XP,YP,PC,R,Rc)
1061
1062      IMPLICIT none
1063      integer IS,NAM
1064      double precision DZ,X,XP,YP,PC,R,Rc
1065      double precision DLAI(4),SAI(4),DES(4)
1066      double precision DLRI(4),DLYI(4),AI(4)
1067      double precision AM(30),QP(30),NPAI
1068      double precision Dc(4),eUm(4)
1069      double precision DYA,W_p
1070c      REAL*4 RNDM4
1071c         REAL*4      RAN_GAUSS
1072      double precision RNDM4
1073      double precision   RAN_GAUSS
1074      CHARACTER*50 PROC              !string that contains the physical process
1075      COMMON /ALAST/DLAI,SAI
1076      COMMON/CRYS/ DLRI,DLYI,AI,DES
1077      common /Proc2/PROC
1078      COMMON/eUc/eUm
1079
1080!------- adding variables --------
1081
1082
1083      double precision cprob_cry(0:5,1:4)
1084      double precision cs(0:5,1:4)
1085      double precision csref_cry(0:5,1:4)
1086      double precision freep(4)
1087      double precision anuc_cry(4)
1088      double precision collnt_cry(4)
1089      double precision bn(4)
1090      double precision bnref_cry(4)
1091
1092      double precision freeco_cry,ppel,ppsd,pptot,pptref_cry,pptco_cry
1093      double precision pperef_cry,sdcoe_cry,pref_cry,ppeco_cry
1094
1095      double precision ecmsq,xln15s,bpp,xm2,bsd,t,teta,va,vb,va2,vb2
1096      double precision tx,tz,r2,zlm,f
1097
1098      integer i,ichoix
1099      double precision aran
1100
1101      common/scat_cry/cprob_cry,csref_cry,collnt_cry,bnref_cry,anuc_cry
1102      common/scat_cry/pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1103      common/scat_cry/pptco_cry,ppeco_cry,freeco_cry
1104
1105      double precision tlcut_cry,hcut_cry(4)
1106      real cgen_cry(200,4),tlow,thigh,ruth_cry
1107      external ruth_cry
1108
1109      common/ruth_scat_cry/tlcut_cry,hcut_cry
1110      common/ruth_scat_cry/cgen_cry,mcurr_cry
1111
1112      integer length_cry,mcurr_cry
1113      real xran_cry(1)
1114
1115
1116
1117!-------daniele--------------
1118! adding new variables for calculation of average density seen
1119!---------------------------
1120
1121      integer Np
1122      double precision aTF,dP,N_am,rho_min,rho_Max,u1,avrrho
1123      double precision x_i,Ueff,pv,Et,Ec,x_min,x_Max,nuc_cl_l
1124      double precision xminU,Umin
1125
1126      common/dech/aTF,dP,u1
1127
1128      double precision rho(4),z(4),Ime(4) 
1129      double precision k,re,me,mp
1130      double precision anuc_cry2(4)
1131      real emr_cry(4)
1132
1133      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1134
1135      double precision csref_tot_rsc,csref_inel_rsc
1136
1137!-------daniele--------------
1138! adding new variables to study the energy loss when the routine is called
1139!---------------------------
1140
1141
1142      double precision PC_in_dan     !daniele     
1143      integer PROC_dan    !daniele
1144      double precision xp_in,yp_in,kxmcs,kymcs
1145
1146
1147      PC_in_dan=PC                   !daniele
1148
1149      xp_in=XP
1150      yp_in=YP
1151
1152
1153!---------------daniele------------
1154! NEW TREATMENT OF SCATTERING ROUTINE BASED ON STANDARD SIXTRACK ROUTINE
1155!----------------------------------
1156
1157
1158!------- useful calculations for cross-section and event topology calculation --------------
1159
1160      ecmsq = 2 * 0.93828d0 * PC
1161      xln15s=log(0.15*ecmsq)
1162! pp(pn) data
1163      pptot = pptref_cry *(PC / pref_cry)** pptco_cry
1164      ppel = pperef_cry *(PC / pref_cry)** ppeco_cry
1165      ppsd = sdcoe_cry * log(0.15d0 * ecmsq)
1166      bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq))
1167
1168!------ distribution for Ruth. scatt.---------
1169
1170      tlow=tlcut_cry
1171      mcurr_cry=IS
1172      thigh=hcut_cry(IS)
1173!      write(*,*)tlow,thigh
1174      call funlxp(ruth_cry,cgen_cry(1,IS),tlow,thigh)
1175
1176!---------- rescale the total and inelastic cross-section accordigly to the average density seen
1177
1178      x_i=X
1179
1180      Np=INT(x_i/dP)       !calculate in which cristalline plane the particle enters
1181      x_i=x_i-Np*dP        !rescale the incoming x at the left crystalline plane
1182      x_i=x_i-(dP/2.d0)    !rescale the incoming x in the middle of crystalline planes
1183
1184      pv=PC*1.d9*PC*1.d9/sqrt(PC*1.d9*PC*1.d9+93828.d6*93828.d6)        !calculate pv=P/E
1185
1186      Ueff=eUm(IS)*(2.d0*x_i/dP)*(2.d0*x_i/dP)+pv*x_i/R   !calculate effective potential
1187
1188      Et=pv*XP*XP/2.+Ueff                                !calculate transverse energy
1189
1190      Ec=eUm(IS)*(1.d0-Rc/R)*(1.d0-Rc/R)                 !calculate critical energy in bent crystals
1191
1192!---- to avoid negative Et
1193
1194      xminU=-dP*dP*PC*1e9/(8.d0*eUm(IS)*R)
1195      Umin=abs(eUm(IS)*(2.d0*xminU/dP)*(2.d0*xminU/dP)+pv*xminU/R)
1196      Et=Et+Umin
1197      Ec=Ec+Umin
1198
1199!-----
1200
1201c      write(*,*)xminU,Umin,Et,Ec,pv,XP,Ueff
1202
1203      x_min=-(dP/2.d0)*Rc/R-(dP/2.d0)*sqrt(Et/Ec);          !calculate min e max of the trajectory between crystalline planes
1204      x_Max=-(dP/2.d0)*Rc/R+(dP/2.d0)*sqrt(Et/Ec);
1205
1206      x_min=x_min-dP/2.d0;                                    !change ref. frame and go back with 0 on the crystalline plane on the left
1207      x_Max=x_Max-dP/2.d0;
1208
1209      N_am=rho(IS)*6.022d23*1.d6/anuc_cry2(IS)           !calculate the "normal density" in m^-3
1210     
1211      rho_Max=N_am*dP/2.d0*(erf(x_Max/sqrt(2.d0*u1*u1))        !calculate atomic density at min and max of the trajectory oscillation
1212     + -erf((dP-x_Max)/sqrt(2*u1*u1)))
1213
1214      rho_min=N_am*dP/2.d0*(erf(x_min/sqrt(2.d0*u1*u1))
1215     + -erf((dP-x_min)/sqrt(2*u1*u1)));
1216
1217      avrrho=(rho_Max-rho_min)/(x_Max-x_min)                        !"zero-approximation" of average nuclear density seen along the trajectory
1218
1219      avrrho=2.d0*avrrho/N_am
1220
1221      csref_tot_rsc=csref_cry(0,IS)*avrrho        !rescaled total ref cs
1222      csref_inel_rsc=csref_cry(1,IS)*avrrho        !rescaled inelastic ref cs
1223
1224c      write(889,*) x_i,pv,Ueff,Et,Ec,N_am,avrrho,
1225c     + csref_tot_rsc,csref_inel_rsc
1226
1227!---------- cross-section calculation -----------
1228
1229!
1230! freep: number of nucleons involved in single scattering
1231        freep(IS) = freeco_cry * anuc_cry(IS)**(1d0/3d0)
1232! compute pp and pn el+single diff contributions to cross-section
1233! (both added : quasi-elastic or qel later)
1234        cs(3,IS) = freep(IS) * ppel
1235        cs(4,IS) = freep(IS) * ppsd
1236!
1237! correct TOT-CSec for energy dependence of qel
1238! TOT CS is here without a Coulomb contribution
1239        cs(0,IS) = csref_tot_rsc + freep(IS) * (pptot - pptref_cry)
1240        bn(IS) = bnref_cry(IS) * cs(0,IS) / csref_tot_rsc               
1241! also correct inel-CS
1242        cs(1,IS) = csref_inel_rsc * cs(0,IS) / csref_tot_rsc
1243!
1244! Nuclear Elastic is TOT-inel-qel ( see definition in RPP)
1245        cs(2,IS) = cs(0,IS) - cs(1,IS) - cs(3,IS) - cs(4,IS)
1246        cs(5,IS) = csref_cry(5,IS)
1247! Now add Coulomb
1248        cs(0,IS) = cs(0,IS) + cs(5,IS)
1249! Interaction length in meter
1250!        xintl(ma) = 0.01d0*anuc(ma)/(fnavo * rho(ma)*cs(0,ma)*1d-24)  !don't need at the moment, take it from pdg
1251
1252
1253! Calculate cumulative probability
1254
1255        do i=1,4
1256          cprob_cry(i,IS)=cprob_cry(i-1,IS)+cs(i,IS)/cs(0,IS)
1257        end do
1258
1259
1260!--------- Multiple Coulomb Scattering ---------
1261
1262      xp=xp*1000
1263      yp=yp*1000
1264
1265!-------- not needed, energy loss by ionization taken into account in the main body
1266
1267
1268c      write(*,*)'xp initial:', xp, 'yp initial', yp
1269C. DEI - dE/dx stoping energy
1270!      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1271
1272C. DYA - rms of coloumb scattering
1273!      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1274c      write(*,*)'dya=',dya
1275
1276!      kxmcs = DYA*RAN_GAUSS(1.)
1277!      kymcs = DYA*RAN_GAUSS(1.)
1278
1279!      XP = xp+kxmcs
1280!      yp = yp+kymcs
1281 
1282cc     XP    = XP+DYA*RAN_GAUSS(1.)
1283cc      YP    = YP+DYA*RAN_GAUSS(1.)
1284
1285      if(NAM .eq. 0) return                    !turn on/off nuclear interactions
1286
1287
1288!--------- Can nuclear interaction happen? -----
1289
1290
1291      nuc_cl_l=collnt_cry(IS)/avrrho        !rescaled nuclear collision length
1292
1293      zlm=-nuc_cl_l*log(dble(rndm4()))
1294
1295c      zlm=0.0 
1296
1297      write(889,*) x_i,pv,Ueff,Et,Ec,N_am,avrrho,
1298     + csref_tot_rsc,csref_inel_rsc,nuc_cl_l
1299
1300      if(zlm.lt.DZ) then
1301
1302!--------- Choose nuclear interaction --------
1303
1304
1305      aran=dble(rndm4())
1306      i=1
1307  10  if ( aran.gt.cprob_cry(i,IS) ) then
1308          i=i+1
1309          goto 10
1310      endif
1311      ichoix=i
1312
1313
1314!-------- Do the interaction ----------
1315
1316c      ichoix=3
1317
1318      if (ichoix.eq.1) then
1319
1320        PROC = 'ch_absorbed'            !deep inelastic, impinging p disappeared
1321
1322      endif
1323
1324      if (ichoix.eq.2) then          ! p-n elastic
1325        PROC = 'ch_pne'
1326        t = -log(dble(rndm4()))/bn(IS)
1327      endif
1328
1329      if ( ichoix .eq. 3 ) then   !p-p elastic
1330        PROC='ch_ppe'
1331        t = -log(dble(rndm4()))/bpp
1332      endif
1333
1334      if ( ichoix .eq. 4 ) then   !single diffractive
1335        PROC='ch_diff'
1336        xm2 = exp( dble(rndm4()) * xln15s )
1337        PC = PC  *(1.d0 - xm2/ecmsq)
1338          if ( xm2 .lt. 2.d0 ) then
1339             bsd = 2.d0 * bpp
1340           elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then
1341                bsd = (106.d0-17.d0*xm2) *  bpp / 26.d0
1342           elseif ( xm2 .gt. 5.d0 ) then
1343                bsd = 7.d0 * bpp / 12.d0
1344          endif
1345        t = -log(dble(rndm4()))/bsd
1346      endif
1347
1348      if ( ichoix .eq. 5 ) then
1349        PROC='ch_ruth'
1350           length_cry=1
1351           call funlux( cgen_cry(1,IS) ,xran_cry,length_cry)
1352           t=xran_cry(1)
1353      endif
1354
1355!---------- calculate the related kick -----------
1356
1357
1358      if ( ichoix .eq. 4) then     
1359        teta = sqrt(t)/PC_in_dan                !DIFF has changed PC!!!
1360      else
1361        teta = sqrt(t)/PC
1362      endif
1363
1364
1365c 100  va  =2d0*rndm4()-1d0 
1366c      vb = dble(rndm4())
1367c      va2 = va*va
1368c      vb2 = vb*vb
1369c      r2 = va2 + vb2
1370c      if ( r2.gt.1.d0) go to 100
1371c      tx = teta * (2.d0*va*vb) / r2
1372c      tz = teta * (va2 - vb2) / r2
1373
1374cc 100  va  =2d0*rndm4()-1d0 
1375cc      vb  =2d0*rndm4()-1d0 
1376cc      va2 = va*va
1377cc      vb2 = vb*vb
1378cc      r2 = va2 + vb2
1379cc      if ( r2.gt.1.d0) go to 100
1380cc      f = SQRT(-2d0*LOG(r2)/r2)
1381cc      tx = teta*f*va
1382cc      tz = teta*f*vb
1383
1384
1385      tx= teta * RAN_GAUSS(1.)
1386      tz= teta * RAN_GAUSS(1.)
1387
1388      tx = tx * 1000
1389      tz = tz * 1000
1390
1391!---------- change p angle --------
1392
1393      XP = XP + tx
1394      YP = YP + tz
1395
1396      end if                !close the if(zlm.lt.DZ)
1397
1398      xp=xp/1000
1399      yp=yp/1000
1400
1401
1402!-----------------------------------------------
1403! print out the energy loss and process experienced
1404!-------------------
1405
1406               if (PROC(1:2).eq.'AM')then
1407                 PROC_dan=1
1408               elseif (PROC(1:2).eq.'VR') then
1409                 PROC_dan=2
1410               elseif (PROC(1:2).eq.'CH')then
1411                 PROC_dan=3
1412               elseif (PROC(1:2).eq.'VC') then
1413                 PROC_dan=4
1414               elseif (PROC(1:3).eq.'out')then
1415                 PROC_dan=-1
1416               elseif (PROC(1:11).eq.'ch_absorbed') then
1417                 PROC_dan=15
1418               elseif (PROC(1:2).eq.'DC')then
1419                 PROC_dan=6
1420               elseif (PROC(1:6).eq.'ch_pne')then               
1421                 PROC_dan=17
1422               elseif (PROC(1:6).eq.'ch_ppe')then               
1423                 PROC_dan=18
1424               elseif (PROC(1:7).eq.'ch_diff')then 
1425                 PROC_dan=19
1426               elseif (PROC(1:7).eq.'ch_ruth')then 
1427                 PROC_dan=20
1428
1429               endif
1430
1431c      write(889,*) PC_in_dan, PC, PROC_dan, tx,tz,
1432c     & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) 
1433
1434!-------- Block Data ----------
1435
1436      RETURN
1437      END
1438C     
1439      BLOCK DATA
1440      implicit none
1441      double precision DLAI(4),SAI(4)
1442      double precision DLRI(4),DLYI(4),AI(4),DES(4)
1443      double precision eUm(4)
1444      COMMON /ALAST/DLAI,SAI
1445      COMMON/CRYS/ DLRI,DLYI,AI,DES
1446      COMMON/eUc/  eUm
1447C-----4 substances: Si(110),W(110),C,Ge----------------------------
1448      DATA DLRI/0.0937,.0035,0.188,.023/    ! radiation  length(m), updated from pdg for Si
1449     +    ,DLYI/.455, .096, .400, .162/      ! nuclear length(m)
1450     +    ,AI /0.96E-7, 0.56E-7, 0.63E-7, 1.E-7/  !Si110 1/2 interplan. dist. mm
1451     +    ,DES/0.56,  3.0,  0.6, 1./         ! energy deposition in subst(GeV/m)
1452     +    ,DLAI/1.6,  0.57, 2.2, 1.0/        ! elastic length(m)
1453     +    ,SAI /42.,  140.,  42., 50./       ! elastic scat. r.m.s(mr)
1454     +    ,eUm/21.34,  21.,   21.,   21./    ! only for Si(110) potent. [eV]
1455
1456
1457!------------- daniele, more data needed---------
1458
1459      double precision cprob_cry(0:5,1:4)
1460!      double precision cs(5,4)
1461      double precision csref_cry(0:5,1:4) 
1462      double precision anuc_cry(4)
1463      double precision collnt_cry(4)
1464      double precision bnref_cry(4)
1465      double precision pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1466      double precision pptco_cry,ppeco_cry,freeco_cry
1467
1468
1469      common/scat_cry/cprob_cry,csref_cry,collnt_cry,bnref_cry,anuc_cry
1470      common/scat_cry/pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1471      common/scat_cry/pptco_cry,ppeco_cry,freeco_cry
1472
1473
1474      integer i
1475
1476      data (cprob_cry(0,i),i=1,4)/4*0.0d0/
1477      data (cprob_cry(5,i),i=1,4)/4*1.0d0/
1478
1479! pp cross-sections and parameters for energy dependence
1480      data pptref_cry,pperef_cry/0.04d0,0.007d0/
1481      data sdcoe_cry,pref_cry/0.00068d0,450.0d0/
1482      data pptco_cry,ppeco_cry,freeco_cry/0.05788d0,0.04792d0,1.618d0/
1483
1484! Atomic mass [g/mole] from pdg
1485
1486      data (anuc_cry(i),i=1,4)/28.08d0,0.0d0,0.0d0,0.0d0/     !implemented only Si
1487
1488! Total and nuclear cross-sections [barn], implemeted only Si
1489      data csref_cry(0,1),csref_cry(1,1)/0.664d0, 0.430d0/  !from pdg
1490c      data csref_cry(0,1),csref_cry(1,1)/0.762d0, 0.504d0/  !with glauber's approx (NIMB 268 (2010) 2655-2659)
1491      data csref_cry(0,2),csref_cry(1,2)/0.0d0, 0.0d0/
1492      data csref_cry(0,3),csref_cry(1,3)/0.0d0, 0.0d0/
1493      data csref_cry(0,4),csref_cry(1,4)/0.0d0, 0.0d0/
1494
1495      data csref_cry(5,1)/0.039d-2/
1496      data csref_cry(5,2)/0.0d0/
1497      data csref_cry(5,3)/0.0d0/
1498      data csref_cry(5,4)/0.0d0/
1499     
1500
1501
1502! Nuclear Collision length [m] from pdg (only for Si for the moment)
1503
1504      data (collnt_cry(i),i=1,4)/0.3016d0,0.0d0,0.0d0,0.0d0/
1505
1506! Nuclear elastic slope from Schiz et al.,PRD 21(3010)1980
1507
1508      data (bnref_cry(i),i=1,4)/123.2d0,0.0d0,0.0d0,0.0d0/      !(only for Si for the moment)
1509
1510! For calculation of dE/dX due to ionization
1511
1512      double precision rho(4),z(4),Ime(4),anuc_cry2(4) !Daniele: material parameters for dE/dX calculation 
1513      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
1514
1515      real emr_cry(4) !nuclear radius for Ruth scatt.
1516
1517      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1518
1519      data (rho(i),i=1,4)/2.33d0,0.0d0,0.0d0,0.0d0/  !material density [g/cm^3]
1520      data (z(i),i=1,4)/14.0d0,0.0d0,0.0d0,0.0d0/    !atomic number
1521      data (Ime(i),i=1,4)/173.0d-6,0.0d0,0.0d0,0.0d0/!mean exitation energy [MeV]
1522      data (anuc_cry2(i),i=1,4)/28.08d0,0.0d0,0.0d0,0.0d0/     !implemented only Si
1523      data (emr_cry(i),i=1,4)/0.306d0,0.0d0,0.0d0,0.0d0/      !nuclear radius take from R_Be*(A_Si/A_Be)^1/3, where R_Be=0.302
1524
1525
1526      data k/0.307075/ !constant in front bethe-bloch [MeV g^-1 cm^2]
1527      data re/2.818d-15/  !electron radius [m]
1528      data me/0.510998910/ !electron mass [MeV/c^2]
1529      data mp/938.272013/ !proton mass [MeV/c^2]
1530
1531! For calculation of dechanneling length
1532
1533      double precision aTF,dP,u1
1534
1535      common/dech/aTF,dP,u1
1536
1537      data aTF/0.194d-10/ !screening function [m]
1538      data dP/1.92d-10/   !distance between planes (110) [m]
1539      data u1/0.075d-10/  !thermal vibrations amplitude
1540
1541      double precision tlcut_cry,hcut_cry(4)   !param. for Ruth scatt.
1542      real cgen_cry(200,4)
1543      integer mcurr_cry
1544
1545      common/ruth_scat_cry/tlcut_cry,hcut_cry
1546      common/ruth_scat_cry/cgen_cry,mcurr_cry
1547
1548      data tlcut_cry/0.0009982d0/
1549      data (hcut_cry(i),i=1,4)/0.02d0,0.0d0,0.0d0,0.0d0/
1550
1551
1552      END
1553
1554!------------------daniele: definition of rutherford scattering formula--------!
1555
1556      function ruth_cry(t_cry)
1557      implicit none
1558      integer nmat_cry,mcurr_cry
1559      parameter(nmat_cry=4)
1560!      double precision z(4),emr_cry(4),tlcut_cry,hcut_cry(4)
1561      double precision tlcut_cry,hcut_cry(4)
1562
1563      double precision rho(4),z(4),Ime(4),anuc_cry2(4) !Daniele: material parameters for dE/dX calculation 
1564      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
1565
1566      real emr_cry(4)
1567
1568      real cgen_cry(200,4)
1569 
1570      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1571
1572!      common/ion/z,emr_cry
1573      common/ruth_scat_cry/tlcut_cry,hcut_cry
1574      common/ruth_scat_cry/cgen_cry,mcurr_cry
1575
1576      real ruth_cry,t_cry
1577      double precision cnorm,cnform
1578      parameter(cnorm=2.607d-4,cnform=0.8561d3)
1579!      parameter(cnorm=2.607d-5,cnform=0.8561d3) !daniele: corrected constant
1580!c      write(6,'('' t,exp'',2e15.8)')t,t*cnform*EMr(mcurr)**2
1581!      write(*,*)mcurr_cry,emr_cry(mcurr_cry),z(mcurr_cry),t_cry     
1582
1583       ruth_cry=cnorm*exp(-t_cry*cnform*emr_cry(mcurr_cry)**2)*
1584     + (z(mcurr_cry)/t_cry)**2
1585
1586      end
1587
1588
1589!------------------------------past treatment commented (c always comment, cc old line)----------------------------------------------------------------------------------
1590
1591
1592
1593cc      xp=xp*1000
1594cc      yp=yp*1000
1595c      write(*,*)'xp initial:', xp, 'yp initial', yp
1596C. DEI - dE/dx stoping energy
1597cc      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1598C. Coulomb scattering
1599C. DYA - rms of coloumb scattering
1600cc      DYA = (14.0/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1601c      write(*,*)'dya=',dya
1602cc      XP    = XP+DYA*RAN_GAUSS(1.)
1603cc      YP    = YP+DYA*RAN_GAUSS(1.)
1604
1605cc      if(NAM .eq. 0) return
1606C.  Elastic scattering
1607cc      IF (rndm4() .LE. DZ/DLAI(IS)) THEN
1608cc        PROC = 'mcs'
1609c         write(*,*)'case 1'
1610c        XP    = XP+SAI(IS)*RAN_GAUSS(1.)/PC
1611c        YP    = YP+SAI(IS)*RAN_GAUSS(1.)/PC
1612cc        xp    = xp+196.*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane
1613cc        yp    = yp+196.*RAN_GAUSS(1.)/PC
1614cc      ENDIF
1615C.  Diffraction interaction
1616cc      IF (rndm4() .LE. DZ/(DLY*6.143)) THEN
1617cc        PROC = 'diff'
1618c         write(*,*)'case 2'
1619
1620!######################################################
1621! daniele: comment old treatement to implement the same as standard SixTrack
1622!######################################################
1623
1624c        XP    = XP+ 257.0*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane
1625c        YP    = YP+ 257.0*RAN_GAUSS(1.)/PC
1626c        W_p = rndm4()
1627c        PC = PC -0.5*(0.3*PC)**W_p             ! m0*c = 1 GeV/c for particle
1628
1629!############## end comment -> "new" treatment of diffractive interactions #########
1630
1631c         ecmsq = 2 * 0.93828d0 * PC
1632c         xln15s=log(0.15*ecmsq) 
1633c         bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq))
1634c         xm2 = exp( dble(rndm4()) * xln15s )
1635c         PC = PC  *(1.d0 - xm2/ecmsq)
1636
1637c         if ( xm2 .lt. 2.d0 ) then
1638c              bsd = 2.d0 * bpp
1639c            elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then
1640c              bsd = (106.d0-17.d0*xm2) *  bpp / 26.d0
1641c            elseif ( xm2 .gt. 5.d0 ) then
1642c              bsd = 7.d0 * bpp / 12.d0
1643c         endif
1644 
1645c         t = -log(dble(rndm4()))/bsd
1646c         teta = sqrt(t)/PC
1647c      10 va  =2d0*rndm4()-1d0 
1648c         vb = dble(rndm4())
1649c         va2 = va*va
1650c         vb2 = vb*vb
1651c         r2 = va2 + vb2
1652c         if ( r2.gt.1.d0) go to 10
1653c         tx = teta * (2.d0*va*vb) / r2
1654c         tz = teta * (va2 - vb2) / r2
1655
1656c         XP = XP + tx
1657c         YP = YP + tz
1658
1659!#################### end of "new" treatment ###############
1660
1661cc      ENDIF
1662C.  Inelastic interaction
1663cc      IF (rndm4() .LE. DZ/DLY) THEN
1664c        PC = 0.                                         !daniele: needed for debugged energy loss in the crystal, otherwise SixTrack get stuck if PC=0 (anyway the particle is "forgot" from now)
1665cc        PROC = 'absorbed'
1666cc      ENDIF
1667c      write(*,*)'xp final:', xp, 'yp final', yp
1668cc      xp=xp/1000
1669cc      yp=yp/1000
1670
1671
1672!---------------DANIELE--------------------------------
1673! print out the energy loss and process experienced
1674!-------------------
1675
1676cc               if (PROC(1:2).eq.'AM')then
1677cc                 PROC_dan=1
1678cc               elseif (PROC(1:2).eq.'VR') then
1679cc                 PROC_dan=2
1680cc               elseif (PROC(1:2).eq.'CH')then
1681cc                 PROC_dan=3
1682cc               elseif (PROC(1:2).eq.'VC') then
1683cc                 PROC_dan=3
1684cc               elseif (PROC(1:3).eq.'out')then
1685cc                 PROC_dan=-1
1686cc              elseif (PROC(1:8).eq.'absorbed') then
1687cc                 PROC_dan=5
1688cc               elseif (PROC(1:2).eq.'DC')then
1689cc                 PROC_dan=6
1690cc               elseif (PROC(1:3).eq.'mcs')then               
1691cc                 PROC_dan=7
1692cc               elseif (PROC(1:4).eq.'diff')then 
1693cc                 PROC_dan=8
1694cc               endif
1695
1696cc      write(889,*) PC_in_dan, PC, PROC_dan 
1697
1698!--------END DANIELE------------------------------------
1699
1700cc      RETURN
1701cc      END
1702C     
1703cc      BLOCK DATA
1704cc      implicit none
1705cc      double precision DLAI(4),SAI(4)
1706cc      double precision DLRI(4),DLYI(4),AI(4),DES(4)
1707cc      double precision eUm(4)
1708cc      COMMON /ALAST/DLAI,SAI
1709cc      COMMON/CRYS/ DLRI,DLYI,AI,DES
1710cc      COMMON/eUc/  eUm
1711C-----4 substances: Si(110),W(110),C,Ge----------------------------
1712cc      DATA DLRI/0.09336,.0035,0.188,.023/    ! radiation  length(m)
1713cc     +    ,DLYI/.455, .096, .400, .162/      ! nuclear length(m)
1714cc     +    ,AI /0.96E-7, 0.56E-7, 0.63E-7, 1.E-7/  !Si110 1/2 interplan. dist. mm
1715cc     +    ,DES/0.56,  3.0,  0.6, 1./         ! energy deposition in subst(GeV/m)
1716cc     +    ,DLAI/1.6,  0.57, 2.2, 1.0/        ! elastic length(m)
1717cc     +    ,SAI /42.,  140.,  42., 50./       ! elastic scat. r.m.s(mr)
1718cc     +    ,eUm/21.34,  21.,   21.,   21./    ! only for Si(110) potent. [eV]
1719cc      END
1720
1721
1722
1723
1724
1725     
1726
1727C                        'MATH.F'   10.07.01
1728C1  RANNOR - Gauss distribution simulation procedure.
1729C2  RF RNDM- Â­Â  [0,1]        *
1730C3  RNDMST - .
1731C4  ...
1732C.**************************************************************************
1733C     subroutine for the generation of random numbers with a gaussian
1734C     distribution
1735C.**************************************************************************
1736C.**************************************************************************
1737      SUBROUTINE RANNOR(A,B) !gaussian with mean 0 and sigma 1
1738C
1739C1    Gauss distribution simulation procedure
1740C
1741      Y = RNDM(1.)
1742      IF (Y .EQ. 0.) Y = RNDM(1.)
1743      Z = RNDM(1.)
1744      X = 6.2831853*Z
1745      A1= SQRT(-2.0*ALOG(Y))
1746      A = A1*SIN(X)
1747      B = A1*COS(X)
1748C
1749      RETURN
1750      END
1751C
1752      real*4 function RNDM(rdummy)
1753      REAL*4          u,c,cd,cm, rrdummy
1754      COMMON/RANDOM/  u(97),c,cd,cm,i,j
1755C
1756      rrdummy = rdummy
1757C
1758      RNDM = u(i)-u(j)
1759      if (RNDM .LT. 0.) RNDM = RNDM+1.
1760      U(i) = RNDM
1761      i = i-1
1762      if ( i .EQ. 0 ) i = 97
1763      j = j-1
1764      if ( j .EQ. 0 ) j = 97
1765      c = c-cd
1766      if ( C .LT. 0.) c = c+cm
1767      RNDM = RNDM - c
1768      if ( RNDM .LT. 0. ) RNDM = RNDM+1.
1769C
1770      return
1771      END
1772C.**************************************************************************
1773C
1774C.**************************************************************************
1775      subroutine RNDMST(NA1,NA2,NA3,NB1)
1776      REAL*4 u,c,cd,cm
1777      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1778C
1779      ma1 = na1
1780      ma2 = na2
1781      ma3 = na3
1782      mb1 = nb1
1783      i   = 97
1784      j   = 33
1785C
1786      do ii2 = 1,97
1787        s = 0.0
1788        t = 0.5
1789        do ii1 = 1,24
1790          mat  = MOD(MOD(ma1*ma2,179)*ma3,179)
1791          ma1  = ma2
1792          ma2  = ma3
1793          ma3  = mat
1794          mb1  = MOD(53*mb1+1,169)
1795          if (MOD(MB1*MAT,64) .GE. 32 ) s = s+t
1796          t = 0.5*t
1797        end do
1798        u(ii2) = s
1799      end do
1800C
1801      c  =   362436./16777216.
1802      cd =  7654321./16777216.
1803      cm = 16777213./16777216.
1804C
1805      return
1806      END
1807C.**************************************************************************
1808C
1809C.**************************************************************************
1810      subroutine RNDMIN(uin,cin,cdin,cmin,iin,jin)
1811      REAL*4 uin(97),cin,cdin,cmin
1812      REAL*4 u,c,cd,cm
1813      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1814C
1815      do kkk = 1,97
1816        u(kkk) = uin(kkk)
1817      end do
1818C
1819      c  = cin
1820      cd = cdin
1821      cm = cmin
1822      i  = iin
1823      j  = jin
1824C
1825      return
1826      END
1827C.**************************************************************************
1828C
1829C.**************************************************************************
1830      subroutine RNDMOU(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
1831      REAL*4 uout(97),cout,cdout,cmout
1832      REAL*4 u,c,cd,cm
1833      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1834C
1835      do kkk = 1,97
1836        uout(kkk) = u(kkk)
1837      end do
1838C
1839      COUT  = C
1840      CDOUT = CD
1841      CMOUT = CM
1842      IOUT  = I
1843      JOUT  = J
1844C
1845      return
1846      END
1847C.**************************************************************************
1848C
1849C.**************************************************************************
1850      subroutine RNDMTE(IO)
1851C
1852C.    *******************************************************************
1853C.    *  SUBROUTINE RNDMTE(IO)                                         
1854C.    *******************************************************************
1855C
1856      REAL*4 uu(97)
1857      REAL*4 u(6),x(6),d(6)
1858      DATA u / 6533892.0 , 14220222.0 ,  7275067.0 ,
1859     +         6172232.0 ,  8354498.0 , 10633180.0 /
1860C
1861      call RNDMOU(UU,CC,CCD,CCM,II,JJ)
1862      call RNDMST(12,34,56,78)
1863C
1864      do ii1 = 1,20000
1865        xx = rndm4()
1866      end do
1867C
1868      sd = 0.0
1869      do ii2 = 1,6
1870        x(II2)  = 4096.*(4096.*rndm4())
1871        d(II2)  = x(II2)-u(II2)
1872        sd = sd+d(II2)
1873      end do
1874C
1875      call RNDMIN(uu,cc,ccd,ccm,ii,jj)
1876      if (IO.EQ.1 .OR. SD .NE. 0.) write(6,10) (u(I),x(I),d(I),i=1,6)
1877C
1878 10   format('  === TEST OF THE RANDOM-GENERATOR ===',/,
1879     +       '    EXPECTED VALUE    CALCULATED VALUE     DIFFERENCE',/,
1880     +       6(F17.1,F20.1,F15.3,/),
1881     +       '  === END OF TEST ;',
1882     +       '  GENERATOR HAS THE SAME STATUS AS BEFORE CALLING RNDMTE')
1883      return
1884      END
1885
1886
Note: See TracBrowser for help on using the repository browser.