*------------------------------------------------------------------------ * HMS tracking of an electron through the HMS * -===- * * HMS monte carlo (including the HMS acceptance) * - reads the start parameters form STDIN and * - writes the results to an PAW/hbook ntuple file * (same format as written by hms.f) * * written by Markus Muehlbauer *------------------------------------------------------------------------ IMPLICIT NONE INTEGER nwpawc,icycle,istat PARAMETER (nwpawc = 1000000) REAL paw(nwpawc) COMMON /pawc/paw CHARACTER*8 tags(11), label(20) REAL event(12),param(20) CHARACTER name*80,title*256 DATA tags/'x0','dx0dz','y0','dy0dz','z0','delta0', + 'x', 'dxdz', 'y', 'dydz', 'z'/ DATA label/'QQ','--','--','E0','dE0', + 'xType','xMin','xMax', + 'yType','yMin','yMax', + 'pType','pMin','pMax', + 'tType','tMin','tMax', + 'dType','dMin','dMax'/ REAL uT(6),u(6),zT,z,th,ph,delta REAL Q2,E,dE INTEGER xn,yn,pn,tn,dn REAL xmin,xmax,ymin,ymax REAL pmin,pmax,tmin,tmax,dmin,dmax INTEGER ntry,nkept,nlost,events LOGICAL lost REAL c PARAMETER (c = 30.) REAL seed,RNDM EXTERNAL RDMIN,RDMOUT,RNDM 10 FORMAT (X,A,$) 20 FORMAT (A) PRINT *,'HMS - monte carlo for Gen' PRINT *,'=== (excludes target magnetic field)' PRINT * PRINT *,'Enter the kinematic parameters:' PRINT 10,'- Momentum transfer Q^2 [GeV/c]: ' READ (5,*) E PRINT 10,'- Beam energy [GeV/c^2]: ' READ (5,*) E PRINT 10,'- Energy loss [GeV/c^2]: ' READ (5,*) dE PRINT * PRINT *,'Enter the simulation parameters:' PRINT *,'(coordinates are HMS target coordinates):' PRINT *,' fixed value: 0,value,0' PRINT *,' uniformly distributed in [min,max] : 1,min,max' PRINT *,' n discrete values in [min,max]: n,min,max' PRINT * PRINT 10,'- out of plane coordinate (x-lab) [cm]: ' READ (5,*) xn,xmin,xmax PRINT 10,'- in plane coordinate (y-lab) [cm]: ' READ (5,*) yn,ymin,ymax PRINT 10,'- out of plane opening of the scattering cone [deg]: ' READ (5,*) pn,pmin,pmax PRINT 10,'- in plane opening of the scattering cone [deg]: ' READ (5,*) tn,tmin,tmax PRINT 10,'- relative momentum displacement [%]: ' READ (5,*) dn,dmin,dmax PRINT * PRINT 10,'Enter the number of events: ' READ (5,*) events PRINT * PRINT 10,'Enter the name of the output file: ' READ (5,20) name PRINT 10,'Enter an output file description: ' READ (5,20) title PRINT * Q2 = 0 param (1) = Q2 param (2) = 0 param (3) = 0 param (4) = E param (5) = dE param (6) = xn param (7) = xmin param (8) = xmax param (9) = yn param (10) = ymin param (11) = ymax param (12) = pn param (13) = pmin param (14) = pmax param (15) = tn param (16) = tmin param (17) = tmax param (18) = dn param (19) = dmin param (20) = dmax seed=2.718 CALL RDMOUT(seed) CALL RDMIN(seed) IF (name .EQ. ' ') name = 'hms' name = name(:INDEX(name,' ')-1)//'.hbook' CALL hlimit(nwpawc) CALL hropen(20,'HMS',name,'NX',4096,ISTAT) CALL htitle(title) CALL hbook1(1,'Parameter',20,0.5,19.5,0.) CALL hlabel(1,20,label,'N') CALL hpak(1,param) CALL hbookn(10,'Lost Events',11,'//HMS',4096,tags) CALL hbookn(20,'Detected Events',11,'//HMS',4096,tags) CALL hmsInitForward ('fwd_maps.dat',9,.FALSE.,E) nkept=0 nlost=0 E = E - dE DO ntry=1,events ! crate an electron (HMS target coordinates) u(1) = xmin IF (xn .EQ. 1) THEN u(1)=xmin+ RNDM(0.) *(xmax-xmin) ELSEIF (xn .GT. 1) THEN u(1)=xmin+INT(xn*RNDM(0.))*(xmax-xmin)/(xn-1) ENDIF u(2) = ymin IF (yn .EQ. 1) THEN u(2)=ymin+ RNDM(0.) *(ymax-ymin) ELSEIF (yn .GT. 1) THEN u(2)=ymin+INT(yn*RNDM(0.))*(ymax-ymin)/(yn-1) ENDIF u(3) = 0 ! moving within a cone along the spectrometer direction ph = pmin IF (tn .EQ. 1) THEN ph=pmin + RNDM(0.)*(pmax-pmin) ELSEIF (pn .GT. 1) THEN ph=pmin + INT(pn*RNDM(0.))*(pmax-pmin)/(pn-1) ENDIF th = tmin IF (tn .EQ. 1) THEN th=tmin + RNDM(0.)*(tmax-tmin) ELSEIF (tn .GT. 1) THEN th=tmin + INT(tn*RNDM(0.))*(tmax-tmin)/(tn-1) ENDIF u(4)= c * sind(ph) ! vx u(5)= c * cosd(ph)*sind(th) ! vy u(6)= c * cosd(ph)*cosd(th) ! vz ! with a momentum displacement delta delta=dmin IF (dn .EQ. 1) THEN delta=dmin + RNDM(0.)*(dmax-dmin) ELSEIF (dn .GT. 1) THEN delta=dmin + INT(dn*RNDM(0.))*(dmax-dmin)/(dn-1) ENDIF delta = 0.01*delta ! store the initial coordinates event(1) = 0.01 * u(1) event(3) = 0.01 * u(2) event(5) = 0.01 * u(3) event(2) = u(4)/u(6) event(4) = u(5)/u(6) event(6) = delta ! run them through HMS uT(1) = 0.01 * u(1) uT(3) = 0.01 * u(2) zT = 0.01 * u(3) uT(2) = u(4)/u(6) uT(4) = u(5)/u(6) uT(6) = delta uT(5) = 0. CALL hmsAccept (uT(1),zT,event(7),z,lost) IF (.not. lost) * CALL hmsForward (uT(1),zT,event(7),z,lost) event(11) = z IF (lost) THEN nlost = nlost+1 CALL hfn (10,event) ELSE nkept = nkept+1 CALL hfn (20,event) ENDIF IF (mod(ntry,1000) .EQ. 0) print *,ntry,nkept,nlost ENDDO CALL hrout(0,icycle,' ') CALL hrend('HMS') STOP END