[12] | 1 | C.************************************************************************** |
---|
| 2 | C SUBROUTINE FOR THE MOVEMENTS OF THE PARTICLES IN THE CRYSTAL |
---|
| 3 | C.************************************************************************** |
---|
[15] | 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 | |
---|
| 12 | c SUBROUTINE CRYST(Length,layerflag,crypara,partpara,crysparaDes, |
---|
| 13 | c + crysparaDlri, |
---|
| 14 | c + crysparaDlyi,crysparaeUmIS,crysparaAIIS,PROC) |
---|
| 15 | C SUBROUTINE CRYST(IS,x,xp,y,yp,PC,Length,p1,p2,p3,p4,p5,p6,p7,PROC,LAYERFLAG,Chann,Dechan,L_chan) |
---|
| 16 | |
---|
[12] | 17 | c |
---|
| 18 | C Simple tranport protons in crystal 2 |
---|
| 19 | C-----------------------------------------------------------C |
---|
| 20 | C J - number of element C |
---|
| 21 | C S - longitudinal coordinate |
---|
| 22 | C IS - number of substance 1-4: Si,W,C,Ge(110) C |
---|
| 23 | C x,xp,y,yp - coordinates at input of crystal C |
---|
| 24 | C PC - momentum of particle*c [GeV] C |
---|
| 25 | C W - weigth of particle C |
---|
| 26 | C-----------------------------------------------------------C |
---|
| 27 | C |
---|
| 28 | C |
---|
| 29 | IMPLICIT none |
---|
| 30 | c |
---|
| 31 | C |
---|
[13] | 32 | C |
---|
[15] | 33 | integer LAYERFLAG |
---|
[13] | 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 |
---|
| 37 | C |
---|
| 38 | double precision p7 !common block: miscut |
---|
| 39 | CHARACTER*50 p8 !common block: PROC |
---|
[15] | 40 | |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | |
---|
[12] | 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 |
---|
| 51 | c 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 |
---|
[15] | 87 | c REAL*4 RNDM !random numbers |
---|
| 88 | c REAL*4 rndm4 |
---|
| 89 | c REAL*4 RAN_GAUSS |
---|
| 90 | double precision RNDM !random numbers |
---|
| 91 | double precision rndm4 |
---|
| 92 | double precision RAN_GAUSS |
---|
[12] | 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 | |
---|
[13] | 105 | C global variables shared by this file and the main file |
---|
| 106 | C "crystaltrack_dan.f". |
---|
| 107 | C common /Par_Cry1/ Cry_length, Rcurv,C_xmax,C_ymax,Alayer,C_orient |
---|
| 108 | C common /miscut/ miscut |
---|
| 109 | C common /Proc2/PROC |
---|
| 110 | |
---|
| 111 | C global variables shared by the subroutines in this file. |
---|
[15] | 112 | C COMMON/NPC/ NAM,ZN |
---|
| 113 | C COMMON/CRYS/ DLRI,DLYI,AI,DES |
---|
| 114 | C COMMON/eUc/ eUm ! |
---|
[12] | 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 | |
---|
[13] | 119 | c global variables shared by this file and the main file |
---|
| 120 | c "crystaltrack_dan.f" |
---|
[15] | 121 | c Cry_length=p1 |
---|
| 122 | c Rcurv=p2 |
---|
| 123 | c C_xmax=p3 |
---|
| 124 | c C_ymax=p4 |
---|
| 125 | c Alayer=p5 |
---|
| 126 | c C_orient=p6 |
---|
| 127 | c miscut=p7 |
---|
[13] | 128 | |
---|
[15] | 129 | IS = IS + 1; |
---|
| 130 | |
---|
[12] | 131 | c common/utils/ counter |
---|
| 132 | C |
---|
| 133 | c |
---|
| 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 | |
---|
| 150 | c DESt=((k*z(IS))/(anuc_cry2(IS)*betar**2))* |
---|
| 151 | c + (0.5*log((2.0d0*me*bgr*bgr*Tmax)/(Ime(IS)*Ime(IS))) |
---|
| 152 | c + -betar**2.0-log(plen/Ime(IS))-log(bgr)+0.5); |
---|
| 153 | |
---|
| 154 | c 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 | |
---|
| 166 | c miscut=0.001000 |
---|
| 167 | c |
---|
[14] | 168 | c write(*,*)"last miscut angle =",miscut |
---|
[12] | 169 | c write(*,*) 'enter crystal subroutine' |
---|
| 170 | c write(*,*) 'particle energy Gev :', PC |
---|
| 171 | c write(*,*) 'x_initial :', x |
---|
| 172 | c write(*,*) 'Length [m]:', Length |
---|
| 173 | c write(*,*) 'Random:', rndm4() |
---|
| 174 | c write(*,*)'xp',xp,'x',x , 's', s |
---|
| 175 | s=0 |
---|
| 176 | s_length=Rcurv*(sin(length/Rcurv)) ! |
---|
| 177 | L_chan=length |
---|
| 178 | |
---|
| 179 | C****************************************** |
---|
| 180 | C Start to call the crystal routine |
---|
| 181 | C the same as: SimCrys::bases() in |
---|
| 182 | C simcrys.cc |
---|
| 183 | C |
---|
[13] | 184 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
[12] | 185 | C |
---|
| 186 | C****************************************** |
---|
| 187 | C |
---|
| 188 | C-----------SimCrys::bases() in simcrys.cc |
---|
| 189 | C |
---|
| 190 | C 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 |
---|
| 199 | C---------------------------- |
---|
| 200 | |
---|
| 201 | |
---|
| 202 | C |
---|
| 203 | C-----------Crystal::Crystal() in crystal.cc |
---|
| 204 | C |
---|
| 205 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 206 | |
---|
| 207 | c FIRST CASE: p don't interact with crystal |
---|
| 208 | ymin = - C_ymax/2 |
---|
| 209 | ymax = C_ymax/2 |
---|
| 210 | C-------------------------------------------- |
---|
| 211 | |
---|
| 212 | C |
---|
| 213 | C-----------SimCrys::cross() in simcrys.cc |
---|
| 214 | C |
---|
| 215 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 216 | |
---|
| 217 | IF (y.LT.ymin .or. y.GT.ymax .or. x.gt.C_xmax) THEN |
---|
[15] | 218 | |
---|
| 219 | write(*,*) "y = ", y |
---|
| 220 | write(*,*) "ymin = ", ymin |
---|
| 221 | write(*,*) "ymax = ", ymax |
---|
| 222 | write(*,*) "x = ", x |
---|
| 223 | write(*,*) "C_xmax = ", C_xmax |
---|
| 224 | |
---|
[12] | 225 | x = x+xp*s_length |
---|
| 226 | y = y+yp*s_length |
---|
| 227 | PROC='out' |
---|
| 228 | GOTO 111 |
---|
[13] | 229 | |
---|
| 230 | C--------------------------- |
---|
| 231 | |
---|
| 232 | |
---|
| 233 | C-------------------------------------------- |
---|
| 234 | |
---|
| 235 | C |
---|
| 236 | C-----------SimCrys::layer() in simcrys.cc |
---|
| 237 | C |
---|
| 238 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 239 | |
---|
[12] | 240 | c 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)) |
---|
[13] | 250 | |
---|
[12] | 251 | if (s .ge. s_length) s=s_length |
---|
[13] | 252 | |
---|
[12] | 253 | x=(xp)*s+x0 |
---|
| 254 | Length_xs=sqrt((x-x0)**2+s**2) |
---|
[13] | 255 | |
---|
[12] | 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) |
---|
[15] | 278 | |
---|
| 279 | LAYERFLAG = 0; |
---|
[12] | 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!' |
---|
[15] | 289 | LAYERFLAG = 0; |
---|
[12] | 290 | GOTO 111 |
---|
| 291 | END IF |
---|
[13] | 292 | |
---|
| 293 | |
---|
| 294 | C-------------------------------------------- |
---|
| 295 | |
---|
| 296 | C |
---|
| 297 | C-----------SimCrys::parameters() in simcrys.cc |
---|
| 298 | C |
---|
| 299 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 300 | C |
---|
[12] | 301 | c |
---|
| 302 | c THIRD CASE: the p interacts with the crystal. |
---|
| 303 | C. Define typical angles/probabilities for orientation 110 |
---|
| 304 | c |
---|
| 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 |
---|
| 311 | c write(*,*) "Critical Radius: ",Rcrit |
---|
[15] | 312 | c write(*,*) "eUm(IS): ", eUm(IS) |
---|
| 313 | c write(*,*) "PC: ", PC |
---|
| 314 | c write(*,*) "xpcrit0: ", xpcrit0 |
---|
[12] | 315 | xpcrit = xpcrit0*(Rcurv-Rcrit)/Rcurv ! critical angle for curved crystal |
---|
| 316 | c----------------valentina approx----------- |
---|
| 317 | c xpcrit = xpcrit0*(1-(Rcrit/Rcurv))**0.5 |
---|
| 318 | c 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 |
---|
| 330 | c 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 |
---|
| 336 | c 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 |
---|
| 340 | cc----------------valentina approx----------- |
---|
| 341 | c Ang_avr=-(xpcrit+xpcrit0) |
---|
| 342 | c Ang_rms=(xpcrit0-xpcrit)/2 |
---|
| 343 | cc-----------end valentina approx-------------- |
---|
| 344 | |
---|
| 345 | c write(*,*) "Rcrit" , Rcrit,"Rcurv",Rcurv, |
---|
| 346 | c c "Ratio: ",ratio,"average VR angle:", ang_avr*1e6,"+-", |
---|
| 347 | c c ang_rms*1e6, "ang crit:", xpcrit0*1e6,xpcrit*1e6 |
---|
| 348 | c |
---|
| 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 |
---|
[13] | 354 | |
---|
| 355 | C--------------------------- |
---|
| 356 | |
---|
| 357 | |
---|
| 358 | |
---|
| 359 | |
---|
| 360 | C-------------------------------------------- |
---|
| 361 | |
---|
| 362 | C |
---|
| 363 | C-----------SimCrys::channel() in simcrys.cc |
---|
| 364 | C |
---|
| 365 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
[12] | 366 | c |
---|
| 367 | C. 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 |
---|
| 374 | c 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 |
---|
| 378 | c 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 | |
---|
| 387 | c 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 |
---|
[13] | 415 | else !channeling |
---|
[12] | 416 | PROC='CH' |
---|
| 417 | xpin=XP |
---|
| 418 | ypin=YP |
---|
| 419 | |
---|
| 420 | c 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 | |
---|
| 424 | c 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 |
---|
| 438 | c xp = xp + Dxp + 2.0*(rndm4()-0.5)*xpcrit |
---|
| 439 | y = y + s_length * yp |
---|
| 440 | c !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 |
---|
| 450 | c Dxp=0.5*(xp_rel/xpcrit+1)*Ang_avr |
---|
| 451 | c xp=xp+Dxp+Ang_rms*RAN_GAUSS(1.) |
---|
| 452 | !next line new from sasha |
---|
[15] | 453 | |
---|
| 454 | c write(*,*) "y = ", y |
---|
[17] | 455 | c write(*,*) "before enter VR: x = ", x, " xp = ",xp |
---|
| 456 | c write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC |
---|
[15] | 457 | |
---|
[12] | 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) |
---|
[15] | 466 | |
---|
[17] | 467 | c write(*,*) "after enter VR: x = ", x, " xp = ",xp |
---|
| 468 | c write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC |
---|
[15] | 469 | |
---|
[12] | 470 | x = x + 0.5*s_length * xp |
---|
| 471 | y = y + 0.5*s_length * yp |
---|
[13] | 472 | ENDIF |
---|
| 473 | |
---|
| 474 | C------------------------------ ! |
---|
| 475 | C |
---|
| 476 | C |
---|
| 477 | C-------------------------------------------- |
---|
| 478 | C |
---|
| 479 | C-----------SimCrys::reflection() in simcrys.cc |
---|
| 480 | C |
---|
| 481 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 482 | |
---|
[12] | 483 | c case 3-2: no good for channeling. check if the can VR |
---|
| 484 | ELSE |
---|
| 485 | Lrefl = (xp_rel)*Rcurv ! distance of refl. point [m] |
---|
| 486 | c 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' |
---|
[15] | 493 | |
---|
[12] | 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) |
---|
[15] | 502 | |
---|
[17] | 503 | c write(*,*) "before enter VR: x = ", x, " xp = ",xp |
---|
| 504 | c write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC |
---|
[15] | 505 | |
---|
[12] | 506 | CALL CALC_ION_LOSS(IS,PC,s_length-Srefl,DESt) |
---|
[15] | 507 | |
---|
[17] | 508 | c write(*,*) "before enter VR: x = ", x, " xp = ",xp |
---|
| 509 | c write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC |
---|
| 510 | c write(*,*) "IS = ", IS, " NAM = ", NAM |
---|
| 511 | c write(*,*) "s_length-Srefl = ",s_length-Srefl |
---|
| 512 | c write(*,*) "DESt = ", DESt, " DLYi(IS) = ", DLYi(IS) |
---|
| 513 | c write(*,*) "DLRi(IS) = ", DLRi(IS) |
---|
[15] | 514 | |
---|
[12] | 515 | CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS), |
---|
| 516 | + DLRi(IS),xp ,yp,PC) |
---|
[15] | 517 | |
---|
[17] | 518 | c write(*,*) "after enter VR: x = ", x, " xp = ",xp |
---|
| 519 | c write(*,*) "y = ",y, "yp = ", yp, "PC = ", PC |
---|
[15] | 520 | |
---|
[12] | 521 | x = x + 0.5 * xp * (s_length - Srefl) |
---|
| 522 | y = y + 0.5 * yp * (s_length - Srefl) |
---|
[13] | 523 | |
---|
| 524 | ELSE !opt 2: VC |
---|
[12] | 525 | x = x + xp * Srefl |
---|
| 526 | y = y + yp * Srefl |
---|
| 527 | c TLdech2= 0.00011*PC**0.25*(1.-1./ratio)**2 ! dechanneling length(m) |
---|
| 528 | c Dechan = log(1.-rndm4()) |
---|
| 529 | c 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 |
---|
| 586 | C. 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 |
---|
[13] | 602 | C--------------------------- |
---|
| 603 | C |
---|
| 604 | C |
---|
[15] | 605 | C print the crystal parameters to an external file |
---|
| 606 | open(unit=833,file='crys_para_check.dat') |
---|
| 607 | |
---|
[12] | 608 | c if (counter .eq. 0) then |
---|
| 609 | 111 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) |
---|
| 617 | c counter=1 |
---|
| 618 | c endif |
---|
| 619 | c |
---|
[17] | 620 | c WRITE(*,*)'Crystal process: ',PROC |
---|
[12] | 621 | c write(*,*)'xp_final :', xp |
---|
| 622 | c WRITE(*,*)'Crystal process: ',PROC,'Chann Angle',Ch_angle/1000, + |
---|
| 623 | c 1'Critical angle: ', xpcrit/1000 |
---|
| 624 | c 2 DLRI(IS),DLYI(IS),AI(IS),DES(IS),eUm(IS),IS,ZN,NAM,C_orient !,W |
---|
[15] | 625 | |
---|
| 626 | close(833) |
---|
[12] | 627 | END |
---|
| 628 | |
---|
| 629 | C.************************************************************************** |
---|
| 630 | C subroutine for the calculazion of the energy loss by ionization |
---|
| 631 | C.************************************************************************** |
---|
| 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 |
---|
[15] | 645 | c REAL*4 RNDM4 |
---|
| 646 | double precision RNDM4 |
---|
[12] | 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 | |
---|
| 683 | c write(*,*)cs_tail,prob_tail,ranc,EnLo*DZ |
---|
| 684 | |
---|
| 685 | RETURN |
---|
| 686 | END |
---|
| 687 | |
---|
| 688 | |
---|
[13] | 689 | |
---|
| 690 | C-------------------------------------------- |
---|
| 691 | |
---|
| 692 | C |
---|
| 693 | C-----------SimCrys::move_am() in simcrys.cc |
---|
| 694 | C |
---|
| 695 | C The two routines are NOT the same!!!!!!! |
---|
| 696 | C |
---|
| 697 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 698 | |
---|
[12] | 699 | C.************************************************************************** |
---|
| 700 | C subroutine for the movement in the amorphous |
---|
| 701 | C.************************************************************************** |
---|
| 702 | SUBROUTINE MOVE_AM_(IS,NAM,DZ,DEI,DLY,DLr, XP,YP,PC) |
---|
| 703 | C. 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 |
---|
[15] | 712 | c REAL*4 RNDM4 |
---|
| 713 | c REAL*4 RAN_GAUSS |
---|
| 714 | double precision RNDM4 |
---|
| 715 | double precision RAN_GAUSS |
---|
[12] | 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 |
---|
[15] | 764 | double precision ran_x, ran_y |
---|
[12] | 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 |
---|
[13] | 828 | C------------------------------ |
---|
[12] | 829 | |
---|
[13] | 830 | C-------------------------------------------- |
---|
[12] | 831 | |
---|
[13] | 832 | C |
---|
| 833 | C-----------SimCrys::cross() in simcrys.cc |
---|
| 834 | C ----The physics are the same. |
---|
| 835 | C Commented by Jianfeng Zhang @ LAL, 29/04/2014 |
---|
| 836 | |
---|
[12] | 837 | !--------- Multiple Coulomb Scattering --------- |
---|
| 838 | |
---|
[15] | 839 | |
---|
[17] | 840 | c write(*,*) "In Move_AM(): " |
---|
[15] | 841 | |
---|
| 842 | |
---|
[12] | 843 | xp=xp*1000 |
---|
| 844 | yp=yp*1000 |
---|
[17] | 845 | c write(*,*)'xp initial:', xp, 'yp initial', yp |
---|
[12] | 846 | C. DEI - dE/dx stoping energy |
---|
| 847 | PC = PC - DEI*DZ ! energy lost because of ionization process[GeV] |
---|
| 848 | |
---|
| 849 | C. DYA - rms of coloumb scattering |
---|
| 850 | DYA = (13.6/PC)*SQRT(DZ/DLr) !MCS (mrad) |
---|
[17] | 851 | c write(*,*)'dya= ',dya |
---|
[15] | 852 | |
---|
| 853 | ran_x = RAN_GAUSS(1.) |
---|
[17] | 854 | c write(*,*) "ran_x = ", ran_x |
---|
[15] | 855 | |
---|
| 856 | ran_y = RAN_GAUSS(1.) |
---|
[17] | 857 | c write(*,*) "ran_y = ", ran_y |
---|
[15] | 858 | |
---|
| 859 | kxmcs = DYA*ran_x |
---|
[16] | 860 | kymcs = DYA*ran_y |
---|
[15] | 861 | |
---|
| 862 | c kxmcs = DYA*RAN_GAUSS(1.) |
---|
[16] | 863 | c kymcs = DYA*RAN_GAUSS(1.) |
---|
[12] | 864 | |
---|
| 865 | XP = xp+kxmcs |
---|
| 866 | yp = yp+kymcs |
---|
[15] | 867 | |
---|
[16] | 868 | c write(*,*) "ran_x = ", ran_x, "ran_y = ", ran_y |
---|
[17] | 869 | c write(*,*) "kxmcs = ", kxmcs, "kymcs = ", kymcs |
---|
| 870 | c write(*,*)'xp + kxmcs:', XP, 'yp + kxmcs', YP |
---|
[15] | 871 | |
---|
[12] | 872 | cc XP = XP+DYA*RAN_GAUSS(1.) |
---|
| 873 | cc 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 | |
---|
| 883 | c 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 | |
---|
| 901 | c 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 | |
---|
[17] | 940 | c write(*,*) "ichoix = ", ichoix |
---|
| 941 | c write(*,*)'xp final:', xp, 'yp final', yp |
---|
[16] | 942 | |
---|
[12] | 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 | |
---|
| 953 | c 100 va =2d0*rndm4()-1d0 |
---|
| 954 | c vb = dble(rndm4()) |
---|
| 955 | c va2 = va*va |
---|
| 956 | c vb2 = vb*vb |
---|
| 957 | c r2 = va2 + vb2 |
---|
| 958 | c if ( r2.gt.1.d0) go to 100 |
---|
| 959 | c tx = teta * (2.d0*va*vb) / r2 |
---|
| 960 | c tz = teta * (va2 - vb2) / r2 |
---|
| 961 | |
---|
| 962 | cc 100 va =2d0*rndm4()-1d0 |
---|
| 963 | cc vb =2d0*rndm4()-1d0 |
---|
| 964 | cc va2 = va*va |
---|
| 965 | cc vb2 = vb*vb |
---|
| 966 | cc r2 = va2 + vb2 |
---|
| 967 | cc if ( r2.gt.1.d0) go to 100 |
---|
| 968 | cc f = SQRT(-2d0*LOG(r2)/r2) |
---|
| 969 | cc tx = teta*f*va |
---|
| 970 | cc 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 | |
---|
[15] | 981 | |
---|
[17] | 982 | c write(*,*) "tx = ", tx, " tz = ", tz |
---|
[15] | 983 | |
---|
[12] | 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 | |
---|
[17] | 992 | c write(*,*)'xp final:', xp, 'yp final', yp |
---|
[12] | 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 | |
---|
| 1022 | c write(889,*) PC_in_dan, PC, PROC_dan, tx,tz, |
---|
| 1023 | c & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) |
---|
| 1024 | |
---|
| 1025 | |
---|
| 1026 | RETURN |
---|
| 1027 | END |
---|
| 1028 | |
---|
| 1029 | |
---|
| 1030 | !--------- Multiple Coulomb Scattering --------- |
---|
| 1031 | |
---|
| 1032 | cc xp=xp*1000 |
---|
| 1033 | cc yp=yp*1000 |
---|
| 1034 | c write(*,*)'xp initial:', xp, 'yp initial', yp |
---|
| 1035 | C. DEI - dE/dx stoping energy |
---|
| 1036 | cc PC = PC - DEI*DZ ! energy lost beacause of ionization process[GeV] |
---|
| 1037 | |
---|
| 1038 | C. DYA - rms of coloumb scattering |
---|
| 1039 | cc DYA = (13.6/PC)*SQRT(DZ/DLr) !MCS (mrad) |
---|
| 1040 | c write(*,*)'dya=',dya |
---|
| 1041 | |
---|
| 1042 | cc kxmcs = DYA*RAN_GAUSS(1.) |
---|
| 1043 | cc kymcs = DYA*RAN_GAUSS(1.) |
---|
| 1044 | |
---|
| 1045 | cc XP = xp+kxmcs |
---|
| 1046 | cc yp = yp+kymcs |
---|
| 1047 | |
---|
| 1048 | c XP = XP+DYA*RAN_GAUSS(1.) |
---|
| 1049 | c YP = YP+DYA*RAN_GAUSS(1.) |
---|
| 1050 | |
---|
| 1051 | cc xp = xp/1000 |
---|
| 1052 | cc yp = yp/1000 |
---|
| 1053 | |
---|
| 1054 | !--------END DANIELE------------------------------------ |
---|
| 1055 | |
---|
| 1056 | |
---|
| 1057 | C.************************************************************************** |
---|
| 1058 | C DANIELE subroutine for check if a nuclear interaction happen while in channeling |
---|
| 1059 | C.************************************************************************** |
---|
| 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 |
---|
[15] | 1070 | c REAL*4 RNDM4 |
---|
| 1071 | c REAL*4 RAN_GAUSS |
---|
| 1072 | double precision RNDM4 |
---|
| 1073 | double precision RAN_GAUSS |
---|
[12] | 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 | |
---|
| 1201 | c 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 | |
---|
| 1224 | c write(889,*) x_i,pv,Ueff,Et,Ec,N_am,avrrho, |
---|
| 1225 | c + 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 | |
---|
| 1268 | c write(*,*)'xp initial:', xp, 'yp initial', yp |
---|
| 1269 | C. DEI - dE/dx stoping energy |
---|
| 1270 | ! PC = PC - DEI*DZ ! energy lost beacause of ionization process[GeV] |
---|
| 1271 | |
---|
| 1272 | C. DYA - rms of coloumb scattering |
---|
| 1273 | ! DYA = (13.6/PC)*SQRT(DZ/DLr) !MCS (mrad) |
---|
| 1274 | c 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 | |
---|
| 1282 | cc XP = XP+DYA*RAN_GAUSS(1.) |
---|
| 1283 | cc 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 | |
---|
| 1295 | c 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 | |
---|
| 1316 | c 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 | |
---|
| 1365 | c 100 va =2d0*rndm4()-1d0 |
---|
| 1366 | c vb = dble(rndm4()) |
---|
| 1367 | c va2 = va*va |
---|
| 1368 | c vb2 = vb*vb |
---|
| 1369 | c r2 = va2 + vb2 |
---|
| 1370 | c if ( r2.gt.1.d0) go to 100 |
---|
| 1371 | c tx = teta * (2.d0*va*vb) / r2 |
---|
| 1372 | c tz = teta * (va2 - vb2) / r2 |
---|
| 1373 | |
---|
| 1374 | cc 100 va =2d0*rndm4()-1d0 |
---|
| 1375 | cc vb =2d0*rndm4()-1d0 |
---|
| 1376 | cc va2 = va*va |
---|
| 1377 | cc vb2 = vb*vb |
---|
| 1378 | cc r2 = va2 + vb2 |
---|
| 1379 | cc if ( r2.gt.1.d0) go to 100 |
---|
| 1380 | cc f = SQRT(-2d0*LOG(r2)/r2) |
---|
| 1381 | cc tx = teta*f*va |
---|
| 1382 | cc 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 | |
---|
| 1431 | c write(889,*) PC_in_dan, PC, PROC_dan, tx,tz, |
---|
| 1432 | c & xp-xp_in-(kxmcs/1000), yp-yp_in-(kymcs/1000) |
---|
| 1433 | |
---|
| 1434 | !-------- Block Data ---------- |
---|
| 1435 | |
---|
| 1436 | RETURN |
---|
| 1437 | END |
---|
| 1438 | C |
---|
| 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 |
---|
| 1447 | C-----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 |
---|
| 1490 | c 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 | |
---|
| 1593 | cc xp=xp*1000 |
---|
| 1594 | cc yp=yp*1000 |
---|
| 1595 | c write(*,*)'xp initial:', xp, 'yp initial', yp |
---|
| 1596 | C. DEI - dE/dx stoping energy |
---|
| 1597 | cc PC = PC - DEI*DZ ! energy lost beacause of ionization process[GeV] |
---|
| 1598 | C. Coulomb scattering |
---|
| 1599 | C. DYA - rms of coloumb scattering |
---|
| 1600 | cc DYA = (14.0/PC)*SQRT(DZ/DLr) !MCS (mrad) |
---|
| 1601 | c write(*,*)'dya=',dya |
---|
| 1602 | cc XP = XP+DYA*RAN_GAUSS(1.) |
---|
| 1603 | cc YP = YP+DYA*RAN_GAUSS(1.) |
---|
| 1604 | |
---|
| 1605 | cc if(NAM .eq. 0) return |
---|
| 1606 | C. Elastic scattering |
---|
| 1607 | cc IF (rndm4() .LE. DZ/DLAI(IS)) THEN |
---|
| 1608 | cc PROC = 'mcs' |
---|
| 1609 | c write(*,*)'case 1' |
---|
| 1610 | c XP = XP+SAI(IS)*RAN_GAUSS(1.)/PC |
---|
| 1611 | c YP = YP+SAI(IS)*RAN_GAUSS(1.)/PC |
---|
| 1612 | cc xp = xp+196.*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane |
---|
| 1613 | cc yp = yp+196.*RAN_GAUSS(1.)/PC |
---|
| 1614 | cc ENDIF |
---|
| 1615 | C. Diffraction interaction |
---|
| 1616 | cc IF (rndm4() .LE. DZ/(DLY*6.143)) THEN |
---|
| 1617 | cc PROC = 'diff' |
---|
| 1618 | c write(*,*)'case 2' |
---|
| 1619 | |
---|
| 1620 | !###################################################### |
---|
| 1621 | ! daniele: comment old treatement to implement the same as standard SixTrack |
---|
| 1622 | !###################################################### |
---|
| 1623 | |
---|
| 1624 | c XP = XP+ 257.0*RAN_GAUSS(1.)/PC ! angle elast. scattering in R plane |
---|
| 1625 | c YP = YP+ 257.0*RAN_GAUSS(1.)/PC |
---|
| 1626 | c W_p = rndm4() |
---|
| 1627 | c 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 | |
---|
| 1631 | c ecmsq = 2 * 0.93828d0 * PC |
---|
| 1632 | c xln15s=log(0.15*ecmsq) |
---|
| 1633 | c bpp = 8.5d0 + 1.086d0 * log(sqrt(ecmsq)) |
---|
| 1634 | c xm2 = exp( dble(rndm4()) * xln15s ) |
---|
| 1635 | c PC = PC *(1.d0 - xm2/ecmsq) |
---|
| 1636 | |
---|
| 1637 | c if ( xm2 .lt. 2.d0 ) then |
---|
| 1638 | c bsd = 2.d0 * bpp |
---|
| 1639 | c elseif (( xm2 .ge. 2.d0 ).and. ( xm2 .le. 5.d0 )) then |
---|
| 1640 | c bsd = (106.d0-17.d0*xm2) * bpp / 26.d0 |
---|
| 1641 | c elseif ( xm2 .gt. 5.d0 ) then |
---|
| 1642 | c bsd = 7.d0 * bpp / 12.d0 |
---|
| 1643 | c endif |
---|
| 1644 | |
---|
| 1645 | c t = -log(dble(rndm4()))/bsd |
---|
| 1646 | c teta = sqrt(t)/PC |
---|
| 1647 | c 10 va =2d0*rndm4()-1d0 |
---|
| 1648 | c vb = dble(rndm4()) |
---|
| 1649 | c va2 = va*va |
---|
| 1650 | c vb2 = vb*vb |
---|
| 1651 | c r2 = va2 + vb2 |
---|
| 1652 | c if ( r2.gt.1.d0) go to 10 |
---|
| 1653 | c tx = teta * (2.d0*va*vb) / r2 |
---|
| 1654 | c tz = teta * (va2 - vb2) / r2 |
---|
| 1655 | |
---|
| 1656 | c XP = XP + tx |
---|
| 1657 | c YP = YP + tz |
---|
| 1658 | |
---|
| 1659 | !#################### end of "new" treatment ############### |
---|
| 1660 | |
---|
| 1661 | cc ENDIF |
---|
| 1662 | C. Inelastic interaction |
---|
| 1663 | cc IF (rndm4() .LE. DZ/DLY) THEN |
---|
| 1664 | c 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) |
---|
| 1665 | cc PROC = 'absorbed' |
---|
| 1666 | cc ENDIF |
---|
| 1667 | c write(*,*)'xp final:', xp, 'yp final', yp |
---|
| 1668 | cc xp=xp/1000 |
---|
| 1669 | cc yp=yp/1000 |
---|
| 1670 | |
---|
| 1671 | |
---|
| 1672 | !---------------DANIELE-------------------------------- |
---|
| 1673 | ! print out the energy loss and process experienced |
---|
| 1674 | !------------------- |
---|
| 1675 | |
---|
| 1676 | cc if (PROC(1:2).eq.'AM')then |
---|
| 1677 | cc PROC_dan=1 |
---|
| 1678 | cc elseif (PROC(1:2).eq.'VR') then |
---|
| 1679 | cc PROC_dan=2 |
---|
| 1680 | cc elseif (PROC(1:2).eq.'CH')then |
---|
| 1681 | cc PROC_dan=3 |
---|
| 1682 | cc elseif (PROC(1:2).eq.'VC') then |
---|
| 1683 | cc PROC_dan=3 |
---|
| 1684 | cc elseif (PROC(1:3).eq.'out')then |
---|
| 1685 | cc PROC_dan=-1 |
---|
| 1686 | cc elseif (PROC(1:8).eq.'absorbed') then |
---|
| 1687 | cc PROC_dan=5 |
---|
| 1688 | cc elseif (PROC(1:2).eq.'DC')then |
---|
| 1689 | cc PROC_dan=6 |
---|
| 1690 | cc elseif (PROC(1:3).eq.'mcs')then |
---|
| 1691 | cc PROC_dan=7 |
---|
| 1692 | cc elseif (PROC(1:4).eq.'diff')then |
---|
| 1693 | cc PROC_dan=8 |
---|
| 1694 | cc endif |
---|
| 1695 | |
---|
| 1696 | cc write(889,*) PC_in_dan, PC, PROC_dan |
---|
| 1697 | |
---|
| 1698 | !--------END DANIELE------------------------------------ |
---|
| 1699 | |
---|
| 1700 | cc RETURN |
---|
| 1701 | cc END |
---|
| 1702 | C |
---|
| 1703 | cc BLOCK DATA |
---|
| 1704 | cc implicit none |
---|
| 1705 | cc double precision DLAI(4),SAI(4) |
---|
| 1706 | cc double precision DLRI(4),DLYI(4),AI(4),DES(4) |
---|
| 1707 | cc double precision eUm(4) |
---|
| 1708 | cc COMMON /ALAST/DLAI,SAI |
---|
| 1709 | cc COMMON/CRYS/ DLRI,DLYI,AI,DES |
---|
| 1710 | cc COMMON/eUc/ eUm |
---|
| 1711 | C-----4 substances: Si(110),W(110),C,Ge---------------------------- |
---|
| 1712 | cc DATA DLRI/0.09336,.0035,0.188,.023/ ! radiation length(m) |
---|
| 1713 | cc + ,DLYI/.455, .096, .400, .162/ ! nuclear length(m) |
---|
| 1714 | cc + ,AI /0.96E-7, 0.56E-7, 0.63E-7, 1.E-7/ !Si110 1/2 interplan. dist. mm |
---|
| 1715 | cc + ,DES/0.56, 3.0, 0.6, 1./ ! energy deposition in subst(GeV/m) |
---|
| 1716 | cc + ,DLAI/1.6, 0.57, 2.2, 1.0/ ! elastic length(m) |
---|
| 1717 | cc + ,SAI /42., 140., 42., 50./ ! elastic scat. r.m.s(mr) |
---|
| 1718 | cc + ,eUm/21.34, 21., 21., 21./ ! only for Si(110) potent. [eV] |
---|
| 1719 | cc END |
---|
| 1720 | |
---|
| 1721 | |
---|
| 1722 | |
---|
| 1723 | |
---|
| 1724 | |
---|
| 1725 | |
---|
| 1726 | |
---|
| 1727 | C 'MATH.F' 10.07.01 |
---|
| 1728 | C1 RANNOR - Gauss distribution simulation procedure. |
---|
| 1729 | C2 RF RNDM- ÂÂ [0,1] * |
---|
| 1730 | C3 RNDMST - . |
---|
| 1731 | C4 ... |
---|
| 1732 | C.************************************************************************** |
---|
| 1733 | C subroutine for the generation of random numbers with a gaussian |
---|
| 1734 | C distribution |
---|
| 1735 | C.************************************************************************** |
---|
| 1736 | C.************************************************************************** |
---|
| 1737 | SUBROUTINE RANNOR(A,B) !gaussian with mean 0 and sigma 1 |
---|
| 1738 | C |
---|
| 1739 | C1 Gauss distribution simulation procedure |
---|
| 1740 | C |
---|
| 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) |
---|
| 1748 | C |
---|
| 1749 | RETURN |
---|
| 1750 | END |
---|
| 1751 | C |
---|
| 1752 | real*4 function RNDM(rdummy) |
---|
| 1753 | REAL*4 u,c,cd,cm, rrdummy |
---|
| 1754 | COMMON/RANDOM/ u(97),c,cd,cm,i,j |
---|
| 1755 | C |
---|
| 1756 | rrdummy = rdummy |
---|
| 1757 | C |
---|
| 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. |
---|
| 1769 | C |
---|
| 1770 | return |
---|
| 1771 | END |
---|
| 1772 | C.************************************************************************** |
---|
| 1773 | C |
---|
| 1774 | C.************************************************************************** |
---|
| 1775 | subroutine RNDMST(NA1,NA2,NA3,NB1) |
---|
| 1776 | REAL*4 u,c,cd,cm |
---|
| 1777 | COMMON/RANDOM/ u(97),c,cd,cm,i,j |
---|
| 1778 | C |
---|
| 1779 | ma1 = na1 |
---|
| 1780 | ma2 = na2 |
---|
| 1781 | ma3 = na3 |
---|
| 1782 | mb1 = nb1 |
---|
| 1783 | i = 97 |
---|
| 1784 | j = 33 |
---|
| 1785 | C |
---|
| 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 |
---|
| 1800 | C |
---|
| 1801 | c = 362436./16777216. |
---|
| 1802 | cd = 7654321./16777216. |
---|
| 1803 | cm = 16777213./16777216. |
---|
| 1804 | C |
---|
| 1805 | return |
---|
| 1806 | END |
---|
| 1807 | C.************************************************************************** |
---|
| 1808 | C |
---|
| 1809 | C.************************************************************************** |
---|
| 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 |
---|
| 1814 | C |
---|
| 1815 | do kkk = 1,97 |
---|
| 1816 | u(kkk) = uin(kkk) |
---|
| 1817 | end do |
---|
| 1818 | C |
---|
| 1819 | c = cin |
---|
| 1820 | cd = cdin |
---|
| 1821 | cm = cmin |
---|
| 1822 | i = iin |
---|
| 1823 | j = jin |
---|
| 1824 | C |
---|
| 1825 | return |
---|
| 1826 | END |
---|
| 1827 | C.************************************************************************** |
---|
| 1828 | C |
---|
| 1829 | C.************************************************************************** |
---|
| 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 |
---|
| 1834 | C |
---|
| 1835 | do kkk = 1,97 |
---|
| 1836 | uout(kkk) = u(kkk) |
---|
| 1837 | end do |
---|
| 1838 | C |
---|
| 1839 | COUT = C |
---|
| 1840 | CDOUT = CD |
---|
| 1841 | CMOUT = CM |
---|
| 1842 | IOUT = I |
---|
| 1843 | JOUT = J |
---|
| 1844 | C |
---|
| 1845 | return |
---|
| 1846 | END |
---|
| 1847 | C.************************************************************************** |
---|
| 1848 | C |
---|
| 1849 | C.************************************************************************** |
---|
| 1850 | subroutine RNDMTE(IO) |
---|
| 1851 | C |
---|
| 1852 | C. ******************************************************************* |
---|
| 1853 | C. * SUBROUTINE RNDMTE(IO) |
---|
| 1854 | C. ******************************************************************* |
---|
| 1855 | C |
---|
| 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 / |
---|
| 1860 | C |
---|
| 1861 | call RNDMOU(UU,CC,CCD,CCM,II,JJ) |
---|
| 1862 | call RNDMST(12,34,56,78) |
---|
| 1863 | C |
---|
| 1864 | do ii1 = 1,20000 |
---|
| 1865 | xx = rndm4() |
---|
| 1866 | end do |
---|
| 1867 | C |
---|
| 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 |
---|
| 1874 | C |
---|
| 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) |
---|
| 1877 | C |
---|
| 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 | |
---|