program simc ! Last modified: ! ! This program represents variations on themes by ! N. C. R. Makins (Argonne National Lab), ! T. G. O'Neill (Argonne National Lab), and ! Seemingly Countless Others (Virtually Everywhere). ! implicit none include 'simulate_init.inc' include 'histograms_init.inc' include 'radc.inc' include 'hbook.inc' include 'sos/struct_sos.inc' include 'hms/struct_hms.inc' include 'hrsr/struct_hrsr.inc' include 'hrsl/struct_hrsl.inc' integer*4 i, ninform real*8 r, sum_sigcc real*8 domega_e, domega_p !populated e/hadron solid angles. logical success character filename*80, genfile*80, histfile*80, timestring1*23 character timestring2*23,genifile*80 record /event/ vertex, orig, recon record /event_main/ main record /contrib/ contrib record /histograms/ H record /event_central/ central real*8 one parameter (one=1.0d0) !double precision 1 for subroutine calls real*8 grnd real*8 ang_targ_earm,ang_targ_parm c ! INITIALIZE !----------- ! ... initialize histogram area for HBOOK call hlimit(PawSize) ! ... read in the data base call dbase_read(H) if (debug(3)) write(6,*) 'Main after dbrd: p,e th', > spec.p.theta,spec.e.theta,using_P_arm_montecarlo,using_E_arm_montecarlo ! ... open diagnostic ntuple file, if requested i = index(base,' ') if (Nntu.gt.0) then filename = 'worksim/'//base(1:i-1)//'.rzdat' call NtupleInit(filename) endif ! ... set up some quantities that the radiative corrections routines will need ! ... N.B. Call this even if we're not using radiative corrections -- the ! ... target thicknesses seen by the incoming and ! ... scattered particles (in radiation lengths) are determined here, and ! ... are needed for the multiple scattering calculations ! ... done in the Monte Carlos. if (debug(2)) write(6,*)'sim: calling radc_init' call radc_init if (debug(2)) write(6,*)'sim: done with radc_init' ! ... compute some quantities for a central event call calculate_central(central) if (debug(2)) write(6,*)'sim: done with calc_central' if (debug(2)) write(6,*)'central.sigcc=',central.sigcc if (debug(4)) write(6,*)'sim: at 1' targetfac=targ.mass_amu/3.75914d+6/(targ.abundancy/100.) > *abs(cos(targ.angle))/(targ.thick*1000.) if (debug(4)) write(6,*)'sim: at 2' ! ... Note: EXPER.charge is in mC and the luminosity comes out in fm^-2 ! ... now luminosity is in ub^-1 luminosity = EXPER.charge/targetfac if (debug(4)) write(6,*)'sim: at 3' ! ... initialize the random number generator and number of attempts r = grnd() ntried = 0 ! GAW - insert calls to initialize target field for both arms ! mkj 10-20-2003 add if statements to determine the angle magnitude ! and sign between target direction and e-arm or p-arm. ! if ( abs(targ_Bphi-spec.e.phi) < .01) then if (targ_Bangle .ge. spec.e.theta) then ang_targ_earm = -1*sin(spec.e.phi)*(targ_Bangle-spec.e.theta) else ang_targ_earm = +1*sin(spec.e.phi)*(spec.e.theta-targ_Bangle) endif elseif ( abs(degrad*(targ_Bphi-spec.e.phi)-180) < .01) then ang_targ_earm = +1*sin(spec.e.phi)*(targ_Bangle+spec.e.theta) else write(*,*) ' Error determining angle between target and e-arm' write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad write(*,*) ' Central e-arm angle =',spec.e.theta*degrad,' e-arm phi =',spec.e.theta*degrad endif c if ( abs(targ_Bphi-spec.p.phi) < .01) then if (targ_Bangle .ge. spec.p.theta) then ang_targ_parm = -1*sin(spec.p.phi)*(targ_Bangle-spec.p.theta) else ang_targ_parm = +1*sin(spec.p.phi)*(targ_Bangle-spec.p.theta) endif elseif ( degrad*abs(targ_Bphi-spec.p.phi)-180 < .01) then ang_targ_parm = +1*sin(spec.p.phi)*(targ_Bangle+spec.p.theta) else write(*,*) ' Error determining angle between target and p-arm' write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad write(*,*) ' Central p-arm angle =',spec.p.theta*degrad,' p-arm phi =',spec.p.theta*degrad endif c write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad write(*,*) ' Angle between target and e-arm',ang_targ_earm*degrad write(*,*) ' Angle between target and p-arm',ang_targ_parm*degrad c call trgInit(tgt_field_file,ang_targ_earm*degrad,0., > ang_targ_parm*degrad,0.) ! LOOP OVER SIMULATED EVENTS !--------------------------- write(6,*) ' ' call date (timestring1) call time (timestring1(11:23)) nevent = 0 ninform = 1000 sum_sigcc = 0.0 !sum of main.sigcc (for generating ave sigcc) do while (nevent.lt.abs(ngen)) ! reset decay distance, initial hadron mass (modified if particle decays), ! and some ntuple variables. Mh2_final = Mh2 ntup.radphot = 0.0 ntup.radarm = 0.0 decdist = 0.0 !decay distance. ntup.resfac = 0.0 ntup.sigcm = 0.0 ntup.krel = 0.0 ntup.mm = 0.0 ntup.mmA = 0.0 ntup.t = 0.0 ! Setup for this event ! ... keep the human interested ntried = ntried + 1 if(debug(2)) then write(6,*)'********************************************' write(6,*)'SIM: ntried =',ntried write(6,*)' ncontribute =',ncontribute endif if (mod(ntried,ninform).eq.1) then write (6,'(1x,a,i9,a,i8,a,e11.4)') 'Generating Event', > ntried, ' ... ', nevent,' successes so far - Monitor:', > wtcontribute*luminosity/ntried if (ntried.ge.5000) ninform = 20000 endif ! ... generate an event call generate(main,vertex,orig,success) if(debug(2)) write(6,*)'sim: after gen, success =',success ! Run the event through various manipulations, checking to see whether ! it made it at each stage ! ... run through spectrometers and check SP cuts if(debug(3)) write(6,*)'sim: before mc: orig.p.E =',orig.p.E if(success)call montecarlo(orig,main,recon,success) if(debug(2)) write(6,*)'sim: after mc, success =',success ! ... calculate everything else about the event if (success) then call complete_recon_ev(recon,success) endif if(debug(2)) write(6,*)'sim: after comp_ev, success =',success if(debug(5)) write(6,*) 'recon.Em,recon.Pm',recon.Em,recon.Pm ! ... calculate remaining pieces of the main structure if (success) call complete_main(.false.,main,vertex,recon,success) ! ... Apply SPedge cuts to success if hard_cuts is set. if (hard_cuts) then if(recon.e.delta .le. (SPedge.e.delta.min+slop.MC.e.delta.used) .or. > recon.e.delta .ge. (SPedge.e.delta.max-slop.MC.e.delta.used) .or. > recon.e.yptar .le. (SPedge.e.yptar.min+slop.MC.e.yptar.used) .or. > recon.e.yptar .ge. (SPedge.e.yptar.max-slop.MC.e.yptar.used) .or. > recon.e.xptar .le. (SPedge.e.xptar.min+slop.MC.e.xptar.used) .or. > recon.e.xptar .ge. (SPedge.e.xptar.max-slop.MC.e.xptar.used) .or. > recon.p.delta .le. (SPedge.p.delta.min+slop.MC.p.delta.used) .or. > recon.p.delta .ge. (SPedge.e.delta.max-slop.MC.p.delta.used) .or. > recon.p.yptar .le. (SPedge.p.yptar.min+slop.MC.p.yptar.used) .or. > recon.p.yptar .ge. (SPedge.p.yptar.max-slop.MC.p.yptar.used) .or. > recon.p.xptar .le. (SPedge.p.xptar.min+slop.MC.p.xptar.used) .or. > recon.p.xptar .ge. (SPedge.p.xptar.max-slop.MC.p.xptar.used) ) > success = .false. if (doing_eep .and. recon.Em .gt. cuts.Em.max) success = .false. endif if (success) sum_sigcc = sum_sigcc + main.sigcc ! Em and Pm histos are NEW. Check that they don't suffer from energy ! shifts (e.g. coulomb distortions - see update_range call in limits_update). if(ntried.gt.0)then call inc(H.geni.e.delta, vertex.e.delta, one) call inc(H.geni.e.yptar, vertex.e.yptar, one) call inc(H.geni.e.xptar, -vertex.e.xptar, one) call inc(H.geni.p.delta, vertex.p.delta, one) call inc(H.geni.p.yptar, vertex.p.yptar, one) call inc(H.geni.p.xptar, -vertex.p.xptar, one) call inc(H.geni.Em, vertex.Em, one) call inc(H.geni.Pm, vertex.Pm, one) endif if (success) then if (debug(2)) write(6,*)'sim: We have a success!' ! ... Output ntuple and histograms if passed all cuts, or if not using ! ... SPedge limits as hard cuts. ! ... increment the "experimental" spectrometer arm and "contribution" ! ... histograms call inc(H.RECON.e.delta, main.RECON.e.delta, main.weight) call inc(H.RECON.e.yptar, main.RECON.e.yptar, main.weight) call inc(H.RECON.e.xptar, main.RECON.e.xptar, main.weight) call inc(H.RECON.p.delta, main.RECON.p.delta, main.weight) call inc(H.RECON.p.yptar, main.RECON.p.yptar, main.weight) call inc(H.RECON.p.xptar, main.RECON.p.xptar, main.weight) call inc(H.RECON.Em, recon.Em, one) call inc(H.RECON.Pm, recon.Pm, one) call inc(H.gen.e.delta, vertex.e.delta, one) call inc(H.gen.e.yptar, vertex.e.yptar, one) call inc(H.gen.e.xptar, -vertex.e.xptar, one) call inc(H.gen.p.delta, vertex.p.delta, one) call inc(H.gen.p.yptar, vertex.p.yptar, one) call inc(H.gen.p.xptar, -vertex.p.xptar, one) call inc(H.gen.Em, vertex.Em, one) ! ... update counters and integrated weights FOR EVENTS INSIDE OF DELTA ! population limits (events are only generated to fill region inside ! of the given delta limits) and below cuts.Em.max. Have to remove slop ! that is added in init.f from the delta limits. if (debug(4)) write(6,*)'sim: cut' ncontribute = ncontribute + 1 if (debug(4)) write(6,*)'sim: ncontribute =',ncontribute if (recon.e.delta.ge.(SPedge.e.delta.min+slop.MC.e.delta.used) .and. > recon.e.delta.le.(SPedge.e.delta.max-slop.MC.e.delta.used) .and. > recon.p.delta.ge.(SPedge.p.delta.min+slop.MC.p.delta.used) .and. > recon.p.delta.le.(SPedge.e.delta.max-slop.MC.p.delta.used) .and. > recon.Em.le.cuts.Em.max) then wtcontribute = wtcontribute + main.weight endif if (.not.rad_proton_this_ev) ncontribute_no_rad_proton = > ncontribute_no_rad_proton + 1 ! ... update the "contribution" and "slop" limits call limits_update(main,vertex,orig,recon,doing_deuterium, > doing_pion,doing_kaon,doing_dvcs,contrib,slop) ! ... write out line to diagnostic ntuple file, if requested endif ! if (Nntu.gt.0) call results_ntu_write(main,vertex,orig,recon,success) ! END of LOOP over events ! Cute thing here, if ngen < 0, then just make |ngen| ATTEMPTS rather than ! insisting on ngen good events ... this is suitable for tests with new ! and uncertain parameters, you don't want the program fishing around all ! day for an event! if (ngen.lt.0) then nevent = nevent+1 else if (success) nevent = nevent+1 endif enddo ! ! END GAME: NORMALIZE, GENERATE OUTPUT !------------------------------------- ! Our last event announcement ... and how long did it all take? write (6,'(1x,''---> Last Event '',i9,'' ...'',i9,'' successes'')') ntried, nevent call date (timestring2) call time (timestring2(11:23)) ! NORMALIZE! ! ... put in the luminosity and efficiency factors normfac = luminosity/ntried*nevent ! ... multiply in the relevant phase spaces (see event.f for description ! ... of generated variables. Electron angles for all cases. ! ... add hadron angles and electron energy for all but H(e,e'p). ! ... add hadron energy for A(e,e'p) (what about phase_space?) ! ... Cross sections are all in microbarn/MeV**i/sr**j (i,j depend on reaction) domega_e=(gen.e.yptar.max-gen.e.yptar.min)*(gen.e.xptar.max-gen.e.xptar.min) domega_p=(gen.p.yptar.max-gen.p.yptar.min)*(gen.p.xptar.max-gen.p.xptar.min) genvol = domega_e ! genvol_inclusive = genvol !may want dOmega, or dE*dOmega ! ... 2-fold to 5-fold. if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon 1 .or.doing_dvcs) then genvol = genvol * domega_p * (gen.e.E.max-gen.e.E.min) endif if (doing_heavy) then !6-fold genvol = genvol * (gen.p.E.max-gen.p.E.min) endif normfac = normfac * genvol if (doing_phsp) normfac = 1.0 wtcontribute = wtcontribute*normfac ! Close diagnostic ntuple file, if used if (Nntu.gt.0) call NtupleClose(filename) ! Produce output files 900 if (ngen.eq.0) goto 1000 ! ... the diagnostic histograms i = index(base,' ') genfile = ' ' write(genfile,'(a,''.gen'')') 'outfiles/'//base(1:i-1) genifile = ' ' write(genifile,'(a,''.geni'')') 'outfiles/'//base(1:i-1) histfile = ' ' write(histfile,'(a,''.hist'')') 'outfiles/'//base(1:i-1) write(6,'(1x,''... writing '',a)') genfile write(6,*) 'Come back soon!' open(unit=7, file=genifile, status='unknown') if (electron_arm.eq.1 .or. hadron_arm.eq.1) then write(7,*) 'HMS Trials: ',hSTOP.trials write(7,*) 'Slit hor/vert/corners ',hSTOP.slit_hor,hSTOP.slit_vert,hSTOP.slit_oct write(7,*) 'Q1 entrance/mid/exit ',hSTOP.Q1_in,hSTOP.Q1_mid,hSTOP.Q1_out write(7,*) 'Q2 entrance/mid/exit ',hSTOP.Q2_in,hSTOP.Q2_mid,hSTOP.Q2_out write(7,*) 'Q3 entrance/mid/exit ',hSTOP.Q3_in,hSTOP.Q3_mid,hSTOP.Q3_out write(7,*) 'Dipole entrance/exit ',hSTOP.D1_in,hSTOP.D1_out write(7,*) 'Events reaching hut ',hSTOP.hut write(7,*) 'DC1, DC2, Scin, Cal ',hSTOP.dc1,hSTOP.dc2,hSTOP.scin,hSTOP.cal write(7,*) 'Successes ',hSTOP.successes write(7,*) endif if (electron_arm.eq.2 .or. hadron_arm.eq.2) then write(7,*) 'SOS Trials: ',sSTOP.trials write(7,*) 'Slit hor/vert/corners ',sSTOP.slit_vert,sSTOP.slit_hor,sSTOP.slit_oct write(7,*) 'Quad entrance/mid/exit',sSTOP.quad_in,sSTOP.quad_mid,sSTOP.quad_out write(7,*) 'D1 entrance/exit ',sSTOP.bm01_in,sSTOP.bm01_out write(7,*) 'D2 entrance/exit ',sSTOP.bm02_in,sSTOP.bm02_out write(7,*) 'Vacuum exit ',sSTOP.exit write(7,*) 'Events reaching hut ',sSTOP.hut write(7,*) 'DC1, DC2, Scin, Cal ',sSTOP.dc1,sSTOP.dc2,sSTOP.scin,sSTOP.cal write(7,*) 'Successes ',sSTOP.successes write(7,*) endif if (electron_arm.eq.3 .or. hadron_arm.eq.3) then write(7,*) 'HRSr Trials: ',rSTOP.trials write(7,*) 'Slit hor/vert ',rSTOP.slit_vert,rSTOP.slit_hor write(7,*) 'Q1 entrance/mid/exit ',rSTOP.Q1_in,rSTOP.Q1_mid,rSTOP.Q1_out write(7,*) 'Q2 entrance/mid/exit ',rSTOP.Q2_in,rSTOP.Q2_mid,rSTOP.Q2_out write(7,*) 'Dipole entrance/exit ',rSTOP.D1_in,rSTOP.D1_out write(7,*) 'Q3 entrance/mid/exit ',rSTOP.Q3_in,rSTOP.Q3_mid,rSTOP.Q3_out write(7,*) 'Events reaching hut ',rSTOP.hut write(7,*) 'VDC1, VDC2 ',rSTOP.dc1,rSTOP.dc2 write(7,*) 'S1, S2, Cal ',rSTOP.s1,rSTOP.s2,rSTOP.cal write(7,*) endif if (electron_arm.eq.4 .or. hadron_arm.eq.4) then write(7,*) 'HRSl Trials: ',lSTOP.trials write(7,*) 'Slit hor/vert ',lSTOP.slit_vert,lSTOP.slit_hor write(7,*) 'Q1 entrance/mid/exit ',lSTOP.Q1_in,lSTOP.Q1_mid,lSTOP.Q1_out write(7,*) 'Q2 entrance/mid/exit ',lSTOP.Q2_in,lSTOP.Q2_mid,lSTOP.Q2_out write(7,*) 'Dipole entrance/exit ',lSTOP.D1_in,lSTOP.D1_out write(7,*) 'Q3 entrance/mid/exit ',lSTOP.Q3_in,lSTOP.Q3_mid,lSTOP.Q3_out write(7,*) 'Events reaching hut ',lSTOP.hut write(7,*) 'VDC1, VDC2 ',lSTOP.dc1,lSTOP.dc2 write(7,*) 'S1, S2, Cal ',lSTOP.s1,lSTOP.s2,lSTOP.cal endif close(7) open(unit=7, file=genfile, status='unknown') write(7,*) 'E arm Experimental Target Distributions:' write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.RECON.e.delta.min+(i-0.5)*H.RECON.e.delta.bin, > H.RECON.e.delta.buf(i), H.RECON.e.yptar.min+(i-0.5)* > H.RECON.e.yptar.bin, H.RECON.e.yptar.buf(i), > H.RECON.e.xptar.min+(i-0.5)*H.RECON.e.xptar.bin, > H.RECON.e.xptar.buf(i) enddo write(7,*) 'P arm Experimental Target Distributions:' write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.RECON.p.delta.min+(i-0.5)*H.RECON.p.delta.bin, > H.RECON.p.delta.buf(i), H.RECON.p.yptar.min+(i-0.5)* > H.RECON.p.yptar.bin, H.RECON.p.yptar.buf(i), > H.RECON.p.xptar.min+(i-0.5)*H.RECON.p.xptar.bin, > H.RECON.p.xptar.buf(i) enddo write(7,*) 'Distributions of Contributing E arm Events:' write(7,'(6a12)') 'delta','CONTRIB','yuptar','CONTRIB','xptar','CONTRIB' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.gen.e.delta.min+(i-0.5)*H.gen.e.delta.bin, > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5)* > H.gen.e.yptar.bin, H.gen.e.yptar.buf(i), > H.gen.e.xptar.min+(i-0.5)*H.gen.e.xptar.bin, H.gen.e.xptar.buf(i) enddo write(7,*) 'Distributions of Contributing P arm Events:' write(7,'(6a12)') 'delta','CONTRIB','yptar','CONTRIB','xptar','CONTRIB' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.gen.p.delta.min+(i-0.5)*H.gen.p.delta.bin, > H.gen.p.delta.buf(i), H.gen.p.yptar.min+(i-0.5)* > H.gen.p.yptar.bin, H.gen.p.yptar.buf(i), > H.gen.p.xptar.min+(i-0.5)*H.gen.p.xptar.bin, H.gen.p.xptar.buf(i) enddo write(7,*) 'Original E arm Events:' write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.gen.e.delta.min+(i-0.5)*H.gen.e.delta.bin, > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5)* > H.gen.e.yptar.bin, H.geni.e.yptar.buf(i), > H.gen.e.xptar.min+(i-0.5)*H.gen.e.xptar.bin, H.geni.e.xptar.buf(i) enddo write(7,*) 'Original P arm Events:' write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN' do i=1,nHbins write(7,'(3(1x,2(e11.4,1x)))') > H.geni.p.delta.min+(i-0.5)*H.geni.p.delta.bin, > H.geni.p.delta.buf(i), H.geni.p.yptar.min+(i-0.5)* > H.geni.p.yptar.bin, H.geni.p.yptar.buf(i), > H.geni.p.xptar.min+(i-0.5)*H.geni.p.xptar.bin, H.geni.p.xptar.buf(i) enddo write(7,*) 'Original Em/Pm distributions:' write(7,'(4a12)') 'Em','ORIGIN','Pm','ORIGIN' do i=1,nHbins write(7,'(3(1x,4(e11.4,1x)))') > H.geni.Em.min+(i-0.5)*H.geni.Em.bin, > H.geni.Em.buf(i), H.geni.Pm.min+(i-0.5)* > H.geni.Pm.bin, H.geni.Pm.buf(i) enddo close(7) ! Counters, etc report to the screen and to the diagnostic histogram file ! call report(6,timestring1,timestring2,central,contrib,sum_sigcc) open(unit=7,file=histfile,status='unknown') call report(7,timestring1,timestring2,central,contrib,sum_sigcc) close(unit=7) 1000 end !------------------------------------------------------------------- subroutine inc(hist,val,weight) implicit none include 'histograms.inc' integer*4 ibin real*8 val, weight record /hist_entry/ hist ibin= nint(0.5+(val-hist.min)/hist.bin) if(ibin.ge.1.and.ibin.le.nHbins)then hist.buf(ibin) = hist.buf(ibin) + weight endif return end !-------------------------------------------------------------------- subroutine report(iun,timestring1,timestring2,central,contrib,sum_sigcc) implicit none include 'simulate.inc' include 'radc.inc' include 'brem.inc' integer*4 iun real*8 sum_sigcc character*23 timestring1, timestring2 record /contrib/ contrib record /event_central/ central ! Output of diagnostics ! Run time write(iun,'(/1x,''BEGIN Time: '',a23)') timestring1 write(iun,'(1x,''END Time: '',a23)') timestring2 ! Kinematics write(iun,*) 'KINEMATICS:' if (doing_eep) then if (doing_hyd_elast) then write(iun,*) ' ****-------- H(e,e''p) --------****' else if (doing_deuterium) then write(iun,*) ' ****-------- D(e,e''p) --------****' else if (doing_heavy) then write(iun,*) ' ****-------- A(e,e''p) --------****' else stop 'I don''t have ANY idea what (e,e''p) we''re doing!!!' endif else if (doing_pion) then if (doing_hydpi) then if (targ.A .eq. 1) then write(iun,*) ' ****-------- H(e,e''pi) --------****' else if (targ.A .ge.3) then write(iun,*) ' ****-------- A(e,e''pi) --------****' endif else if (doing_deutpi) then write(iun,*) ' ****-------- D(e,e''pi) --------****' else if (doing_hepi) then write(iun,*) ' ****-------- A(e,e''pi) --------****' else stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!' endif if (which_pion .eq. 0) then write(iun,*) ' ****---- Default Final State ----****' else if (which_pion .eq. 10) then write(iun,*) ' ****---- Final State is A+pi ----****' endif else if (doing_dvcs) then if (doing_hyddvcs) then if (targ.A .eq. 1) then write(iun,*) ' ****-------- H(e,e''gamma) --------****' else if (targ.A .ge.3) then write(iun,*) ' ****-------- A(e,e''gamma) --------****' endif else if (doing_deutdvcs) then write(iun,*) ' ****-------- D(e,e''gamma) --------****' else if (doing_hedvcs) then write(iun,*) ' ****-------- A(e,e''gamma) --------****' else stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!' endif if (which_dvcs .eq. 0) then write(iun,*) ' ****---- Default Final State ----****' else if (which_dvcs .eq. 10) then write(iun,*) ' ****---- Final State is A+gamma ----****' endif else if (doing_kaon) then if (doing_hydkaon) then if (targ.A .eq. 1) then write(iun,*) ' ****-------- H(e,e''K) --------****' else if (targ.A .ge.3) then write(iun,*) ' ****-------- A(e,e''K) --------****' endif else if (doing_deutkaon) then write(iun,*) ' ****-------- D(e,e''K) --------****' else if (doing_hekaon) then write(iun,*) ' ****-------- A(e,e''K) --------****' else stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!' endif if (which_kaon.eq.0) then write(iun,*) ' ****---- producing a LAMBDA ----****' else if (which_kaon.eq.1) then write(iun,*) ' ****---- producing a SIGMA0 ----****' else if (which_kaon.eq.2) then write(iun,*) ' ****---- producing a SIGMA- ----****' else if (which_kaon.eq.10) then write(iun,*) ' ****---- WITH BOUND LAMBDA ----****' else if (which_kaon.eq.11) then write(iun,*) ' ****---- WITH BOUND SIGMA0 ----****' else if (which_kaon.eq.12) then write(iun,*) ' ****---- WITH BOUND SIGMA- ----****' else stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!' endif else if (doing_phsp) then write(iun,*) ' ****--- PHASE SPACE - NO physics, NO radiation (may not work)---****' else stop 'I don''t have ANY idea what we''re doing!!!' endif write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'Ebeam', Ebeam, 'MeV' write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') > '(dE/E)beam', dEbeam/Ebeam, '(full wid)' write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'x-width', gen.xwid, 'cm' write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'y-width', gen.ywid, 'cm' write(iun,'(9x,a12,'' = '',i15,2x,a16)') > 'fr_pattern',targ.fr_pattern, '1=square,2=circ' write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr1', targ.fr1, 'cm' write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr2', targ.fr2, 'cm' write(iun,*) ' ' write(iun,'(9x,18x,2(a15,2x))') '____E arm____','____P arm____' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'angle', > spec.e.theta*degrad, spec.p.theta*degrad, 'deg' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'momentum', > spec.e.p, spec.p.p, 'MeV/c' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'x offset', > spec.e.offset.x, spec.p.offset.x, 'cm' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'y offset', > spec.e.offset.y, spec.p.offset.y, 'cm' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'z offset', > spec.e.offset.z, spec.p.offset.z, 'cm' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'xptar offset', > spec.e.offset.xptar, spec.p.offset.xptar, 'mr' write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'yptar offset', > spec.e.offset.yptar, spec.p.offset.yptar, 'mr' write(iun,'(6x,''Central Values: '',a5,'' = '',f15.4,2x,a9)') > 'Q2', central.Q2/1.d6, '(GeV/c)^2' write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'q', central.q, 'MeV/c' write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'omega', > central.omega, 'MeV' write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'Em', central.Em, 'MeV' write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'Pm', central.Pm, 'MeV/c' ! Target write(iun,*) 'TARGET specs:' 9911 format(2x,2(5x,a10,' = ',e12.6,1x,a5)) write(iun,9911) 'A', targ.A, ' ', 'Z', targ.Z, ' ' write(iun,9911) 'mass', targ.mass_amu, 'amu', 'mass', targ.M,'MeV' write(iun,9911) 'Mrec', targ.mrec_amu, 'amu', 'Mrec', targ.Mrec,'MeV' write(iun,9911) 'Mtar_struck',targ.Mtar_struck,'MeV', > 'Mrec_struck',targ.Mrec_struck,'MeV' write(iun,9911) 'rho', targ.rho, 'g/cm3', 'thick', targ.thick,'g/cm2' write(iun,9911) 'angle', targ.angle*degrad, 'deg', 'abundancy', > targ.abundancy, '%' write(iun,9911) 'X0', targ.X0, 'g/cm2', 'X0_cm', targ.X0_cm,'cm' write(iun,9911) 'length',targ.length,'cm','zoffset',targ.zoffset,'cm' write(iun,9911) 'xoffset',targ.xoffset,'cm','yoffset',targ.yoffset,'cm' write(iun,'(t12,3a15)') '__ave__', '__lo__', '__hi__' 9912 format(1x,a15,3f15.5,2x,a6) write(iun,9912) 'Coulomb', targ.Coulomb.ave, targ.Coulomb.min, > targ.Coulomb.max, 'MeV' write(iun,9912) 'Eloss_beam', targ.Eloss(1).ave, > targ.Eloss(1).min, targ.Eloss(1).max, 'MeV' write(iun,9912) 'Eloss_e', targ.Eloss(2).ave, targ.Eloss(2).min, > targ.Eloss(2).max, 'MeV' write(iun,9912) 'Eloss_p', targ.Eloss(3).ave, targ.Eloss(3).min, > targ.Eloss(3).max, 'MeV' write(iun,9912) 'teff_beam', targ.teff(1).ave, targ.teff(1).min, > targ.teff(1).max, 'radlen' write(iun,9912) 'teff_e', targ.teff(2).ave, targ.teff(2).min, > targ.teff(2).max, 'radlen' write(iun,9912) 'teff_p', targ.teff(3).ave, targ.teff(3).min, > targ.teff(3).max, 'radlen' 9913 format(1x,a15,t25,f15.5,2x,a6) write(iun,9913) 'musc_nsig_max', targ.musc_nsig_max, ' ' write(iun,9913) 'musc_max_beam', targ.musc_max(1)*1000., 'mr' write(iun,9913) 'musc_max_e', targ.musc_max(2)*1000., 'mr' write(iun,9913) 'musc_max_p', targ.musc_max(3)*1000., 'mr' ! Flags write(iun,*) 'FLAGS:' write(iun,'(5x,2(2x,a19,''='',l2))') 'doing_dvcs', doing_dvcs, > 'which_dvcs', which_dvcs write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_eep', doing_eep, > 'doing_kaon', doing_kaon, 'doing_pion', doing_pion write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_phsp', doing_phsp, > 'which_pion', which_pion, 'which_kaon', which_kaon write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hyd_elast', doing_hyd_elast, > 'doing_deuterium', doing_deuterium, 'doing_heavy', doing_heavy write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hyddvcs', doing_hyddvcs, > 'doing_deutdvcs', doing_deutdvcs, 'doing_hedvcs', doing_hedvcs write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydpi', doing_hydpi, > 'doing_deutpi', doing_deutpi, 'doing_hepi', doing_hepi write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydkaon', doing_hydkaon, > 'doing_deutkaon', doing_deutkaon, 'doing_hekaon', doing_hekaon write(iun,'(5x,3(2x,a19,''='',l2))') 'mc_smear', mc_smear, > 'electron_arm', electron_arm, 'hadron_arm', hadron_arm write(iun,'(5x,3(2x,a19,''='',l2)))') 'using_Eloss', using_Eloss, > 'using_Coulomb',using_Coulomb,'deForest_flag',deForest_flag write(iun,'(5x,3(2x,a19,''='',l2)))') 'correct_Eloss', correct_Eloss, > 'correct_raster',correct_raster, 'doing_decay', doing_decay write(iun,'(5x,2(2x,a19,''='',l2))') > 'using_E_arm_montecarlo', using_E_arm_montecarlo, > 'using_P_arm_montecarlo', using_P_arm_montecarlo write(iun,'(7x,a11,''='',f10.3,a4)') 'ctau',ctau,'cm' ! Counters write(iun,*) 'COUNTERS:' write(iun,'(12x,''Ngen (request) = '',i10)') ngen write(iun,'(12x,''Ntried = '',i10)') ntried write(iun,'(12x,''Ncontribute = '',i10)') ncontribute write(iun,'(12x,''Nco_no_rad_prot= '',i10)') ncontribute_no_rad_proton write(iun,'(12x,''-> %no_rad_prot= '',f10.3)') > (100.*ncontribute_no_rad_proton/max(dfloat(ncontribute),0.1)) write(iun,'(/1x,''INTEGRATED WEIGHTS (number of counts in delta/Em cuts!):'')') write(iun,'(1x,'' MeV: wtcontr= '',e16.8)') wtcontribute/nevent ! Radiative corrections write(iun,*) 'RADIATIVE CORRECTIONS:' if (.not.using_rad) then write(iun,'(x,a14,''='',l3)') 'using_rad', using_rad else write(iun,'(4(x,a14,''='',l3))') 'use_expon',use_expon, > 'include_hard',include_hard,'calc_spence',calculate_spence write(iun,'(2(x,a14,''='',l3))') 'using_rad', using_rad, > 'use_offshell_rad', use_offshell_rad write(iun,'(4(x,a14,''='',i3))') 'rad_flag',rad_flag, > 'extrad_flag', extrad_flag, 'one_tail', one_tail, > 'intcor_mode', intcor_mode write(iun,'(x,a14,''='',f11.3)') 'dE_edge_test',dE_edge_test write(iun,'(x,a14,''='',f11.3)') 'Egamma_max', Egamma_tot_max 9914 format(1x,a18,' = ',4f11.3) write(iun,*) 'Central Values:' write(iun,9914) 'hardcorfac', central.rad.hardcorfac write(iun,9914) 'etatzai', central.rad.etatzai write(iun,9914) 'frac(1:3)', central.rad.frac write(iun,9914) 'lambda(1:3)', central.rad.lambda write(iun,9914) 'bt(1:2)', central.rad.bt write(iun,9914) 'c_int(0:3)', central.rad.c_int write(iun,9914) 'c_ext(0:3)', central.rad.c_ext write(iun,9914) 'c(0:3)', central.rad.c write(iun,9914) 'g_int', central.rad.g_int write(iun,9914) 'g_ext', central.rad.g_ext write(iun,9914) 'g(0:3)', central.rad.g endif ! Miscellaneous write(iun,*) 'MISCELLANEOUS:' 9915 format(12x,a14,' = ',e16.6,1x,a6) write(iun,*) 'Note that central.sigcc is for central delta,theta,phi in both spectrometers' write(iun,*) ' and may give non-physical kinematics, esp. for Hydrogen' write(iun,*) 'Note also that AVE.sigcc is really AVER.weight (the two arenot exactly equal)' write(iun,9915) 'CENTRAL.sigcc', central.sigcc, ' ' write(iun,9915) 'AVERAGE.sigcc', sum_sigcc/nevent, ' ' write(iun,9915) 'charge', EXPER.charge, 'mC' write(iun,9915) 'targetfac', targetfac, ' ' write(iun,9915) 'luminosity', luminosity, 'ub^-1' write(iun,9915) 'luminosity', luminosity*(hbarc/100000.)**2, 'GeV^2' write(iun,9915) 'genvol', genvol, ' ' write(iun,9915) 'normfac', normfac, ' ' 9916 format(12x,a14,' = ',f6.1,' to ',f6.1,' in ',i4,' bins of ',f7.2,1x,a5) if (doing_heavy) then write(iun,'(12x,''Theory file: '',a)') > theory_file(1:index(theory_file,' ')-1) endif ! Input spectrometer limits. write(iun,*) 'Input Spectrometer Limits:' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.delta.min/max', > SPedge.e.delta.min,SPedge.e.delta.max,'%' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.yptar.min/max', > SPedge.e.yptar.min,SPedge.e.yptar.max,'rad' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.xptar.min/max', > SPedge.e.xptar.min,SPedge.e.xptar.max,'rad' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.delta.min/max', > SPedge.p.delta.min,SPedge.p.delta.max,'%' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.yptar.min/max', > SPedge.p.yptar.min,SPedge.p.yptar.max,'rad' write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.xptar.min/max', > SPedge.p.xptar.min,SPedge.p.xptar.max,'rad' ! Edges used in generation and checking, as well as range of events found ! to contribute within those edges ! ... on GENERATED qties write(iun,*) 'Limiting VERTEX values (vertex.e/p.*,Em,Pm,Trec)' write(iun,*) ' USED limits are gen.e/p.*, and VERTEXedge.Em,Pm,Trec' write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') '______used______', > '_____found______' write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 9917 format(1x,a18,t21,2f12.3,t50,2f10.3,2x,a5) write(iun,9917) 'E arm delta', gen.e.delta.min, gen.e.delta.max, > contrib.gen.e.delta.lo, contrib.gen.e.delta.hi, '%' write(iun,9917) 'E arm yptar', gen.e.yptar.min*1000., > gen.e.yptar.max*1000., > contrib.gen.e.yptar.lo*1000., contrib.gen.e.yptar.hi*1000.,'mr' write(iun,9917) 'E arm xptar', gen.e.xptar.min*1000., > gen.e.xptar.max*1000., > contrib.gen.e.xptar.lo*1000., contrib.gen.e.xptar.hi*1000.,'mr' write(iun,9917) 'P arm delta', gen.p.delta.min, gen.p.delta.max, > contrib.gen.p.delta.lo, contrib.gen.p.delta.hi, '%' write(iun,9917) 'P arm yptar', gen.p.yptar.min*1000., > gen.p.yptar.max*1000., > contrib.gen.p.yptar.lo*1000., contrib.gen.p.yptar.hi*1000.,'mr' write(iun,9917) 'P arm xptar', gen.p.xptar.min*1000., > gen.p.xptar.max*1000., > contrib.gen.p.xptar.lo*1000., contrib.gen.p.xptar.hi*1000.,'mr' write(iun,9917) 'sumEgen', gen.sumEgen.min, gen.sumEgen.max, > contrib.gen.sumEgen.lo, contrib.gen.sumEgen.hi, 'MeV' ! ... on VERTEX qties ! write(iun,*) 'Limiting VERTEX values:' ! write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') ! > '______used______', '_____found______' ! write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' write(iun,9917) 'Trec', VERTEXedge.Trec.min, VERTEXedge.Trec.max, > contrib.vertex.Trec.lo, contrib.vertex.Trec.hi, 'MeV' write(iun,9917) 'Em', VERTEXedge.Em.min, VERTEXedge.Em.max, > contrib.vertex.Em.lo, contrib.vertex.Em.hi, 'MeV' write(iun,9917) 'Pm', VERTEXedge.Pm.min, VERTEXedge.Pm.max, > contrib.vertex.Pm.lo, contrib.vertex.Pm.hi, 'MeV/c' if ((doing_deuterium .or. doing_pion .or. doing_kaon .or. doing_dvcs) .and. using_rad) > write(iun,*) ' *** NOTE: sumEgen.min only used in GENERATE_RAD' ! ... on TRUE qties write(iun,*) 'Limiting ORIGINAL values: orig.e/p.*,Em,Pm,Trec (no edge.* limits for Pm,Trec)' write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') > '______used______', '_____found______' write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' write(iun,9917) 'E arm E', edge.e.E.min, edge.e.E.max, > contrib.tru.e.E.lo, contrib.tru.e.E.hi, 'MeV' write(iun,9917) 'E arm yptar', edge.e.yptar.min*1000., > edge.e.yptar.max*1000., > contrib.tru.e.yptar.lo*1000., contrib.tru.e.yptar.hi*1000.,'mr' write(iun,9917) 'E arm xptar', edge.e.xptar.min*1000., edge.e.xptar.max*1000., > contrib.tru.e.xptar.lo*1000., contrib.tru.e.xptar.hi*1000., 'mr' write(iun,9917) 'P arm E', edge.p.E.min, edge.p.E.max, > contrib.tru.p.E.lo, contrib.tru.p.E.hi, 'MeV' write(iun,9917) 'P arm yptar', edge.p.yptar.min*1000., > edge.p.yptar.max*1000., > contrib.tru.p.yptar.lo*1000., contrib.tru.p.yptar.hi*1000.,'mr' write(iun,9917) 'P arm xptar', edge.p.xptar.min*1000., edge.p.xptar.max*1000., > contrib.tru.p.xptar.lo*1000., contrib.tru.p.xptar.hi*1000., 'mr' write(iun,9917) 'Em', max(-999999.999d0,edge.Em.min), min(999999.999d0,edge.Em.max), > contrib.tru.Em.lo, contrib.tru.Em.hi, 'MeV' write(iun,9917) 'Pm', 0., 0., contrib.tru.Pm.lo, > contrib.tru.Pm.hi, 'MeV' write(iun,9917) 'Trec', 0., 0., contrib.tru.Trec.lo, > contrib.tru.Trec.hi, 'MeV' ! ... on SPECTROMETER qties write(iun,*) 'Limiting SPECTROMETER values' write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') > '______used______', '_____found______' write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' write(iun,9917) 'E arm delta', SPedge.e.delta.min, > SPedge.e.delta.max, contrib.SP.e.delta.lo, > contrib.SP.e.delta.hi, '%' write(iun,9917) 'E arm yptar', SPedge.e.yptar.min*1000., > SPedge.e.yptar.max*1000., > contrib.SP.e.yptar.lo*1000., contrib.SP.e.yptar.hi*1000., 'mr' write(iun,9917) 'E arm xptar', SPedge.e.xptar.min*1000., > SPedge.e.xptar.max*1000., > contrib.SP.e.xptar.lo*1000., contrib.SP.e.xptar.hi*1000., 'mr' write(iun,9917) 'P arm delta', SPedge.p.delta.min, > SPedge.p.delta.max, contrib.SP.p.delta.lo, > contrib.SP.p.delta.hi, '%' write(iun,9917) 'P arm yptar', SPedge.p.yptar.min*1000., > SPedge.p.yptar.max*1000., > contrib.SP.p.yptar.lo*1000., contrib.SP.p.yptar.hi*1000., 'mr' write(iun,9917) 'P arm xptar', SPedge.p.xptar.min*1000., > SPedge.p.xptar.max*1000., > contrib.SP.p.xptar.lo*1000., contrib.SP.p.xptar.hi*1000., 'mr' ! ... on RADIATION qties if (using_rad) then write(iun,*) 'Limiting RADIATION values CONTRIBUTING to the (Em,Pm) distributions:' write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') > '______used______', '_____found______' write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' write(iun,9917) 'Egamma(1)', 0., Egamma1_max, > contrib.rad.Egamma(1).lo, contrib.rad.Egamma(1).hi, 'MeV' write(iun,9917) 'Egamma(2)', 0., Egamma2_max, contrib.rad.Egamma(2).lo, > contrib.rad.Egamma(2).hi, 'MeV' write(iun,9917) 'Egamma(3)', 0., Egamma3_max, contrib.rad.Egamma(3).lo, > contrib.rad.Egamma(3).hi, 'MeV' write(iun,9917) 'Egamma_total', 0., Egamma_tot_max, > contrib.rad.Egamma_total.lo, contrib.rad.Egamma_total.hi,'MeV' endif ! ... on slops write(iun,*) 'ACTUAL and LIMITING SLOP values used/obtained:' write(iun,'(t25,3a10)') '__used__', '__min__', '__max__' 9918 format(1x,a10,a12,t25,3f10.3,2x,a5) write(iun,9918) 'slop.MC ', 'E arm delta', slop.MC.e.delta.used, > slop.MC.e.delta.lo, slop.MC.e.delta.hi, '%' write(iun,9918) ' ', 'E arm yptar', slop.MC.e.yptar.used*1000., > slop.MC.e.yptar.lo*1000., slop.MC.e.yptar.hi*1000., 'mr' write(iun,9918) ' ', 'E arm xptar', slop.MC.e.xptar.used*1000., > slop.MC.e.xptar.lo*1000., slop.MC.e.xptar.hi*1000., 'mr' write(iun,9918) ' ', 'P arm delta', slop.MC.p.delta.used, > slop.MC.p.delta.lo, slop.MC.p.delta.hi, '%' write(iun,9918) ' ', 'P arm yptar', slop.MC.p.yptar.used*1000., > slop.MC.p.yptar.lo*1000., slop.MC.p.yptar.hi*1000., 'mr' write(iun,9918) ' ', 'P arm xptar', slop.MC.p.xptar.used*1000., > slop.MC.p.xptar.lo*1000., slop.MC.p.xptar.hi*1000., 'mr' write(iun,9918) 'slop.total', 'Em', slop.total.Em.used, > slop.total.Em.lo, slop.total.Em.hi, 'MeV' write(iun,9918) ' ', 'Pm', 0., slop.total.Pm.lo, > slop.total.Pm.hi, 'MeV/c' write(iun,'(/)') return end !----------------------------------------------------------------------- subroutine calculate_central(central) implicit none include 'simulate.inc' include 'radc.inc' integer i logical success record /event_main/ main record /event/ vertex0, recon0 record /event_central/ central ! Pop in generation values for central event if (debug(2)) write(6,*)'calc_cent: entering...' main.target.x = 0.0 main.target.y = 0.0 main.target.z = 0.0 vertex0.Ein = Ebeam_vertex_ave main.target.Coulomb = targ.Coulomb.ave main.target.Eloss(1) = targ.Eloss(1).ave main.target.teff(1) = targ.teff(1).ave vertex0.e.delta = 0.0 vertex0.e.yptar = 0.0 vertex0.e.xptar = 0.0 vertex0.e.theta = spec.e.theta vertex0.e.phi = spec.e.phi vertex0.e.P = spec.e.P*(1.+vertex0.e.delta/100.) vertex0.e.E = vertex0.e.P vertex0.p.delta = 0.0 vertex0.p.yptar = 0.0 vertex0.p.xptar = 0.0 vertex0.p.theta = spec.p.theta vertex0.p.phi = spec.p.phi vertex0.p.P = spec.p.P*(1.+vertex0.p.delta/100.) vertex0.p.E = sqrt(vertex0.p.P**2 + Mh2) ! Complete_recon_ev for vertex0. Note: complete_recon_ev doesn't ! call radc_init_ev or calculate teff(2,3). ! JRA: Do we want complete_ev or complete_recon_ev? Do we want to calculate ! and/or dump other central values (for pion/kaon production, for example). call complete_recon_ev(vertex0,success) if (debug(2)) write(6,*)'calc_cent: done with comp_ev' if (.not.success) stop 'COMPLETE_EV failed trying to complete a CENTRAL event!' central.Q2 = vertex0.Q2 central.q = vertex0.q central.omega = vertex0.omega central.Em = vertex0.Em central.Pm = vertex0.Pm main.target.teff(2) = targ.teff(2).ave main.target.teff(3) = targ.teff(3).ave if (debug(2)) write(6,*)'calc_cent: calling radc_init_ev' if (debug(4)) write(6,*)'calc_cent: Ein =',vertex0.Ein call radc_init_ev(main,vertex0) if (debug(2)) write(6,*)'calc_cent: done with radc_init_ev' if (using_rad) then central.rad.hardcorfac = hardcorfac central.rad.etatzai= etatzai central.rad.g_int = g_int central.rad.g_ext = g_ext do i = 0, 3 if (i.gt.0) central.rad.frac(i) = frac(i) if (i.gt.0) central.rad.lambda(i) = lambda(i) if (i.gt.0.and.i.lt.3) central.rad.bt(i) = bt(i) central.rad.c_int(i) = c_int(i) central.rad.c_ext(i) = c_ext(i) central.rad.c(i) = c(i) central.rad.g(i) = g(i) enddo endif if (debug(4)) write(6,*)'calc_cent: at 1' do i = 1, neventfields recon0.all(i) = vertex0.all(i) enddo call complete_main(.true.,main,vertex0,recon0,success) central.sigcc = main.sigcc if (debug(2)) write(6,*)'calc_cent: ending...' return end !------------------------------------------------------------------------- subroutine montecarlo(orig,main,recon,success) implicit none include 'simulate.inc' * Track-coordinate and spectrometer common blocks real*8 x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm,delta_E_arm real*8 x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm,delta_P_arm real*8 xfp, yfp, dxfp, dyfp real*8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2) logical success logical ok_E_arm, ok_P_arm record /event/ orig, recon record /event_main/ main real*8 beta_electron parameter (beta_electron = 1.) real*8 tmpfact real*8 fry !fast raster y position. real*8 frx !fast raster x position - left as looking downstream real*8 m2 !mass for call to mc_hms(sos). Changes if decay real*8 pathlen real*8 ztemp real*8 zero parameter (zero=0.0d0) !double precision zero for subroutine calls ! Prepare the event for the Monte Carlo's and/or spectrometer cuts success = .false. ntup.resfac = 0.0 !resfac (see simulate.inc) c if (correct_raster) then fry = -main.target.rastery frx = main.target.rasterx c else c fry = 0.0 c frx = 0.0 c endif !BEAM MULTIPLE SCATTERING: NOT YET IMPLEMENTED, JUST GETTING IT READY! ! ... multiple scattering of the beam. Generate angles for beam deflection, ! ... and then apply those angles to the electron and proton. ! ... The yptar offset goes directly to yptar of both arms (small angle approx). ! ... xptar is multiplied by cos(theta) to get the xptar offsets. if (mc_smear) then call target_musc(orig.Ein,beta_electron,main.target.teff(1),dang_in) else dang_in(1)=0.0 !yp offset, goes directly to yp of both arms dang_in(2)=0.0 !xp offset, mult. by cos_th for xp of both arms endif if (debug(3)) write(6,*) 'mc: using p,e arm mc =', > using_P_arm_montecarlo,using_E_arm_montecarlo !____ P arm ________________ ! Go from TRUE to SPECTROMETER quantities by computing target distortions ! ... ionization loss correction (if requested) if (using_Eloss) then if (debug(3)) write(6,*)'mc: p arm stuff0 =', > orig.p.E,main.target.Eloss(3),spec.p.P main.SP.p.delta = (sqrt(abs((orig.p.E-main.target.Eloss(3))**2 > -Mh2))-spec.p.P) / spec.p.P*100. else if (debug(3)) write(6,*)'mc: p arm stuff1 =',orig.p.delta main.SP.p.delta = orig.p.delta endif ! ... multiple scattering if (mc_smear) then beta = orig.p.p/orig.p.E call target_musc(orig.p.p, beta, main.target.teff(3), dangles) else dangles(1)=0.0 dangles(2)=0.0 endif main.SP.p.yptar = orig.p.yptar + dangles(1) + dang_in(1) main.SP.p.xptar = orig.p.xptar + dangles(2) + dang_in(2)*spec.p.cos_th ! CASE 1: Using the spectrometer Monte Carlo if (using_P_arm_montecarlo) then ! ... change to P arm spectrometer coordinates (TRANSPORT system), if (abs(cos(spec.p.phi)).gt.0.0001) then !phi not at +/- pi/2 write(6,*) 'y_P_arm, z_P_arm will be incorrect if spec.p.phi <> pi/2 or 3*pi/2' write(6,*) 'spec.p.phi=',spec.p.phi,'=',spec.p.phi*180/pi,'degrees' endif delta_P_arm = main.SP.p.delta x_P_arm = -main.target.y y_P_arm = -main.target.x*spec.p.cos_th - main.target.z*spec.p.sin_th*sin(spec.p.phi) z_P_arm = main.target.z*spec.p.cos_th + main.target.x*spec.p.sin_th*sin(spec.p.phi) ! ... Apply spectrometer offset (using spectrometer coordinate system). x_P_arm = x_P_arm - spec.p.offset.x y_P_arm = y_P_arm - spec.p.offset.y z_P_arm = z_P_arm - spec.p.offset.z dx_P_arm = main.SP.p.xptar - spec.p.offset.xptar dy_P_arm = main.SP.p.yptar - spec.p.offset.yptar c write(6,*) 'pion BEFORE target field:' c write(6,*) 'xp, yp:', dx_P_arm, dy_P_arm ! GAW -project to z=0 to compare with reconstructed target positions c in_P_arm.z = y_P_arm - z_P_arm*dy_P_arm c write(6,*) 'using_tgt_field',using_tgt_field c c c write(*,*) 'sign_hms_part =' ,sign_hms_part if (using_tgt_field) then ! do target field tracking - GAW c if (debug(6)) then c write(*,*) '------------------------------' c write(*,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_P_arm c endif call track_from_tgt(x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm, > sign_hms_part*spec.p.P*(1+delta_P_arm/100.),Mh,1,ok_P_arm) endif ! GAW - end 99/11/3 ! ........ drift this position back to z=0, the plane through the target center x_P_arm = x_P_arm - z_P_arm*dx_P_arm y_P_arm = y_P_arm - z_P_arm*dy_P_arm z_P_arm = 0.0 c call print_coord2('virtual track z=0: ',x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm) ! GAW ! ... Monte Carlo through the P arm! if we succeed, continue ... ! ... Here's what's passed: ! spectrometer central momentum ! spectrometer central theta angle ! momentum delta ! x, y, z, theta, and phi in TRANSPORT coordinates ! x, y, theta, and phi at the focal plane ! radiation length of target (NOT USED!) ! particle mass (modified if the hadron decays) ! multiple scattering flag ! wcs smearing flag ! decay flag ! resmult=resolution factor for DCs (see simulate.inc) ! vertical position at target (for reconstruction w/raster MEs). ! ok flag m2 = Mh2 pathlen = 0.0 if (hadron_arm.eq.1) then call mc_hms(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > m2, mc_smear, mc_smear, doing_decay, > ntup.resfac, frx, fry, ok_P_arm, pathlen, using_tgt_field > ,sign_hms_part) else if (hadron_arm.eq.2) then call mc_sos(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > m2, mc_smear, mc_smear, doing_decay, > ntup.resfac, fry, ok_P_arm, pathlen) else if (hadron_arm.eq.3) then call mc_hrsr(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > m2, mc_smear, mc_smear, doing_decay, > ntup.resfac, fry, ok_P_arm, pathlen) else if (hadron_arm.eq.4) then call mc_hrsl(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > m2, mc_smear, mc_smear, doing_decay, > ntup.resfac, fry, ok_P_arm, pathlen) else if (hadron_arm.eq.5) then call mc_calo(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > m2, mc_smear, mc_smear, doing_decay, > ntup.resfac, frx, fry, ok_P_arm, pathlen, using_tgt_field, > main.target.z*cos(spec.p.theta),ztemp,drift_to_cal > ) endif if (.not.ok_P_arm) then if (debug(3)) write(6,*)'mc: ok_P_arm =',ok_P_arm return endif c call print_coord2('recon tgt: ',x_P_arm,y_P_arm,0.,dx_P_arm,dy_P_arm) ! GAW main.RECON.p.delta = delta_P_arm main.RECON.p.yptar = dy_P_arm main.RECON.p.xptar = dx_P_arm main.RECON.p.z = y_P_arm main.FP.p.x = xfp main.FP.p.dx = dxfp main.FP.p.y = yfp main.FP.p.dy = dyfp main.FP.p.path = pathlen ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the ! OUT record else main.RECON.p.delta = main.SP.p.delta main.RECON.p.yptar = main.SP.p.yptar main.RECON.p.xptar = main.SP.p.xptar endif ! Fill recon quantities. recon.p.delta = main.RECON.p.delta recon.p.yptar = main.RECON.p.yptar recon.p.xptar = main.RECON.p.xptar recon.p.z = main.RECON.p.z recon.p.P = spec.p.P*(1.+recon.p.delta/100.) recon.p.E = sqrt(recon.p.P**2 + Mh2) call physics_angles(spec.p.theta,spec.p.phi,recon.p.xptar, > recon.p.yptar,recon.p.theta,recon.p.phi) ! ... correct for energy loss - use most probable (last flag = 4) if (correct_Eloss) then call trip_thru_target (3, zero, recon.p.E, > recon.p.theta, eloss_P_arm, r,Mh,4) recon.p.E = recon.p.E + eloss_P_arm recon.p.E = max(recon.p.E,sqrt(Mh2+0.000001)) !can get P~0 when calculating hadron momentum-->P<0 after eloss recon.p.P = sqrt(recon.p.E**2-Mh2) recon.p.delta = (recon.p.P-spec.p.P)/spec.p.P*100. endif !___ E arm ______________ ! Go from TRUE to SPECTROMETER quantities by computing target distortions ! ... ionization loss correction (if requested) and Coulomb deceleration if (debug(3)) write(6,*)'mc: e arm stuff0 =', > orig.e.E,main.target.Eloss(2),spec.e.P main.SP.e.delta=100.* (orig.e.E - main.target.Eloss(2) > - main.target.Coulomb - spec.e.P) / spec.e.P ! ... multiple scattering if (mc_smear) then call target_musc(orig.e.p, beta_electron, main.target.teff(2), dangles) else dangles(1)=0.0 dangles(2)=0.0 endif main.SP.e.yptar = orig.e.yptar + dangles(1) + dang_in(1) main.SP.e.xptar = orig.e.xptar + dangles(2) + dang_in(2)*spec.p.cos_th ! CASE 1: Using the spectrometer Monte Carlo if (using_E_arm_montecarlo) then ! ... change to E arm spectrometer coordinates (TRANSPORT system), if (abs(cos(spec.e.phi)).gt.0.0001) then !phi not at +/- pi/2 write(6,*) 'y_E_arm, z_E_arm will be incorrect if spec.e.phi <> pi/2 or 3*pi/2' write(6,*) 'spec.e.phi=',spec.e.phi,'=',spec.e.phi*180/pi,'degrees' endif delta_E_arm = main.SP.e.delta x_E_arm = -main.target.y y_E_arm = -main.target.x*spec.e.cos_th - main.target.z*spec.e.sin_th*sin(spec.e.phi) z_E_arm = main.target.z*spec.e.cos_th + main.target.x*spec.e.sin_th*sin(spec.e.phi) ! ... Apply spectrometer offset (using spectrometer coordinate system). x_E_arm = x_E_arm - spec.e.offset.x y_E_arm = y_E_arm - spec.e.offset.y z_E_arm = z_E_arm - spec.e.offset.z dx_E_arm = main.SP.e.xptar - spec.e.offset.xptar dy_E_arm = main.SP.e.yptar - spec.e.offset.yptar ! GAW -project to z=0 to compare with reconstructed target positions c in_E_arm.z = y_E_arm - z_E_arm*dy_E_arm if (using_tgt_field) then ! do target field tracking - GAW c if (debug(6)) then c write(*,*) '------------------------------' c write(*,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_E_arm c endif c write(6,*) 'electron BEFORE target field:' c write(6,*) 'xp, yp:', dx_E_arm, dy_E_arm call track_from_tgt(x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm, > -spec.e.P*(1+delta_E_arm/100.),Me,-1,ok_E_arm) c write(6,*) 'electron AFTER target field:' c write(6,*) 'xp, yp:', dx_E_arm, dy_E_arm endif ! GAW - end 99/11/3 ! ........ drift this position back to z=0, the plane through the target center x_E_arm = x_E_arm - z_E_arm*dx_E_arm y_E_arm = y_E_arm - z_E_arm*dy_E_arm z_E_arm = 0.0 c call print_coord2('virtual track z=0: ',x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm) ! GAW main.SP.e.z=y_E_arm ! ... Monte Carlo through the E arm! if we succeed, continue ... ! ... Here's what's passed: ! spectrometer central momentum ! spectrometer central theta angle ! momentum delta ! x, y, z, theta, and phi in TRANSPORT coordinates ! x, y, theta, and phi at the focal plane ! radiation length of target (NOT USED!) ! particle mass ! multiple scattering flag ! wcs smearing flag ! resmult=resolution factor for DCs (see simulate.inc) ! decay flag (hardwired to .false. for electron arm). ! ok flag m2 = me2 pathlen = 0.0 if (electron_arm.eq.1) then call mc_hms(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > me2, mc_smear, mc_smear, .false., > tmpfact, frx, fry, ok_E_arm, pathlen,using_tgt_field) else if (electron_arm.eq.2) then call mc_sos(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > me2, mc_smear, mc_smear, .false., > tmpfact, fry, ok_E_arm, pathlen) else if (electron_arm.eq.3) then call mc_hrsr(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > me2, mc_smear, mc_smear, .false., > tmpfact, fry, ok_E_arm, pathlen) else if (electron_arm.eq.4) then call mc_hrsl(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > me2, mc_smear, mc_smear, .false., > tmpfact, fry, ok_E_arm, pathlen) else if (electron_arm.eq.5) then ztemp = (recon.p.z/sin(spec.p.theta)) call mc_calo(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > me2, mc_smear, mc_smear, .false., > tmpfact, frx,fry, ok_E_arm, pathlen,using_tgt_field, > main.target.z,ztemp,drift_to_cal) endif ntup.resfac=ntup.resfac+tmpfact if (.not.ok_E_arm) then if (debug(3)) write(6,*)'mc: ok_E_arm =',ok_E_arm return endif c call print_coord2('recon tgt: ',x_E_arm,y_E_arm,0.,dx_E_arm,dy_E_arm) ! GAW main.RECON.e.delta = delta_E_arm main.RECON.e.yptar = dy_E_arm main.RECON.e.xptar = dx_E_arm main.RECON.e.z = y_E_arm main.FP.e.x = xfp main.FP.e.dx = dxfp main.FP.e.y = yfp main.FP.e.dy = dyfp main.FP.e.path = pathlen ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the ! OUT record else main.RECON.e.delta = main.SP.e.delta main.RECON.e.yptar = main.SP.e.yptar main.RECON.e.xptar = main.SP.e.xptar endif ! Fill recon quantities. recon.e.delta = main.RECON.e.delta recon.e.yptar = main.RECON.e.yptar recon.e.xptar = main.RECON.e.xptar recon.e.z = main.RECON.e.z recon.e.P = spec.e.P*(1.+recon.e.delta/100.) recon.e.E = recon.e.P call physics_angles(spec.e.theta,spec.e.phi,recon.e.xptar, 1 recon.e.yptar,recon.e.theta,recon.e.phi) ! ... correct for energy loss and correct for Coulomb deceleration if (correct_Eloss) then call trip_thru_target (2, zero, recon.e.E, recon.e.theta, & eloss_E_arm, r, Me, 4) recon.e.E = recon.e.E + eloss_E_arm endif recon.e.E = recon.e.E + targ.Coulomb.ave recon.e.P = recon.e.E recon.e.delta = (recon.e.P-spec.e.P)/spec.e.P*100. ! Made it! success = .true. return end