program main_cwn implicit none include 'const.inc' integer i, j, jj, k, kk, kk2, m, n, ll, trigger, trig28 parameter(m=90, n=500000, ll=20) integer ranno,counter,count28, eventno real xlam, HRNDM1, vect(4) double precision Rx, Rx2, Rz, Rad, aa(m) parameter(Rz=14.0) parameter(Rad=1000.0) double precision xmax, ymax, zmax parameter(xmax=5.5,ymax=4.5,zmax=28.) ! in cm double precision xmin, ymin, zmin parameter(xmin=0.,ymin=0.,zmin=0.) double precision thindex, mum(m), sqm(m), th0, ph0 double precision x0, y0, z0, thx0, thy0 double precision x, y, z, thx, thxp, thy double precision xtest, ytest, ztest, yztest, xztest double precision xx(ll),yy(ll),zz(ll),thxx(ll),thyy(ll),thc(ll) double precision corupt, dx(ll) c----------------------------------- c fill the ntuple 'sm.rz' //// ntuple variables .... integer max_paw, iquest integer istat, icycle parameter(max_paw = 10000000) ! PAW common block size real quest, p common/quest/iquest(100) common/pawc/p(max_paw) real nlam, nx0, ny0, nz0, nthx0, nthy0 real nx, ny, nz, nthx, nthy common /smblk/ i, j, nlam, nx0, ny0, nz0, nthx0, nthy0, & nx, ny, nz, nthx, nthy c HLIMIT (must be called before any other HBOOK routine) call hlimit(max_paw) iquest(10) = 6500000 c HROPEN - open a unit for hist output: call hropen(1,'Myntuple','sm.rz','NQ',8190,istat) c HBNT (10 - identifer of ntuple call hbnt(10,'smnt',' ') c HBNAME - describe varialbles to be stores in CWN, c CALL HBNAME (ID, CHBLOK, VARIABLE, CHFORM), c where: c ID -Identifier of the Ntuple as in the call to HBNT. c CHBLOK - Character variable of maximum length 8 characters c specifying the name by which the block of variables described c VARIABLE - The first variable that is described in CHFORM. Variables c must be in common blocks but may not be in a ZEBRA bank. For c example, given the common block CEXAM described below, one c would call HBNAME with the argument IEVENT. In the case of c character variables, the routine HBNAMC must be used. In all c other cases one should use HBNAME. c CHFORM - Can be either a character string describing the variables to c be stored in block CHBLOK call hbname(10,'sm',i,'i:I,j:I,nlam:R,nx0:R,ny0:R,nz0:R,'// & 'nthx0:R,nthy0:R,nx:R,ny:R,nz:R,nthx:R,nthy:R') c----------------------------------------------------------------------- c CALL HRGET (ID,CHFILE,CHOPT) - c Read a histogram from a given direct access file. This routine cannot c be used for Ntuples. c ID - Histogram identifier. ID=0 read all histograms into the current c directory. c CHFILE - Character variable or constant defining the input filename. C If CHFILE=' ' the histogram is read from the current working c directory on disk. c CHOPT - Character option specifying the desired option. c call hrget(100,'lam.hbook','') c----------------------------------------------------------------------- print *,"asking for seed" read(*,*) ranno; c ranno=7811 c CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ C!!! one 32-bit integer INT and sets Luxury Level LUX ++ C!!! which is integer between zero and MAXLEV, or if ++ C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ C!!! should be set to zero unless restarting at a break++ C!!! point given by output of RLUXAT (see RLUXAT). ++ call rluxgo(3,ranno,0,0) ! RANLUX luxury slope, slopeatc, inter, level 2 c----------------------------------------------------------------------- 11 format(i3,2x,i3,2x,f5.2) c open(15,file="tmpdat",status="new") open(19,file="hits.dat",status="unknown"); c----------------------------------------------------------------------- c see explanation is expl_sm_exec.txt: print *,"asking for Event number" read(*,*) eventno c eventno=464 counter=0 do 100 i = 1, n c if(i.ne.eventno) goto 100 c----------------------------------------------------------------------- c HRNDM1(ID) - returns a random number accordind to the contents c of 1-dimensional ID xlam = HRNDM1(100) ! Randomize in the range of 100 hist c xlam = 5. ! Angstrom c----------------------------------------------------------------------- c CALL RANLUX(RVEC,LEN) - returns a vector RVEC of LEN 32-bit random c floating point numbers in the interval (0,1), c the end points excluded. RVEC is an array of c type REAL and of length LEN at least. call ranlux(vect,4) c write(*,*) (vect(j),j=1,4); if(i.eq.11) stop c c - So this is used to izortopicaly randomize th and phi c----------------------------------------------------------------------- thindex = 0.03*xlam ! theta bigger than 0.03 ! index of refraction for SM C This is very convinient to start from polar angle, bacause of index c of refraction. c ! th0 = thindex*vect(1) ! polar angle ph0 = 2.00*pi*vect(2) ! azimuthal angle (radians) x0 = (xmax-xmin)*vect(3) + xmin ! initial position for x y0 = (ymax-ymin)*vect(4) + ymin ! initial position for y z0 = zmin ! inital position for z c c The all th0,phi0,x0,y0,z0 are independent. Starting position of frame is c x0,y0,z0 c Changing from polar to Cartesian frame: thx0 = datan(dtan(th0)*dcos(ph0)) ! double prec. arctan thy0 = datan(dtan(th0)*dsin(ph0)) if(i.eq.eventno) write(19,*) z0,x0,thx0 c write(22,*) z0,x0,thx0 c c----------------------------------------------------------------------- c Moving particle in y position - to the end of SM (returns: y, thy) call refy(i,y0,thy0,y,thy) c c----------------------------------------------------------------------- c xx(1) = x0 zz(1) = z0 thxx(1) = thx0 do j = 1, ll ! ll=20 (maximal number of hits) jj=j+1 call refx(i,j,xx(j),zz(j),thxx(j),Rad,Rz,zz(jj),xx(jj), & thc(jj),thxx(jj)) if(zz(jj).eq.-100) goto 100 ! hits the edge of the SM if(zz(jj).lt.zz(j)) then counter = counter + 1 c--- goto 100 endif if(i.eq.eventno) write(19,*) zz(jj),xx(jj),thxx(jj) dx(j) = abs(xx(jj)-xx(j)) if(zz(jj).eq.28) then x = xx(jj) z = zz(jj) thxp = thxx(jj) goto 101 elseif(28-zz(jj).lt.zz(jj)-zz(j)) then thxp = thxx(jj) z = 28. x = dtan(thxp)*z + (xx(jj)-dtan(thxp)*zz(jj)) goto 101 endif enddo 101 trigger = 1 ! for end : if(i.eq.eventno) write(*,*) jj if(thxp.gt.pi) then thx = thxp - 2.*pi elseif(thxp.lt.-1*pi) then thx = thxp + 2.*pi else thx = thxp endif if(x.lt.xmin.or.x.gt.xmax) goto 100 write(35,*) i,j,x,z c c c nlam = xlam nx0 = x0 ny0 = y0 nz0 = z0 nthx0 = thx0 nthy0 = thy0 nx = x ny = y nz = z nthx = thx nthy = thy c Write to file : ? call hfnt(10) c c----------------------------------------------------------------------- c 100 continue close(19) c c ------------------------- call hrout(10,icycle,' ') call hrend('Myntuple') c ------------------------- c corupt =100.*float(counter)/float(n) write(*,*) "No. of corrupted events = ",counter ," / ",corupt,"%" c return end