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

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

(1) Modify the bug to connect to the Fortran files. (2) Fix the bug of the definition type of the random functions used in Fortran file crystal_dan_FINAL....f.
File size: 63.1 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
455           write(*,*) "before enter VR: x = ", x, "  xp = ",xp
456           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     
467          write(*,*) "after enter VR: x = ", x, "  xp = ",xp
468          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
503            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
504            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
505           
506            CALL CALC_ION_LOSS(IS,PC,s_length-Srefl,DESt)
507           
508            write(*,*) "before enter VR: x = ", x, "  xp = ",xp
509            write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC
510            write(*,*) "IS = ", IS, " NAM = ", NAM
511            write(*,*)  "s_length-Srefl = ",s_length-Srefl
512            write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS)
513            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     
518            write(*,*) "after enter VR: x = ", x, "  xp = ",xp
519            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     
620      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     
840      write(*,*) "In Move_AM(): "
841
842
843      xp=xp*1000
844      yp=yp*1000
845      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)
851      write(*,*)'dya= ',dya
852     
853      ran_x = RAN_GAUSS(1.)
854      write(*,*) "ran_x = ", ran_x
855     
856      ran_y = RAN_GAUSS(1.)
857      write(*,*) "ran_y = ", ran_y
858     
859      kxmcs = DYA*ran_x
860      kxmcs = DYA*ran_x
861     
862c      kxmcs = DYA*RAN_GAUSS(1.)
863c      kymcs = DYA*RAN_GAUSS(1.)
864
865      XP = xp+kxmcs
866      yp = yp+kymcs
867     
868       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
940!---------- calculate the related kick -----------
941
942
943      if ( ichoix .eq. 4) then     
944        teta = sqrt(t)/PC_in_dan                !DIFF has changed PC!!!
945      else
946        teta = sqrt(t)/PC
947      endif
948
949
950c 100  va  =2d0*rndm4()-1d0 
951c      vb = dble(rndm4())
952c      va2 = va*va
953c      vb2 = vb*vb
954c      r2 = va2 + vb2
955c      if ( r2.gt.1.d0) go to 100
956c      tx = teta * (2.d0*va*vb) / r2
957c      tz = teta * (va2 - vb2) / r2
958
959cc 100  va  =2d0*rndm4()-1d0 
960cc      vb  =2d0*rndm4()-1d0 
961cc      va2 = va*va
962cc      vb2 = vb*vb
963cc      r2 = va2 + vb2
964cc      if ( r2.gt.1.d0) go to 100
965cc      f = SQRT(-2d0*LOG(r2)/r2)
966cc      tx = teta*f*va
967cc      tz = teta*f*vb
968
969
970      tx= teta * RAN_GAUSS(1.)
971      tz= teta * RAN_GAUSS(1.)
972
973      tx = tx * 1000
974      tz = tz * 1000
975
976!---------- change p angle --------
977
978
979      write(*,*) "tx = ", tx, " tz = ", tz
980     
981      XP = XP + tx
982      YP = YP + tz
983
984      end if                !close the if(zlm.lt.DZ)
985
986      xp=xp/1000
987      yp=yp/1000
988
989      write(*,*)'xp final:', xp, 'yp final', yp
990!-----------------------------------------------
991! print out the energy loss and process experienced
992!-------------------
993
994               if (PROC(1:2).eq.'AM')then
995                 PROC_dan=1
996               elseif (PROC(1:2).eq.'VR') then
997                 PROC_dan=2
998               elseif (PROC(1:2).eq.'CH')then
999                 PROC_dan=3
1000               elseif (PROC(1:2).eq.'VC') then
1001                 PROC_dan=4
1002               elseif (PROC(1:3).eq.'out')then
1003                 PROC_dan=-1
1004               elseif (PROC(1:8).eq.'absorbed') then
1005                 PROC_dan=5
1006               elseif (PROC(1:2).eq.'DC')then
1007                 PROC_dan=6
1008               elseif (PROC(1:3).eq.'pne')then               
1009                 PROC_dan=7
1010               elseif (PROC(1:3).eq.'ppe')then               
1011                 PROC_dan=8
1012               elseif (PROC(1:4).eq.'diff')then 
1013                 PROC_dan=9
1014               elseif (PROC(1:4).eq.'ruth')then 
1015                 PROC_dan=10
1016
1017               endif
1018
1019c      write(889,*) PC_in_dan, PC, PROC_dan, tx,tz,
1020c     & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) 
1021
1022
1023      RETURN
1024      END
1025
1026
1027!--------- Multiple Coulomb Scattering ---------
1028
1029cc      xp=xp*1000
1030cc      yp=yp*1000
1031c      write(*,*)'xp initial:', xp, 'yp initial', yp
1032C. DEI - dE/dx stoping energy
1033cc      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1034
1035C. DYA - rms of coloumb scattering
1036cc      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1037c      write(*,*)'dya=',dya
1038
1039cc      kxmcs = DYA*RAN_GAUSS(1.)
1040cc      kymcs = DYA*RAN_GAUSS(1.)
1041
1042cc      XP = xp+kxmcs
1043cc      yp = yp+kymcs
1044 
1045c      XP    = XP+DYA*RAN_GAUSS(1.)
1046c      YP    = YP+DYA*RAN_GAUSS(1.)
1047
1048cc      xp = xp/1000
1049cc      yp = yp/1000 
1050
1051!--------END DANIELE------------------------------------
1052
1053
1054C.**************************************************************************
1055C     DANIELE subroutine for check if a nuclear interaction happen while in channeling
1056C.**************************************************************************
1057      SUBROUTINE MOVE_CH_(IS,NAM,DZ,X,XP,YP,PC,R,Rc)
1058
1059      IMPLICIT none
1060      integer IS,NAM
1061      double precision DZ,X,XP,YP,PC,R,Rc
1062      double precision DLAI(4),SAI(4),DES(4)
1063      double precision DLRI(4),DLYI(4),AI(4)
1064      double precision AM(30),QP(30),NPAI
1065      double precision Dc(4),eUm(4)
1066      double precision DYA,W_p
1067c      REAL*4 RNDM4
1068c         REAL*4      RAN_GAUSS
1069      double precision RNDM4
1070      double precision   RAN_GAUSS
1071      CHARACTER*50 PROC              !string that contains the physical process
1072      COMMON /ALAST/DLAI,SAI
1073      COMMON/CRYS/ DLRI,DLYI,AI,DES
1074      common /Proc2/PROC
1075      COMMON/eUc/eUm
1076
1077!------- adding variables --------
1078
1079
1080      double precision cprob_cry(0:5,1:4)
1081      double precision cs(0:5,1:4)
1082      double precision csref_cry(0:5,1:4)
1083      double precision freep(4)
1084      double precision anuc_cry(4)
1085      double precision collnt_cry(4)
1086      double precision bn(4)
1087      double precision bnref_cry(4)
1088
1089      double precision freeco_cry,ppel,ppsd,pptot,pptref_cry,pptco_cry
1090      double precision pperef_cry,sdcoe_cry,pref_cry,ppeco_cry
1091
1092      double precision ecmsq,xln15s,bpp,xm2,bsd,t,teta,va,vb,va2,vb2
1093      double precision tx,tz,r2,zlm,f
1094
1095      integer i,ichoix
1096      double precision aran
1097
1098      common/scat_cry/cprob_cry,csref_cry,collnt_cry,bnref_cry,anuc_cry
1099      common/scat_cry/pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1100      common/scat_cry/pptco_cry,ppeco_cry,freeco_cry
1101
1102      double precision tlcut_cry,hcut_cry(4)
1103      real cgen_cry(200,4),tlow,thigh,ruth_cry
1104      external ruth_cry
1105
1106      common/ruth_scat_cry/tlcut_cry,hcut_cry
1107      common/ruth_scat_cry/cgen_cry,mcurr_cry
1108
1109      integer length_cry,mcurr_cry
1110      real xran_cry(1)
1111
1112
1113
1114!-------daniele--------------
1115! adding new variables for calculation of average density seen
1116!---------------------------
1117
1118      integer Np
1119      double precision aTF,dP,N_am,rho_min,rho_Max,u1,avrrho
1120      double precision x_i,Ueff,pv,Et,Ec,x_min,x_Max,nuc_cl_l
1121      double precision xminU,Umin
1122
1123      common/dech/aTF,dP,u1
1124
1125      double precision rho(4),z(4),Ime(4) 
1126      double precision k,re,me,mp
1127      double precision anuc_cry2(4)
1128      real emr_cry(4)
1129
1130      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1131
1132      double precision csref_tot_rsc,csref_inel_rsc
1133
1134!-------daniele--------------
1135! adding new variables to study the energy loss when the routine is called
1136!---------------------------
1137
1138
1139      double precision PC_in_dan     !daniele     
1140      integer PROC_dan    !daniele
1141      double precision xp_in,yp_in,kxmcs,kymcs
1142
1143
1144      PC_in_dan=PC                   !daniele
1145
1146      xp_in=XP
1147      yp_in=YP
1148
1149
1150!---------------daniele------------
1151! NEW TREATMENT OF SCATTERING ROUTINE BASED ON STANDARD SIXTRACK ROUTINE
1152!----------------------------------
1153
1154
1155!------- useful calculations for cross-section and event topology calculation --------------
1156
1157      ecmsq = 2 * 0.93828d0 * PC
1158      xln15s=log(0.15*ecmsq)
1159! pp(pn) data
1160      pptot = pptref_cry *(PC / pref_cry)** pptco_cry
1161      ppel = pperef_cry *(PC / pref_cry)** ppeco_cry
1162      ppsd = sdcoe_cry * log(0.15d0 * ecmsq)
1163      bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq))
1164
1165!------ distribution for Ruth. scatt.---------
1166
1167      tlow=tlcut_cry
1168      mcurr_cry=IS
1169      thigh=hcut_cry(IS)
1170!      write(*,*)tlow,thigh
1171      call funlxp(ruth_cry,cgen_cry(1,IS),tlow,thigh)
1172
1173!---------- rescale the total and inelastic cross-section accordigly to the average density seen
1174
1175      x_i=X
1176
1177      Np=INT(x_i/dP)       !calculate in which cristalline plane the particle enters
1178      x_i=x_i-Np*dP        !rescale the incoming x at the left crystalline plane
1179      x_i=x_i-(dP/2.d0)    !rescale the incoming x in the middle of crystalline planes
1180
1181      pv=PC*1.d9*PC*1.d9/sqrt(PC*1.d9*PC*1.d9+93828.d6*93828.d6)        !calculate pv=P/E
1182
1183      Ueff=eUm(IS)*(2.d0*x_i/dP)*(2.d0*x_i/dP)+pv*x_i/R   !calculate effective potential
1184
1185      Et=pv*XP*XP/2.+Ueff                                !calculate transverse energy
1186
1187      Ec=eUm(IS)*(1.d0-Rc/R)*(1.d0-Rc/R)                 !calculate critical energy in bent crystals
1188
1189!---- to avoid negative Et
1190
1191      xminU=-dP*dP*PC*1e9/(8.d0*eUm(IS)*R)
1192      Umin=abs(eUm(IS)*(2.d0*xminU/dP)*(2.d0*xminU/dP)+pv*xminU/R)
1193      Et=Et+Umin
1194      Ec=Ec+Umin
1195
1196!-----
1197
1198c      write(*,*)xminU,Umin,Et,Ec,pv,XP,Ueff
1199
1200      x_min=-(dP/2.d0)*Rc/R-(dP/2.d0)*sqrt(Et/Ec);          !calculate min e max of the trajectory between crystalline planes
1201      x_Max=-(dP/2.d0)*Rc/R+(dP/2.d0)*sqrt(Et/Ec);
1202
1203      x_min=x_min-dP/2.d0;                                    !change ref. frame and go back with 0 on the crystalline plane on the left
1204      x_Max=x_Max-dP/2.d0;
1205
1206      N_am=rho(IS)*6.022d23*1.d6/anuc_cry2(IS)           !calculate the "normal density" in m^-3
1207     
1208      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
1209     + -erf((dP-x_Max)/sqrt(2*u1*u1)))
1210
1211      rho_min=N_am*dP/2.d0*(erf(x_min/sqrt(2.d0*u1*u1))
1212     + -erf((dP-x_min)/sqrt(2*u1*u1)));
1213
1214      avrrho=(rho_Max-rho_min)/(x_Max-x_min)                        !"zero-approximation" of average nuclear density seen along the trajectory
1215
1216      avrrho=2.d0*avrrho/N_am
1217
1218      csref_tot_rsc=csref_cry(0,IS)*avrrho        !rescaled total ref cs
1219      csref_inel_rsc=csref_cry(1,IS)*avrrho        !rescaled inelastic ref cs
1220
1221c      write(889,*) x_i,pv,Ueff,Et,Ec,N_am,avrrho,
1222c     + csref_tot_rsc,csref_inel_rsc
1223
1224!---------- cross-section calculation -----------
1225
1226!
1227! freep: number of nucleons involved in single scattering
1228        freep(IS) = freeco_cry * anuc_cry(IS)**(1d0/3d0)
1229! compute pp and pn el+single diff contributions to cross-section
1230! (both added : quasi-elastic or qel later)
1231        cs(3,IS) = freep(IS) * ppel
1232        cs(4,IS) = freep(IS) * ppsd
1233!
1234! correct TOT-CSec for energy dependence of qel
1235! TOT CS is here without a Coulomb contribution
1236        cs(0,IS) = csref_tot_rsc + freep(IS) * (pptot - pptref_cry)
1237        bn(IS) = bnref_cry(IS) * cs(0,IS) / csref_tot_rsc               
1238! also correct inel-CS
1239        cs(1,IS) = csref_inel_rsc * cs(0,IS) / csref_tot_rsc
1240!
1241! Nuclear Elastic is TOT-inel-qel ( see definition in RPP)
1242        cs(2,IS) = cs(0,IS) - cs(1,IS) - cs(3,IS) - cs(4,IS)
1243        cs(5,IS) = csref_cry(5,IS)
1244! Now add Coulomb
1245        cs(0,IS) = cs(0,IS) + cs(5,IS)
1246! Interaction length in meter
1247!        xintl(ma) = 0.01d0*anuc(ma)/(fnavo * rho(ma)*cs(0,ma)*1d-24)  !don't need at the moment, take it from pdg
1248
1249
1250! Calculate cumulative probability
1251
1252        do i=1,4
1253          cprob_cry(i,IS)=cprob_cry(i-1,IS)+cs(i,IS)/cs(0,IS)
1254        end do
1255
1256
1257!--------- Multiple Coulomb Scattering ---------
1258
1259      xp=xp*1000
1260      yp=yp*1000
1261
1262!-------- not needed, energy loss by ionization taken into account in the main body
1263
1264
1265c      write(*,*)'xp initial:', xp, 'yp initial', yp
1266C. DEI - dE/dx stoping energy
1267!      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1268
1269C. DYA - rms of coloumb scattering
1270!      DYA = (13.6/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1271c      write(*,*)'dya=',dya
1272
1273!      kxmcs = DYA*RAN_GAUSS(1.)
1274!      kymcs = DYA*RAN_GAUSS(1.)
1275
1276!      XP = xp+kxmcs
1277!      yp = yp+kymcs
1278 
1279cc     XP    = XP+DYA*RAN_GAUSS(1.)
1280cc      YP    = YP+DYA*RAN_GAUSS(1.)
1281
1282      if(NAM .eq. 0) return                    !turn on/off nuclear interactions
1283
1284
1285!--------- Can nuclear interaction happen? -----
1286
1287
1288      nuc_cl_l=collnt_cry(IS)/avrrho        !rescaled nuclear collision length
1289
1290      zlm=-nuc_cl_l*log(dble(rndm4()))
1291
1292c      zlm=0.0 
1293
1294      write(889,*) x_i,pv,Ueff,Et,Ec,N_am,avrrho,
1295     + csref_tot_rsc,csref_inel_rsc,nuc_cl_l
1296
1297      if(zlm.lt.DZ) then
1298
1299!--------- Choose nuclear interaction --------
1300
1301
1302      aran=dble(rndm4())
1303      i=1
1304  10  if ( aran.gt.cprob_cry(i,IS) ) then
1305          i=i+1
1306          goto 10
1307      endif
1308      ichoix=i
1309
1310
1311!-------- Do the interaction ----------
1312
1313c      ichoix=3
1314
1315      if (ichoix.eq.1) then
1316
1317        PROC = 'ch_absorbed'            !deep inelastic, impinging p disappeared
1318
1319      endif
1320
1321      if (ichoix.eq.2) then          ! p-n elastic
1322        PROC = 'ch_pne'
1323        t = -log(dble(rndm4()))/bn(IS)
1324      endif
1325
1326      if ( ichoix .eq. 3 ) then   !p-p elastic
1327        PROC='ch_ppe'
1328        t = -log(dble(rndm4()))/bpp
1329      endif
1330
1331      if ( ichoix .eq. 4 ) then   !single diffractive
1332        PROC='ch_diff'
1333        xm2 = exp( dble(rndm4()) * xln15s )
1334        PC = PC  *(1.d0 - xm2/ecmsq)
1335          if ( xm2 .lt. 2.d0 ) then
1336             bsd = 2.d0 * bpp
1337           elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then
1338                bsd = (106.d0-17.d0*xm2) *  bpp / 26.d0
1339           elseif ( xm2 .gt. 5.d0 ) then
1340                bsd = 7.d0 * bpp / 12.d0
1341          endif
1342        t = -log(dble(rndm4()))/bsd
1343      endif
1344
1345      if ( ichoix .eq. 5 ) then
1346        PROC='ch_ruth'
1347           length_cry=1
1348           call funlux( cgen_cry(1,IS) ,xran_cry,length_cry)
1349           t=xran_cry(1)
1350      endif
1351
1352!---------- calculate the related kick -----------
1353
1354
1355      if ( ichoix .eq. 4) then     
1356        teta = sqrt(t)/PC_in_dan                !DIFF has changed PC!!!
1357      else
1358        teta = sqrt(t)/PC
1359      endif
1360
1361
1362c 100  va  =2d0*rndm4()-1d0 
1363c      vb = dble(rndm4())
1364c      va2 = va*va
1365c      vb2 = vb*vb
1366c      r2 = va2 + vb2
1367c      if ( r2.gt.1.d0) go to 100
1368c      tx = teta * (2.d0*va*vb) / r2
1369c      tz = teta * (va2 - vb2) / r2
1370
1371cc 100  va  =2d0*rndm4()-1d0 
1372cc      vb  =2d0*rndm4()-1d0 
1373cc      va2 = va*va
1374cc      vb2 = vb*vb
1375cc      r2 = va2 + vb2
1376cc      if ( r2.gt.1.d0) go to 100
1377cc      f = SQRT(-2d0*LOG(r2)/r2)
1378cc      tx = teta*f*va
1379cc      tz = teta*f*vb
1380
1381
1382      tx= teta * RAN_GAUSS(1.)
1383      tz= teta * RAN_GAUSS(1.)
1384
1385      tx = tx * 1000
1386      tz = tz * 1000
1387
1388!---------- change p angle --------
1389
1390      XP = XP + tx
1391      YP = YP + tz
1392
1393      end if                !close the if(zlm.lt.DZ)
1394
1395      xp=xp/1000
1396      yp=yp/1000
1397
1398
1399!-----------------------------------------------
1400! print out the energy loss and process experienced
1401!-------------------
1402
1403               if (PROC(1:2).eq.'AM')then
1404                 PROC_dan=1
1405               elseif (PROC(1:2).eq.'VR') then
1406                 PROC_dan=2
1407               elseif (PROC(1:2).eq.'CH')then
1408                 PROC_dan=3
1409               elseif (PROC(1:2).eq.'VC') then
1410                 PROC_dan=4
1411               elseif (PROC(1:3).eq.'out')then
1412                 PROC_dan=-1
1413               elseif (PROC(1:11).eq.'ch_absorbed') then
1414                 PROC_dan=15
1415               elseif (PROC(1:2).eq.'DC')then
1416                 PROC_dan=6
1417               elseif (PROC(1:6).eq.'ch_pne')then               
1418                 PROC_dan=17
1419               elseif (PROC(1:6).eq.'ch_ppe')then               
1420                 PROC_dan=18
1421               elseif (PROC(1:7).eq.'ch_diff')then 
1422                 PROC_dan=19
1423               elseif (PROC(1:7).eq.'ch_ruth')then 
1424                 PROC_dan=20
1425
1426               endif
1427
1428c      write(889,*) PC_in_dan, PC, PROC_dan, tx,tz,
1429c     & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) 
1430
1431!-------- Block Data ----------
1432
1433      RETURN
1434      END
1435C     
1436      BLOCK DATA
1437      implicit none
1438      double precision DLAI(4),SAI(4)
1439      double precision DLRI(4),DLYI(4),AI(4),DES(4)
1440      double precision eUm(4)
1441      COMMON /ALAST/DLAI,SAI
1442      COMMON/CRYS/ DLRI,DLYI,AI,DES
1443      COMMON/eUc/  eUm
1444C-----4 substances: Si(110),W(110),C,Ge----------------------------
1445      DATA DLRI/0.0937,.0035,0.188,.023/    ! radiation  length(m), updated from pdg for Si
1446     +    ,DLYI/.455, .096, .400, .162/      ! nuclear length(m)
1447     +    ,AI /0.96E-7, 0.56E-7, 0.63E-7, 1.E-7/  !Si110 1/2 interplan. dist. mm
1448     +    ,DES/0.56,  3.0,  0.6, 1./         ! energy deposition in subst(GeV/m)
1449     +    ,DLAI/1.6,  0.57, 2.2, 1.0/        ! elastic length(m)
1450     +    ,SAI /42.,  140.,  42., 50./       ! elastic scat. r.m.s(mr)
1451     +    ,eUm/21.34,  21.,   21.,   21./    ! only for Si(110) potent. [eV]
1452
1453
1454!------------- daniele, more data needed---------
1455
1456      double precision cprob_cry(0:5,1:4)
1457!      double precision cs(5,4)
1458      double precision csref_cry(0:5,1:4) 
1459      double precision anuc_cry(4)
1460      double precision collnt_cry(4)
1461      double precision bnref_cry(4)
1462      double precision pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1463      double precision pptco_cry,ppeco_cry,freeco_cry
1464
1465
1466      common/scat_cry/cprob_cry,csref_cry,collnt_cry,bnref_cry,anuc_cry
1467      common/scat_cry/pptref_cry,pperef_cry,sdcoe_cry,pref_cry
1468      common/scat_cry/pptco_cry,ppeco_cry,freeco_cry
1469
1470
1471      integer i
1472
1473      data (cprob_cry(0,i),i=1,4)/4*0.0d0/
1474      data (cprob_cry(5,i),i=1,4)/4*1.0d0/
1475
1476! pp cross-sections and parameters for energy dependence
1477      data pptref_cry,pperef_cry/0.04d0,0.007d0/
1478      data sdcoe_cry,pref_cry/0.00068d0,450.0d0/
1479      data pptco_cry,ppeco_cry,freeco_cry/0.05788d0,0.04792d0,1.618d0/
1480
1481! Atomic mass [g/mole] from pdg
1482
1483      data (anuc_cry(i),i=1,4)/28.08d0,0.0d0,0.0d0,0.0d0/     !implemented only Si
1484
1485! Total and nuclear cross-sections [barn], implemeted only Si
1486      data csref_cry(0,1),csref_cry(1,1)/0.664d0, 0.430d0/  !from pdg
1487c      data csref_cry(0,1),csref_cry(1,1)/0.762d0, 0.504d0/  !with glauber's approx (NIMB 268 (2010) 2655-2659)
1488      data csref_cry(0,2),csref_cry(1,2)/0.0d0, 0.0d0/
1489      data csref_cry(0,3),csref_cry(1,3)/0.0d0, 0.0d0/
1490      data csref_cry(0,4),csref_cry(1,4)/0.0d0, 0.0d0/
1491
1492      data csref_cry(5,1)/0.039d-2/
1493      data csref_cry(5,2)/0.0d0/
1494      data csref_cry(5,3)/0.0d0/
1495      data csref_cry(5,4)/0.0d0/
1496     
1497
1498
1499! Nuclear Collision length [m] from pdg (only for Si for the moment)
1500
1501      data (collnt_cry(i),i=1,4)/0.3016d0,0.0d0,0.0d0,0.0d0/
1502
1503! Nuclear elastic slope from Schiz et al.,PRD 21(3010)1980
1504
1505      data (bnref_cry(i),i=1,4)/123.2d0,0.0d0,0.0d0,0.0d0/      !(only for Si for the moment)
1506
1507! For calculation of dE/dX due to ionization
1508
1509      double precision rho(4),z(4),Ime(4),anuc_cry2(4) !Daniele: material parameters for dE/dX calculation 
1510      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
1511
1512      real emr_cry(4) !nuclear radius for Ruth scatt.
1513
1514      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1515
1516      data (rho(i),i=1,4)/2.33d0,0.0d0,0.0d0,0.0d0/  !material density [g/cm^3]
1517      data (z(i),i=1,4)/14.0d0,0.0d0,0.0d0,0.0d0/    !atomic number
1518      data (Ime(i),i=1,4)/173.0d-6,0.0d0,0.0d0,0.0d0/!mean exitation energy [MeV]
1519      data (anuc_cry2(i),i=1,4)/28.08d0,0.0d0,0.0d0,0.0d0/     !implemented only Si
1520      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
1521
1522
1523      data k/0.307075/ !constant in front bethe-bloch [MeV g^-1 cm^2]
1524      data re/2.818d-15/  !electron radius [m]
1525      data me/0.510998910/ !electron mass [MeV/c^2]
1526      data mp/938.272013/ !proton mass [MeV/c^2]
1527
1528! For calculation of dechanneling length
1529
1530      double precision aTF,dP,u1
1531
1532      common/dech/aTF,dP,u1
1533
1534      data aTF/0.194d-10/ !screening function [m]
1535      data dP/1.92d-10/   !distance between planes (110) [m]
1536      data u1/0.075d-10/  !thermal vibrations amplitude
1537
1538      double precision tlcut_cry,hcut_cry(4)   !param. for Ruth scatt.
1539      real cgen_cry(200,4)
1540      integer mcurr_cry
1541
1542      common/ruth_scat_cry/tlcut_cry,hcut_cry
1543      common/ruth_scat_cry/cgen_cry,mcurr_cry
1544
1545      data tlcut_cry/0.0009982d0/
1546      data (hcut_cry(i),i=1,4)/0.02d0,0.0d0,0.0d0,0.0d0/
1547
1548
1549      END
1550
1551!------------------daniele: definition of rutherford scattering formula--------!
1552
1553      function ruth_cry(t_cry)
1554      implicit none
1555      integer nmat_cry,mcurr_cry
1556      parameter(nmat_cry=4)
1557!      double precision z(4),emr_cry(4),tlcut_cry,hcut_cry(4)
1558      double precision tlcut_cry,hcut_cry(4)
1559
1560      double precision rho(4),z(4),Ime(4),anuc_cry2(4) !Daniele: material parameters for dE/dX calculation 
1561      double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
1562
1563      real emr_cry(4)
1564
1565      real cgen_cry(200,4)
1566 
1567      common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
1568
1569!      common/ion/z,emr_cry
1570      common/ruth_scat_cry/tlcut_cry,hcut_cry
1571      common/ruth_scat_cry/cgen_cry,mcurr_cry
1572
1573      real ruth_cry,t_cry
1574      double precision cnorm,cnform
1575      parameter(cnorm=2.607d-4,cnform=0.8561d3)
1576!      parameter(cnorm=2.607d-5,cnform=0.8561d3) !daniele: corrected constant
1577!c      write(6,'('' t,exp'',2e15.8)')t,t*cnform*EMr(mcurr)**2
1578!      write(*,*)mcurr_cry,emr_cry(mcurr_cry),z(mcurr_cry),t_cry     
1579
1580       ruth_cry=cnorm*exp(-t_cry*cnform*emr_cry(mcurr_cry)**2)*
1581     + (z(mcurr_cry)/t_cry)**2
1582
1583      end
1584
1585
1586!------------------------------past treatment commented (c always comment, cc old line)----------------------------------------------------------------------------------
1587
1588
1589
1590cc      xp=xp*1000
1591cc      yp=yp*1000
1592c      write(*,*)'xp initial:', xp, 'yp initial', yp
1593C. DEI - dE/dx stoping energy
1594cc      PC    = PC - DEI*DZ    ! energy lost beacause of ionization process[GeV]
1595C. Coulomb scattering
1596C. DYA - rms of coloumb scattering
1597cc      DYA = (14.0/PC)*SQRT(DZ/DLr)             !MCS (mrad)
1598c      write(*,*)'dya=',dya
1599cc      XP    = XP+DYA*RAN_GAUSS(1.)
1600cc      YP    = YP+DYA*RAN_GAUSS(1.)
1601
1602cc      if(NAM .eq. 0) return
1603C.  Elastic scattering
1604cc      IF (rndm4() .LE. DZ/DLAI(IS)) THEN
1605cc        PROC = 'mcs'
1606c         write(*,*)'case 1'
1607c        XP    = XP+SAI(IS)*RAN_GAUSS(1.)/PC
1608c        YP    = YP+SAI(IS)*RAN_GAUSS(1.)/PC
1609cc        xp    = xp+196.*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane
1610cc        yp    = yp+196.*RAN_GAUSS(1.)/PC
1611cc      ENDIF
1612C.  Diffraction interaction
1613cc      IF (rndm4() .LE. DZ/(DLY*6.143)) THEN
1614cc        PROC = 'diff'
1615c         write(*,*)'case 2'
1616
1617!######################################################
1618! daniele: comment old treatement to implement the same as standard SixTrack
1619!######################################################
1620
1621c        XP    = XP+ 257.0*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane
1622c        YP    = YP+ 257.0*RAN_GAUSS(1.)/PC
1623c        W_p = rndm4()
1624c        PC = PC -0.5*(0.3*PC)**W_p             ! m0*c = 1 GeV/c for particle
1625
1626!############## end comment -> "new" treatment of diffractive interactions #########
1627
1628c         ecmsq = 2 * 0.93828d0 * PC
1629c         xln15s=log(0.15*ecmsq) 
1630c         bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq))
1631c         xm2 = exp( dble(rndm4()) * xln15s )
1632c         PC = PC  *(1.d0 - xm2/ecmsq)
1633
1634c         if ( xm2 .lt. 2.d0 ) then
1635c              bsd = 2.d0 * bpp
1636c            elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then
1637c              bsd = (106.d0-17.d0*xm2) *  bpp / 26.d0
1638c            elseif ( xm2 .gt. 5.d0 ) then
1639c              bsd = 7.d0 * bpp / 12.d0
1640c         endif
1641 
1642c         t = -log(dble(rndm4()))/bsd
1643c         teta = sqrt(t)/PC
1644c      10 va  =2d0*rndm4()-1d0 
1645c         vb = dble(rndm4())
1646c         va2 = va*va
1647c         vb2 = vb*vb
1648c         r2 = va2 + vb2
1649c         if ( r2.gt.1.d0) go to 10
1650c         tx = teta * (2.d0*va*vb) / r2
1651c         tz = teta * (va2 - vb2) / r2
1652
1653c         XP = XP + tx
1654c         YP = YP + tz
1655
1656!#################### end of "new" treatment ###############
1657
1658cc      ENDIF
1659C.  Inelastic interaction
1660cc      IF (rndm4() .LE. DZ/DLY) THEN
1661c        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)
1662cc        PROC = 'absorbed'
1663cc      ENDIF
1664c      write(*,*)'xp final:', xp, 'yp final', yp
1665cc      xp=xp/1000
1666cc      yp=yp/1000
1667
1668
1669!---------------DANIELE--------------------------------
1670! print out the energy loss and process experienced
1671!-------------------
1672
1673cc               if (PROC(1:2).eq.'AM')then
1674cc                 PROC_dan=1
1675cc               elseif (PROC(1:2).eq.'VR') then
1676cc                 PROC_dan=2
1677cc               elseif (PROC(1:2).eq.'CH')then
1678cc                 PROC_dan=3
1679cc               elseif (PROC(1:2).eq.'VC') then
1680cc                 PROC_dan=3
1681cc               elseif (PROC(1:3).eq.'out')then
1682cc                 PROC_dan=-1
1683cc              elseif (PROC(1:8).eq.'absorbed') then
1684cc                 PROC_dan=5
1685cc               elseif (PROC(1:2).eq.'DC')then
1686cc                 PROC_dan=6
1687cc               elseif (PROC(1:3).eq.'mcs')then               
1688cc                 PROC_dan=7
1689cc               elseif (PROC(1:4).eq.'diff')then 
1690cc                 PROC_dan=8
1691cc               endif
1692
1693cc      write(889,*) PC_in_dan, PC, PROC_dan 
1694
1695!--------END DANIELE------------------------------------
1696
1697cc      RETURN
1698cc      END
1699C     
1700cc      BLOCK DATA
1701cc      implicit none
1702cc      double precision DLAI(4),SAI(4)
1703cc      double precision DLRI(4),DLYI(4),AI(4),DES(4)
1704cc      double precision eUm(4)
1705cc      COMMON /ALAST/DLAI,SAI
1706cc      COMMON/CRYS/ DLRI,DLYI,AI,DES
1707cc      COMMON/eUc/  eUm
1708C-----4 substances: Si(110),W(110),C,Ge----------------------------
1709cc      DATA DLRI/0.09336,.0035,0.188,.023/    ! radiation  length(m)
1710cc     +    ,DLYI/.455, .096, .400, .162/      ! nuclear length(m)
1711cc     +    ,AI /0.96E-7, 0.56E-7, 0.63E-7, 1.E-7/  !Si110 1/2 interplan. dist. mm
1712cc     +    ,DES/0.56,  3.0,  0.6, 1./         ! energy deposition in subst(GeV/m)
1713cc     +    ,DLAI/1.6,  0.57, 2.2, 1.0/        ! elastic length(m)
1714cc     +    ,SAI /42.,  140.,  42., 50./       ! elastic scat. r.m.s(mr)
1715cc     +    ,eUm/21.34,  21.,   21.,   21./    ! only for Si(110) potent. [eV]
1716cc      END
1717
1718
1719
1720
1721
1722     
1723
1724C                        'MATH.F'   10.07.01
1725C1  RANNOR - Gauss distribution simulation procedure.
1726C2  RF RNDM- Â­Â  [0,1]        *
1727C3  RNDMST - .
1728C4  ...
1729C.**************************************************************************
1730C     subroutine for the generation of random numbers with a gaussian
1731C     distribution
1732C.**************************************************************************
1733C.**************************************************************************
1734      SUBROUTINE RANNOR(A,B) !gaussian with mean 0 and sigma 1
1735C
1736C1    Gauss distribution simulation procedure
1737C
1738      Y = RNDM(1.)
1739      IF (Y .EQ. 0.) Y = RNDM(1.)
1740      Z = RNDM(1.)
1741      X = 6.2831853*Z
1742      A1= SQRT(-2.0*ALOG(Y))
1743      A = A1*SIN(X)
1744      B = A1*COS(X)
1745C
1746      RETURN
1747      END
1748C
1749      real*4 function RNDM(rdummy)
1750      REAL*4          u,c,cd,cm, rrdummy
1751      COMMON/RANDOM/  u(97),c,cd,cm,i,j
1752C
1753      rrdummy = rdummy
1754C
1755      RNDM = u(i)-u(j)
1756      if (RNDM .LT. 0.) RNDM = RNDM+1.
1757      U(i) = RNDM
1758      i = i-1
1759      if ( i .EQ. 0 ) i = 97
1760      j = j-1
1761      if ( j .EQ. 0 ) j = 97
1762      c = c-cd
1763      if ( C .LT. 0.) c = c+cm
1764      RNDM = RNDM - c
1765      if ( RNDM .LT. 0. ) RNDM = RNDM+1.
1766C
1767      return
1768      END
1769C.**************************************************************************
1770C
1771C.**************************************************************************
1772      subroutine RNDMST(NA1,NA2,NA3,NB1)
1773      REAL*4 u,c,cd,cm
1774      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1775C
1776      ma1 = na1
1777      ma2 = na2
1778      ma3 = na3
1779      mb1 = nb1
1780      i   = 97
1781      j   = 33
1782C
1783      do ii2 = 1,97
1784        s = 0.0
1785        t = 0.5
1786        do ii1 = 1,24
1787          mat  = MOD(MOD(ma1*ma2,179)*ma3,179)
1788          ma1  = ma2
1789          ma2  = ma3
1790          ma3  = mat
1791          mb1  = MOD(53*mb1+1,169)
1792          if (MOD(MB1*MAT,64) .GE. 32 ) s = s+t
1793          t = 0.5*t
1794        end do
1795        u(ii2) = s
1796      end do
1797C
1798      c  =   362436./16777216.
1799      cd =  7654321./16777216.
1800      cm = 16777213./16777216.
1801C
1802      return
1803      END
1804C.**************************************************************************
1805C
1806C.**************************************************************************
1807      subroutine RNDMIN(uin,cin,cdin,cmin,iin,jin)
1808      REAL*4 uin(97),cin,cdin,cmin
1809      REAL*4 u,c,cd,cm
1810      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1811C
1812      do kkk = 1,97
1813        u(kkk) = uin(kkk)
1814      end do
1815C
1816      c  = cin
1817      cd = cdin
1818      cm = cmin
1819      i  = iin
1820      j  = jin
1821C
1822      return
1823      END
1824C.**************************************************************************
1825C
1826C.**************************************************************************
1827      subroutine RNDMOU(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
1828      REAL*4 uout(97),cout,cdout,cmout
1829      REAL*4 u,c,cd,cm
1830      COMMON/RANDOM/ u(97),c,cd,cm,i,j
1831C
1832      do kkk = 1,97
1833        uout(kkk) = u(kkk)
1834      end do
1835C
1836      COUT  = C
1837      CDOUT = CD
1838      CMOUT = CM
1839      IOUT  = I
1840      JOUT  = J
1841C
1842      return
1843      END
1844C.**************************************************************************
1845C
1846C.**************************************************************************
1847      subroutine RNDMTE(IO)
1848C
1849C.    *******************************************************************
1850C.    *  SUBROUTINE RNDMTE(IO)                                         
1851C.    *******************************************************************
1852C
1853      REAL*4 uu(97)
1854      REAL*4 u(6),x(6),d(6)
1855      DATA u / 6533892.0 , 14220222.0 ,  7275067.0 ,
1856     +         6172232.0 ,  8354498.0 , 10633180.0 /
1857C
1858      call RNDMOU(UU,CC,CCD,CCM,II,JJ)
1859      call RNDMST(12,34,56,78)
1860C
1861      do ii1 = 1,20000
1862        xx = rndm4()
1863      end do
1864C
1865      sd = 0.0
1866      do ii2 = 1,6
1867        x(II2)  = 4096.*(4096.*rndm4())
1868        d(II2)  = x(II2)-u(II2)
1869        sd = sd+d(II2)
1870      end do
1871C
1872      call RNDMIN(uu,cc,ccd,ccm,ii,jj)
1873      if (IO.EQ.1 .OR. SD .NE. 0.) write(6,10) (u(I),x(I),d(I),i=1,6)
1874C
1875 10   format('  === TEST OF THE RANDOM-GENERATOR ===',/,
1876     +       '    EXPECTED VALUE    CALCULATED VALUE     DIFFERENCE',/,
1877     +       6(F17.1,F20.1,F15.3,/),
1878     +       '  === END OF TEST ;',
1879     +       '  GENERATOR HAS THE SAME STATUS AS BEFORE CALLING RNDMTE')
1880      return
1881      END
1882
1883
Note: See TracBrowser for help on using the repository browser.