*------------------------------------------------------------------------ * VIRTUAL calculating the virtual image of electrons accepted by HMS * -=======- * * Combined Target and HMS monte carlo calculating the virtual image * of the electrons accepted by HMS * - reads the start parameters form STDIN and * - writes the results to an PAW/hbook ntuple file * (same format as written by hms.f,electron.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*256,title*256 DATA tags/'x0','dx0dz','y0','dy0dz','z0','delta0', + 'x', 'dxdz', 'y', 'dydz', 'z'/ DATA label/'QQ','theta','thetaB','E0','dE0', + 'xType','xMin','xMax', + 'yType','yMin','yMax', + 'zType','zMin','zMax', + 'tType','tMin','tMax', + 'dType','dMin','dMax'/ REAL uT(6),u(6),zT,z,SINth,COSth,phi,delta,dzdt REAL tht_e,tht_B,E,dE,Q2 INTEGER un(3),tn,dn REAL umin(3),umax(3),tmin,tmax,dmin,dmax INTEGER ntry,nkept,nlost,events,i LOGICAL lost REAL c,pi180 PARAMETER (c = 30.) PARAMETER (pi180 = 3.141592653/180.) REAL seed,RNDM EXTERNAL RDMIN,RDMOUT,RNDM 10 FORMAT (X,A,$) 20 FORMAT (A) PRINT *,'Virtual: - image monte carlo for Gen' PRINT *,'======== includes target magnetic field' PRINT * PRINT *,'Enter the kinematic parameters:' PRINT *,' (a field orientation of 0 will neglect the field)' PRINT 10,'- Momentum transfer Q^2 [GeV/c]: ' READ (5,*) Q2 PRINT 10,'- Scattering angle [deg]: ' READ (5,*) tht_e PRINT 10,'- Field orientation [deg]: ' READ (5,*) tht_b 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 *,' 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,*) un(1),umin(1),umax(1) PRINT 10,'- in plane coordinate (y-lab) [cm]: ' READ (5,*) un(2),umin(2),umax(2) PRINT 10,'- in beam coordinate (z-lab) [cm]: ' READ (5,*) un(3),umin(3),umax(3) PRINT 10,'- 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 * param (1) = Q2 param (2) = tht_e param (3) = tht_B param (4) = E param (5) = dE param (6) = un(1) param (7) = umin(1) param (8) = umax(1) param (9) = un(2) param (10) = umin(2) param (11) = umax(2) param (12) = un(3) param (13) = umin(3) param (14) = umax(3) 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 = 'virtual' name = name(:INDEX(name,' ')-1)//'.hbook' CALL hlimit(nwpawc) CALL hropen(20,'Electron',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,'//Electron',4096,tags) CALL hbookn(20,'Detected Events',11,'//Electron',4096,tags) IF (ABS(tht_B) .GT. 0.0001) > CALL trgInit ('field_map.dat',tht_B,0.) CALL hmsInitForward ('fwd_maps.dat',9,.FALSE.,E) nkept=0 nlost=0 E = E - dE tmin = COS(tmin*pi180) tmax = COS(tmax*pi180) tht_e = tht_e*pi180 DO ntry=1,events ! crate an electron at the target coordinates DO i=1,3 u(i)=umin(i) IF (un(i) .EQ. 1) THEN u(i)=u(i)+ RNDM(0.) *(umax(i)-u(i)) ELSEIF (un(i) .GT. 1) THEN u(i)=u(i)+INT(un(i)*RNDM(0.))*(umax(i)-u(i))/(un(i)-1) ENDIF ENDDO ! moving within a cone along the spectrometer direction COSth = tmin IF (tn .EQ. 1) THEN COSth=tmin + RNDM(0.)*(tmax-tmin) ELSEIF (tn .GT. 1) THEN COSth=tmin + INT(tn*RNDM(0.))*(tmax-tmin)/(tn-1) ENDIF SINth = SQRT (1.-COSth**2) phi = 2.*3.141592653589*RNDM(0.) u(4)= c * SINth*SIN(phi) ! vx u(5)= c * (-COSth*SIN(tht_e) + SINth*COS(phi)*COS(tht_e)) ! vy u(6)= c * ( COSth*COS(tht_e) - SINth*COS(phi)*SIN(tht_e)) ! 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 ! convert the initial coordinates into spectrometer coordinates event(5) = 0.01*(-SIN(tht_e)*u(2) + COS(tht_e)*u(3)) event(3) = 0.01*( COS(tht_e)*u(2) + SIN(tht_e)*u(3)) event(1) = 0.01* u(1) dzdt = -SIN(tht_e)*u(5) + COS(tht_e)*u(6) event(4) = (COS(tht_e)*u(5) + SIN(tht_e)*u(6))/dzdt event(2) = u(4)/dzdt event(6) = delta ! run the ray tracing IF (ABS(tht_B) .GT. 0.0001) > CALL trgTrack (u,-E*(1+delta)*1000.,1.,100.) ! convert the resulting coordinates into spectrometer coordinates zT = 0.01*(-SIN(tht_e)*u(2) + COS(tht_e)*u(3)) uT(3) = 0.01*( COS(tht_e)*u(2) + SIN(tht_e)*u(3)) uT(1) = 0.01* u(1) dzdt = -SIN(tht_e)*u(5) + COS(tht_e)*u(6) uT(4) = (COS(tht_e)*u(5) + SIN(tht_e)*u(6))/dzdt uT(2) = u(4)/dzdt uT(6) = delta uT(5) = 0. ! run them through HMS CALL hmsAccept (uT(1),zT,event(7),z,lost) event(7) = uT(1) - uT(2)*zT event(8) = uT(2) event(9) = uT(3) - uT(4)*zT event(10) = uT(4) event(11) = zT 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('Electron') STOP END