c c SENECA 1.01 c c COSMIC RAY AIR SHOWER SIMULATION WITH CASCADE EQUATIONS c main author: Hans-Joachim Drescher drescher@th.physik.uni-frankfurt.de c c $Id: conex.F,v 1.29 2003/03/11 10:47:17 drescher Exp $ c c $Log: conex.F,v $ c Revision 1.29 2003/03/11 10:47:17 drescher c inclined density precision fixed c c Revision 1.28 2003/02/25 10:28:19 drescher c interaction test c c Revision 1.27 2003/02/24 17:49:35 hjd1 c hadron test c c Revision 1.26 2003/02/18 16:36:09 hjd1 c hadronic models update c ---------------------------------------------------------------------- c c Revision 1.25 2003/02/01 17:12:24 hjd1 c auger detectors/hadronic model check c c Revision 1.24 2002/12/19 20:12:20 hjd1 c NA c c Revision 1.23 2002/12/17 21:23:22 hjd1 c removed files c c Revision 1.22 2002/12/09 17:08:54 hjd1 c divers c c Revision 1.21 2002/12/05 15:23:40 hjd1 c variable names updated c c Revision 1.20 2002/11/15 20:14:58 hjd1 c atmospheric depth issues c c Revision 1.19 2002/11/12 23:10:06 hjd1 c f in caskip/emcip c c Revision 1.18 2002/11/06 03:04:47 hjd1 c MeV/GeV Bug in coemc - kshort/klong in coneu c c Revision 1.17 2002/11/05 18:27:43 hajo c curved option c c Revision 1.15 2002/11/01 15:38:23 hjd1 c routine c c Revision 1.14 2002/10/25 21:55:10 hjd1 c weekend update c c Revision 1.13 2002/10/24 00:47:25 hjd1 c minor c c Revision 1.12 2002/10/24 00:31:27 hjd1 c platform dependent changes c c Revision 1.11 2002/10/23 16:09:04 hjd1 c det area/glauber c c Revision 1.10 2002/10/22 14:32:01 hjd1 c consistent analysis of cascade equations c c Revision 1.9 2002/10/21 21:54:18 hjd1 c display/user_var c c Revision 1.8 2002/10/20 19:53:55 hjd1 c cxana c c Revision 1.7 2002/10/19 14:56:24 hjd1 c cana example/getvalue/gdb-batch mode c c Revision 1.6 2002/10/15 14:30:09 hjd1 c egcalor in Geant/cana-option c c Revision 1.5 2002/10/14 14:13:28 hjd1 c fbase/#ifndef c c Revision 1.4 2002/10/07 19:05:22 hjd1 c detector particle output option c c Revision 1.3 2002/09/27 21:28:05 hjd1 c bugfix ares c c Revision 1.2 2002/09/27 20:22:13 hjd1 c area problem with detectors c c Revision 1.1 2002/09/26 19:13:13 hjd1 c start of CVS c c project startet: 2000-11-01 hans-joachim drescher #include "config.h" program seneca PARAMETER (TWOPI=6.28318530717958648,PI=3.14159265358979324) COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax parameter (mxblk=12) parameter (mxstk = mxblk*50*2) COMMON /DEBUG/ DEBUG common /cdebug/icd,icdsrc common /cstack/stack(mxstk,2),istack(2),irec(2),ifst(2) common /cstene/estack(2) common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /ceshow/ e1,e2,egamma,eran common /ccegs4/ ioegs4,eecrit,icesrc common /cmshow/ i0show,ishow,Xfirst,Xfel logical lproj common /clshow/lproj,ngener,ecut2,iogmag common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /cbase2/hinj,haxmax,hax0,haxc,dinj common /cbaseb/ theta1,theta2,tran,phi1,phi2 common /ccore/xcore,ycore,hcore,tcore common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis common /cprim/ hprim,xprim,yprim,tprim,ioprim common /cthin/ ilong,iothin,thinH,thinE $ ,ethin,hthin,wtmaxH,wtmaxE,rthmax,gthinH,gthinE common /ccask/e0cmin,ecmin,eclow,ecmax,fcask,iocask,thin2,athin2 $ ,iscmin common /chaint/iohadr,ifrmod,ihadt character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan COMMON/RANDOM/IXX double precision seed0,seed1,seed2,seed3 common /ccsee/ seed0,seed1,seed2,seed3 common /ccsee2/ Seed2Engy,fSeed2Engy,Stack2Engy,fStack2Engy data imxstk/mxstk/ parameter (mxz=3000,mxie=151) common /cslant/ mz common /elost/ elost(mxz) COMMON /ccel0/ eelow,eemin,eemax,fecask,isemin,themin COMMON /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask common /cnecsk/ necask logical kill integer iutime(5) double precision starttime,atime common /cptlty/ imode,icoll,elowhi integer*8 nptlc1,nptlc2,nptlcm,nptle1,nptle2,nptlem,nptlog $ ,nptllc,nptlle COMMON /CNPTLC/ nptlc1,nptlc2,nptlcm,nptle1,nptle2,nptlem,nptlog $ ,nptllc,nptlle character lzeros*10 data lzeros/"0000000000"/ common /cnptcl/nptcl common /cddsum/ddsum parameter (nmlast=1000) common /clast/ zlast(0:nmlast),elast(0:nmlast),idlast(0:nmlast) $ ,modlast(0:nmlast) character arg*80 common /ccread2/ ifread,echo logical echo common /curve1/rearth c.....initialize general global variables call iconex1 c.....read user instructions 10 call cread call OpenStack() if(Stack2Engy.lt.e0.or.Stack2Engy.lt.e2.or.fStack2Engy.gt.0. $ .and.fStack2Engy.lt.1.0)then else if(ifst(2).ne.0)then close(unit=ifst(2),status='delete') endif endif c.....other user-dependend intializations call iconex2 c.....initialize high energy hadronic model call inithh if(iohadr.eq.1) call Add2ToCheckfile('! hadonic model is neXus&') if(iohadr.eq.2) $ call Add2ToCheckfile(' ! hadonic model is QGSJET98&') if(iohadr.eq.3) $ call Add2ToCheckfile(' ! hadonic model is QGSJET01&') if(iohadr.eq.4) $ call Add2ToCheckfile(' ! hadonic model is SIBYLL21&') #ifdef GEANT c.....initialize GEANT call CGEANT #else c.....initialize gheisha call gheini #endif c.....initialize EGS4 if (ioegs4.ne.0) call iegs4 c.....initialize cask if (iocask.ne.0.or.iohadr.ge.10) call icask c.....initialize ecask if (iecask.ne.0) call iemc write(*,'("e0:",3(1pg8.2,2x))') e0 write(*,'("costet,phi:",3(1pg10.4,2x))') costet,aphi write(*,'("observation level at ",1pg13.6," meters")') hobs write(*,'("nshow:",i5,", seed: ",1d26.15)') nshow,seed0 c..... if(icdsrc.eq.2) call casksourceVerif(0) print *,"rhoair",rhoair(height(1023.0))/100. print *,rhoair(1000.0),depth(1000.0),rhoair(height(41.6930034)) iaxis=0 c.....calculate showers.................................... do ishow=1+i0show,nshow+i0show if(seed1.ne.0d0) then if(seed2.eq.0d0) call ranfgt(seed2) print *,"activating seed1",seed1 call ranfst(seed1) endif if(imodan.ne.0.and.mod(ishow,imodan).eq.0) then Xfirst=0. Xfel=0. endif call push xprim=dptl(6) yprim=dptl(7) hprim=dptl(8) tprim=dptl(9) etot=estack(1)+estack(2) hthin=thinH*(etot/e1)**gthinH*etot ethin=thinE*(etot/e1)**gthinE*etot*1e3 if(gthinH.ne.0..and.thinH.ne.0..and.mod(ishow,mshow).eq.0)then print *,"Energy dependent had. thinning level:",hthin/etot endif if(gthinE.ne.0..and.thinE.ne.0..and.mod(ishow,mshow).eq.0)then print *,"Energy dependent em. thinning level:",ethin/etot/1e3 endif ffcask=fcask if(fStack2Engy.gt.0.) Stack2Engy=etot*fStack2Engy if(fSeed2Engy .gt.0.) Seed2Engy =etot*fSeed2Engy if(iaxis.eq.0)then call initaxis hax0=0. hax0=-haxis(xprim,yprim,hprim) print *,"be zero",haxis(xprim,yprim,hprim) haxc=haxis(xcore,ycore,hcore) print *,"haxc:",haxc r=sqrt(dptl(6)**2+dptl(7)**2) h=dptl(8) atheta=acos(costet) call curini(atheta,aphi,r,h) iaxis=1 endif #ifdef USER_ANALYSIS call cini() #endif igo=0 ncask=0 necask=0 icesrc=0 iccsrc=0 if(iecask.ne.0)eecrit=min(eemax,e0*fecask)*1e3 ifi0=ifirst nptcl=0 ddsum=0.0 call stget(igo) do while (igo.eq.0) if(dptl(4).lt.Seed2Engy)then if(seed2.ne.0d0)then print *,"activating seed2",seed2 call ranfst(seed2) seed2=0d0 endif endif id=int(dptl(10)) if( iocask.ne.0 .and. e0.gt.e0cmin .and. dptl(4).le.ecmax $ .and. dptl(4).le.ffcask*e0 $ .and. dptl(4).gt.ecmin*iscmin.and.iccsrc.eq.0 $ .and. itype(id).ne.0 $ .and.gettheta(dptl(1),dptl(2),dptl(3)).lt.themin $ .and.ddepth().lt.zcemax $ ) then call caskip ncask=ncask+1 elseif(iecask.ne.0.and.(id.eq.10.or.abs(id).eq.12) $ .and.dptl(4).le.eemax.and.dptl(4).le.e0*fecask $ .and.(dptl(4).gt.eemin*isemin.and.icesrc.eq.0) $ .and.gettheta(dptl(1),dptl(2),dptl(3)).lt.themin $ .and.ddepth().lt.zcemax $ )then call emcip necask=necask+1 else if(icdsrc.eq.2.and.dptl(4).le.ecmin) call casksourceVerif(1) call cprop ida=iabs(id) zz=ddepth() if(dptl(6).eq.0..and.dptl(7).eq.0..and.dptl(8).eq.0.)then print *,"nuller" endif if( dptl(4) .gt. 0.9*e0 ) then if (nptcl.eq.0) then X0first=ddepth() X0fel=ddepth() else if (ddepth().gt.X0fel) X0fel=ddepth() endif endif if(imode.eq.3.or.imode.eq.4.or.imode.eq.6)then if (imode.eq.3) call decay if (imode.eq.4.or.imode.eq.6) call ccoll igener=dptl(12) if(igener.le.nmlast)then elast(igener)=dptl(4) zlast(igener)=ddepth() idlast(igener)=dptl(10) modlast(igener)=imode endif dptl(12)=dptl(12)+1 if(imode.ne.1) call ccuts endif if(imode.gt.1) call cptl2stack endif if(imode.ne.1) then call stget(igo) ddsum=0.0 nptcl=nptcl+1 endif if(igo.eq.1.and.ncask.ge.1.and.(iocask.eq.1 $ .or.(iocask.eq.2.and.ishow.eq.nshow)))then call cask call caskana call casksource if (ioegs4.eq.0.and.iecask.eq.0) then call casknkg else call caskemc endif if (iecask.eq.1) necask=necask+1 call stget(igo) ncask=0 iccsrc=1 endif if(igo.eq.1.and.necask.ne.0.and.(iecask.eq.1 $ .or.(iecask.eq.2.and.ishow.eq.nshow)))then call emc call emcsource(1) call emcana call stget(igo) necask=0 icesrc=1 eecrit=0. endif nptlc1=nptlc1+1 if(mod(nptlc1,nptlcm).eq.0.and.icd.ge.1)then nptlc1=0 nptlc2=nptlc2+nptllc ii=int(log10(float(nptlcm/nptllc))) write (*,'(" h",i10,a,g13.6)') $ nptlc2,lzeros(11-ii:10),estack(1) inquire(file='conex.icd',exist=kill) if(kill)then open(2,file="conex.icd") read (2,*)icd close(2) endif if(nptlog.eq.1.and.nptlc2.ne.1)then if(float(nint(log10(float(nptlc2)))) $ .eq.log10(float(nptlc2)))then nptlcm=nptlcm*10 nptllc=nptllc*10 endif endif endif enddo c one shower is finished here if(mod(ishow,mshow).eq.0)then call ranfgt(seed0) inquire(file='conex.stop.'//host,exist=kill) if(kill)stop'conex.stop host ' inquire(file='conex.stop',exist=kill) if(kill)stop'conex.stop ' write(*,'("shower:",1i6,1x,1i10,1x,1d26.15,1x,1i11)') $ ishow,nptcl,seed0,IXX ! call flush(6) endif Xfirst=Xfirst+X0first Xfel=Xfel+X0fel if(imodan.ne.0) then if(mod(ishow,imodan).eq.0) then call cfinal endif endif #ifdef USER_ANALYSIS call cfin() #endif call creset2 enddo write(*,*) 'all done',imodan if (icdsrc.eq.2) call casksourceVerif(2) if (imodan.eq.0) call cfinal goto 10 end c------------------------------------------------------------------------ subroutine OpenStack() c------------------------------------------------------------------------ character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan parameter (mxblk=12) parameter (mxstk = mxblk*50*2) common /cstack/stack(mxstk,2),istack(2),irec(2),ifst(2) common /cstene/estack(2) data imxstk/mxstk/ save iopen logical iopen data iopen/.false./ if(iopen)return c.....initialize stack 1 istack(1)=0 irec(1)=0 ifst(1)=53 estack(1)=0. if(index(ftmp,' ').eq.1)then open(unit=ifst(1),status='scratch', * form='unformatted',access='direct',recl=imxstk*8) else open(unit=ifst(1),file=ftmp(1:index(ftmp,' ')-1)//".stack1", * form='unformatted',access='direct',recl=imxstk*8) endif c.....stack 2 .... istack(2)=0 irec(2)=0 estack(2)=0. ifst(2)=54 if(index(ftmp,' ').eq.1)then open(unit=ifst(2),status='scratch', * form='unformatted',access='direct',recl=imxstk*8) else open(unit=ifst(2),file=ftmp(1:index(ftmp,' ')-1)//".stack2", * form='unformatted',access='direct',recl=imxstk*8) endif iopen=.true. end c------------------------------------------------------------------------ subroutine CloseStack() c------------------------------------------------------------------------ parameter (mxblk=12) parameter (mxstk = mxblk*50*2) common /cstack/stack(mxstk,2),istack(2),irec(2),ifst(2) common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 c.....close stack close(unit=ifst(1),status='delete') if(ifst(2).ne.0)then close(unit=ifst(2),status='delete') endif if(checkfile)close(52) end c---------------------------------------------------------------------- function zint(x1,y1,h1,x2,y2,h2) common/curve1/rearth c integrates density from point 1 to point b c for checking purposes c c---------------------------------------------------------------------- dl=sqrt((x2-x1)**2+(y2-y1)**2+(h2-h1)**2) r1=sqrt(x1*x1+y1*y1) ha1=sqrt((rearth+h1)**2+r1**2)-rearth ! altitude over sea level ddl=1.0 ddl=dl/int(dl/ddl) z=0.0 do a=0.0,dl+ddl/2,ddl x=x1+(x2-x1)/dl*a y=y1+(y2-y1)/dl*a h=h1+(h2-h1)/dl*a r=sqrt(x*x+y*y) ha=sqrt((rearth+h)**2+r**2)-rearth ! altitude over sea level if(a.gt.0)then z=z+(rhoair(ha)+rhoalt)/2*ddl endif rhoalt=rhoair(ha) enddo zint=z end c---------------------------------------------------------------------- subroutine curini(theta,phi,r,h) c---------------------------------------------------------------------- c c prepares data for thickness of atmosphere c akr,bkr = rho parameterization akr+X*bkr c c X=dX*(iz-1) c iz=int(X/dX)+1 c c ha(iz)=height over sea-level (=rearth) c hl(iz)=length along shower-axis from (xprim,yprim,hprim) c c 2002/11 Hans-Joachim Drescher c common/cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common/cbase2/hinj,haxmax,hax0,haxc,dinj common/curdef/arho(20000),brho(20000),ha(20000),hl(20000) common/curve1/rearth common /cslant/ mz COMMON /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax Xmax=37000.0 pz= cos(theta) pr=-sin(theta) p=1.0 dX=2.5 h1=h r1=r #ifdef CURVED ha1=sqrt((rearth+h1)**2+r1**2)-rearth #else ha1=h1 #endif rho1=rhoair(ha1) print *,"rho1",rho1 X1=0 iz=0 h0=h1 ! save initial coordinates r0=r1 ! dx1=-cos(phi)*r1 dy1=-sin(phi)*r1 dh1=h1 htop=500000.0 ii0=100 ii=ii0 print *,"phi",phi do while(ha1.gt.hground.and.ha1.lt.htop.and.X1.lt.Xmax) ii=ii-1 if(ii.eq.0)then iz=iz+1 endif 110 ha1c=sqrt((rearth+h1)**2+r1**2)-rearth #ifdef CURVED sina=r1/(rearth+ha1) cosa=(rearth+h1)/(rearth+ha1) pza= cosa*pz-sina*pr pra= sina*pz+cosa*pr cta=pza/sqrt(pza**2+pra**2) #else cta=pz/sqrt(pz**2+pr**2) #endif dXv=dX/ii0*cta ha2=height(depth(ha1)+dXv) dha=(ha1-ha2) if(dha.lt.1.0)then dl= dX/ii0 / rhoair(ha1) else dl= dha / cta endif h2=h1-dl*pz/p r2=r1+dl*pr/p X2=X1+dX/ii0 c check ha2=sqrt((rearth+h2)**2+r2**2)-rearth if(ii.eq.ii0-1)then dha1=ha1 dx1=-cos(phi)*r1 dy1=-sin(phi)*r1 dh1=h1 c call toaxis(dx1,dy1,dh1) endif c if(ii.eq.0)then cc print *,"ZINT:",X2,x1,y1,h1,x2,y2,h2,ha1,ha2 c endif c if(ii.eq.0)then ha(iz)=dha1 hl(iz)=haxis(dx1,dy1,dh1) cc PRINT*,X1,iz,hl(iz),dx1,dy1,dh1 dx2=-cos(phi)*r2 dy2=-sin(phi)*r2 dh2=h2 cc print *,"ZINT:",X2,dx1,dy1,dh1,dx2,dy2,dh2 cc $ ," dz:",zint(dx1,dy1,dh1,dx2,dy2,dh2) endif if(ii.eq.0)then rho2=rhoair(ha2) print *,iz,rho1,rho2,X1,X2,X2-dX,dha1,rhoair(dha1) brho(iz)=(rho2-rho1)/dX arho(iz)=rho1-brho(iz)*(X2-dX) rho1=rho2 c brho(iz)=X1/rho1 12 format("-->",3(1x,g12.6),"--",3(1x,g12.6),"--",6(1x,g12.6)) c if(iz.ge.2) c $ write (*,12) h1,r1,ha1, h2,r2,ha2, X1,hl(iz),hl(iz-1),dl endif h1=h2 r1=r2 X1=X2 ha1=ha2 if(ii.eq.0)then ii0=1 if(X2.lt.100)ii0=10 ii=ii0 endif enddo print *,"zmax:",X1,iz mz=iz izmax=min(int(zcemax/2.5),mz) iezmax=izmax zmax=X1 dZZ=2.5 c return do iz=1,mz X=dX*(iz-1) write (51,'(20(g12.6,1x))') $ X,ha(iz),hl(iz),arho(iz),brho(iz),height(X*costet) enddo do X=0.0,1000.0-5,0.1 iz=int(X/2.5)+1 dh=ha(iz)+mod(X,2.5)*(ha(iz+1)-ha(iz))/2.5 rho=arho(iz)+brho(iz)*X write(50,*) X,rho,rhoair(dh),dh enddo return end c----------------------------------------------------------------------------- subroutine si c----------------------------------------------------------------------------- c h.drescher c implicit double precision (a-h,o-z) real airz,aira,airw,airavz,airava,sighh common /compair/airz(3),aira(3),airw(3),airavz,airava COMMON /SIGM/ SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO DOUBLE PRECISION SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO do ee=2,12,0.5 e=10.0**ee cn=14.54*1660.44 / rlam(1,e,0d0) cp=14.54*1660.44 / rlam(2,e,0d0) ck=14.54*1660.44 / rlam(3,e,0d0) s1=0. s2=0. s3=0. do i=1,3 s1=s1+sighh(sngl(e),1,aira(i))*airw(i) s2=s2+sighh(sngl(e),2,aira(i))*airw(i) s3=s3+sighh(sngl(e),3,aira(i))*airw(i) enddo write (91,'(17(g12.6,1x))') e,s1,s2,s3,cn,cp,ck enddo stop end c----------------------------------------------------------------------- integer function itype(id) c----------------------------------------------------------------------- c returns type for a particle dptl which is good for cask c 0 otherwise c h.drescher c c ida=abs(id) if(ida.eq.1120.or.ida.eq.1220)then itype=1 elseif(ida.eq.120)then itype=2 elseif(ida.eq.130)then itype=3 elseif(id.eq.-20) then itype=4 elseif(id.eq. 20) then itype=5 c elseif(id.eq. 110) then c itype=6 c elseif(id.eq. 10) then c itype=7 else itype=0 endif end c----------------------------------------------------------------------- subroutine iconex1 c----------------------------------------------------------------------- c h.drescher c COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax parameter (mxblk=12) common /cdebug/icd,icdsrc common /compair/airz(3),aira(3),airw(3),airavz,airava common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cecut/ ehcut,eecut,emcut,epcut logical lproj common /ceshow/e1,e2,egamma,eran common /clshow/lproj,ngener,ecut2,iogmag common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /cbase2/hinj,haxmax,hax0,haxc,dinj common /cbaseb/ theta1,theta2,tran,phi1,phi2 common /ccegs4/ioegs4,eecrit,icesrc double precision seed0,seed1,seed2,seed3 common /ccsee/ seed0,seed1,seed2,seed3 common /ccsee2/ Seed2Engy,fSeed2Engy,Stack2Engy,fStack2Engy common /chaint/iohadr,ifrmod,ihadt common /ccask/e0cmin,ecmin,eclow,ecmax,fcask,iocask,thin2,athin2 $ ,iscmin common /ccask2/thin3,athin3 common /ccel8/thin4,athin4 common /cthin/ ilong,iothin,thinH,thinE $ ,ethin,hthin,wtmaxH,wtmaxE,rthmax,gthinH,gthinE common /cmshow/ i0show,ishow,Xfirst,Xfel common /ccel0/ eelow,eemin,eemax,fecask,isemin,themin common /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask integer*8 nptlc1,nptlc2,nptlcm,nptle1,nptle2,nptlem,nptlog $ ,nptllc,nptlle COMMON /CNPTLC/ nptlc1,nptlc2,nptlcm,nptle1,nptle2,nptlem,nptlog $ ,nptllc,nptlle c.....data character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan parameter (mxz=3000,mxie=151) common /londat/ hz(6,mxie,mxz),del(mxz),dpo(mxz),dmuon(mxz) common /elost/ elost(mxz) parameter (mxla=300,mrla=60) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis common /latda3/ dlat(4,mxla,mxla), $ xlmin,xlmax,ylmin,ylmax,nxlbin,nylbin common /alldat3/ SaveGrndPtcls,SaveLongPtcls logical SaveGrndPtcls,SaveLongPtcls common /cptlty/ imode,icoll,elowhi common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac common /cgeom/ zbound common /cprim/ hprim,xprim,yprim,tprim,ioprim common /ccore/xcore,ycore,hcore,tcore common/curve1/rearth common /ccana2/ dzana #ifdef NEXUS2 character*80 nexhome common /ccnex/nexhome #endif #ifdef GEANT c geant stuff common/cihad/IHAD,egcalor #endif #ifdef DETECTORS c.....detectors #include "detectors.inc" #endif #ifdef DETRESP #include "detresp.inc" #endif c............... common /cgetxor/ getxor,stsort logical getxor common/cfakesig/ fakesig common /cgsig/ iosign common /CINEX/inex(50) !inex(i)=nexus-id for geant-id (i) common /QNEX/iqnex(-10:10) common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) data inex/10,-12,12,11,-14,14,110,120,-120,-20 $ ,130,-130,1220,1120,-1120,20,220,2130,1130,1230 $ ,2230,1330,2330,3331,-1220,-2130,-2230,-1230,-2230,-1330 $ ,-2330,-3331,-16,16,240,-240,140,-140,340,-340 $ ,2140,42,43,44,45,46,47,48,0,0/ data iqnex/ * 0, 0, 0, 0, -2130, -20, -130, -1220, -1120, -120, * 110, * 120, 1120, 1220, 130, 20, 2130, 0, 0, 0, 220/ DATA HATM / -6e3 , 4e3 , 1e4 ,4e4, 1e5 / DATA AATM / -186.5562, -94.9199, 0.61289,0e0,.01128292 / DATA BATM / 1222.6562,1144.9069,1305.5948,540.1778,1.0 / DATA CATM / 9941.8638,8781.5355,6361.4304,7721.7016,1e7/ common /cgeoB/ geoBx,geoBy,geoBz,geoB c.....data files........................... character fload*80 common/cload/nload,fload(10) #ifdef QGSJET98 parameter (nloadi=3) data nload/nloadi/ data (fload(i),i=1,nloadi)/ c $ "ppjghe.dat" $ "ppjqgs.dat" $ ,"ppk.dat" $ ,"ppem3.dat" $ / #endif #ifdef NEXUS2 parameter (nloadi=3) data nload/nloadi/ data (fload(i),i=1,nloadi)/ c $ "ppjghe.dat" $ "ppjnxs.dat" $ ,"ppk.dat" $ ,"ppem3.dat" $ / #endif common/ccread2/ ifread,echo logical echo common /cupkill/ upkill logical upkill c.....basic variables ifread=5 ! std input echo=.true. i0show=0 ! shower-number offset ioegs4=1 ! EM-shower EGS4 (1) or NGK (0) eecrit=0 ! ipmax=1 ! absolute highest energy index ipcmax=1 ! highest energy index in cask izmin=mxz icd=1 ! debug - output (as ish) icdsrc=0 costet=1e0 ! cos(theta) theta-incident angle phi1=0. phi2=0. theta1=0. ! if theta2>theta1 theta2=0. ! ! randomly between theta1 and theta2 tran=2. ! aphi=0. ! incident azimuth angle emin=1e0 ! emax=1e12 ! nde=10 ! number of energy bins per decade ioprim=0 ! xprim=0 ! yprim=0 ! hprim=0 ! tprim=0 ! #ifdef CURVED rearth=6371.01e3 ! radius of earth #else rearth=0.0 #endif c.....primary energy e0=1e5 ! primary energy e1=1e5 ! if e2>e1 -> e0 is chosen randomly e2=1e5 ! according to egamma=1 ! dn/de = e^gamma eran=0.0 ! discretize e0 z0=.0 ! thinE=0e0 ! thinning fraction thinH=0e0 ! thinning fraction gthinE=0.0 ! thinning energy dependence gthinH=0.0 ! thinning energy dependence wtmaxH=1e30 ! max weight hadrons wtmaxE=1e30 ! max weight em c.....energy cut-offs ehcut= 0.05 ! hadrons kinetic energy cutoff emcut= 0.05 ! muons kinetic energy cutoff eecut= 1e-3 ! electrons kinetic energy cutoff epcut= 1e-3 ! photons kinetic energy cutoff c.....initial particle (chosen by random) id0(1)=1120 ! initial particle 1 do i=2,16 id0(i)=0 ! initial particle 2..9 if not 0 enddo nshow=1 ! number of showers mshow=1 ! show every mshow-th shower ifirst=4 ! force first interaction 4=collision 3=decay zmin=0 ! min depth zmax=1 ! max depth, will be adapted zcemax=3000.0 ! hground=100.0 ! hground=height(zmax) hobs=100.0 ! observation level in meters hcore=hobs xcore=0. ycore=0. tcore=0. hinj=100000. ! injection height dinj=0.0 ! distance of injection from hcore haxmax=0.0 dz=2.5 ! step size in slant-depth dzana=10.0 ! for analysis inexan=0 icasan=2 ! analysis mode 1=least 10=most imodan=0 ! if not 0 every imodan-th shower is analysed iohadr=0 ! choice for hadronic interaction model ! is set by initialization ! ifrmod=2 ! fragmentation of target (only QGS-Jet so far) ! 1=single nucleons ! 2=realistic fragmentation ! 3=single nucleus lproj=.false. ! project particles below ehcut to ground ecut2=0. ! iogmag=0 ! consider geomagnetic field ngener=-1 ! maximal number of generations (-1=off) ilong=0 ! if 1, no lateral developpment is done ! (for comparing shower and cask-mode) iothin=1 ! thinning-modus Seed2Engy=0 ! threshold for activating seed2 fSeed2Engy=0 ! Stack2Engy=1e30 ! upper threshold for stack 2 fStack2Engy=0 ! seed0=0d0 ! seed for begin of run seed1=0d0 ! seed for begin of each event seed2=0d0 ! seed for first particle below Seed2Engy ! to use in combination with Stack2Engy seed3=1234 ! seed for EGS4 upkill=.false. c.....el-ph cascade equations zbound=0.0 eelow=1e-3 iecask=1 fecask=1e-3 eemax=1e30 eemin=10 isemin=0 ! put lower energies into source function themin=0.1 call creset2 elowhi=89.125093 ! threshold for low-energy hadronic #ifdef GEANT ihad=1 #endif c.....cask fcask=1e-3 ! energy-fraction below which cask-mode is done ecmax=1e30 ! max energy for cask-algorithm ecmin=1e4 ! min energy iscmin=0 ! 0=put lower energies into source function eclow=1e0 ! lower limit source function e0cmin=0e0 ! min primary energy to use cask-mode iocask=1 ! cask modus ! = 0 : off ! = 1 : do all particles of one shower at once ! = 2 : do all particles of all showers at once ! thin2=1.0 ! fraction of energy from source-function athin2=1e6 ! maximal energy in source-function thin3=1.0 ! fraction of energy from source-function pi0 athin3=1e30 ! maximal energy in source-function pi0 thin4=1.0 ! fraction of energy from source-function elph athin4=1e5 ! maximal energy in source-function elph c.....lateral data-sampling SaveGrndPtcls=.false. ! SaveLongPtcls=.false. ! xlmin=-3000 ! analyse in grid xlmax= 3000 ! ylmin=-3000 ! ylmax= 3000 ! nxlbin=300 ! number of bins nylbin=300 ! rlmin=10.0**1.5 ! analyse as radial function rlmax=10.0**3.5 ! nrlbin=20 ! number of bins home='data/ ' ! where are all the dat-files #ifdef NEXUS2 nexhome=home #endif host=' ' fbase='casc ' ! longitudinal data ftmp=' ' ! iraxis=1 ! c.....initialize air aira(1)=14. airz(1)=7. airw(1)=0.781 aira(2)=16. airz(2)=8. airw(2)=.21 aira(3)=40. airz(3)=18. airw(3)=0.009 airava=aira(1)*airw(1)+aira(2)*airw(2)+aira(3)*airw(3) airavz=airz(1)*airw(1)+airz(2)*airw(2)+airz(3)*airw(3) geoBx= 20.0e-6 * 0.299792458 geoBy= 0. geoBz=-40.0e-6 * 0.299792458 geoB=sqrt(geoBx**2+geoBy**2+geoBz**2) #ifdef DETECTORS c..... idetmod=0 ndet=0 zdetmin=1e30 zdetmax=0 xdetshift=0 ydetshift=0 detangle=-1 detdrad=1.05 rdetmax=1e4 hdetfac=1 ! adetabs=0.0 c ahdetabs=0.0 tdetpoiss=0.0 t1res=500.0 ! intdet=.false. detresc=.false. c tdetmin=-40012.5 ! nano seconds c tdetmax= 40012.5 ! c ntdetbin=3201 ! dtdet=25.0 ntdetbin=400 iadet=0111 tshiftdet=.false. c surfdet=.true. #endif #ifdef DETRESP depemn=0 !7.9432823e-05 ! incident energy range depemx=0 !125.89254 ndepe =0 !31 depxmn=0 !8.9125094e-05 ! output value range depxmx=0 !112.20184 ndepx =0 !61 iodepg=.true. !.true. ! true/false is for log depamn=0 !0.0 ! angle range in sec or rad depamx=0 !1.5707963267949 ndepa =0 !18 iodepa=.false. !.false. ! true if angle is sec fdepa=1.0 ! detver='AGASA ' aratio=0.0 ares=0.0 vip=0.0110 #endif iosign=0 fakesig=0.0 call deletearrays getxor=.false. stsort=0.0 call cxiniall print *,"###############################################" print *,"# SENECA 1.0 beta #" print *,"# main reference: #" print *,"# H.J.Drescher and G.Farrar astro/ph-0212018 #" print *,"###############################################" end c----------------------------------------------------------------------- subroutine iconex2 c----------------------------------------------------------------------- c h.drescher c PARAMETER (TWOPI=6.28318530717958648,PI=3.14159265358979324) parameter (mxblk=12) common /cdebug/icd,icdsrc common /compair/airz(3),aira(3),airw(3),airavz,airava common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground parameter (mxz=3000,mxie=151) common /cslant/ mz double precision seed0,seed1,seed2,seed3 common /ccsee/ seed0,seed1,seed2,seed3 common /ccsee2/ Seed2Engy,fSeed2Engy,Stack2Engy,fStack2Engy common /cmshow/ i0show,ishow,Xfirst,Xfel common /ccask/e0cmin,ecmin,eclow,ecmax,fcask,iocask,thin2,athin2 $ ,iscmin common /ccask2/thin3,athin3 common /clshow/lproj,ngener,ecut2,iogmag common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) c.....data character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan character*80 fname common /londat/ hz(6,mxie,mxz),del(mxz),dpo(mxz),dmuon(mxz) parameter (mxla=300,mrla=60) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis common /latda3/ dlat(4,mxla,mxla), $ xlmin,xlmax,ylmin,ylmax,nxlbin,nylbin common /elost/ elost(mxz) c.................................... common/cjprob/Q(10),JPROB,IMUMOD PARAMETER (EMMU=0.105658389) PARAMETER (EMASS=0.00051099906) c.................................... include 'cptl.inc' common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis #ifdef DETECTORS c.....detectors #include "detectors.inc" #endif c.....agasa detectors.... #ifdef DETRESP #include "detresp.inc" logical da #endif character filename*120 call cdecin(.false.) if(i0show.le.-1) call i0ini() cc mz=int((zmax/costet-zmin)/dz)+1 if(imodan.eq.0)then oosho=1e0/float(nshow) else oosho=1e0/float(imodan) endif if(iogmag.ge.2) call SETIAUSFL(6) nptl=0 c.....some muon stuff for GEANT SE=SQRT(2.71828) do i=1,3 C1=airz(i)**0.333333 JPROB=i*3 Q(JPROB-2)=LOG(189.*EMMU/(EMASS*C1)) IF(airz(i).GT.10)Q(JPROB-2)=Q(JPROB-2)+LOG(0.666666/C1) Q(JPROB-1)=189.*SE*EMMU*EMMU/(2.*EMASS*C1) Q(JPROB )=0.75*SE*EMMU*C1 enddo #ifdef DETRESP c.....check for tables..... ! detver now defined in cread h.drescher ! c=normal b=BIRKS law included ! th = with theta dependence filename=home(1:index(home,' ')-1) $ //'/edep_'//detver(1:index(detver,' ')-1)//'.conf' inquire(file=filename,exist=da) if(da)then open(2,file=filename,status='old') read(2,*) depemn,depemx,ndepe print *,"detector - response file: " write(*,*) depemn,depemx,ndepe if(ndepx.eq.0)then read(2,*) depxmn,depxmx,ndepx,iodepg else read(2,*) adummy,adummy,idummy write(*,*) "override by pre-defined values:" endif write(*,*) depxmn,depxmx,ndepx,iodepg if(ndepa.eq.0)then read(2,*) depamn,depamx,ndepa,iodepa else write(*,*) "override by pre-defined values:" read(2,*) adummy,adummy,idummy endif write(*,*) depamn,depamx,ndepa,iodepa read(2,*) fdepa write(*,*) fdepa close(2) else print *,filename stop 'no detector response configuration file !!!!!!!!!!' endif do ld=1,6 if(ld.eq.1) filename=home(1:index(home,' ')-1) $ //'/pho_'//detver(1:index(detver,' ')-1)//'.dat' if(ld.eq.2) filename=home(1:index(home,' ')-1) $ //'/elec_'//detver(1:index(detver,' ')-1)//'.dat' if(ld.eq.3) filename=home(1:index(home,' ')-1) $ //'/pos_'//detver(1:index(detver,' ')-1)//'.dat' if(ld.eq.4) filename=home(1:index(home,' ')-1) $ //'/muon_'//detver(1:index(detver,' ')-1)//'.dat' if(ld.eq.5) filename=home(1:index(home,' ')-1) $ //'/neu_'//detver(1:index(detver,' ')-1)//'.dat' if(ld.eq.6) filename=home(1:index(home,' ')-1) $ //'/pro_'//detver(1:index(detver,' ')-1)//'.dat' inquire(file=filename,exist=da) if(da)then print *,"reading from file:",filename open(2,file=filename,status='old') c depneu(ix,iy,is,ld) = energy deposit probability c ix = 1 .. 31 c = primary particle kinetic energy ey=10.0**( (float(ix)-0.5)*(0.2)-4.10) c iy = 1 .. 61 c = energy deposit ex=10.0**( (float(ix)-0.5)*(0.1)-4.05) c ix=0 means undeflow below lowest bin (including zero) c is = 1 .. 6 c = secans angle sec=1.0+(is-1)*0.2 c do is=1,ndepa do ix=1,ndepx c print *,is,ix,ndepe read (2,*) (depneu(iy,ix,is,ld),iy=1,ndepe) enddo do iy=1,ndepe c if(depneu(iy,0,is,ld).eq.0.) then sum=0.0 do ix=1,ndepx sum=sum+depneu(iy,ix,is,ld) enddo c print *,"setting",is,iy,ld," to ",1.0-sum depneu(iy,0,is,ld)=1.-sum c endif enddo enddo close(2) else print *,"NO such file !!!!!!!!!!! ",filename do is=1,ndepa do iy=1,ndepe depneu(iy,0,is,ld)=1. enddo enddo endif enddo #endif #ifdef DETECTORS if(mod(iadet/1000,10).ne.0.and.ndet.gt.0) then idetout=mod(iadet/1000,10) fname=fbase(1:index(fbase,' ')) write (fname(index(fbase,' '):),'("-det4.dat ")') call chkname(fname,fname) open(86,file=fname(1:index(fname,' ')-1)) write (86,*) "# particle list for detectors " write (86,*) "# columns:" write (86,*) "# 1: number of detector" write (86,*) "# 2: particle id" write (86,*) "# 3: arrival time" write (86,*) "# 4: weight" write (86,*) "# 5: px" write (86,*) "# 6: py" write (86,*) "# 7: pz" write (86,*) "# 8: kinetic energy" #ifdef DETRESP write (86,*) "# 9: energy deposit" #endif cc write(86,186) i,id,ttt,wtnew,px,py,pz,ek,edep else idetout=0 endif #endif do il=1,5 DATM(il)=depth(HATM(il)) print *,"layer",il,HATM(il),DATM(il) enddo return end c----------------------------------------------------------------------- subroutine initaxis c----------------------------------------------------------------------- c h.drescher c PARAMETER (TWOPI=6.28318530717958648,PI=3.14159265358979324) common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /ccore/xcore,ycore,hcore,tcore COMMON /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask parameter (mxla=300,mrla=60) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis c.....detectors...................... parameter (mxz=3000,mxie=151) #ifdef DETECTORS #include "detectors.inc" #endif c.................................... common /cslant/ mz common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 atheta=acos(costet) xsaxis=-cos(aphi)*sin(atheta) ysaxis=-sin(aphi)*sin(atheta) zsaxis=costet print *,"shower-angles:",atheta,aphi,sin(atheta),sqrt(1-costet**2) c predefine rotation to shower-axis-system aphi2=aphi tsaxis(1)= cos(atheta) * cos(-aphi2) tsaxis(2)= cos(atheta) *(-sin(-aphi2)) tsaxis(3)= sin(atheta) tsaxis(4)= sin(-aphi2) tsaxis(5)= cos(-aphi2) tsaxis(6)= 0. tsaxis(7)=-sin(atheta) * cos(-aphi2) tsaxis(8)=-sin(atheta) *(-sin(-aphi2)) tsaxis(9)= cos(atheta) c predefine rotation back to normal system fsaxis(1)= cos(aphi2)*cos(-atheta) fsaxis(2)=-sin(aphi2) fsaxis(3)= cos(aphi2)*sin(-atheta) fsaxis(4)= sin(aphi2)*cos(-atheta) fsaxis(5)= cos(aphi2) fsaxis(6)= sin(aphi2)*sin(-atheta) fsaxis(7)= -sin(-atheta) fsaxis(8)= 0. fsaxis(9)= cos(-atheta) c check that the product is unity do j=1,3 write (*,'(" ( ",$)') do i=1,3 c=0. do k=1,3 i1=(i-1)*3+k-1+1 i2=(k-1)*3+j-1+1 c=c+tsaxis(i1)*fsaxis(i2) enddo write (*,'(f6.4," ",$)') c enddo write (*,'(" )")') enddo xx=xsaxis yy=ysaxis zz=zsaxis print *,xx,yy,zz call toaxis(xx,yy,zz) print *,xx,yy,zz #ifdef DETECTORS c.....detectors...................................................... if(checkfile)then write (52,*) 'Detectors' endif if(idetmod.eq.2)then open(91,file='det.check') write(91,*) "GroundDetectors" if(detresc)then print *,"detector radius rescaling" endif do i=1,ndet rr=sqrt(xdet(i)**2+ydet(i)**2) if(detresc)then xr=xdet(i) yr=ydet(i) hr=0.0 call toaxis(xr,yr,hr) f=sqrt(xr**2+yr**2)/rr xdet(i)=xdet(i)/f ydet(i)=ydet(i)/f endif xr=xdet(i) yr=ydet(i) hr=0.0 if(iraxis.eq.1) call toaxis(xr,yr,hr) rr=sqrt(xr**2+yr**2) r1det(i)=rr/(detdrad) r2det(i)=rr*(detdrad) hdet(i)=hobs h1det(i)=hobs h2det(i)=hobs aSdet(i)=sqrt( (r1det(i)-r2det(i))**2 $ + (h1det(i)-h2det(i))**2 )*(r1det(i)+r2det(i))/2.*twopi if(iraxis.eq.1) then aSdet(i)=aSdet(i)/costet endif if(detangle.ge.0.) aSdet(i)=aSdet(i)*detangle/twopi if(adet(i).le.0.) adet(i)=aSdet(i) angdet(i)=atan2(yr,xr) dangdet(i)=detangle/2 if(r1det(i).gt.rdetmax)then ! switch off actdet(i)=.false. else actdet(i)=.true. hdetmin=hobs hdetmax=hobs endif it0det(i)=0 tdetmn(i)=1e10 if(checkfile)then write (52,*) detname(i),xdet(i),ydet(i),hdet(i),adet(i) endif write (91,9) actdet(i),r1det(i),r2det(i),h1det(i),h2det(i) $ ,aSdet(i),adet(i),aSdet(i)/adet(i) enddo close(91) endif 9 format(L1,1x,30(g13.6,1x)) if(checkfile)then write (52,*) 'endDetectors' endif c if(zmax.lt.depth(hdetmin))then c hground=hdetmin c write (*,*) 'adapting hground to lowest detector height:' c $ ,hground,"g/cm^2" c endif #endif !DETECTORS c.................................... end c------------------------------------------------------------------------ function haxis(x,y,z) common /cbase2/hinj,haxmax,hax0,haxc,dinj c h.drescher c common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis haxis=tsaxis(7)*x+tsaxis(8)*y+tsaxis(9)*z-hax0 haxis=-haxis end c------------------------------------------------------------------------ subroutine toaxis(x,y,z) c h.drescher c common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis xx=tsaxis(1)*x+tsaxis(2)*y+tsaxis(3)*z yy=tsaxis(4)*x+tsaxis(5)*y+tsaxis(6)*z zz=tsaxis(7)*x+tsaxis(8)*y+tsaxis(9)*z x=xx y=yy z=zz end subroutine fromaxis(x,y,z) c h.drescher c common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis xx=fsaxis(1)*x+fsaxis(2)*y+fsaxis(3)*z yy=fsaxis(4)*x+fsaxis(5)*y+fsaxis(6)*z zz=fsaxis(7)*x+fsaxis(8)*y+fsaxis(9)*z x=xx y=yy z=zz end function gettheta(px,py,pz) c h.drescher c x=-px y=-py z=pz call toaxis(x,y,z) gettheta=atan2(sqrt(x*x+y*y),z) end c----------------------------------------------------------------------- subroutine cread c----------------------------------------------------------------------- c h.drescher c parameter (mxblk=12) common /cdebug/icd,icdsrc common /compair/airz(3),aira(3),airw(3),airavz,airava common /ccegs4/ioegs4,eecrit,icesrc common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /ccore/xcore,ycore,hcore,tcore common /cbase2/hinj,haxmax,hax0,haxc,dinj common /cbaseb/ theta1,theta2,tran,phi1,phi2 COMMON /PRNTFL/INBCD,NEWBCD,INBIN,NEWBIN,NPEVT,NEVTP,LPRT,NPRT(10) LOGICAL LPRT,NPRT common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cecut/ ehcut,eecut,emcut,epcut common /ceshow/e1,e2,egamma,eran logical lproj common /cgeoB/ geoBx,geoBy,geoBz,geoB common /clshow/lproj,ngener,ecut2,iogmag common /chaint/iohadr,ifrmod,ihadt common /ccask/e0cmin,ecmin,eclow,ecmax,fcask,iocask,thin2,athin2 $ ,iscmin common /ccask2/thin3,athin3 common /ccel8/thin4,athin4 common /cthin/ ilong,iothin,thinH,thinE $ ,ethin,hthin,wtmaxH,wtmaxE,rthmax,gthinH,gthinE common /cptlty/ imode,icoll,elowhi include 'cptl.inc' common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax double precision seed0,seed1,seed2,seed3 common /ccsee/ seed0,seed1,seed2,seed3 common /ccsee2/ Seed2Engy,fSeed2Engy,Stack2Engy,fStack2Engy character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan parameter (mxla=300,mrla=60) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis common /latda3/ dlat(4,mxla,mxla), $ xlmin,xlmax,ylmin,ylmax,nxlbin,nylbin common /cgetxor/ getxor,stsort logical getxor common /cmshow/ i0show,ishow,Xfirst,Xfel common /cthpar/ thpar(10) COMMON /ccel0/ eelow,eemin,eemax,fecask,isemin,themin COMMON /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask common /alldat3/ SaveGrndPtcls,SaveLongPtcls logical SaveGrndPtcls,SaveLongPtcls parameter (mxz=3000,mxie=151) common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac common /cgeom/ zbound common /ccana2/ dzana #ifdef NEXUS2 character*80 nexhome common /ccnex/nexhome #endif C.....detectors............... #ifdef DETECTORS #include "detectors.inc" #endif #ifdef DETRESP #include "detresp.inc" #endif c............... character fload*80 common/cload/nload,fload(10) common/cfakesig/ fakesig common /cgsig/ iosign common/ccread2/ ifread,echo logical haveword,checkfile,echo common/ccread/haveword,word,iunit common/ccread1/ checkfile,checkfilename character word*80,word2*80,checkfilename*80 double precision getvalue dimension ia(8) common /cupkill/ upkill logical upkill #ifdef GEANT c geant stuff common/cihad/IHAD,egcalor #endif #ifdef USER_VAR #include "user_var.inc" #endif checkfile=.false. haveword=.false. call getw('init ') do while(.true.) if(.not.haveword)call getw(word) haveword=.false. if(word.eq.'input')then ifread=ifread+1 if(ifread.eq.6)ifread=10 call getw(word) open(unit=ifread,file=word,status='old') call getw('init ') print *,"input:",ifread,word elseif(word.eq.'unload')then nload=0 elseif(word.eq.'load')then nload=nload+1 call getw(fload(nload)) elseif(word.eq.'set')then call getw(word) if(word.eq.'nload') then nload=nint(getvalue()) elseif(word.eq.'chkfile')then call getw(checkfilename) checkfile=.true. open(52,file=checkfilename(1:index(checkfilename,' ')-1)) elseif(word.eq.'i0show') then i0show=getvalue() elseif(word.eq.'Seed2Energy') then Seed2Engy=getvalue() elseif(word.eq.'Stack2Energy') then Stack2Engy=getvalue() elseif(word.eq.'Seed2Engy') then Seed2Engy=getvalue() elseif(word.eq.'fSeed2Engy') then fSeed2Engy=getvalue() elseif(word.eq.'Stack2Engy') then Stack2Engy=getvalue() elseif(word.eq.'fStack2Engy') then fStack2Engy=getvalue() elseif(word.eq.'seed'.or.word.eq.'seed0') then seed0 =getvalue() call ranfst(seed0) elseif(word.eq.'seed1') then seed1 =getvalue() elseif(word.eq.'seed2') then seed2=getvalue() elseif(word.eq.'seed3') then seed3=getvalue() elseif(word.eq.'nshow') then nshow =getvalue() elseif(word.eq.'mshow') then mshow =getvalue() elseif(word.eq.'ifirst') then ifirst=getvalue() elseif(word.eq.'zmin') then zmin=getvalue() elseif(word.eq.'zcemax') then zcemax=getvalue() elseif(word.eq.'zmax') then zmax=getvalue() hground=height(zmax) hobs=hground hcore=hobs elseif(word.eq.'hground')then hground=getvalue() hobs=hground hcore=hobs elseif(word.eq.'hgrnd') then hobs=getvalue() hcore=hobs elseif(word.eq.'hobs') then hobs=getvalue() hcore=hobs elseif(word.eq.'hcore') then hcore=getvalue() elseif(word.eq.'xcore') then xcore=getvalue() elseif(word.eq.'ycore') then ycore=getvalue() elseif(word.eq.'tcore') then tcore=getvalue() elseif(word.eq.'haxmax') then haxmax=getvalue() elseif(word.eq.'hinj') then hinj=getvalue() if(iunit.eq.2)then hinj=height(hinj) endif elseif(word.eq.'zinj') then zinj=getvalue() elseif(word.eq.'dinj') then dinj=getvalue() elseif(word.eq.'dzana') then dzana=getvalue() elseif(word.eq.'e0') then e0=getvalue() elseif(word.eq.'sqs') then sqs=getvalue() amproj=.93891897 amnuc =.93891897 c sqs=sqrt( 2.0*elab*amnuc+amnuc**2+amproj**2 ) e0=(sqs**2-amproj**2-amnuc**2)/2/amnuc elseif(word.eq.'e1') then e1=getvalue() elseif(word.eq.'e2') then e2=getvalue() elseif(word.eq.'egamma')then egamma=getvalue() elseif(word.eq.'eran') then eran=getvalue() elseif(word.eq.'z0') then z0=getvalue() elseif(word.eq.'id0') then id0(1)=getvalue() elseif(word.eq.'id1') then id0(2)=getvalue() elseif(word.eq.'id2') then id0(3)=getvalue() elseif(word.eq.'id3') then id0(4)=getvalue() elseif(word.eq.'id4') then id0(5)=getvalue() elseif(word.eq.'id5') then id0(6)=getvalue() elseif(word.eq.'id6') then id0(7)=getvalue() elseif(word.eq.'id7') then id0(8)=getvalue() elseif(word.eq.'id9') then id0(9)=getvalue() elseif(word.eq.'id10') then id0(10)=getvalue() elseif(word.eq.'id11') then id0(11)=getvalue() elseif(word.eq.'id12') then id0(12)=getvalue() elseif(word.eq.'id13') then id0(14)=getvalue() elseif(word.eq.'id14') then id0(15)=getvalue() elseif(word.eq.'id15') then id0(16)=getvalue() elseif(word.eq.'idproj') then id0(1)=getvalue() elseif(word.eq.'xprim') then xprim=getvalue() ioprim=1 elseif(word.eq.'yprim') then yprim=getvalue() ioprim=1 elseif(word.eq.'hprim') then hprim=getvalue() ioprim=1 elseif(word.eq.'tprim') then tprim=getvalue() ioprim=1 elseif(word.eq.'gthinH') then !-------------- gthinH=getvalue() elseif(word.eq.'gthinE') then gthinE=getvalue() elseif(word.eq.'gthin') then gthinH=getvalue() gthinE=gthinH elseif(word.eq.'thinH') then thinH=getvalue() elseif(word.eq.'thinE') then thinE=getvalue() elseif(word.eq.'thin') then thinH=getvalue() thinE=thinH elseif(word.eq.'wtmaxE') then wtmaxE=getvalue() elseif(word.eq.'wtmaxH') then wtmaxH=getvalue() elseif(word.eq.'wtmax') then wtmaxH=getvalue() wtmaxE=wtmaxH elseif(word.eq.'rthmax') then rthmax=getvalue() elseif(word.eq.'thin2'.or.word.eq.'fElowH') then thin2=getvalue() elseif(word.eq.'athin2'.or.word.eq.'ElowH') then athin2=getvalue() elseif(word.eq.'thin3') then thin3=getvalue() elseif(word.eq.'athin3') then athin3=getvalue() elseif(word.eq.'thin4'.or.word.eq.'fElowE') then thin4=getvalue() elseif(word.eq.'athin4'.or.word.eq.'ElowE') then athin4=getvalue() elseif(word.eq.'ehcut'.or.word.eq.'EcutH') then ehcut=getvalue() elseif(word.eq.'emcut'.or.word.eq.'EcutM') then emcut=getvalue() elseif(word.eq.'eecut'.or.word.eq.'EcutE') then eecut=getvalue() elseif(word.eq.'epcut'.or.word.eq.'EcutP') then epcut=getvalue() elseif(word.eq.'icasan') then icasan=getvalue() elseif(word.eq.'imodan') then imodan=getvalue() elseif(word.eq.'icd') then icd =getvalue() elseif(word.eq.'icdsrc') then icdsrc =getvalue() elseif(word.eq.'host') then call getw(host) elseif(word.eq.'home') then call getw(home) #ifdef NEXUS2 nexhome=home elseif(word.eq.'nexhome') then call getw(nexhome) #endif elseif(word.eq.'fbase') then call getw(fbase) elseif(word.eq.'ftmp') then call getw(ftmp) elseif(word.eq.'costet') then costet=getvalue() elseif(word.eq.'theta') then costet=cos(getvalue()) elseif(word.eq.'phi1') then phi1=getvalue() elseif(word.eq.'phi2') then phi2=getvalue() elseif(word.eq.'theta1') then theta1=getvalue() elseif(word.eq.'theta2') then theta2=getvalue() elseif(word.eq.'tran') then tran=getvalue() elseif(word.eq.'phi') then aphi=getvalue() elseif(word.eq.'fcask'.or.word.eq.'fEmaxH') then fcask=getvalue() elseif(word.eq.'fEmax') then fcask=getvalue() fecask=fcask elseif(word.eq.'ecmax'.or.word.eq.'EmaxH') then ecmax=getvalue() elseif(word.eq.'ecmin'.or.word.eq.'EminH') then ecmin=getvalue() elseif(word.eq.'Emax') then ecmax=getvalue() eemax=ecmax elseif(word.eq.'iscmin') then iscmin=getvalue() elseif(word.eq.'eclow') then eclow=getvalue() elseif(word.eq.'e0cmin') then e0cmin=getvalue() elseif(word.eq.'iecask') then iecask=getvalue() elseif(word.eq.'fecask'.or.word.eq.'fEmaxE') then fecask=getvalue() elseif(word.eq.'eelow') then eelow=getvalue() elseif(word.eq.'themin') then themin=getvalue() elseif(word.eq.'eemin'.or.word.eq.'EminE') then eemin=getvalue() elseif(word.eq.'eemax'.or.word.eq.'EmaxE') then eemax=getvalue() elseif(word.eq.'iocask') then iocask=getvalue() elseif(word.eq.'ifrmod') then ifrmod=getvalue() elseif(word.eq.'ioegs4') then ioegs4=getvalue() elseif(word.eq.'eecrit') then eecrit=getvalue()*1000 ! in MeV elseif(word.eq.'ngener') then ngener=getvalue() elseif(word.eq.'ecut2') then ecut2=getvalue() elseif(word.eq.'iogmag') then iogmag=getvalue() elseif(word.eq.'geoBx') then geoBx=getvalue() * 0.299792458 geoB=sqrt(geoBx**2+geoBy**2+geoBz**2) elseif(word.eq.'geoBy') then geoBy=getvalue() * 0.299792458 geoB=sqrt(geoBx**2+geoBy**2+geoBz**2) elseif(word.eq.'geoBz') then geoBz=getvalue() * 0.299792458 * (-1.0) geoB=sqrt(geoBx**2+geoBy**2+geoBz**2) elseif(word.eq.'xlmin') then xlmin=getvalue() elseif(word.eq.'xlmax') then xlmax=getvalue() elseif(word.eq.'ylmin') then ylmin=getvalue() elseif(word.eq.'ylmax') then ylmax=getvalue() elseif(word.eq.'nxlbin') then nxlbin=getvalue() elseif(word.eq.'nylbin') then nylbin=getvalue() elseif(word.eq.'rlmin') then rlmin=getvalue() elseif(word.eq.'rlmax') then rlmax=getvalue() elseif(word.eq.'nrlbin') then nrlbin=getvalue() elseif(word.eq.'iraxis') then iraxis=getvalue() elseif(word.eq.'ilong') then ilong=getvalue() elseif(word.eq.'iothin') then iothin=getvalue() elseif(word.eq.'thinpar1') then thpar(1)=getvalue() elseif(word.eq.'thinpar2') then thpar(2)=getvalue() elseif(word.eq.'thinpar3') then thpar(3)=getvalue() elseif(word.eq.'thinpar4') then thpar(4)=getvalue() elseif(word.eq.'upkill') then value=getvalue() if(value.eq.1.) upkill=.true. if(value.eq.0.) upkill=.false. #ifdef DETECTORS elseif(word.eq.'tshiftdet') then value=getvalue() if(value.eq.1.) tshiftdet=.true. if(value.eq.0.) tshiftdet=.false. c elseif(word.eq.'surfdet') then c value=getvalue() c if(value.eq.1.) surfdet=.true. c if(value.eq.0.) surfdet=.false. elseif(word.eq.'intdet') then value=getvalue() if(value.eq.1.) intdet=.true. if(value.eq.0.) intdet=.false. elseif(word.eq.'detresc') then value=getvalue() if(value.eq.1.) detresc=.true. if(value.eq.0.) detresc=.false. elseif(word.eq.'detangle') then detangle=getvalue() elseif(word.eq.'detdrad') then detdrad=getvalue() elseif(word.eq.'rdetmax') then rdetmax=getvalue() elseif(word.eq.'hdetfac') then hdetfac=getvalue() elseif(word.eq.'adetabs') then adetabs=getvalue() c elseif(word.eq.'ahdetabs') then c ahdetabs=getvalue() elseif(word.eq.'tdetpoiss') then tdetpoiss=getvalue() elseif(word.eq.'t1res')then t1res=getvalue() elseif(word.eq.'iadet') then iadet=getvalue() elseif(word.eq.'tdetmin') then tdetmin=getvalue() elseif(word.eq.'dtdet') then dtdet=getvalue() elseif(word.eq.'ntdetbin') then ntdetbin=getvalue() if(ntdetbin.gt.mxti)then write (*,*) "too many time bins" write (*,*) "change ntdetbin or mxti in detectors.inc" stop endif #endif elseif(word.eq.'stsort') then stsort=getvalue() elseif(word.eq.'elowhi') then elowhi=getvalue() elseif(word.eq.'ihad') then ihad=getvalue() elseif(word.eq.'zbound') then zbound=getvalue() #ifdef GEANT #ifdef GCALOR elseif(word.eq.'egcalor') then egcalor=getvalue() #endif #endif #ifdef DETRESP elseif(word.eq.'vip') then vip=getvalue() elseif(word.eq.'aratio') then aratio=getvalue() elseif(word.eq.'ares') then ares=getvalue() elseif(word.eq.'depx') then depxmn=getvalue() depxmx=getvalue() ndepx=getvalue() if(.not.haveword)call getw(word) haveword=.false. if(word.eq.'T ') then iodepg=.true. elseif(word.eq.'F ') then iodepg=.false. else stop 'check depx optns !!! ' endif elseif(word.eq.'depa') then depamn=getvalue() depamx=getvalue() ndepa=getvalue() if(.not.haveword)call getw(word) haveword=.false. if(word.eq.'T ') then iodepa=.true. elseif(word.eq.'F ') then iodepa=.false. else print *,word stop 'check depa optns !!! ' endif elseif(word.eq.'detver') then call getw(word) detver=word #endif elseif(word.eq.'iosign')then iosign=getvalue() elseif(word.eq.'PrimaryEnergy') then e0=getvalue() elseif(word.eq.'DebugLevel') then icd=getvalue() elseif(word.eq.'AnalysisLevel') then icasan=int(getvalue()) elseif(word.eq.'PrimaryType') then ii=0 111 ii=ii+1 id0(ii)=getvalue() if(word(1:1).ge.'0'.and.word(1:1).le.'9') goto 111 #ifdef USER_VAR #include "user_var.F" #endif else print *,"unknown variable:",word stop endif elseif(word.eq.'LowEnergyHadronicModel')then call getw(word) #ifdef GEANT if(word.eq.'GHEISHA')then ihad=1 elseif(word.eq.'G-FLUKA')then ihad=2 elseif(word(1:1).eq.'G'.and.word(2:).eq.'CALOR')then ihad=3 else #else if(word.eq.'GHEISHA')then ihad=1 else #endif print *,"possible choices:" print *,"GHEISHA" #ifdef GEANT print *,"G-FLUKA" print *,"G","CALOR" #endif stop endif elseif(word.eq.'HighEnergyHadronicModel')then call getw(word) #ifdef NEXUS2 if(word.ne.'NEXUS2')then print *,"This version is compiled for " print *,"NEXUS2" print *,"only." stop endif #endif #ifdef QGSJET98 if(word.ne.'QGSJET98')then print *,"This version is compiled for " print *,"QGSJET98" print *,"only." stop endif #endif elseif(word.eq.'SaveGroundParticles')then SaveGrndPtcls=.true. elseif(word.eq.'SaveLongParticles')then SaveLongPtcls=.true. elseif(word.eq.'upkill')then upkill=.true. elseif(word.eq.'noupkill')then upkill=.false. #ifdef DETECTORS elseif(word.eq.'DetResc')then detresc=.true. elseif(word.eq.'NoDetResc')then detresc=.false. elseif(word.eq.'IntDet')then intdet=.true. elseif(word.eq.'NoIntDet')then intdet=.false. elseif(word.eq.'GroundDetectors')then idetmod=2 elseif(word.eq.'Detectors')then echo=.false. if(idetmod.eq.0)idetmod=2 120 if(.not.haveword)call getw(word) haveword=.false. if(word.ne.'endDetectors') then ndet=ndet+1 ! number of detectors detname(ndet)=word ! xdet(ndet)=getvalue() + xdetshift ! ydet(ndet)=getvalue() + ydetshift ! hdet(ndet)=getvalue() * hdetfac ! adet(ndet)=getvalue() ! c if(.not.surfdet)then c ahdet(ndet)=getvalue() ! c endif c if(ahdetabs.ne.0.)ahdet(ndet)=ahdetabs ! if(adetabs.ne.0.) adet(ndet)=adetabs ! dangdet(ndet)=detangle ! dxdet(ndet)=sqrt(adet(ndet))/2.0 goto 120 endif haveword=.false. print *,"detectors:",ndet,zdetmax,zdetmin if(checkfile)then write (52,'(/,"endDetectors")') endif echo=.true. elseif(word.eq.'level')then call getw(word) if(word.eq.'Detectors')then print*,"attention: DETECTOR RELEVELING!!!!!!!!!!!!!!!!!!" value=getvalue() do i=1,ndet hdet(i)=value enddo endif elseif(word.eq.'AdjustObservationLevel')then rmin=1e30 do i=1,ndet r=sqrt((xdet(i))**2+(ydet(i))**2) if(r.lt.rmin)then rmin=r hobs=hdet(i) hcore=hobs endif enddo if(checkfile)then write (52,'("Observation level at:",1(g12.6,1x)," meters")') $ hobs endif print *,"Ground level adjusted to ",hobs,"meters" elseif(word.eq.'RandomCorePositionXY')then xmin=getvalue() xmax=getvalue() ymin=getvalue() ymax=getvalue() xshift=-(xmin+(xmax-xmin)*rangen()) yshift=-(ymin+(ymax-ymin)*rangen()) do i=1,ndet xdet(i)= xdet(i) + xshift ydet(i)= ydet(i) + yshift enddo print *,"random core position at",-xshift,-yshift xdetshift=xdetshift+xshift ydetshift=ydetshift+yshift if(checkfile)then write (52,'(/,"core at:",2(g12.6,1x))') -xdetshift,-ydetshift endif elseif(word.eq.'RandomCorePosition')then radius=getvalue() xmin=1e30 ymin=1e30 xmax=-1e30 ymax=-1e30 do i=1,ndet if(xdet(i).lt.xmin) xmin=xdet(i) if(xdet(i).gt.xmax) xmax=xdet(i) if(ydet(i).lt.ymin) ymin=ydet(i) if(ydet(i).gt.ymax) ymax=ydet(i) enddo 100 xshift=-(xmin+(xmax-xmin)*rangen()) yshift=-(ymin+(ymax-ymin)*rangen()) ni=0 iamax=4 do ii=1,iamax ia(ii)=0 enddo do i=1,ndet r=sqrt((xdet(i)+xshift)**2+(ydet(i)+yshift)**2) if(r.gt.radius) then ni=ni+1 angle=atan2(ydet(i)+yshift,xdet(i)+xshift) if(angle.lt.0.)angle=angle+2.0*3.14159265 ii=angle/2.0/3.14159265*float(iamax)+1 ia(ii)=1 endif enddo do ii=1,iamax if(ia(ii).eq.0)goto 100 enddo do i=1,ndet xdet(i)= xdet(i) + xshift ydet(i)= ydet(i) + yshift enddo print *,"random core position at",-xshift,-yshift xdetshift=xdetshift+xshift ydetshift=ydetshift+yshift if(checkfile)then write (52,'(/,"core at:",2(g12.6,1x))') -xdetshift,-ydetshift endif elseif(word.eq.'shift')then call getw(word) if(word.eq.'Detectors')then xshift=getvalue() yshift=getvalue() do i=1,ndet xdet(i)= xdet(i) + xshift ydet(i)= ydet(i) + yshift enddo xdetshift=xdetshift+xshift ydetshift=ydetshift+yshift endif #endif elseif(word.eq.'fakesig')then print*,"attention: altered CROSS-SECTIONS !!!!!!!!!!!!!!!!!!" fakesig=getvalue() elseif(word.eq.'switch')then call getw(word) call getw(word2) if(word.eq.'lateral') then if(word2.eq.'on') ilong=0 if(word2.eq.'off') ilong=1 elseif(word.eq.'projection') then if(word2.eq.'on') lproj=.true. if(word2.eq.'off') lproj=.false. elseif(word.eq.'HadronicCascadeEquations' $ .or.word.eq.'hce') then if(word2.eq.'on') iocask=1 if(word2.eq.'off') iocask=0 elseif(word.eq.'ElecCascadeEquations'.or.word.eq.'ece') then if(word2.eq.'on') iecask=1 if(word2.eq.'off') iecask=0 elseif(word.eq.'hsdj') then if(word2.eq.'on') iecask=1 if(word2.eq.'off') iecask=0 elseif(word.eq.'EGS4') then if(word2.eq.'on') ioegs4=1 if(word2.eq.'off') ioegs4=0 elseif(word.eq.'projection') then if(word2.eq.'on') lproj=.true. if(word2.eq.'off') lproj=.false. else print *,"unknown switch",word stop endif elseif(word.eq.'air')then airz(1)=getvalue() aira(1)=getvalue() airw(1)=getvalue() airz(2)=getvalue() aira(2)=getvalue() airw(2)=getvalue() airz(3)=getvalue() aira(3)=getvalue() airw(3)=getvalue() airava=aira(1)*airw(1)+aira(2)*airw(2)+aira(3)*airw(3) airavz=airz(1)*airw(1)+airz(2)*airw(2)+airz(3)*airw(3) elseif(word.eq.'atmosphere')then do i=1,5 hatm(i)=getvalue() aatm(i)=getvalue() batm(i)=getvalue() catm(i)=getvalue() enddo elseif(word.eq.'analysis')then inexan=1 elseif(word.eq.'analysis2')then inexan=11 elseif(word.eq.'push')then call push elseif(word.eq.'return')then if(ifread.eq.5)then if(checkfile)then write(52,*) endif return else word='ret ' call getw(word) haveword=.true. endif elseif(word.eq.'run')then if(checkfile)then write(52,*) endif return elseif(word.eq.'write')then call getw(word) elseif(word.eq.'writearray')then nco=getvalue() do k=1,nrbins if(nco.eq.2)write(ifhi,'(3(e14.7,1x))')ar(k,1),ar(k,3) if(nco.eq.3)write(ifhi,'(3(e14.7,1x))')ar(k,1),ar(k,3) $ ,ar(k,4) enddo elseif(word.eq.'beginhisto')then call cxini elseif(word.eq.'stop')then call CloseStack() stop else print *,"unknown command:",word stop endif enddo end c----------------------------------------------------------------------- double precision function getvalue() c----------------------------------------------------------------------- c h.drescher c common/ccread/haveword,word,iunit logical haveword character word*80,word2*80 common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 common/ccread2/ ifread,echo logical echo character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan save iunit=0 if(.not.haveword)call getw(word) haveword=.false. i2=index(word,' ')-1 if(index(word(1:i2),'%%').ge.1)then c$$$ write (*,*) fbase c$$$ write (*,*) "--->",word(1:index(word(1:i2),'%%')-1),"<---" c$$$ write (*,*) "--->",index(fbase,word(1:index(word(1:i2),'%%')-1)) c$$$ write (*,*) "--->",word(index(word(1:i2),'%%')+2:i2),"<---" c$$$ write (*,*) "--->",index(word(1:i2),'%%') c$$$ write (*,*) "--->",i2," (i2)" c$$$ write (*,*) "--->",index(fbase c$$$ $ ,word(1:index(word(1:i2),'%%')-1))+index(word(1:i2),'%%')-1 c$$$ write (*,*) "--->",index(fbase c$$$ $ ,word(index(word(1:i2),'%%')+2:i2))-1 c$$$ write (*,*) "--->",word(index(word(1:i2),'%%')+2:i2),"<--" c$$$ write (*,*) "--->" j1=index(fbase,word(1:index(word(1:i2),'%%')-1))+ $ index(word(1:i2),'%%')-1 j2=index(fbase,word(index(word(1:i2),'%%')+2:i2))-1 if(j1.ge.1.and.j2.eq.0)then write (*,*) j1,j2 word=fbase(j1:) i2=index(word,' ')-1 write(*,'("---->",a,"<---",$)') word(1:i2) if(echo.and.checkfile)then write(52,'(a," ",$)') word(1:i2) endif elseif(j1.ge.1.and.j2.ge.j1)then write (*,*) j1,j2 word(1:j2-j1+1)=fbase(j1:j2) word(j2-j1+2:j2-j1+2)=' ' i2=index(word,' ')-1 write(*,'("---->",a,"<---",$)') word(1:i2) if(echo.and.checkfile)then write(52,'(a," ",$)') word(1:i2) endif else print *,"no match in %% expanding" stop endif endif if(index(word(1:i2),'e').ne.0.or.index(word(1:i2),'E').ne.0)then i1=max(index(word(1:i2),'e'),index(word(1:i2),'E')) if(index(word(i1:i2),'.').gt.0)then read(word(1:i1-1),*) getvalue read(word(i1+1:i2),*) exponent c print *,"exp:",word(1:i1-1),"<->",word(i1+1:i2) c $ ,"<->",getvalue,exponent,getvalue*10.0**exponent getvalue=getvalue*10.0**exponent else read(word,*) getvalue endif elseif(index(word(1:i2),'/').gt.0)then i1=index(word(1:i2),'/') read(word(1:i1-1),*) getvalue read(word(i1+1:i2),*) v2 getvalue=getvalue/v2 elseif(index(word(1:i2),'*').gt.0)then i1=index(word(1:i2),'*') read(word(1:i1-1),*) getvalue read(word(i1+1:i2),*) v2 getvalue=getvalue*v2 elseif(index(word(1:i2),'^').gt.0)then i1=index(word(1:i2),'^') read(word(1:i1-1),*) getvalue read(word(i1+1:i2),*) v2 getvalue=getvalue**v2 else read(word,*) getvalue endif call getw(word) if(word.eq.'eV ')then getvalue=getvalue / 1e9 elseif(word.eq.'MeV ')then getvalue=getvalue / 1e3 elseif(word.eq.'GeV ')then elseif(word.eq.'TeV ')then getvalue=getvalue * 1e3 elseif(word.eq.'PeV ')then getvalue=getvalue * 1e6 elseif(word.eq.'EeV ')then getvalue=getvalue * 1e9 elseif(word.eq.'m ')then elseif(word.eq.'g/cm^2 ')then iunit=2 elseif(word.eq.'km ')then getvalue=getvalue * 1e3 elseif(word.eq.'cm ')then getvalue=getvalue / 1e2 elseif(word.eq.'rad ')then elseif(word.eq.'deg ')then getvalue=getvalue/180*3.14159265 elseif(word.eq.'T ')then elseif(word.eq.'uT ')then getvalue=getvalue*1e-6 elseif(word.eq.'nT ')then getvalue=getvalue*1e-9 elseif(word.eq.'mT ')then getvalue=getvalue*1e-3 else haveword=.true. endif return end c----------------------------------------------------------------------- subroutine getw(word) c----------------------------------------------------------------------- c h.drescher c character line*120,word*80 save i,j,line common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 common/ccread2/ ifread,echo logical echo if(word(1:5).eq.'init ')then j=120-1 line(120:120)=' ' return endif if(word(1:4).eq.'ret ') goto 9999 10 i=j+1 do while(i.le.120.and.line(i:i).eq.' ') i=i+1 if(i.eq.121)then read (ifread,'(a)',end=9999) line i=1 if(echo.and.checkfile)write(52,'()') endif enddo j=i i0=i iquote=0 if(line(j:j).eq."'") then iquote=1 elseif(line(j:j+1).eq.'"{') then iquote=2 elseif(line(j:j).eq.'"') then iquote=3 endif do while((line(j+1:j+1).ne.' '.or.iquote.ne.0).and.j+1.le.120) if((line(j+1:j+1).eq."'" $ .or.line(j+1:j+1).eq.'"') $ .and.(iquote.eq.1.or.iquote.eq.3)) then i=i+1 iquote=0 endif if(line(j+1:j+2).eq.'}"'.and.iquote.eq.2) then i=i+2 iquote=0 endif j=j+1 enddo if(line(i:i).eq.'!'.or.line(i:i).eq.'#') then j=120 do while(line(j:j).eq.' '.and.j.gt.i) j=j-1 enddo j=120-1 goto 10 endif word=line(i:j-(i-i0))//' ' c write(*,'("-->",a,"<--")') word if(echo.and.checkfile.and.index(line(i:j),'%%').eq.0)then write(52,'(a," ",$)') line(i0:j) endif return 9999 close(ifread) print *,"end:",ifread if(ifread.eq.9) then word='stop ' return endif ifread=ifread-1 if(ifread.eq.9)ifread=5 if(ifread.eq.4) stop j=120-1 line(120:120)=' ' goto 10 end c............................................................................ subroutine restofcheckfile(starttime,atime,iutime) c........................................................................... c h.drescher c double precision starttime,atime double precision dtime integer iutime(5) common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 if(checkfile)then write(52,*) write(52,'("!",2(1pg13.6,1x)," user/sys sec total")') $ iutime(1)/100.,iutime(2)/100. dtime=(atime-starttime) write(52,'("!",1pg13.6," real sec total ")') dtime write(52,'("!",5x,"=",i3,"days ",i2,":",i2.2,":",i2.2,".",i3.3)') $ int(dtime/86400),int(mod(dtime/3600,24d0)) $ ,int(mod(dtime/60,60d0)),int(mod(dtime,60d0)) $ ,int(mod(dtime/.001,1000d0)) endif end c---------------------------------------------------------------------------- subroutine AddToCheckfile(text,val) c---------------------------------------------------------------------------- c h.drescher c common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80,text*80 if(checkfile)then if(index(text,'$').eq.0)then write (52,'(a,1pg12.6)') text(1:index(text,'&')-1),val else write (52,'(a,1pg12.6,$)') text(1:index(text,'&')-1),val endif endif end c--------------------------------------------------------------------------- subroutine Add2ToCheckfile(text) c---------------------------------------------------------------------------- c h.drescher c common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80,text*80 if(checkfile)then if(index(text,'$').eq.0)then write (52,'(a)') text(1:index(text,'&')-1) else write (52,'(a,$)') text(1:index(text,'&')-1) endif endif end c$$$c----------------------------------------------------------------------- c$$$ subroutine cout c$$$c----------------------------------------------------------------------- c$$$ open(2,file="out") c$$$ if(e2.gt.e1)then c$$$ write(2,'("PrimaryEnergyRange:",2(g13.6,1x))') e1,e2 c$$$ write(2,'("Exponent of Spectrum:",2(g13.6,1x))') egamma c$$$ else c$$$ write(2,'("PrimaryEnergy:",1g13.6)'),e0 c$$$ endif c$$$ write(2,'("Number of showers:"),i6') nshow c$$$ imax=16 c$$$ do while (id0(imax).eq.0) c$$$ imax=imax-1 c$$$ enddo c$$$ write(2,'("PrimaryType:",12i6)') (id0(i),i=0,imax) c$$$ write(2,'("IncidentAngle:",2(g13.6,1x))') acos(costet),phi c$$$ write(2,'("EnergyCutoffs:",4(g13.6,1x))') ehcut,eecut,emcut,epcut c$$$ close(2) c$$$ end c----------------------------------------------------------------------- subroutine stput c----------------------------------------------------------------------- c puts particle i from dptl common block on stack c write to disc if full c common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cdebug/icd,icdsrc double precision seed0,seed1,seed2,seed3 real Seed2Engy,Stack2Engy common /ccsee/ seed0,seed1,seed2,seed3 common /ccsee2/ Seed2Engy,fSeed2Engy,Stack2Engy,fStack2Engy parameter (mxblk=12) parameter (mxstk = mxblk*50*2) common /cstack/stack(mxstk,2),istack(2),irec(2),ifst(2) common /cstene/estack(2) data imxstk/mxstk/ common /coptl/dptl(mxblk) include 'cptl.inc' save maxst k=1 if(dptl(4).ge.Stack2Engy) k=2 if(istack(k).eq.mxstk)then write(ifst(k),rec=irec(k)+1) (stack(i,k),i=1 ,imxstk/2) write(ifst(k),rec=irec(k)+2) (stack(i,k),i=imxstk/2+1,imxstk ) irec(k)=irec(k)+2 istack(k)=0 endif do i=1,12 stack(istack(k)+i,k)=dptl(i) enddo c$$$ if(istack(k)+irec(k)*mxstk.gt.maxst) then c$$$ maxst=(istack(k)+irec(k)*mxstk) c$$$ print *,maxst,irec(k) c$$$ endif istack(k)=istack(k)+mxblk estack(k)=estack(k)+dptl(4) if(icd.ge.5)write(*,100) k,istack(k),dptl 100 format('stput:',i2,1x,i5,1x,12(g10.4,1x)) end c----------------------------------------------------------------------- subroutine stget(iret) c----------------------------------------------------------------------- c h.drescher c c gets particle from stack c or from a source function and put it in dptl c read from disk if empty, iret=1 if no more particles common /cdebug/icd,icdsrc parameter (mxblk=12) parameter (mxstk = mxblk*50*2) common /cstack/stack(mxstk,2),istack(2),irec(2),ifst(2) common /cstene/estack(2) logical isemc common /cstemc/isemc data imxstk/mxstk/ common /coptl/dptl(mxblk) common /cptlty/ imode,icoll,elowhi save ini data ini/1/ double precision s k=2 10 if(istack(k).eq.0)then if(irec(k).eq.0)then if(k.eq.2)then if(ini.eq.0)then call ranfgt(s) print *,"stack 2 now empty, seed=",s ini=1 endif k=1 goto 10 endif if(isemc)then call emcsource(2) iret=0 if(.not.isemc) iret=1 else iret=1 endif return endif read(ifst(k),rec=irec(k)) (stack(i,k),i=1,imxstk/2) irec(k)=irec(k)-1 istack(k)=imxstk/2 else if(k.eq.2) then if(ini.eq.1) then call ranfgt(s) print *,"stack 2 not empty, seed=",s endif ini=0 endif endif iret=0 istack(k)=istack(k)-mxblk do i=1,12 dptl(i)=stack(istack(k)+i,k) enddo estack(k)=estack(k)-dptl(4) if(icd.ge.5)write(*,100) k,istack(k),dptl 100 format('stget:',i2,1x,i5,1x,12(g10.4,1x)) end c----------------------------------------------------------------------- subroutine cptl2stack c----------------------------------------------------------------------- c h.drescher c include 'cptl.inc' parameter (mxblk=12) common /coptl/dptl(mxblk) common /cptlty/ imode,icoll,elowhi common /cgetxor/ getxor,stsort logical getxor ida0=int(dptl(10)) if(stsort.ne.0.0)then njor=0 do i=1,nptl if(istptl(i).eq.0) then njor=njor+1 jorptl(njor)=i do j=1,njor-1 if(pptl(4,jorptl(njor))*stsort $ .gt.pptl(4,jorptl(j))*stsort)then jj=jorptl(njor) jorptl(njor)=jorptl(j) jorptl(j)=jj endif enddo endif enddo do j=1,njor i=jorptl(j) do k=1,5 ! dptl(k)=pptl(k,i) ! copy particle enddo ! to dptl dptl(10)=float(idptl(i)) dptl(11)=qsqptl(i) ! weight if(getxor)then dptl(6)=xorptl(1,i) ! dptl(7)=xorptl(2,i) ! dptl(8)=xorptl(3,i) ! dptl(9)=xorptl(4,i) ! endif call stput ! put dptl on stack enddo else do i=1,nptl if(istptl(i).eq.0)then do j=1,5 ! dptl(j)=pptl(j,i) ! copy particle enddo ! to dptl dptl(10)=float(idptl(i)) dptl(11)=qsqptl(i) ! weight if(getxor)then dptl(6)=xorptl(1,i) ! dptl(7)=xorptl(2,i) ! dptl(8)=xorptl(3,i) ! dptl(9)=xorptl(4,i) ! endif call stput ! put dptl on stack endif enddo endif nptl=0 getxor=.false. end c----------------------------------------------------------------------- subroutine cmuon c----------------------------------------------------------------------- c h.drescher c common/cjprob/Q(10),JPROB,IMUMOD common /cptlty/ imode,icoll,elowhi common /cgetxor/ getxor,stsort logical getxor if(IMUMOD.eq.1) then call gpairm elseif(IMUMOD.eq.2) then call gbremm endif imode=1 ! go back in propagate mode ! not to destroy the actual muon getxor=.true. ! take position from xor not dptl end c----------------------------------------------------------------------- subroutine cegs4 c----------------------------------------------------------------------- c h.drescher c c common /cdebug/icd,icdsrc parameter (mxblk=12) common /coptl/dptl(mxblk) parameter (mxz=3000,mxie=151) common /cslant/ mz common /elost/ elost(mxz) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /chaint/iohadr,ifrmod,ihadt common /compair/airz(3),aira(3),airw(3),airavz,airava common /clshow/lproj,ngener,ecut2,iogmag c ra=rhoair((dptl(8)*costet)) ra=rhoair(dptl(8)) ZLOW=0. if(int(dptl(10)).eq. 12) IQIN=-1 if(int(dptl(10)).eq.-12) IQIN= 1 if(int(dptl(10)).eq. 10) IQIN= 0 EIN=dptl(4)*1000. p=sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) XMIN=dptl(6) YMIN=dptl(7) HMIN=dptl(8) XIN= 0.0 YIN= 0.0 ZIN= depth(dptl(8)) !dptl(8)*costet UIN=dptl(1)/p ! 0.0 VIN=dptl(2)/p ! 0.0 WIN=dptl(3)/p ! 1.0 IRIN=2 ! we start in medium air WTIN=dptl(11) TMIN=dptl(9) c if(ein.gt.10000)then c print *,"auuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu" c $ ,gettheta(dptl(1),dptl(2),dptl(3)) c $ ,ddepth() c endif c print *,IQIN,EIN,XMIN,YMIN,ZIN,TMIN,UIN,VIN,WIN,IRIN,WTIN CALL SHOW(IQIN,EIN,XMIN,YMIN,HMIN,TMIN,UIN,VIN,WIN,IRIN,WTIN) call egs2pptl(0,0,0) ! if anything in egs-list, put it into pptl return end c----------------------------------------------------------------------- subroutine cgheisha c----------------------------------------------------------------------- c h.drescher c c common /cdebug/icd,icdsrc parameter (mxblk=12) common /coptl/dptl(mxblk) parameter (mxz=3000,mxie=151) common /cslant/ mz common /elost/ elost(mxz) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /chaint/iohadr,ifrmod,ihadt common /compair/airz(3),aira(3),airw(3),airavz,airava include 'cptl.inc' nptl=nptl+1 do i=1,5 pptl(i,nptl)=dptl(i) enddo idptl(nptl)=int(dptl(10)) istptl(nptl)=0 iptl=nptl if(icd.ge.4) $ write(*,'(a,i5,1x,3g10.5,$)') 'calling gheisha for:' $ ,idptl(nptl),dptl(4) $ ,dptl(11) call GHEISH(iptl) iptl=1 do while(iptl.le.nptl) nptlo=nptl id=idptl(iptl) ida=abs(id) if(ida.ne.1120.and.ida.ne.1220.and.ida.ne.120.and.ida.ne.130 $ .and.ida.ne.20.and.id.ne.110) then call cdecay(iptl,iret) if( nptlo.ne.nptl ) then istptl(iptl)=1 do i=nptlo+1,nptl istptl(i)=0 enddo ifrptl(1,iptl)=nptlo+1 ifrptl(2,iptl)=nptl endif endif iptl=iptl+1 enddo istptl(1)=1 iorptl(1)=-1 maproj=1 matarg=0 if(icd.ge.7) call clist('list for gheisha&',1,nptl) if(icd.ge.4) write(*,*) ' nptl=',nptl end c c----------------------------------------------------------------------- subroutine ccoll c----------------------------------------------------------------------- c h.drescher c c interactions for dptl particle iptl with air c calls gheisha for low energies c include 'cptl.inc' common/cgenant/siggeant,ngeant common /cdebug/icd,icdsrc parameter (mxblk=12) common /coptl/dptl(mxblk) parameter (mxz=3000,mxie=151) common /cslant/ mz common /elost/ elost(mxz) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /chaint/iohadr,ifrmod,ihadt common /compair/airz(3),aira(3),airw(3),airavz,airava common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /cptlty/ imode,icoll,elowhi double precision ww dimension www(1057) common /chsig/ sighad(3),sigtot common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac c id=int(dptl(10)) ida=abs(id) if(id.eq.10.or.ida.eq.12)then call cegs4 elseif(ida.eq.14)then call cmuon() elseif(icoll.eq.1) then #ifdef GEANT nptl=ngeant ngeant=0 if(icd.ge.7) call clist('list for GEANT&',1,nptl) #else call cgheisha #endif elseif(iohadr.le.9)then !--- high energy hadronic model ------------- call collhh elseif(iohadr.eq.9999)then c$$$ iappl = 1 c$$$ nevent = 1 c$$$ engy = -1 c$$$ pnll = -1 c$$$ elab = dptl(4) c$$$ iframe = 12 ! lab frame c$$$ maproj = 1 ! proton c$$$ laproj = -1 ! -1 for special projectile -> idproj c$$$ idproj = int(dptl(10)) c$$$ if(abs(idproj).eq.20) idproj=230 c$$$c if(abs(idproj).eq.1120.or.abs(idproj).eq.1220) idproj=1120 c$$$c idproj=iabs(idproj) c$$$ if(idproj.ge.100.and.mod(idproj,100).eq.0)then c$$$ maproj= idproj/100 c$$$ laproj= float(maproj) / 2.15 + 0.7 c$$$ idproj=1220 c$$$ elab=elab/maproj c$$$ endif c$$$ i=0 ! choose air-molecule c$$$ r=rangen()*sigtot ! c$$$ do while(r.gt.0.) ! c$$$ i=i+1 ! c$$$ r=r-sighad(i) ! c$$$ enddo ! c$$$ latarg = airz(i) ! c$$$ matarg = aira(i) ! c$$$ c$$$ if(icd.ge.4) c$$$ $ write(*,'(a,2g10.5,i5,5x,1g10.5$)') 'calling nexus for e=' c$$$ $ ,dptl(4) c$$$ $ ,dptl(11),idproj,elost(mz-1)/e0 c$$$ call ainit c$$$ ! call aseed(2) c$$$ call anexus c$$$ call afinal c$$$ if(icd.ge.4)then c$$$ write(*,*) ' nptl=',nptl c$$$ endif c$$$ ifch=6 c$$$ if(icd.ge.7)call clist('list for neXus&',1,nptl) elseif(iohadr.eq.11.or.iohadr.eq.12)then id=int(dptl(10)) ida=abs(id) if(ida.eq.1120.or.ida.eq.1220)then n1=1 elseif(ida.eq.120)then n1=2 elseif(ida.eq.130)then n1=3 elseif(ida.eq.20)then n1=4 else stop 'I don''t know what to do' endif ie=1+nint(log10( dptl(4)/emin )*nde) sum=0. do ii=1,ie do n2=1,7 ee= emin*10.**(float(ii-1)/nde) w=ww(ie,ii,n1,n2) sum=sum+w www(n2+ii*7)=w enddo enddo etot=dptl(4) nfail=0 nptl=0 do while (etot.gt.0.) 123 r=rangen() s=sum*r ! choose n2,ii ii=1 n2=0 do while(s.gt.0) n2=n2+1 if(n2.gt.7)then n2=1 ii=ii+1 endif s=s-www(n2+ii*7) enddo ! chosen ee= emin*10.**(float(ii-1)/nde) if(etot-ee.lt.-1.0)then sum=r*sum nfail=nfail+1 goto 123 endif if(n2.eq.1)id=(1120+100*nint(rangen()))*(2*nint(rangen())-1) if(n2.eq.2)id=120*(2*nint(rangen())-1) if(n2.eq.3)id=130*(2*nint(rangen())-1) if(n2.eq.4)id=-20 if(n2.eq.5)id= 20 if(n2.eq.6)id=110 if(n2.eq.7)id= 10 call icmass(id,amass) pz=sqrt(ee**2-amass**2) nptl=nptl+1 pptl(1,nptl)=0. pptl(2,nptl)=0. pptl(3,nptl)=pz pptl(4,nptl)=ee pptl(5,nptl)=amass istptl(nptl)=0 idptl(nptl)=id iorptl(nptl)=0 iorptl(nptl)=0 ifrptl(1,nptl)=0 ifrptl(2,nptl)=0 etot=etot-ee enddo print *,"mult:",nptl,nfail if(icd.ge.7) call clist('list for W-matrices&',1,nptl) endif if(inexan.ge.1) then call cxana if(inexan.ge.10) call conex2oscar endif return end c----------------------------------------------------------------------- subroutine decay c----------------------------------------------------------------------- c h.drescher c c decays particles in dptl, results in nexus common cptl c common /cdebug/icd,icdsrc parameter (mxblk=12) common /coptl/dptl(mxblk) common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan include 'cptl.inc' common /cgetxor/ getxor,stsort logical getxor id=int(dptl(10)) ida=abs(id) nptl=nptl+1 c pptl(1,nptl)=0 c pptl(2,nptl)=0 c pptl(3,nptl)=sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) do i=1,5 pptl(i,nptl)=dptl(i) enddo idptl(nptl)=id istptl(nptl)=0 iptl=nptl if(ida.eq.1120.or.ida.eq.1220) then print *,dptl stop'ayyyyyyyyyyyy' endif if(id.eq.110) then call cdecay(iptl,iret) ! elseif(ida.eq.120) then call cdecay(iptl,iret) elseif(ida.eq.130) then call cdecay(iptl,iret) elseif(ida.eq.20) then ! etc... call cdecay(iptl,iret) elseif(ida.eq.220) then ! etc... call cdecay(iptl,iret) elseif(ida.eq.14) then ! etc... call cdecay(iptl,iret) else ! do nothing, particle is lost endif do i=iptl+1,nptl if(getxor)then xorptl(1,i)=dptl(6) xorptl(2,i)=dptl(7) xorptl(3,i)=dptl(8) xorptl(4,i)=dptl(9) endif istptl(i)=0 enddo istptl(iptl)=1 ! decayed particle if(icd.ge.6)then istptl(1)=1 iorptl(1)=-1 maproj=1 matarg=0 call clist('list for decay&',1,nptl) c write(*,*) 'decay products:' c do i=1,nptl c write(*,'("nptl ",5(g12.6,1x),3i4)') c $ pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i) c $ ,pptl(5,i),idptl(i),istptl(i) c enddo endif if(inexan.ge.1) then call cxana endif return end c----------------------------------------------------------------------- subroutine cprop c----------------------------------------------------------------------- c h.drescher c parameter (pi=3.14159265) common /cdebug/icd,icdsrc COMMON /GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP + ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD parameter (mxblk=12) common /coptl/dptl(mxblk) common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /cbase2/hinj,haxmax,hax0,haxc,dinj common /ccegs4/ioegs4,eecrit,icesrc common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /CINEX/inex(50) !inex(i)=nexus-id for geant-id (i) common /compair/airz(3),aira(3),airw(3),airavz,airava double precision rlam logical lproj common /clshow/lproj,ngener,ecut2,iogmag parameter (mxz=3000,mxie=151) COMMON /SIGM/ SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO DOUBLE PRECISION SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO common/cross1/xcross(3,10,56),ecross(10),icross common/cjprob/Q(10),JPROB,IMUMOD dimension sigMB(3),sigMP(3) common /cptlty/ imode,icoll,elowhi common /ptlcor/ x1,y1,h1,x2,y2,h2,z1,z2,t1,t2 common /cecut/ ehcut,eecut,emcut,epcut common /smuon/ imuon,dptlm(12) common /chsig/ sighad(3),sigtot common /cgsig/ iosign common /cgeoB/ geoBx,geoBy,geoBz,geoB common/cgenant/siggeant,ngeant common/curve1/rearth logical cnew common /ccore/xcore,ycore,hcore,tcore common /ccta/cta common /cddsum/ddsum common/cnptcl/nptcl common /cmshow/ i0show,ishow,Xfirst,Xfel common /cprim/ hprim,xprim,yprim,tprim,ioprim p=sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) #ifdef CURVED x1=dptl(6) y1=dptl(7) h1=dptl(8) px=dptl(1) py=dptl(2) pz=dptl(3) r1=sqrt(dptl(6)*dptl(6)+dptl(7)*dptl(7)) if(r1.gt.0.0)then cosb=x1/r1 sinb=y1/r1 else cosb=0.0 sinb=0.0 endif pxb= cosb*px+sinb*py c pyb=-sinb*px+cosb*py pzb= pz ha1=sqrt((rearth+h1)**2+r1**2)-rearth ! altitude over sea level sina=r1/(rearth+ha1) cosa=(rearth+h1)/(rearth+ha1) c pxa= cosa*pxb+sina*pzb c pya= pyb pza=-sina*pxb+cosa*pzb cta=pza/p z1=depth(ha1) ! vertical depth Xv t1=dptl(9) #else cta=dptl(3)/sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) x1=dptl(6) y1=dptl(7) h1=dptl(8) z1=depth(h1) ! vertical depth Xv t1=dptl(9) ha1=h1 #endif id=int(dptl(10)) wt=dptl(11) EK=dptl(4)-dptl(5) c if ((ngener.ge.0.and.ngener.le.int(dptl(12))) c $ .or.dptl(4).lt.ecut2) then c imode=5 c dd=0 c if (lproj)then c dd = (zmax - depth(dptl(8)))/cta c imode=2 c endif c goto 9999 c endif c.....search for geant code for id ...................... ipart=0 i=1 if(id.eq.110.or.id.eq.220)id=120 ! just temporarily do while(i.lt.48.and.inex(i).ne.id) i=i+1 enddo ipart=i DENS=0. CORR=0. NLM=3 id=int(dptl(10)) ida=abs(id) np=0 rzd=1e30 ! initialize decay length rzi=1e30 ! initialize interaction length sigtot=0. icoll=0 if(ida.ne.14) imuon=0 if(ida.eq.1120.or.ida.eq.1220)then np=1 B = -1 ! negative means stable elseif(ida.eq.120)then np=2 B = 114.0 elseif(ida.eq.130)then np=3 B = 852.0 elseif(id.eq.-20)then np=4 B = 205.0 elseif(id.eq.20)then np=5 B = 1.19069e5 elseif(id.eq.220)then ! eta imode=3 ! decay return elseif(id.eq.110)then ! pi0 np=2 B=3.5e10 imode=3 ! decay ig=dptl(12) call cana2(x1,y1,h1,z1,t1,x1,y1,h1,z1,t1 $ ,dptl(1),dptl(2),dptl(3),dptl(4),dptl(5),wt,id,ig) goto 9999 elseif(id.eq.10.or.ida.eq.12)then ! photon or electron if(ioegs4.ne.0)then dd=0. !DZEGS/ct np=6 imode=6 goto 9999 ! jump back all is done in EGS4 else stop 'do not set ioegs4 to zero' endif elseif(ida.eq.14) then !muons imuon=imuon+1 if(imuon.eq.1)then do i=1,12 dptlm(i)=dptl(i) enddo endif B=1.02 imode=1 rzi=10 sigtot=-1.0 sigMBT=0. sigMPT=0. do i=1,3 sigMB(i)=gbrsgm(airz(i),dptl(4))*airw(i) sigMP(i)=gprsgm(airz(i),dptl(4))*airw(i) sigMBT=sigMBT+sigMB(i) sigMPT=sigMPT+sigMP(i) enddo sigtot=(sigMBT+sigMPT)/1000. rziB=-log(rangen())*airava*1.660/max(sigMBT,1e-10) rziP=-log(rangen())*airava*1.660/max(sigMPT,1e-10) rzi=min(rziP,rziB) if(icd.ge.8)then write (*,*) "sigMBT:",sigMBT,"sigMPT:",sigMPT write (*,*) "i-length:",rziB, rziP,10. endif if(rzi.gt.1e1)then ! just propagate imode=1 rzi=1e1 ! muon-step sigtot=-1.0 elseif(rzi.eq.rziB)then ! bremsstrahlung r=rangen() if(r.lt.sigMB(1)/sigMBT)then JPROB=1 elseif(r.lt.(sigMB(1)+sigMB(2))/sigMBT)then JPROB=2 else JPROB=3 endif IMUMOD=1 imode=4 else ! pair-production r=rangen() if(r.lt.sigMP(1)/sigMPT)then JPROB=1 elseif(r.lt.(sigMP(1)+sigMP(2))/sigMPT)then JPROB=2 else JPROB=3 endif imode=4 IMUMOD=2 endif 320 continue elseif(id.ge.100.and.mod(id,100).eq.0)then B=-1. np=10 #ifdef GEANT elseif(id.eq.45)then ! 47 = d B=-1.0 elseif(id.eq.46)then ! 47 = t B=-1.0 elseif(id.eq.47)then ! 47 = alpha B=-1.0 #endif else print *,"discard",id,ek,dptl(5) imode = 5 ! discard (or do something with it ) goto 9999 endif if(icd.ge.5)then write (*,*) 'height (g/cm^2)',depth(dptl(8)) ,h1,rhoair(h1)/100. endif c.....decay ???......... if(B.gt.0.)then r=rangen() dl= -log(r)*p*6400./B if( ha1-dl*cta .lt. 0. ) then rzd=1e30 elseif(dl*cta.gt.10)then rzd= ( depth(ha1-dl*cta)-depth(ha1) ) /cta else rzd=dl*rhoair(ha1) endif endif c.....or collision ???.... if(sigtot.eq.0.) then if( dptl(4) .lt. elowhi ) then #ifdef GEANT icoll=1 call GRUN call greset dzi=0.001205/siggeant sigtot=airava*1660.44/dzi if(icd.ge.5) write(*,*) 'GEANTSIG:',sigtot rzi = -log(rangen()) * dzi #else icoll=1 sigtot=GHESIG(P,EK,airava,airA,airZ,airW,NLM,DENS,CORR,IPART) if(icd.ge.5) write(*,*) 'GHEISHASIG:',sigtot dzi = airava*1660.44/sigtot ! interaction length rzi = -log(rangen()) * dzi #endif else icoll=2 if(np.lt.10)then sigtot=0. do i=1,3 sighad(i)=sighh(dptl(4),np,aira(i)) * airw(i) sigtot=sigtot+sighad(i) enddo dzi = airava*1660.44 / sigtot rzi = -log(rangen()) * dzi c print *,sigtot,dzi,rzi elseif(id.ge.100.and.mod(id,100).eq.0)then if(iosign.eq.0)then ii=0 ! search cross section in data base elab=dptl(4)/float(int(id/100)) do i=1,10 if(elab.eq.ecross(i)) ii=i enddo if(ii.eq.0) then icross=mod(icross,10)+1 ii=icross ecross(ii)=elab do j=1,56 do i=1,3 xcross(i,ii,j)=0. enddo enddo endif !------------------------------------ cnew=.true. if(xcross(1,ii,id/100).le.0.)then c----------------------------new sigtot=0. do i=1,3 sigma=sighh(dptl(4),1,1.0) sighad(i)=crtot(float(id/100),aira(i),sigma)*airw(i) xcross(i,ii,id/100)=sighad(i) sigtot=sigtot+sighad(i) if(icd.ge.1)print *,'sigma ',aira(i),'-(',id/100,',' $ ,int(float(id)/100/ 2.15 + 0.7),')' $ ,sighad(i)/max(airw(i),1e-9) enddo if(icd.ge.1)print *,'sigma Air-(',id/100,',' $ ,int(float(id)/100/ 2.15 + 0.7),')',sigtot c---------------------------- else c----------------------------new sigtot=0. do i=1,3 sighad(i)=xcross(i,ii,id/100) sigtot=sigtot+sighad(i) enddo endif else sigtot=0.0 do i=1,3 sigma=sighh(dptl(4),id,aira(i)) sighad(i)=sigma*airw(i) sigtot=sigtot+sighad(i) if(icd.ge.1)print *,'sigma ',aira(i),'-(',id/100,',' $ ,int(float(id)/100/ 2.15 + 0.7),')' $ ,sighad(i)/max(airw(i),1e-9) enddo if(icd.ge.1)print *,'sigma Air-(',id/100,',' $ ,int(float(id)/100/ 2.15 + 0.7),')',sigtot endif dzi = airava*1660.44 / sigtot rzi = -log(rangen()) * dzi endif endif endif c.....choose one ............ if ( (rzd .gt. rzi .and. ifi0.ne.3).or.ifi0.eq.4 ) then if(sigtot.gt.0)imode=4 ! interaction dd = rzi else icoll=0 ! imode=3 ! decay dd = rzd if(ifi0.eq.3)dd=1 endif cc ddmax=2.5 cc dd=10.0 cc if(dd.gt.ddmax)then cc imode=1 cc dd=ddmax cc ddsum=ddsum+dd cc if(nptcl.eq.0) write(60,*) ddsum cc else cc ddsum=ddsum+dd cc if(nptcl.eq.0) write(60,*) ddsum cc if(nptcl.eq.0)then cc print *,ddsum,nptcl,ishow cc write(59,*) ishow,ddsum cc endif cc endif ifi0=0 dhl=dd*sqrt(1-cta**2)/rhoair(ha1) if(dhl.gt.1000.0)then dd=dd*1000.0/dhl imode=1 endif cc dhl=dd*sqrt(1-cta**2)/rhoair(ha1) 123 z2=z1+dd*cta ! second vertical depth if ( z2.le.0.) then ! particle leaves atmosphere ???? imode=5 goto 9999 endif ha2=height(z2) dha=(ha1-ha2) if(abs(dha).lt.1.0)then dl= dd / rhoair(ha1) else dl= dha / cta endif call icchrg(id,chrg) if(iogmag.ge.1.and.chrg.ne.0.)then c.......propagate one half step................... dl=dl/2 x2=x1 + dl * dptl(1) / p y2=y1 + dl * dptl(2) / p h2=h1 - dl * dptl(3) / p t2=t1 + dl * dptl(4) / p dha=dl*cta ha2=ha1-dha z2=depth(ha2) dptl(6)= x2 dptl(7)= y2 dptl(8)= h2 ! depth(h2)/costet dptl(9)= t2 ig=dptl(12) call cana2(x1,y1,h1,z1,t1,x2,y2,h2,z2,t2 $ ,dptl(1),dptl(2),dptl(3),dptl(4),dptl(5),wt,id,ig) x1=x2 y1=y2 h1=h2 ha1=ha2 z1=z2 t1=t2 c.......adjust momentum........................... c tesla * e = 0.299792458 GeV/c /m c the times 2 in the end is to correct for the half step width dl/2 dpx=dptl(1) + ( dptl(2)*geoBz-dptl(3)*geoBy ) * chrg * dl / p*2. dpy=dptl(2) + (-dptl(1)*geoBz+dptl(3)*geoBx ) * chrg * dl / p*2. dpz=dptl(3) + ( dptl(1)*geoBy-dptl(2)*geoBx ) * chrg * dl / p*2. c print *,(dptl(i),i=1,3),dpx,dpy,dpz pn=sqrt(dpx**2+dpy**2+dpz**2) dptl(1) = dpx * p/pn dptl(2) = dpy * p/pn dptl(3) = dpz * p/pn c.......moliere scattering for muons ...................... if (ida.eq.14) then beta2=(p/dptl(4))**2 sigma=dd*(0.021/(dptl(4)*beta2))**2/37.7 theta=sigma*sin(rangen()*2.*pi)*sqrt(-2.0*log(rangen())) phi=rangen()*2.0*pi ddx=sin(theta)*cos(phi) ddy=sin(theta)*sin(phi) ddz=sqrt(1.0-ddx*ddx+ddy*ddy) call utrotc(-1,ddx,ddy,ddz,dptl(1),dptl(2),dptl(3)) endif c.......propagate one more half step.............. c print *,(dptl(i),i=1,3) #ifdef CURVED r1=sqrt(x1*x1+y1*y1) if(r1.gt.0.0)then cosb=x1/r1 sinb=y1/r1 else sinb=0.0 cosb=0.0 endif pxb= cosb*px+sinb*py c pyb=-sinb*px+cosb*py pzb= pz ha1=sqrt((rearth+h1)**2+r1**2)-rearth ! altitude over sea level sina=r1/(rearth+ha1) cosa=(rearth+h1)/(rearth+ha1) c pxa= cosa*pxb+sina*pzb c pya= pyb pza=-sina*pxb+cosa*pzb ! p=sqrt(pza*pza+pxa*pxa+pya*pya) cta=pza/p #else cta=dptl(3)/p #endif dha=dl * cta endif c.....energy loss ......................................... if (ida.ge.14.and.chrg.ne.0.) then gamma2=(dptl(4)/dptl(5))**2 if(gamma2.gt.1.001)then beta2=(gamma2-1e0)/gamma2 dEdx = gamma2/(gamma2-1.0)*0.153287 $ * (log(gamma2-1.0)-beta2+9.386 ) dptl(4) = max(dptl(4) - dEdx * 1e-3 * dd,0.) pn=sqrt(max(dptl(4)**2-dptl(5)**2,0.)) dptl(1)=dptl(1)*pn/p dptl(2)=dptl(2)*pn/p dptl(3)=dptl(3)*pn/p p=pn if(ida.eq.14)then if (dptl(4)-dptl(5).le.emcut) then imode=5 dptl(4)=dptl(5) goto 9999 endif else if (dptl(4)-dptl(5).le.ehcut) then imode=5 dptl(4)=dptl(5) goto 9999 endif endif else imode=5 dptl(4)=dptl(5) goto 9999 endif endif if( ha1-dha .le. hground )then imode=2 ! particle reaches ground endif x2=x1 + dl * dptl(1) / p y2=y1 + dl * dptl(2) / p h2=h1 - dl * dptl(3) / p t2=t1 + dl * dptl(4) / p dha=dl*cta ha2=ha1-dha z2=depth(ha2) if(icd.ge.2)then d1=haxis(x1,y1,h1) d2=haxis(x2,y2,h2) dd2=38.0 dl2=alength(dd2) endif dptl(6) = x2 dptl(7) = y2 dptl(8) = h2 ! z2/costet ! slant depth dptl(9) = t2 ig=dptl(12) call cana2(x1,y1,h1,z1,t1,x2,y2,h2,z2,t2 $ ,dptl(1),dptl(2),dptl(3),dptl(4),dptl(5),wt,id,ig) c dzint=zint(xprim,yprim,hprim,x2,y2,h2) if(nptcl.eq.0) write (61,*) dzint,ddsum,z2,z2/cta c call zint(x1,y1,h1,x2,y2,h2) c if(ida.ne.14) c $ write (99,'(i5,10(1x,g12.6))') id,x1,y1,h1,t1,x2,y2,h2,t2 if(haxmax.ne.0.0)then if(haxis(x2,y2,h2).gt.haxmax)imode=2 endif 9999 if(icd.ge.5)then write (*,*) 'cprop:',imode write (*,*) 'decay,collide length:',rzd,rzi if(imode.eq.1) write(*,*) 'propagate' if(imode.eq.2) write(*,*) 'ground' if(imode.eq.3) write(*,*) 'decay' if(imode.eq.4) write(*,*) 'collide' if(imode.eq.5) write(*,*) 'discard' if(imode.eq.6) write(*,*) 'EGS4' endif end c----------------------------------------------------------------------- function depth(h) c----------------------------------------------------------------------- c c common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) if ( h .lt. HATM(2) ) then depth = aatm(1) + batm(1) * exp ( -h / catm(1) ) elseif ( h .lt. HATM(3) ) then depth = aatm(2) + batm(2) * exp ( -h / catm(2) ) elseif ( h .lt. HATM(4) ) then depth = aatm(3) + batm(3) * exp ( -h / catm(3) ) elseif ( h .lt. HATM(5) ) then depth = aatm(4) + batm(4) * exp ( -h / catm(4) ) else depth = aatm(5) - h * batm(5)/catm(5) endif end c----------------------------------------------------------------------- function height(t) c----------------------------------------------------------------------- c c common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) if ( t .gt. DATM(2) ) then height = catm(1) * log ( batm(1) / (t - aatm(1)) ) elseif ( t .gt. DATM(3) ) then height = catm(2) * log ( batm(2) / (t - aatm(2)) ) elseif ( t .gt. DATM(4) ) then height = catm(3) * log ( batm(3) / (t - aatm(3)) ) elseif ( t .gt. DATM(5) ) then height = catm(4) * log ( batm(4) / (t - aatm(4)) ) else height = (aatm(5) - t) * catm(5)/batm(5) endif end c----------------------------------------------------------------------- function rhoair(h) c----------------------------------------------------------------------- c c attention units are g/(cm^2 m) c common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) if ( h .lt. HATM(2) ) then rhoair = batm(1) / catm(1) * exp ( -h / catm(1) ) elseif ( h .lt. HATM(3) ) then rhoair = batm(2) / catm(2) * exp ( -h / catm(2) ) elseif ( h .lt. HATM(4) ) then rhoair = batm(3) / catm(3) * exp ( -h / catm(3) ) elseif ( h .lt. HATM(5) ) then rhoair = batm(4) / catm(4) * exp ( -h / catm(4) ) else rhoair = batm(5) / catm(5) endif end c----------------------------------------------------------------------- function rhoaiX(X) c----------------------------------------------------------------------- c c same but argument is depth (g/cm^2) c attention units are g/(cm^2 m) c common /atmos/ AATM(5),BATM(5),CATM(5),HATM(5),DATM(5) if ( X .gt. DATM(2) ) then rhoaiX = (X - aatm(1)) / catm(1) elseif ( X .gt. DATM(3) ) then rhoaiX = (X - aatm(2)) / catm(2) elseif ( X .gt. DATM(4) ) then rhoaiX =(X - aatm(3)) / catm(3) elseif ( X .gt. DATM(5) ) then rhoaiX =(X - aatm(4)) / catm(4) else rhoaiX = batm(5)/catm(5) endif end subroutine nptlreset2 common/cnptcl/nptcl include 'cptl.inc' nptcl=nptcl+nptl nptl=0 end subroutine nptlreset include 'cptl.inc' nptl=0 end c----------------------------------------------------------------------- subroutine cana2(x1,y1,h1,z1,t1,x2,y2,h2,z2,t2,px,py,pz,E,am $ ,wt,id,ig) c----------------------------------------------------------------------- c h.drescher c c analyzes a particle going from (x1,y1,h1,z1,t1) to (x2,y2,h2,z2,t2) c c id = nexus/isajet particle id c x1,y1,h1,x2,y2,h2 = position in meter above sea-level c z1,z2 = slant depth in g/cm^2 c t1,t2 = time in meter c divide by 0.299792458 to get nano-seconds c px,py,pz,E,am = momentum,mass in GeV,GeV/c^2 c wt = weight c c in case of abs(id)<=1 we have EGS4 which means: c 0 = photon c -1 = electron c 1 = positron c px,py,pz are the directional vectors c (not necessarily normalized to 1, but almost) c E and am are in MeV units c parameter (mxblk=12) common /coptl/dptl(mxblk) PARAMETER (TWOPI=6.28318530717958648,PI=3.14159265358979324) PARAMETER (EMASS2=.00000026112003932088) common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /ccore/xcore,ycore,hcore,tcore parameter (mxz=3000,mxie=151) common /cslant/ mz common /londat/ hz(6,mxie,mxz),del(mxz),dpo(mxz),dmuon(mxz) common /loelph/ hez(2,-30:mxie,mxz) parameter (mxla=300,mrla=60) common /latda2/ drad2(10,mrla),drad3(10,mrla),drad4(10,mrla) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis common /latda3/ dlat(4,mxla,mxla), $ xlmin,xlmax,ylmin,ylmax,nxlbin,nylbin common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis common /cprim/ hprim,xprim,yprim,tprim,ioprim common /ccana2/ dzana common /cthin/ ilong,iothin,thinH,thinE $ ,ethin,hthin,wtmaxH,wtmaxE,rthmax,gthinH,gthinE parameter (nmlast=1000) common /clast/ zlast(0:nmlast),elast(0:nmlast),idlast(0:nmlast) $ ,modlast(0:nmlast) common /alldat3/ SaveGrndPtcls,SaveLongPtcls logical SaveGrndPtcls,SaveLongPtcls #ifdef DETECTORS c.....detectors............................ #include "detectors.inc" #endif #ifdef DETRESP c.....detectors...................... #include "detresp.inc" #endif c.................. common /smuon/ imuon,dptlm(12) c.................. double precision a,b,c,ss logical igo data c/0.299792458/ save c common /ccta/cta c write(*,'(i5,20(1x,g12.6))') c $ id,z1/cta,z2/cta,x1,y1,h1,x2,y2,h2,cta ee=E amm=am ida=abs(id) if(ida.le.1)then ee=ee/1000. ! if (ida.le.1) means EGS4 which has MeV units amm=am/1000. endif #ifdef USER_ANALYSIS if(ida.le.1)then if(ida.eq.1)then p=sqrt(ee*ee-EMASS2) idisa=-sign(12,id) else p=ee idisa=10 endif ppx=px*p ppy=py*p ppz=pz*p else idisa=id ppx=px ppy=py ppz=pz endif call cana(x1,y1,h1,z1,t1,x2,y2,h2,z2,t2,ppx,ppy,ppz,ee,amm $ ,wt,idisa,ig) #endif il1=int(max(0.,min(z1,z2))/abs(costet)/dzana) !?> -1 il2=int(min(max(z1,z2),zmax)/abs(costet)/dzana) !?> -1 IF (il1.lt.il2.and.id.ne.110) THEN ie=nint(log10(ee/emin)*nde) + 1 do ll=il1+1,il2 if (id.eq.0) THEN ie=max(ie,-30) hez(1,ie,ll)=hez(1,ie,ll)+wt elseif(ida.eq.1)then ie=max(ie,-30) if(id.eq.1)then dpo(ll)=dpo(ll)+wt else del(ll)=del(ll)+wt endif hez(2,ie,ll)=hez(2,ie,ll)+wt elseif(id.eq.14)then ie=max(ie,1) dmuon(ll)=dmuon(ll) + wt else ie=max(ie,1) in=0 if(ida.eq.1120.or.ida.eq.1220) then in=1 elseif(ida.eq.120) then in=2 elseif(ida.eq.130) then in=3 elseif(id.eq.-20) then in=4 elseif(id.eq. 20) then in=5 endif if(in.gt.0) hz(in,ie,ll)=hz(in,ie,ll)+wt endif enddo elseif(id.eq.110)then il1=nint(z1/costet/dzana) if(il1.gt.0)then ie=nint(log10(ee/emin)*nde) + 1 ie=max(ie,1) in=6 hz(in,ie,il1)=hz(in,ie,il1)+wt endif return endif c............................................................................. if(SaveLongPtcls)then xref=x2 yref=y2 href=h2-hcore call toaxis(xref,yref,href) r2=(xref)*(xref)+(yref)*(yref) if(r2.gt.rthmax*rthmax.or.rangen().lt.0.001.or.h1.gt.20000)then write (99,'(i5,10(1x,g12.6))') id,x1,y1,h1,t1,x2,y2,h2,t2 c if(ida.eq.1)then c write (99,'(i5,10(1x,g12.6))') id,x1,y1,z1,t1,px,py,pz,ee,wt c endif endif endif c.....check radial stuff................................................ if(h2.le.hobs.and.h1.gt.hobs)then dh=hobs-h1 if(abs(dh).gt.0.001)then t=dh/(h2-h1) xx=x1+t*(x2-x1) yy=y1+t*(y2-y1) tt=t1+t*(t2-t1) else xx=x1 yy=y1 tt=t1 endif rg=sqrt(xx**2+yy**2) ! distance along the ground if(iraxis.eq.1)then hr=0.0 xr=xx yr=yy ! call toaxis(xr,yr,hr) r=sqrt(xr**2+yr**2) ! distance to shower-axis else xr=xx yr=yy r=rg endif IF (id.eq.0) THEN ! here check for particles ll=3 ld=1 elseif(id.eq.-1)then ll=1 ld=2 elseif(id.eq.1)then ll=1 ld=3 elseif(ida.eq.14)then ll=2 ld=4 elseif(id.eq.1120)then ll=4 ld=6 elseif(id.eq.-1120)then ld=0 ll=5 elseif(id.eq.1220)then ld=5 ll=6 elseif(id.eq.-1220)then ld=0 ll=7 elseif(ida.eq.120)then ld=0 ll=8 elseif(ida.eq.130)then ld=0 ll=9 else ld=0 ll=0 ENDIF if(SaveGrndPtcls)then xref=x2 yref=y2 href=h2-hcore call toaxis(xref,yref,href) r2=(xref)*(xref)+(yref)*(yref) if(r2.gt.10000.or.rangen().lt.0.001.or.h1.gt.20000)then write (99,'(i5,10(1x,g12.6))') id,x1,y1,h1,t1,x2,y2,h2,t2 c if(ida.eq.1)then c write (99,'(i5,10(1x,g12.6))') id,x1,y1,z1,t1,px,py,pz,ee,wt c endif endif endif IF (r.gt.rlmin.and.r.lt.rlmax) THEN ir=int(log(r/rlmin)/log(rlmax/rlmin)*nrlbin)+1 if(ll.ne.0)then drad(ll,ir)=drad(ll,ir)+wt drad2(ll,ir)=drad2(ll,ir)+wt**2 drad3(ll,ir)=drad3(ll,ir)+1.0 #ifdef DETRESP c...........edep in rad ................................................... if(ld.ne.0)then sec= sqrt(px*px+py*py+pz*pz)/pz ccc write (99,*) i,sec,dphi,tt,id,xx,yy,hh ccc write (99,*) i,xxx,yyy,hhh,xx,yy,hh ek=ee-amm ! kinetic Energy in GeV edep=getedep(ek,sec,ld,fa) drad4(ll,ir)=drad4(ll,ir)+wt*edep*fa else edep=0.0 fa=1.0 endif #endif c.......................................................................... endif endif ix=int((xx-xlmin)/(xlmax-xlmin)*nxlbin)+1 iy=int((yy-ylmin)/(ylmax-ylmin)*nylbin)+1 if(ix.ge.1.and.ix.le.nxlbin.and.iy.ge.1.and.iy.le.nylbin)then ll=0 if(ida.eq.1)then ll=1 elseif(ida.eq.14)then ll=2 endif if(ll.ne.0)then dlat(ll,ix,iy)=dlat(ll,ix,iy)+wt if(dlat(ll+2,ix,iy).ne.0.)then if(tt.lt.dlat(ll+2,ix,iy)) dlat(ll+2,ix,iy)=tt else dlat(ll+2,ix,iy)=tt endif endif endif #ifdef DETECTORS c.....check GroundDetectors............................................. if(idetmod.eq.2) then phi=atan2(yr,xr) if (id.eq.0) then ! here check for particles ll=3 ld=1 elseif(id.eq.-1)then ll=1 ld=2 elseif(id.eq.1)then ll=1 ld=3 elseif(ida.eq.14)then ll=2 ld=4 elseif(id.eq.1220)then ll=4 ld=5 elseif(id.eq.1120)then ll=5 ld=6 else ll=0 ld=0 endif do i=1,ndet dphi=angdet(i)-phi if (dphi.gt. pi) dphi= dphi-twopi if (dphi.lt.-pi) dphi= dphi+twopi igo=.true. if(ll.eq.0) igo=.false. if(dangdet(i).ge.0.)then if (abs(dphi).ge.dangdet(i)) igo=.false. endif if(igo)then c-c exchange r by rg if(r.gt.r1det(i).and.r.lt.r2det(i))then sec= sqrt(px*px+py*py+pz*pz)/pz ccc write (99,*) i,sec,dphi,tt,id,xx,yy,hh ccc write (99,*) i,xxx,yyy,hhh,xx,yy,hh ek=ee-amm ! kinetic Energy in GeV ie=nint(log10(max(ek,0.0008))*10.0) + 1 if(ek.ge.depemn)then #ifdef DETRESP edep=getedep(ek,sec,ld,fa) #else fa=1.0 #endif c if(surfdet)then f=wt/aSdet(i)*adet(i)*fa c else c ttan=sqrt(px*px+py*py)/pz c f=wt/aSdet(i)*(adet(i)+ttan*ahdet(i))*fa/sec c endif if(f.lt.tdetpoiss)then expmean=exp(-f) !mean of poisson pir = 1 N = -1 do while(.true.) N=N+1 pir = pir*rangen() if(pir.le.expmean)goto 9 enddo 9 wtnew=float(N) else if(intdet)then wtnew=float(int(f+rangen())) else wtnew=f endif endif if(tshiftdet)then xr=xx-xdet(i) yr=yy-ydet(i) hr=0.0 call toaxis(xr,yr,hr) ttt=(tt+hr)/0.299792458 else ttt=tt/0.299792458 endif itr=int((ttt-tdetmn(i))/dtdet) if (ttt.lt.tdetmn(i)) itr=itr-1 it=mod(itr+it0det(i),ntdetbin) itr0=itr if(itr.lt.0)then itr=max(itr,-ntdetbin) do itt=-1,itr,-1 ii=mod(itt+it0det(i)+ntdetbin,ntdetbin) do l2=1,5 #ifdef DETRESP edepdet(l2,i,ntdetbin)=edepdet(l2,i,ntdetbin) $ +edepdet(l2,i,ii) edepdet(l2,i,ii)=0.0 #endif tidet(l2,i,ntdetbin)=tidet(l2,i,ntdetbin) $ +tidet(l2,i,ii) tidet(l2,i,ii)=0.0 enddo enddo it0det(i)=mod(it0det(i)+ntdetbin+itr,ntdetbin) tdetmn(i)=tdetmn(i)+dtdet*itr0 it=it0det(i) elseif(itr.ge.ntdetbin)then it=ntdetbin else !print *,"normal",i,it endif #ifdef DETRESP edepdet(ll,i,it)=edepdet(ll,i,it)+edep*wtnew if(idetout.gt.0.and.wtnew.gt.0.)then if(idetout.eq.1)then write(86,186) i,id,ttt,wtnew,px,py,pz,ek $ ,edep elseif(idetout.eq.2)then write(86,186) i,id,ttt,wtnew,xx,yy,hh,px,py,pz,ek $ ,edep endif endif #else if(idetout.gt.0.and.wtnew.gt.0.)then if(idetout.eq.1)then write(86,186) i,id,ttt,wtnew,px,py,pz,ek else(idetout.eq.2)then write(86,186) i,id,ttt,wtnew,xx,yy,hh,px,py,pz,ek endif endif #endif 186 format(i4,1x,i5,18(1x,1pg13.6)) tidet(ll,i,it)=tidet(ll,i,it)+wtnew endet(ll,i,ie)=endet(ll,i,ie)+wtnew ! wtnew or f endif endif endif enddo endif #endif ! detectors endif 9999 return end c----------------------------------------------------------------------------- function getedep(ek,sec,ld,fa) c----------------------------------------------------------------------------- c h.drescher c c gets energy deposit (or equivalent) for detectors c as a function of angle and Ekin and particle type c 2002/09/20 h.drescher c 2002/09/26 interpolation added h.drescher include 'detresp.inc' save theta fa=0.0 if(ld.ne.0)then iek=int(log(ek/depemn)/log(depemx/depemn)*ndepe)+1 if(iek.gt.ndepe)iek=ndepe if(iodepa)then isec=int((sec-depamn)/(depamx-depamn)*ndepa)+1 else theta=acos(1.0/sec) isec=int((theta-depamn)/(depamx-depamn)*ndepa)+1 endif if(isec.gt.ndepa)isec=ndepa if(isec.ge.1.and.iek.ge.1)then ix=0 fa=max(1.0,fdepa*(1.0-depneu(iek,ix,isec,ld))) edep=0. ran=1.0-fa/fdepa*rangen() do while(ix.le.ndepx.and.ran.gt.depneu(iek,ix,isec,ld)) ran=ran-depneu(iek,ix,isec,ld) ix=ix+1 enddo if(ix.eq.0)then edep=0.0 elseif(ix.gt.ndepx)then ix=ix-1 do while(depneu(iek,ix,isec,ld).eq.0.0) ix=ix-1 enddo if(iodepg)then edep=depxmn*(depxmx/depxmn)**((float(ix)-.5)/ndepx) else edep=depxmn+(depxmx-depxmn)*((float(ix)-.5)/ndepx) endif else if(ix.ge.2.and.ix.le.ndepx-1)then y1=depneu(iek,ix-1,isec,ld) y2=depneu(iek,ix ,isec,ld) y3=depneu(iek,ix+1,isec,ld) A1=(y1+y2)/2 A2=(y2+y3)/2 r=A1/(A1+A2) if(rangen().lt.r)then r=rangen() rr=(y1**2-2.0*A1*r*y1+2.0*r*A1*y2) if(rr.lt.-1e-6)then stop'rr<0' elseif(rr.lt.0.0)then rr=0.0 endif if(y1.ne.y2) then x1=(y1+sqrt(rr))/(y1-y2) x2=(y1-sqrt(rr))/(y1-y2) else x2=r endif c print *,A1,A2,ix,x,"left",rr if(x2.lt.0.or.x2.gt.1.0) then 10 continue goto 10 endif x=float(ix)-1.0+x2 else r=rangen() rr=(y2**2-2.0*A2*r*y2+2.0*r*A2*y3) if(rr.lt.-1e-6)then stop'rr<0' elseif(rr.lt.0.0)then rr=0.0 endif if(y2.ne.y3) then x1=(y2+sqrt(rr))/(y2-y3) x2=(y2-sqrt(rr))/(y2-y3) else x2=r endif c print *,A1,A2,ix,x,"right",rr if(x2.lt.0.or.x2.gt.1.0) then 20 continue goto 20 endif x=float(ix)+x2 endif else x=float(ix) endif if(iodepg)then edep=depxmn*(depxmx/depxmn)**((x-.5)/ndepx) else edep=depxmn+(depxmx-depxmn)*((x-.5)/ndepx) endif endif c print *,ld,theta,sec,isec,ek,iek,ix,edep,iodepa c print *,depamn,depamx,ndepa,sec,isec c print *,depxmn,depxmx,ndepx,edep,ix,iodepg endif if(ares.ne.0.0)then r=(rangen()+rangen()+rangen()+rangen()-2.0)/2.0 edep=(edep)*10.0**(r*ares) endif else edep=0.0 endif getedep=edep end c----------------------------------------------------- subroutine chkname(in,out) c h.drescher c common /cmshow/ i0show,ishow,Xfirst,Xfel character*80 in,out,nu j=index(in,'%n') if( j.ne.0) then ii=log(float(ishow))/log(10.) if(ii.eq.0)write(nu,'(i1)') ishow if(ii.eq.1)write(nu,'(i2)') ishow if(ii.eq.2)write(nu,'(i3)') ishow if(ii.eq.3)write(nu,'(i4)') ishow if(j.eq.1)then out=nu(1:1+ii)//in(j+2:index(in,' ')) else out=in(1:j-1)//nu(1:1+ii)//in(j+2:index(in,' ')) endif else out=in endif end c----------------------------------------------------- subroutine i0ini c h.drescher c common /cmshow/ i0show,ishow,Xfirst,Xfel logical schonda character*80 fname character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 schonda=.true. if(index(fbase,'%n').eq.0)then print *,'no %n' return endif i0show=-i0show-1 do while(schonda) i0show=i0show+1 ishow=1+i0show call chkname(fbase,fname) inquire(file=fname,exist=schonda) print *,fname(1:20),schonda enddo end c----------------------------------------------------------------------- subroutine cfinal c----------------------------------------------------------------------- c c 10/10/2002 description added h.drescher c c h.drescher c COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax common /cdebug/icd,icdsrc parameter (mxblk=12) common /coptl/dptl(mxblk) common /compair/airz(3),aira(3),airw(3),airavz,airava common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground parameter (mxz=3000,mxie=151) common /cslant/ mz character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan common /londat/ hz(6,mxie,mxz),del(mxz),dpo(mxz),dmuon(mxz) common /loelph/ hez(2,-30:mxie,mxz) parameter (mxla=300,mrla=60) common /latda2/ drad2(10,mrla),drad3(10,mrla),drad4(10,mrla) common /latdat/ drad(10,mrla),rlmin,rlmax,nrlbin,iraxis common /latda3/ dlat(4,mxla,mxla), $ xlmin,xlmax,ylmin,ylmax,nxlbin,nylbin common /elost/ elost(mxz) common /ccask/e0cmin,ecmin,eclow,ecmax,fcask,iocask,thin2,athin2 $ ,iscmin c.....casmod common character*80 fname dimension a(4),sig(mxz),ia(4),covar(4,4),alpha(4,4),dyda(4) $ ,yfit(mxz),xrad(mrla),zfit(mxz) external fnel,frad,f2rad,f3rad,f4rad common /cmshow/ i0show,ishow,Xfirst,Xfel COMMON /ccel1/ dZZ,iezmax,ieemax,ieemin,ieelow,iecask common /ccegs4/ioegs4,eecrit,icesrc common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac common /ccana2/ dzana logical iaa(8) #ifdef DETECTORS c.....detectors....... #include "detectors.inc" #endif c..... common /cclel/ hes(2,-150:601,3000) common /cppj4/ mneie,mxeie,mxen1,mxen2,iende c..... maxZana=zmax/dzana c........................................................................ c radial normalization pi=acos(-1.0) do ir=1,nrlbin r1=rlmin*(rlmax/rlmin)**((float(ir)-1.0)/nrlbin) r2=rlmin*(rlmax/rlmin)**((float(ir+1)-1.0)/nrlbin) fnorm=(r2**2-r1**2)*pi/oosho if(iraxis.eq.1.and.costet.gt.0.0) fnorm=fnorm/costet r=rlmin*(rlmax/rlmin)**((float(ir)-0.5)/nrlbin) do j=1,10 drad(j,ir)=drad(j,ir)/fnorm drad2(j,ir)=sqrt(drad2(j,ir))/fnorm drad3(j,ir)=drad3(j,ir)/fnorm drad4(j,ir)=drad4(j,ir)/fnorm enddo enddo c longitudinal normalization do j=1,ipmax do k=1,mxz do i=1,6 hz(i,j,k)=hz(i,j,k)*oosho enddo enddo enddo do j=-30,ipmax do k=1,mxz do i=1,2 hez(i,j,k)=hez(i,j,k)*oosho enddo enddo enddo do i=1,mxz del(i)=del(i)*oosho dpo(i)=dpo(i)*oosho dmuon(i)=dmuon(i)*oosho c elost(i)=elost(i)*oosho enddo c.......................................................................... if(icasan.ge.1)then c.......longitudinal....................................................... call chkname(fbase(1:index(fbase,' ')-1)//"-l.dat",fname) if(icasan.eq.1. $ .or.(imodan.ne.0.and.index(fbase,'%n').eq.0)) then open(2,file=fname,status='unknown',access='append') print *,'append to l-output file' else open(2,file=fname,status='unknown') write(2,'(a)') '# ' write(2,'(a)') '# table of particles in atmospheric depth' write(2,'(a)') '# ' write(2,'(a)') '# fit-paramter in #max line' write(2,'(2a)') '# prim energy,prim type,' $ ,'Nmax,X0,Xmax,Lambda,Xfirst,Xfel' write(2,'(a)') '# Xfirst: depth of first interaction' write(2,'(2a)') '# Xfel: depth of last interaction' $ ,' with e>0.9 e0' write(2,'(a)') '# ' write(2,'(a,a)') $ '# depth nucleons pions kaons+- Klong', $ ' Kshort pi0 photon e+- mu+-' print *,'l-output file opened' endif if(icasan.ge.2)then do j1=1,maxZana sn=0e0 sp=0e0 skz=0e0 skl=0e0 sks=0e0 spi0=0e0 sph=0e0 sel=0e0 do i=-30,0 sph = sph + ( hez(1,i,j1) ) !/ ee sel = sel + ( hez(2,i,j1) ) !/ ee enddo do i=1,ipmax ee= emin*10.**(float(i-1)/nde) sn = sn + ( hz(1,i,j1) ) ! / ee sp = sp + ( hz(2,i,j1) ) ! / ee skz = skz + ( hz(3,i,j1) ) ! / ee skl = skl + ( hz(4,i,j1) ) ! / ee sks = sks + ( hz(5,i,j1) ) ! / ee spi0 = spi0 + ( hz(6,i,j1) ) ! / ee sph = sph + ( hez(1,i,j1) ) !/ ee sel = sel + ( hez(2,i,j1) ) !/ ee enddo stot=sn+sp+skz+skl+sks write(2,124) zmin+dzana*(j1-1),sn,sp,skz,skl,sks,spi0,sph $ ,sel,dmuon(j1) !,del(j1),dpo(j1) enddo endif c.......fitting del(i) to formula c.......Nmax*((x-x0)/(Xmax-x0))**((Xmax-X0)/labda)*exp((Xmax-x)/lambda) c comment out if you don't want the fitting c goto 512 a(1)=0. a(2)=-100. a(4)=60. do i=1,maxZana sig(i)=0.1 yfit(i)=(del(i)+dpo(i))/e0 zfit(i)=zmin+dzana*(i-1) if( yfit(i).gt.a(1)) then imax=i a(1)=yfit(i) a(3)=zmin+dzana*(i-1) endif enddo do i=1,4 ia(i)=1 enddo alamda=-1. ma=4 ndata=maxZana chisq=1e30 ochisq=2e30 ntry=0 do while (ochisq.gt.chisq.or.ntry.le.10) ochisq=chisq call mrqmin(zfit,yfit,sig,ndata,a,ia,ma $ ,covar,alpha,ma,chisq,fnel,alamda,iret) if(iret.eq.1) then print *,'error in longitudinal fitting' goto 512 endif ntry=ntry+1 enddo c$$$ do i=1,maxZana c$$$ call fnel(z(i),a,y,dyda,ma) c$$$ write (21,*) z(i),yfit(i),y c$$$ enddo c$$$ call flush(21) c$$$ rewind (21) x=max( 0. , a(2) ) do while ( x. lt. a(3) ) y=a(1)*((x-a(2))/(a(3)-a(2)))**((a(3)-a(2))/a(4)) $ * exp((a(3)-x)/a(4)) if( y.lt.a(1)/100. ) then x10=x else x=a(3) endif x=x+0.1 enddo print *,'--------------------------------------------------' print *,"fit to Gaisser-Hillas function:" print *,"f(x)=Nmax*((x-X0)/(Xmax-X0))**((Xmax-X0)/Lambda)"// $ "*exp((Xmax-x)/Lambda)" print *,'Nmax= ',a(1)*e0,' # (',del(imax),')' print *,'Xmax= ',a(3),' # (',zmin+dzana*imax,')' print *,'X0= ',a(2) print *,'Lambda= ',a(4) write (2,'("#max: ",(g12.6,1x),i5,1x,10(g12.6,1x))') $ e0,id0(0), $ a(1)*e0,a(2),a(3),a(4),chisq/ndata $ ,Xfirst,Xfel goto 512 512 write (2,'(/)') close(2) write(*,*) 'done:',fname(1:index(fname,' ')-1) endif c.....radial.............................................................. if(icasan.ge.2)then call chkname(fbase(1:index(fbase,' ')-1)//"-r.dat",fname) if(icasan.eq.1. $ .or.(imodan.ne.0.and.index(fbase,'%n').eq.0)) then open(2,file=fname,status='unknown',access='append') print *,'append to r-output file' else open(2,file=fname,status='unknown') write(2,*) '# ' write(2,*) '# table of ground particles ' write(2,*) '# column 1: r in m' write(2,*) '# 2: log10(r)' write(2,*) $ '# 3-11: density of e+-, mu+-, gamma, ' $ ,'p,ap,n,an, pi+-, K+-' write(2,'(a,a)') '# 13-21: sqrt(-**2)', $ ' of same particles' #ifdef DETRESP write(2,'(a,a)') '# 23-31: detector-response', $ ' of same particles' #endif print *,'r-output file opened' endif c c call chkname(fbase(1:index(fbase,' ')-1),fname) c c open(2,file=fname,status='unknown') do ir=1,nrlbin r=rlmin*(rlmax/rlmin)**((float(ir)-0.5)/nrlbin) write(2,123)r,log10(r),(drad(j,ir),j=1,10) $ ,((drad2(j,ir)),j=1,10) #ifdef DETRESP $ ,(drad4(j,ir),j=1,10) #endif enddo write (2,'(/)') close(2) 123 FORMAT(80(1pe12.4)) write(*,*) 'done:',fname(1:index(fname,' ')-1) endif c........................................................................... if(inexan.ge.1) then call cfnex endif #ifdef DETECTORS c.....output detector data.................................................... if(ndet.gt.0)then do i=1,8 if(mod(iadet/10**(i-1),10).eq.1)then iaa(i)=.true. else iaa(i)=.false. endif enddo vipe=vip*exp(-1.0) dt=dtdet if (iaa(1)) then fname=fbase(1:index(fbase,' ')) write (fname(index(fbase,' '):),'("-det1.dat ")') call chkname(fname,fname) open(11,file=fname(1:index(fname,' ')-1)) endif if (iaa(2)) then fname=fbase(1:index(fbase,' ')) write (fname(index(fbase,' '):),'("-det2.dat ")') call chkname(fname,fname) open(12,file=fname(1:index(fname,' ')-1)) endif if (iaa(3)) then fname=fbase(1:index(fbase,' ')) write (fname(index(fbase,' '):),'("-det3.dat ")') call chkname(fname,fname) open(13,file=fname(1:index(fname,' ')-1)) endif #ifdef DETRESP #endif if (iaa(2)) then write (12,'("# number densities in detectors:")') write (12,'("# columns are:")') write (12,'("# 1 rr distance from axis")') write (12,'("# 2 xdet x-coordinate of detector")') write (12,'("# 3 ydet y-coordinate of detector")') write (12,'("# 4 hdet height of detector")') write (12,'("# 5 adet area of detector")') write (12,'("# 6 nep number of e+-")') write (12,'("# 7 nmu number of muons")') write (12,'("# 8 npho number of photons")') write (12,'("# 9 nneu number of neutrons")') write (12,'("# 10 npro number of protons")') #ifdef DETRESP write (12,'("# 11-15 energy deposit of same ptcls")') #endif endif if (iaa(3)) then write (13,'("# arrival times in detectors: ")') write (13,'("# columns:")') write (13,'("# 1 time [ns]")') write (13,'("# 2 number of e+-")') write (13,'("# 3 number of muons")') write (13,'("# 4 number of photons")') write (13,'("# 5 number of neutrons")') write (13,'("# 6 number of protons")') #ifdef DETRESP write (13,'("# 7-11 energy deposit of same ptcls")') #endif endif if (iaa(1)) then write (11,'("# kinetic energy spectrum of electrons" $ ,",muons,photons,neutrons,protons ")') write (11,'("# colums: ")') write (11,'("# 1 : energy (10 bins per decade) ")') write (11,'("# 2 : number of e+- in bin ")') write (11,'("# 3 : number of muons ")') write (11,'("# 4 : number of photons ")') write (11,'("# 5 : number of neutrons ")') write (11,'("# 6 : number of protons ")') write (11,'("# ")') endif do i=1,ndet xr=xdet(i) yr=ydet(i) hr=0.0 call toaxis(xr,yr,hr) rr=sqrt(xr**2+yr**2) iemx=mxie do while( $ endet(1,i,iemx).eq.0.and. $ endet(2,i,iemx).eq.0.and. $ endet(3,i,iemx).eq.0.and. $ endet(4,i,iemx).eq.0.and. $ endet(5,i,iemx).eq.0.and. $ iemx.gt.-30) iemx=iemx-1 enddo if (iaa(1)) $ write(11,500) i,detname(i),depth(hdet(i)),hdet(i),rr 500 format('# detector ',i3,1x,a10 $ ,'at ',f5.1,"g/cm^2=",f5.1,"m" $ ," r_axis=",g13.6,"m ") sum1=0 sum2=0 sum3=0 sum4=0 sum5=0 do ie=-30,iemx ee=emin*10.**(float(ie-1)/nde) write(11,'(10(g13.6,1x))') $ ee,(endet(j,i,ie),j=1,5) sum1=sum1+endet(1,i,ie) sum2=sum2+endet(2,i,ie) sum3=sum3+endet(3,i,ie) sum4=sum4+endet(4,i,ie) sum5=sum5+endet(5,i,ie) enddo if (iaa(1)) then write(11,'("#",40a)') ("-",j=1,40) write(11,'("#",14x,10(g13.6,1x))') sum1,sum2,sum3,sum4,sum5 write(11,'("#",40a,/,/)') ("-",j=1,40) endif #ifdef DETRESP c.........energy deposit in detectors edep1=0. edep2=0. edep3=0. edep4=0. edep5=0. itmin=1 tim1=0. tim2=0. tim3=0. tim4=0. itmin=0 c itmax=ntdetbin-1 c do while (edpndet(i,itmax).eq.0..and.itmax.ge.1) c itmax=itmax-1 c enddo itmax=ntdetbin-1! is ok if we print only != 0 values if (iaa(3)) then write(13,'("# detector:",i3,1x,a20,1g13.6)') i,detname(i),rr endif do itr=itmin,itmax it=mod(itr+it0det(i),ntdetbin) t=tdetmn(i)+dtdet*(float(itr)+0.5) sum=0. do j=1,5 sum=sum+tidet(j,i,it) enddo if(sum.ne.0.)then if (iaa(3)) then write(13,'(15(1x,g13.6))') t $ ,(tidet(j,i,it) ,j=1,5) #ifdef DETRESP $ ,(edepdet(j,i,it),j=1,5) #endif endif edep1=edep1+edepdet(1,i,it) edep2=edep2+edepdet(2,i,it) edep3=edep3+edepdet(3,i,it) edep4=edep4+edepdet(4,i,it) edep5=edep5+edepdet(5,i,it) endif enddo if (iaa(3)) then write(13,'(14x,15(1x,g13.6))') (tidet(j,i,mxti),j=1,5) endif edep1=edep1+edepdet(1,i,ntdetbin) edep2=edep2+edepdet(2,i,ntdetbin) edep3=edep3+edepdet(3,i,ntdetbin) edep4=edep4+edepdet(4,i,ntdetbin) edep5=edep5+edepdet(5,i,ntdetbin) if (iaa(3)) then write(13,'(/)') endif #endif if (iaa(2)) then write(12,'(15(1x,g13.6))') $ rr $ ,xdet(i),ydet(i) $ ,hdet(i),adet(i) #ifdef DETRESP $ ,sum1,sum2,sum3,sum4,sum5,edep1,edep2,edep3,edep4,edep5 #else $ ,sum1,sum2,sum3,sum4,sum5 #endif ! DETRESP endif enddo ! end detector-loop print *,"iadet",iadet if(iaa(1))close(11) if(iaa(2))close(12) if(iaa(3))close(13) endif #endif ! DETECTORS #ifdef DETECTORS if(idetout.gt.0) close(86) #endif c....................................................... entry deletearrays write (*,'("deleting arrays .... ",$)') do j=1,mxie do k=1,mxz do i=1,6 hz(i,j,k)=0 enddo enddo enddo do j=-30,mxie do k=1,mxz do i=1,2 hez(i,j,k)=0 enddo enddo enddo do i=1,mxz del(i)=0. dpo(i)=0. dmuon(i)=0. elost(i)=0. enddo do i=1,4 do j=1,mxla do k=1,mxla dlat(i,j,k)=0. enddo enddo enddo do i=1,10 do j=1,mrla drad(i,j)=0. drad2(i,j)=0. drad3(i,j)=0. drad4(i,j)=0. enddo enddo #ifdef DETECTORS do i=1,ndet do ie=-30,mxie do n=1,mxdp endet(n,i,ie)=0. enddo enddo do it=0,mxti #ifdef DETRESP do n=1,mxdp edepdet(n,i,it)=0. enddo #endif do n=1,mxdp tidet(n,i,it)=0. enddo enddo enddo #endif !DETECTORS ipmax=1 izmin=mz write (*,'("done")') 124 FORMAT(F7.1,20(1pe12.3)) 125 FORMAT(F7.1,4(1pE12.3)) 126 FORMAT(10(1pE12.3)) end c----------------------------------------------------------------------- subroutine cfnex c h.drescher c----------------------------------------------------------------------- c.....nexus include 'cptl.inc' c..... common/ccread1/checkfile,checkfilename logical checkfile character checkfilename*80 character word*80 common/cana4/ar(751,5),nrbins,ifhi,inexan,hisfac character*80 home,host,fbase,ftmp common /cdata/home,host,fbase,ftmp,icasan,imodan double precision getvalue common/ccread2/ ifread,echo logical echo,echold write(6,'(a)')'rewind copy-file'//checkfilename//'<--' ifold=ifread echold=echo ifread=9 echo=.false. close(52) open(ifread,file=checkfilename,status='old') ifhi=31 open(ifhi,file=fbase(1:index(fbase,' ')-1)//'.histo') n=0 word=' ' hisfac=1.0 call getw('init ') do while(word.ne.'stop ') call getw(word) if(word.eq.'write ')then call getw(word) write (ifhi,*) word elseif(word.eq.'set ')then call getw(word) if(word.eq.'hisfac ')then call getw(word) read(word,*) hisfac endif elseif(word.eq.'endhisto ')then n=n+1 print *,"histo:",n call cxhis(n) elseif(word.eq.'writearray')then call getw(word) read(word,*)nco write(ifhi,'(a,i3)')'array',nco do k=1,nrbins if(nco.eq.2)write(ifhi,'(3(e14.7,1x))')ar(k,1),ar(k,3) if(nco.eq.3)write(ifhi,'(3(e14.7,1x))')ar(k,1),ar(k,3) $ ,ar(k,4) enddo write(ifhi,'(a)')' endarray' endif enddo close(ifhi) ifread=ifold echo=echold open(52,file=checkfilename(1:index(checkfilename,' ')-1) $ ,access='append') end c------------------------------------------------------------------------ function alength(z) c------------------------------------------------------------------------ c c returns length along shower axis for particle at slant depth z c common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground parameter (mxblk=12) common /coptl/dptl(mxblk) common /cslant/ mz common/curdef/arho(20000),brho(20000),ha(20000),hl(20000) iz1=int((z-zmin)/dz)+1 iz2=int((z-zmin)/dz)+2 h1=hl(iz1) h2=hl(iz2) X1=max(dz*(iz1-1),1e-5) X2=dz*(iz2-1) h0=(h2-h1)/log(X1/X2) X0=X2/exp(-h2/h0) alength=-h0*log(z/X0) return end c------------------------------------------------------------------------ function sdepth(xx,yy,hh) c------------------------------------------------------------------------ c c returns slant depth along shower axis for particle (x,y,h) c common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground parameter (mxblk=12) common /coptl/dptl(mxblk) common /cslant/ mz common/curdef/arho(20000),brho(20000),ha(20000),hl(20000) save iz data iz/1/ hax=haxis(xx,yy,hh) iz=1 do while(hax.gt.hl(iz).and.iz.le.mz) iz=iz+1 enddo do while(hax.lt.hl(iz-1).and.iz.gt.1) iz=iz-1 enddo if(iz.eq.1.and.hax.lt.hl(iz))then sdepth=-2.5 return endif h2=hl(iz) h1=hl(iz-1) X2=dz*(iz-1) X1=max(dz*(iz-2),1e-5) c ddepth2=X1+dz*(hax-h1)/(h2-h1) dh=(h2-h1) if(abs(dh).gt.10.0)then h0=(h2-h1)/log(X1/X2) X0=X2/exp(-h2/h0) X=X0*exp(-hax/h0) else X=X1+dz*(hax-h1)/dh endif c print *,X0*exp(-hax/h0),X0,h0,X1/exp(-h1/h0) c ddepth=depth(dptl(8))/costet c print *,ddepth,ddepth2,depth(dptl(8)),X sdepth=X !sdepth=depth(h)/costet return end c------------------------------------------------------------------------ function ddepth() c------------------------------------------------------------------------ c c returns slant depth along shower axis for particle in dptl c common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground parameter (mxblk=12) common /coptl/dptl(mxblk) common /cslant/ mz common/curdef/arho(20000),brho(20000),ha(20000),hl(20000) save iz data iz/1/ c c prim ------------>-------->---------->------- core c hl(iz-1) -- hax -- hl(iz) c ddepth=depth(dptl(8))/costet c return hax=haxis(dptl(6),dptl(7),dptl(8)) do while(iz.lt.mz.and.hax.ge.hl(iz)) iz=iz+1 enddo do while(iz.gt.1.and.hax.lt.hl(iz-1)) iz=iz-1 enddo if(iz.ge.mz.and.hax.gt.hl(iz))then ddepth=dz*(mz+1) return endif if(iz.eq.1.and.hax.le.hl(iz))then ddepth=-2.5 return endif h2=hl(iz) h1=hl(iz-1) X2=dz*(iz-1) X1=max(dz*(iz-2),1e-5) ddepth2=X1+dz*(hax-h1)/(h2-h1) dh=(h2-h1) if(abs(dh).gt.10.0)then h0=(h2-h1)/log(max(1e-5,X1)/X2) X0=X2/exp(-h2/h0) X=X0*exp(-hax/h0) else X=X1+dz*(hax-h1)/dh endif c-c print *,X0*exp(-hax/h0),X0,h0,X1/exp(-h1/h0) c-c print *,ddepth,ddepth2,depth(dptl(8)),X ddepth=X c-c ssl=alength(X) c-c print *,ssl return end c----------------------------------------------------------------------- subroutine ccuts c----------------------------------------------------------------------- c h.drescher c c apply cuts, thin sampling c c include 'cptl.inc' common /cdebug/icd,icdsrc common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /ccore/xcore,ycore,hcore,tcore parameter (mxz=3000,mxie=151) common /cslant/ mz parameter (mxblk=12) common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cecut/ ehcut,eecut,emcut,epcut common /cthin/ ilong,iothin,thinH,thinE $ ,ethin,hthin,wtmaxH,wtmaxE,rthmax,gthinH,gthinE common /elost/ elost(mxz) dimension etype(0:10) common /cthpar/ thpar(10) common /chaint/iohadr,ifrmod,ihadt common /cgetxor/ getxor,stsort logical getxor logical lproj common /clshow/lproj,ngener,ecut2,iogmag c.....EGS4 common/ccegs4/ioegs4,eecrit,icesrc COMMON/BOUNDS/ECUT(3),PCUT(3),VACDST common /cptlty/ imode,icoll,elowhi common /cupkill/ upkill logical upkill common /ccana2/ dzana e1=0. e2=0. nbaryon=0 ida0=iabs(int(dptl(10))) njor=0 if(ngener.eq.int(dptl(12))) then call nptlreset2 imode=0 return endif do i=1,nptl if(istptl(i).eq.0)then c initialize weight if(.not.getxor) qsqptl(i)=dptl(11) if(ilong.ne.0)then if(mod(ilong,10).eq.1.and.(ida0.ne.10.and.ida0.ne.12) $ .or.mod(ilong/10,10).eq.1.and.(ida0.eq.10.or.ida0.eq.12) )then c......... if you want to test cask....................... if(ida.ge.100.or.ida.eq.20.and.pptl(4,i).gt.0.1)then ee=pptl(4,i) ie=1+nint(log10( pptl(4,i) /emin )*nde) if(ie.ge.1) pptl(4,i) = emin*10.**(float(ie-1)/nde) qsqptl(i)= qsqptl(i) * ee/pptl(4,i) endif pptl(3,i)=sqrt(max(0.,pptl(4,i)**2-pptl(5,i)**2)) c pptl(3,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2) pptl(1,i)=0. pptl(2,i)=0. endif endif if(icoll.eq.2) then cc rotate to initial particle cc sca=pptl(3,i)*sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) call utrotc(-1,dptl(1),dptl(2),dptl(3) $ ,pptl(1,i),pptl(2,i),pptl(3,i) ) cc scn=0. cc do k=1,3 cc scn=scn+dptl(k)*pptl(k,i) cc enddo cc print *,sca,scn,dptl(3),sqrt(dptl(1)**2+dptl(2)**2) cc $ /sqrt(dptl(1)**2+dptl(2)**2+dptl(3)**2) endif c if ( abs(pptl(2,i)).gt.3000.0)then c print *,"stop here" c istop=1 c 10 continue c if(istop.eq.1) goto 10 c endif c.........main cut ...................................... id=idptl(i) ida=abs(id) if (ida.eq.11.or.ida.eq.13) istptl(i)=1 if(ida.ge.20)then ! ------ hadrons ------------------ if (ida.ge.1120) nbaryon=nbaryon+sign(1,id) if( pptl(4,i)-pptl(5,i) .lt. ehcut ) then istptl(i)=1 endif c kill upwards if( upkill.and.pptl(3,i) .lt.0 ) istptl(i)=1 endif ! ------------------------------------------ if(ida.eq.12.or.ida.eq.10)then ! --- leptons + photons ----------- if(ida.eq.12.and. pptl(4,i)-0.511e-3.lt.eecut ) istptl(i)=1 if(ida.eq.10.and. pptl(4,i) .lt.epcut ) istptl(i)=1 endif c........................................................ c if(istptl(i).eq.1)then c elost(iz)=elost(iz)+pptl(4,i)*dptl(11) c e2=e2+pptl(4,i) c else c e1=e1+pptl(4,i) c endif c........................................................ endif enddo if(imode.eq.4)then if(icd.ge.7) print *,'baryons:',nbaryon if (dptl(4).lt.elowhi .or. iohadr.eq.1 ) then if(abs(dptl(10)).eq.1120.or.abs(dptl(10)).eq.1220) $ nbaryon=nbaryon-sign(1,int(abs(dptl(10)))) if(int(dptl(10)).ge.100.and.mod(int(dptl(10)),100).eq.0)then if(iohadr.eq.1)then nbaryon=nbaryon-dptl(10)/100 endif endif c elost(iz)=elost(iz) - 0.93828*nbaryon*dptl(11)*f endif endif c if(icd.ge.7) print *, 'elost',elost(mz),e1,e2,icoll c.....thin sampling... c using iorptl() for index c if( dptl(4) . lt. hthin ) then c if( thin.gt.0. ) then do i=0,10 etype(i)=0. enddo etot=0. ntot=0 do i=1,nptl if(istptl(i).eq.0.and.pptl(4,i).lt.hthin) then ! ntot=ntot+1 ! number iorptl(ntot)=i ! index for valid particles etot=etot+pptl(4,i) ! sum up energy ity=itype(idptl(i)) ! etype(ity)=etype(ity)+pptl(4,i) c print *,'thin1:',pptl(4,i),nptl endif enddo enthin=0. if(iothin.eq.1)then c.........classical thinning................................ ernd=rangen()*etot i=1 do while(i.le.ntot) iptl=iorptl(i) ernd=ernd-pptl(4,iptl) if(ernd.lt.0.)then istptl(iptl)=0 ernd=2.0*etot ! so all others get istptl=1 enthin = pptl(4,iptl) qsqptl(iptl)=qsqptl(iptl) * etot / enthin else istptl(iptl)=1 endif i=i+1 enddo c.......second option max weight=wtmaxH..................... elseif(iothin.eq.2)then in=1 do while(in.le.ntot) iptl=iorptl(in) ity=itype(idptl(iptl)) prb=pptl(4,iptl)/etot if(qsqptl(iptl)/prb .gt. wtmaxH )prb=1. prb=min(prb,1.0) if( rangen() .lt. prb )then !accepted istptl(iptl)=0 weight=1./prb qsqptl(iptl) = qsqptl(iptl) * weight else istptl(iptl)=1 endif in=in+1 enddo elseif(iothin.eq.4)then if(dptl(11).eq.1.0)then in=1 do while(in.le.ntot) iptl=iorptl(in) ity=itype(idptl(iptl)) prb=1.0/wtmaxH if( rangen() .lt. prb )then !accepted istptl(iptl)=0 weight=wtmaxH qsqptl(iptl) = qsqptl(iptl) * weight else istptl(iptl)=1 endif in=in+1 enddo endif c.......no lateral thinning (experimental).................... elseif(iothin.eq.3)then xx=dptl(6) yy=dptl(7) hh=dptl(8)-hcore call toaxis(xx,yy,hh) if(xx*xx+yy*yy.lt.rthmax*rthmax)then in=1 do while(in.le.ntot) iptl=iorptl(in) ity=itype(idptl(iptl)) prb=pptl(4,iptl)/etot if(qsqptl(iptl)/prb .gt. wtmaxH )prb=1. prb=min(prb,1.0) if( rangen() .lt. prb )then !accepted istptl(iptl)=0 weight=1./prb qsqptl(iptl) = qsqptl(iptl) * weight else istptl(iptl)=1 endif in=in+1 enddo endif endif c............................................................ endif c.....end thinning... end c----------------------------------------------------------------------- subroutine iegs4 c----------------------------------------------------------------------- c calls initialization routines c c h.drescher c COMMON/BOUNDS/ECUT(3),PCUT(3),VACDST COMMON/MEDIA/ RLC(1),RLDU(1),RHO(1),MSGE(1),MGE(1),MSEKE(1),MEKE(1 *),MLEKE(1),MCMFP(1),MRANGE(1),IRAYLM(1),NMED COMMON/MEDIAC/MEDIA(24,1) CHARACTER*4 MEDIA COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(3),RHOR(3),IRAYLR(3) COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(1),AE(1),UP(1),UE(1),TE(1),THMOLL *(1) CHARACTER*4 MEDARR(24) COMMON/RANDOM/IXX DATA MEDARR /'A','I','R',21*' '/ common /cecut/ ehcut,eecut,emcut,epcut double precision seed0,seed1,seed2,seed3 common /ccsee/ seed0,seed1,seed2,seed3 common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /cppj4/ mneie,mxeie,mxen1,mxen2,iende iende=50 DO I=1,24 MEDIA(I,1)=MEDARR(I) ENDDO MED(1)=0 MED(2)=1 MED(3)=0 DUNIT=-2 ECUT(2)=eecut*1000.+0.511 ! cut-off for electrons in MeV PCUT(2)=epcut*1000. ! cut-off for photons in MeV print *,epcut,eecut WRITE(6,1020) 1020 FORMAT('1 START '//' CALL HATCH TO GET CROSS-SECTION DATA'/ *) CALL HATCH WRITE(6,1030)AE(1)-0.511, AP(1) 1030 FORMAT('0KNOCK-ON ELECTRONS CAN BE CREATED AND ANY ELECTRON FOLLOW *ED DOWN TO' /T40,F8.3,' MeV KINETIC ENERGY'/ ' BREM PHOTONS CAN *BE CREATED AND ANY PHOTON FOLLOWED DOWN TO ', /T40,F8.3,' MeV * ') IXX=int(seed3) CALL HATCH2 return end c----------------------------------------------------------------------- subroutine push c----------------------------------------------------------------------- c h.drescher c parameter (pi=3.14159265) COMMON /ccask1/ ipmax,izmin,ipcmax,izmax,zcemax parameter (mxblk=12) common /coptl/dptl(mxblk) common /cnshow/nshow,oosho,e0,z0,id0(0:16),mshow,ifirst,ifi0 common /cbase/emin,emax,nde,zmin,zmax,dz,costet,aphi,hobs,hground common /ceshow/ e1,e2,egamma,eran common /cbaseb/ theta1,theta2,tran,phi1,phi2 common /cprim/ hprim,xprim,yprim,tprim,ioprim common /curve1/rearth common /cbase2/hinj,haxmax,hax0,haxc,dinj common /caxis/ xsaxis,ysaxis,zsaxis,tsaxis(9),fsaxis(9),iaxis common /ccore/ xcore,ycore,hcore,tcore call OpenStack() if(e2.gt.e1)then !--------- choose e0 --------------- if(eran.ge.0.)then rengy=rangen() if(eran.gt.0.) rengy=eran*int(rengy/eran) elseif(eran.lt.0.)then rengy=float(ishow-1)/max(nshow-1,1) endif r=rengy if(egamma.ne.1.)then g=1.0-egamma e0=(r*(e2**g-e1**g)+e1**g)**(1.0/g) print *,"energy chosen:",e0,"(",e1,e2,")" else e0=e1*(e2/e1)**r print *,"energy chosen:",e0,"(",e1,e2,")" endif call AddToCheckfile("random energy chosen:&",e0) endif if(theta2.gt.theta1) then !---------- choose costet -------------- if(tran.eq.1.0)then r=rangen() c atheta=acos(sqrt(sin(theta1)+r*(sin(theta2)-sin(theta1)))) atheta=acos((cos(theta1)+r*(cos(theta2)-cos(theta1)))) costet=cos(atheta) elseif(tran.eq.2.0)then r=rangen() atheta=acos(sqrt(cos(theta1)**2 $ +r*(cos(theta2)**2-cos(theta1)**2))) costet=cos(atheta) elseif(tran.eq.3.0)then #include "theta.F" endif print *,"theta chosen:",atheta,"(",theta1,theta2,")" print *,"phi chosen:",aphi,"(",phi1,phi2,")" call AddToCheckfile("random theta chosen:&",atheta/pi*180) call AddToCheckfile("costet :&",costet) iaxis=0 endif if(phi2.gt.phi1) then !---------- choose phi -------------- r=rangen() aphi=rangen()*(phi2-phi1)+phi1 call AddToCheckfile("random phi chosen :&",aphi/pi*180) iaxis=0 endif i=1+NINT(LOG10(e0/emin)*nde) if ( i.gt.ipmax ) ipmax=i ! save maximal primary energy j=int((z0-zmin)/dz)+1 if ( j.lt.izmin ) izmin=j ! save minimal depth i=16 do while(id0(i).eq.0) ! choose primary i=i-1 enddo i=int(rangen()*float(i))+1 id=id0(i) id0(0)=id if(id.eq.230)id=-20 if(id.ge.100.and.mod(id,100).eq.0)then dptl(5)=0.938*int(id/100) else call icmass(id,dptl(5)) endif p=sqrt(e0**2-dptl(5)**2) dptl(1)=p*cos(aphi)*sqrt(1.0-costet**2) dptl(2)=p*sin(aphi)*sqrt(1.0-costet**2) dptl(3)=p*costet dptl(4)=e0 c dptl(8)=max(0.0013,z0)/costet ! height in slant depth g/cm^2 c dh=height( dptl(8)*costet ) - hobs ! observation level atheta=acos(costet) if(dinj.eq.0)then ! ioprim==0 #ifdef CURVED if(abs(costet).gt.1e-6)then A=1.0+tan(atheta)**2 B=2.0*rearth-2.0*hcore*tan(atheta)**2 C=-hinj**2-2*rearth*hinj+hcore**2*tan(atheta)**2 ss=(B**2-4.0*A*C) h=(-B+sqrt(ss))/2./A r=h*tan(atheta) else r=sqrt((rearth+hinj)**2-(rearth+hcore)**2) h=hcore endif dh=(h-hcore) dl=sqrt(dh**2+r**2) dptl(6)=xcore-cos(aphi)*r dptl(7)=ycore-sin(aphi)*r dptl(8)=h dptl(9)=-dl #else if(atheta.gt.60.0/180.0*pi)then print *,"warning: flat earth with large angle !!!!" endif if(costet.eq.0.0)then print*, "don''t use theta=90 deg with flat earth" stop endif dh=hinj-hcore dptl(6)=xcore-dh*cos(aphi)*sqrt(1.0-costet**2)/costet ! x dptl(7)=ycore-dh*sin(aphi)*sqrt(1.0-costet**2)/costet ! y dptl(8)=hcore+dh dptl(9)=tcore-dh/costet ! timing #endif else dptl(6)=xcore-dinj*cos(aphi)*sqrt(1.0-costet**2) dptl(7)=ycore-dinj*sin(aphi)*sqrt(1.0-costet**2) dptl(8)=hcore+dinj*costet dptl(9)=tprim-dinj endif dptl(10)=id ! id dptl(11)=1. ! weight dptl(12)=0. ! generation call stput end #ifdef NEXUS2 #else c----------------------------------------------------------------------- real function ranf() c----------------------------------------------------------------------- c uniform random number generator from cern library c----------------------------------------------------------------------- double precision dranf, g900gt, g900st double precision ds(2), dm(2), dseed double precision dx24, dx48 double precision dl, dc, du, dr logical single data ds / 1665 1885.d0, 286 8876.d0 / data dm / 1518 4245.d0, 265 1554.d0 / data dx24 / 1677 7216.d0 / data dx48 / 281 4749 7671 0656.d0 / single = .true. goto 10 entry dranf() single = .false. 10 dl = ds(1) * dm(1) dc = dint(dl/dx24) dl = dl - dc*dx24 du = ds(1)*dm(2) + ds(2)*dm(1) + dc ds(2) = du - dint(du/dx24)*dx24 ds(1) = dl dr = (ds(2)*dx24 + ds(1)) / dx48 if(single) then ranf = sngl(dr) else dranf = dr endif return entry g900gt() g900gt = ds(2)*dx24 + ds(1) return entry g900st(dseed) ds(2) = dint(dseed/dx24) ds(1) = dseed - ds(2)*dx24 g900st = ds(1) return end c----------------------------------------------------------------------- subroutine ranfgt(seed) c----------------------------------------------------------------------- double precision seed, g900gt, g900st, dummy seed = g900gt() return entry ranfst(seed) dummy = g900st(seed) return end c----------------------------------------------------------------------- function rangen() c----------------------------------------------------------------------- c generates a random number c----------------------------------------------------------------------- 1 rangen=ranf() if(rangen.le.0.)goto1 if(rangen.ge.1.)goto1 return entry rangel(xmax) 2 rangel=ranf()*xmax if(rangel.le.0.)goto2 if(rangel.ge.xmax)goto2 return end #endif c----------------------------------------------------------------------- subroutine utrotc(isig,ax,ay,az,x,y,z) c K.Werner c c----------------------------------------------------------------------- c performs a rotation c----------------------------------------------------------------------- c if(az.ge.0.)then rx=ax ry=ay rz=az c else c rx=-ax c ry=-ay c rz=-az c endif if(rz.eq.0..and.ry.eq.0.)then alp=0. stop else c alp=abs(acos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry) alp=atan2(ry,rz) endif c bet= c *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx) bet=atan2(rx,sqrt(rz**2+ry**2)) cosa=cos(alp) sina=sin(alp) cosb=cos(bet) sinb=sin(bet) if(isig.gt.0)then xs=x*cosb-y*sina*sinb-z*cosa*sinb ys= y*cosa -z*sina zs=x*sinb+y*sina*cosb+z*cosa*cosb elseif(isig.lt.0)then xs= x*cosb +z*sinb ys=-x*sinb*sina+y*cosa+z*cosb*sina zs=-x*sinb*cosa-y*sina+z*cosb*cosa endif x=xs y=ys z=zs return end subroutine icchrg(id,chrg) c h.drescher c ida=iabs(id) if(ida.eq.0.or.ida.eq.10)then chrg=0. elseif(ida.eq.1)then chrg=float(id) elseif(ida.eq.12)then chrg=float(sign(1,-id)) elseif(ida.eq.14)then chrg=float(sign(1,-id)) elseif(ida.eq.120)then chrg=float(sign(1,id)) elseif(ida.eq.220)then chrg=0. elseif(ida.eq.110)then chrg=0. elseif(ida.eq.130)then chrg=float(sign(1,id)) elseif(ida.eq.20)then chrg=0. elseif(ida.eq.1120)then chrg=float(sign(1,id)) elseif(ida.eq.1220)then chrg=0. elseif(ida.eq.1130)then chrg=float(sign(1,id)) elseif(ida.eq.2130)then chrg=0. elseif(ida.eq.1230)then chrg=0. elseif(ida.eq.2230)then chrg=-float(sign(1,id)) elseif(id.eq.45)then ! geant code for deut chrg=1.0 elseif(id.eq.46)then ! geant code for ???? chrg=1.0 elseif(id.eq.47)then ! geant code for alpha chrg=2.0 elseif(id.eq.49)then ! geant code for ???? chrg=0.0 print *,"geancode 49 " elseif(mod(id,100).eq.0)then chrg=(id/100)/2 else print *,"did't find charge for",id stop endif end integer function idgeant(idisa) c h.drescher c idgeant=id end integer function idisa(idgeant) ik=idgeant if(ik.eq.1)then id=10 elseif(ik.eq.2)then id=12 elseif(ik.eq.3)then id=-12 elseif(ik.eq.4)then id=13 elseif(ik.eq.5)then id=14 elseif(ik.eq.6)then id=-14 elseif(ik.eq.7)then id=110 elseif(ik.eq.8)then id=120 elseif(ik.eq.9)then id=-120 elseif(ik.eq.11)then id=130 elseif(ik.eq.12)then id=-130 elseif(ik.eq.10)then id=-20 ! KLONG elseif(ik.eq.16)then id=20 ! KSHORT elseif(ik.eq.14)then id=1120 elseif(ik.eq.15)then id=-1120 elseif(ik.eq.13)then id=1220 elseif(ik.eq.25)then id=-1220 elseif(ik.eq.17)then ! eta id=220 elseif(ik.eq.18)then id=2130 elseif(ik.eq.19)then id=1130 elseif(ik.eq.20)then id=1230 elseif(ik.eq.21)then id=2230 elseif(ik.eq.22)then id=1330 elseif(ik.eq.23)then id=2330 elseif(ik.eq.24)then id=3331 elseif(ik.eq.26)then id=-2130 elseif(ik.eq.27)then id=-2230 elseif(ik.eq.28)then id=-1230 elseif(ik.eq.29)then id=--2230 elseif(ik.eq.30)then id=-1330 elseif(ik.eq.31)then id=-2330 elseif(ik.eq.32)then id=-3331 else id=0 endif idisa=id end subroutine icmass(id,amass) c h.drescher c ida=iabs(id) if(ida.eq.0.or.ida.eq.10)then amass=0. elseif(ida.eq.1.or.ida.eq.12)then amass=0.00051099906 elseif(ida.eq.11.or.ida.eq.13)then amass=0. elseif(ida.eq.14)then amass=0.105658389 elseif(ida.eq.120)then amass=.13956995 elseif(ida.eq.110)then amass=.1349764 elseif(ida.eq.220)then amass=.54730 elseif(ida.eq.130)then amass=0.493677 elseif(ida.eq.20)then amass=0.497672 elseif(ida.eq.1120)then amass=0.93827231 elseif(ida.eq.1220)then amass=0.93956563 elseif(ida.eq.1130)then amass=.11894E+01 elseif(ida.eq.1230)then amass=.11925E+01 elseif(ida.eq.2130)then amass=1.115683 elseif(ida.eq.2230)then amass=.11974E+01 elseif(ida.eq.1330)then amass=.13149E+01 elseif(ida.eq.2330)then amass=.13213E+01 elseif(ida.eq.3331)then amass=.16722E+01 else print *,"did''t find mass for",id stop endif end c----------------------------------------------------------------------- subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4) c----------------------------------------------------------------------- c performs a lorentz boost c----------------------------------------------------------------------- real beta(4),z(4) if(p5.le.0.)then write(*,*)'***** mass <= 0.' write(*,*)'p(5): ',p1,p2,p3,p4,p5 stop endif z(1)=x1 z(2)=x2 z(3)=x3 z(4)=x4 beta(1)=-p1/p5 beta(2)=-p2/p5 beta(3)=-p3/p5 beta(4)= p4/p5 bp=0. do 220 k=1,3 220 bp=bp+z(k)*isig*beta(k) do 230 k=1,3 230 z(k)=z(k)+isig*beta(k)*z(4) *+isig*beta(k)*bp/(beta(4)+1.) z(4)=beta(4)*z(4)+bp x1=z(1) x2=z(2) x3=z(3) x4=z(4) return end c----------------------------------------------------------------------- subroutine utlob5(yboost,x) c----------------------------------------------------------------------- c K. Werner c dimension x(5) amt=sqrt(x(5)**2+x(1)**2+x(2)**2) y=sign(1.,x(3))*alog((x(4)+abs(x(3)))/amt) y=y-yboost x(4)=amt*cosh(y) x(3)=amt*sinh(y) return end