1 jones 1.1 program simc
2
3 ! Last modified:
4 !
5 ! This program represents variations on themes by
6 ! N. C. R. Makins (Argonne National Lab),
7 ! T. G. O'Neill (Argonne National Lab), and
8 ! Seemingly Countless Others (Virtually Everywhere).
9 !
10 implicit none
11 include 'simulate_init.inc'
12 include 'histograms_init.inc'
13 include 'radc.inc'
14 include 'hbook.inc'
15 include 'sos/struct_sos.inc'
16 include 'hms/struct_hms.inc'
17 include 'hrsr/struct_hrsr.inc'
18 include 'hrsl/struct_hrsl.inc'
19
20 integer*4 i, ninform
21 real*8 r, sum_sigcc
22 jones 1.1 real*8 domega_e, domega_p !populated e/hadron solid angles.
23 logical success
24 character filename*80, genfile*80, histfile*80, timestring1*23
25 character timestring2*23,genifile*80
26 record /event/ vertex, orig, recon
27 record /event_main/ main
28 record /contrib/ contrib
29 record /histograms/ H
30 record /event_central/ central
31
32 real*8 one
33 parameter (one=1.0d0) !double precision 1 for subroutine calls
34
35 real*8 grnd
36 real*8 ang_targ_earm,ang_targ_parm
37 c
38
39 ! INITIALIZE
40 !-----------
41
42 ! ... initialize histogram area for HBOOK
43 jones 1.1
44 call hlimit(PawSize)
45
46 ! ... read in the data base
47
48 call dbase_read(H)
49
50 if (debug(3)) write(6,*) 'Main after dbrd: p,e th',
51 > spec.p.theta,spec.e.theta,using_P_arm_montecarlo,using_E_arm_montecarlo
52
53 ! ... open diagnostic ntuple file, if requested
54
55 i = index(base,' ')
56 if (Nntu.gt.0) then
57 filename = 'worksim/'//base(1:i-1)//'.rzdat'
58 call NtupleInit(filename)
59 endif
60
61 ! ... set up some quantities that the radiative corrections routines will need
62 ! ... N.B. Call this even if we're not using radiative corrections -- the
63 ! ... target thicknesses seen by the incoming and
64 jones 1.1 ! ... scattered particles (in radiation lengths) are determined here, and
65 ! ... are needed for the multiple scattering calculations
66 ! ... done in the Monte Carlos.
67
68 if (debug(2)) write(6,*)'sim: calling radc_init'
69 call radc_init
70 if (debug(2)) write(6,*)'sim: done with radc_init'
71
72 ! ... compute some quantities for a central event
73
74 call calculate_central(central)
75 if (debug(2)) write(6,*)'sim: done with calc_central'
76 if (debug(2)) write(6,*)'central.sigcc=',central.sigcc
77 if (debug(4)) write(6,*)'sim: at 1'
78
79 targetfac=targ.mass_amu/3.75914d+6/(targ.abundancy/100.)
80 > *abs(cos(targ.angle))/(targ.thick*1000.)
81 if (debug(4)) write(6,*)'sim: at 2'
82
83 ! ... Note: EXPER.charge is in mC and the luminosity comes out in fm^-2
84 ! ... now luminosity is in ub^-1
85 jones 1.1
86 luminosity = EXPER.charge/targetfac
87 if (debug(4)) write(6,*)'sim: at 3'
88
89 ! ... initialize the random number generator and number of attempts
90
91 r = grnd()
92 ntried = 0
93
94 ! GAW - insert calls to initialize target field for both arms
95 ! mkj 10-20-2003 add if statements to determine the angle magnitude
96 ! and sign between target direction and e-arm or p-arm.
97 !
98 if ( abs(targ_Bphi-spec.e.phi) < .01) then
99 if (targ_Bangle .ge. spec.e.theta) then
100 ang_targ_earm = -1*sin(spec.e.phi)*(targ_Bangle-spec.e.theta)
101 else
102 ang_targ_earm = +1*sin(spec.e.phi)*(spec.e.theta-targ_Bangle)
103 endif
104 elseif ( abs(degrad*(targ_Bphi-spec.e.phi)-180) < .01) then
105 ang_targ_earm = +1*sin(spec.e.phi)*(targ_Bangle+spec.e.theta)
106 jones 1.1 else
107 write(*,*) ' Error determining angle between target and e-arm'
108 write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad
109 write(*,*) ' Central e-arm angle =',spec.e.theta*degrad,' e-arm phi =',spec.e.theta*degrad
110 endif
111 c
112 if ( abs(targ_Bphi-spec.p.phi) < .01) then
113 if (targ_Bangle .ge. spec.p.theta) then
114 ang_targ_parm = -1*sin(spec.p.phi)*(targ_Bangle-spec.p.theta)
115 else
116 ang_targ_parm = +1*sin(spec.p.phi)*(targ_Bangle-spec.p.theta)
117 endif
118 elseif ( degrad*abs(targ_Bphi-spec.p.phi)-180 < .01) then
119 ang_targ_parm = +1*sin(spec.p.phi)*(targ_Bangle+spec.p.theta)
120 else
121 write(*,*) ' Error determining angle between target and p-arm'
122 write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad
123 write(*,*) ' Central p-arm angle =',spec.p.theta*degrad,' p-arm phi =',spec.p.theta*degrad
124 endif
125 c
126 write(*,*) ' targ_Bangle =',targ_Bangle*degrad,' targ_Bphi =',targ_Bphi*degrad
127 jones 1.1 write(*,*) ' Angle between target and e-arm',ang_targ_earm*degrad
128 write(*,*) ' Angle between target and p-arm',ang_targ_parm*degrad
129 c
130
131 call trgInit(tgt_field_file,ang_targ_earm*degrad,0.,
132 > ang_targ_parm*degrad,0.)
133
134
135 ! LOOP OVER SIMULATED EVENTS
136 !---------------------------
137 write(6,*) ' '
138 call date (timestring1)
139 call time (timestring1(11:23))
140 nevent = 0
141 ninform = 1000
142 sum_sigcc = 0.0 !sum of main.sigcc (for generating ave sigcc)
143
144 do while (nevent.lt.abs(ngen))
145
146 ! reset decay distance, initial hadron mass (modified if particle decays),
147 ! and some ntuple variables.
148 jones 1.1
149 Mh2_final = Mh2
150 ntup.radphot = 0.0
151 ntup.radarm = 0.0
152 decdist = 0.0 !decay distance.
153 ntup.resfac = 0.0
154 ntup.sigcm = 0.0
155 ntup.krel = 0.0
156 ntup.mm = 0.0
157 ntup.mmA = 0.0
158 ntup.t = 0.0
159
160 ! Setup for this event
161
162 ! ... keep the human interested
163
164 ntried = ntried + 1
165
166 if(debug(2)) then
167 write(6,*)'********************************************'
168 write(6,*)'SIM: ntried =',ntried
169 jones 1.1 write(6,*)' ncontribute =',ncontribute
170 endif
171
172 if (mod(ntried,ninform).eq.1) then
173 write (6,'(1x,a,i9,a,i8,a,e11.4)') 'Generating Event',
174 > ntried, ' ... ', nevent,' successes so far - Monitor:',
175 > wtcontribute*luminosity/ntried
176 if (ntried.ge.5000) ninform = 20000
177 endif
178
179 ! ... generate an event
180
181 call generate(main,vertex,orig,success)
182 if(debug(2)) write(6,*)'sim: after gen, success =',success
183
184 ! Run the event through various manipulations, checking to see whether
185 ! it made it at each stage
186
187 ! ... run through spectrometers and check SP cuts
188
189 if(debug(3)) write(6,*)'sim: before mc: orig.p.E =',orig.p.E
190 jones 1.1 if(success)call montecarlo(orig,main,recon,success)
191 if(debug(2)) write(6,*)'sim: after mc, success =',success
192
193 ! ... calculate everything else about the event
194
195 if (success) then
196 call complete_recon_ev(recon,success)
197 endif
198 if(debug(2)) write(6,*)'sim: after comp_ev, success =',success
199 if(debug(5)) write(6,*) 'recon.Em,recon.Pm',recon.Em,recon.Pm
200
201 ! ... calculate remaining pieces of the main structure
202
203 if (success) call complete_main(.false.,main,vertex,recon,success)
204
205 ! ... Apply SPedge cuts to success if hard_cuts is set.
206
207 if (hard_cuts) then
208 if(recon.e.delta .le. (SPedge.e.delta.min+slop.MC.e.delta.used) .or.
209 > recon.e.delta .ge. (SPedge.e.delta.max-slop.MC.e.delta.used) .or.
210 > recon.e.yptar .le. (SPedge.e.yptar.min+slop.MC.e.yptar.used) .or.
211 jones 1.1 > recon.e.yptar .ge. (SPedge.e.yptar.max-slop.MC.e.yptar.used) .or.
212 > recon.e.xptar .le. (SPedge.e.xptar.min+slop.MC.e.xptar.used) .or.
213 > recon.e.xptar .ge. (SPedge.e.xptar.max-slop.MC.e.xptar.used) .or.
214 > recon.p.delta .le. (SPedge.p.delta.min+slop.MC.p.delta.used) .or.
215 > recon.p.delta .ge. (SPedge.e.delta.max-slop.MC.p.delta.used) .or.
216 > recon.p.yptar .le. (SPedge.p.yptar.min+slop.MC.p.yptar.used) .or.
217 > recon.p.yptar .ge. (SPedge.p.yptar.max-slop.MC.p.yptar.used) .or.
218 > recon.p.xptar .le. (SPedge.p.xptar.min+slop.MC.p.xptar.used) .or.
219 > recon.p.xptar .ge. (SPedge.p.xptar.max-slop.MC.p.xptar.used) )
220 > success = .false.
221 if (doing_eep .and. recon.Em .gt. cuts.Em.max) success = .false.
222 endif
223
224 if (success) sum_sigcc = sum_sigcc + main.sigcc
225
226 ! Em and Pm histos are NEW. Check that they don't suffer from energy
227 ! shifts (e.g. coulomb distortions - see update_range call in limits_update).
228
229 if(ntried.gt.0)then
230 call inc(H.geni.e.delta, vertex.e.delta, one)
231 call inc(H.geni.e.yptar, vertex.e.yptar, one)
232 jones 1.1 call inc(H.geni.e.xptar, -vertex.e.xptar, one)
233 call inc(H.geni.p.delta, vertex.p.delta, one)
234 call inc(H.geni.p.yptar, vertex.p.yptar, one)
235 call inc(H.geni.p.xptar, -vertex.p.xptar, one)
236 call inc(H.geni.Em, vertex.Em, one)
237 call inc(H.geni.Pm, vertex.Pm, one)
238 endif
239
240 if (success) then
241 if (debug(2)) write(6,*)'sim: We have a success!'
242
243 ! ... Output ntuple and histograms if passed all cuts, or if not using
244 ! ... SPedge limits as hard cuts.
245
246 ! ... increment the "experimental" spectrometer arm and "contribution"
247 ! ... histograms
248 call inc(H.RECON.e.delta, main.RECON.e.delta, main.weight)
249 call inc(H.RECON.e.yptar, main.RECON.e.yptar, main.weight)
250 call inc(H.RECON.e.xptar, main.RECON.e.xptar, main.weight)
251 call inc(H.RECON.p.delta, main.RECON.p.delta, main.weight)
252 call inc(H.RECON.p.yptar, main.RECON.p.yptar, main.weight)
253 jones 1.1 call inc(H.RECON.p.xptar, main.RECON.p.xptar, main.weight)
254 call inc(H.RECON.Em, recon.Em, one)
255 call inc(H.RECON.Pm, recon.Pm, one)
256 call inc(H.gen.e.delta, vertex.e.delta, one)
257 call inc(H.gen.e.yptar, vertex.e.yptar, one)
258 call inc(H.gen.e.xptar, -vertex.e.xptar, one)
259 call inc(H.gen.p.delta, vertex.p.delta, one)
260 call inc(H.gen.p.yptar, vertex.p.yptar, one)
261 call inc(H.gen.p.xptar, -vertex.p.xptar, one)
262 call inc(H.gen.Em, vertex.Em, one)
263
264 ! ... update counters and integrated weights FOR EVENTS INSIDE OF DELTA
265 ! population limits (events are only generated to fill region inside
266 ! of the given delta limits) and below cuts.Em.max. Have to remove slop
267 ! that is added in init.f from the delta limits.
268 if (debug(4)) write(6,*)'sim: cut'
269
270 ncontribute = ncontribute + 1
271 if (debug(4)) write(6,*)'sim: ncontribute =',ncontribute
272 if (recon.e.delta.ge.(SPedge.e.delta.min+slop.MC.e.delta.used) .and.
273 > recon.e.delta.le.(SPedge.e.delta.max-slop.MC.e.delta.used) .and.
274 jones 1.1 > recon.p.delta.ge.(SPedge.p.delta.min+slop.MC.p.delta.used) .and.
275 > recon.p.delta.le.(SPedge.e.delta.max-slop.MC.p.delta.used) .and.
276 > recon.Em.le.cuts.Em.max) then
277 wtcontribute = wtcontribute + main.weight
278 endif
279 if (.not.rad_proton_this_ev) ncontribute_no_rad_proton =
280 > ncontribute_no_rad_proton + 1
281
282 ! ... update the "contribution" and "slop" limits
283
284 call limits_update(main,vertex,orig,recon,doing_deuterium,
285 > doing_pion,doing_kaon,doing_dvcs,contrib,slop)
286
287 ! ... write out line to diagnostic ntuple file, if requested
288
289 endif ! <success>
290
291 if (Nntu.gt.0) call results_ntu_write(main,vertex,orig,recon,success)
292
293 ! END of LOOP over events
294 ! Cute thing here, if ngen < 0, then just make |ngen| ATTEMPTS rather than
295 jones 1.1 ! insisting on ngen good events ... this is suitable for tests with new
296 ! and uncertain parameters, you don't want the program fishing around all
297 ! day for an event!
298
299 if (ngen.lt.0) then
300 nevent = nevent+1
301 else
302 if (success) nevent = nevent+1
303 endif
304
305 enddo ! <loop over ntried>
306
307
308 ! END GAME: NORMALIZE, GENERATE OUTPUT
309 !-------------------------------------
310
311 ! Our last event announcement ... and how long did it all take?
312 write (6,'(1x,''---> Last Event '',i9,'' ...'',i9,'' successes'')') ntried, nevent
313 call date (timestring2)
314 call time (timestring2(11:23))
315
316 jones 1.1 ! NORMALIZE!
317
318 ! ... put in the luminosity and efficiency factors
319
320 normfac = luminosity/ntried*nevent
321
322 ! ... multiply in the relevant phase spaces (see event.f for description
323 ! ... of generated variables. Electron angles for all cases.
324 ! ... add hadron angles and electron energy for all but H(e,e'p).
325 ! ... add hadron energy for A(e,e'p) (what about phase_space?)
326 ! ... Cross sections are all in microbarn/MeV**i/sr**j (i,j depend on reaction)
327
328 domega_e=(gen.e.yptar.max-gen.e.yptar.min)*(gen.e.xptar.max-gen.e.xptar.min)
329 domega_p=(gen.p.yptar.max-gen.p.yptar.min)*(gen.p.xptar.max-gen.p.xptar.min)
330
331 genvol = domega_e
332 ! genvol_inclusive = genvol !may want dOmega, or dE*dOmega
333
334 ! ... 2-fold to 5-fold.
335 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon
336 1 .or.doing_dvcs) then
337 jones 1.1 genvol = genvol * domega_p * (gen.e.E.max-gen.e.E.min)
338 endif
339
340 if (doing_heavy) then !6-fold
341 genvol = genvol * (gen.p.E.max-gen.p.E.min)
342 endif
343
344 normfac = normfac * genvol
345 if (doing_phsp) normfac = 1.0
346 wtcontribute = wtcontribute*normfac
347
348 ! Close diagnostic ntuple file, if used
349
350 if (Nntu.gt.0) call NtupleClose(filename)
351
352 ! Produce output files
353
354 900 if (ngen.eq.0) goto 1000
355
356 ! ... the diagnostic histograms
357
358 jones 1.1 i = index(base,' ')
359 genfile = ' '
360 write(genfile,'(a,''.gen'')') 'outfiles/'//base(1:i-1)
361 genifile = ' '
362 write(genifile,'(a,''.geni'')') 'outfiles/'//base(1:i-1)
363 histfile = ' '
364 write(histfile,'(a,''.hist'')') 'outfiles/'//base(1:i-1)
365 write(6,'(1x,''... writing '',a)') genfile
366
367 write(6,*) 'Come back soon!'
368
369 open(unit=7, file=genifile, status='unknown')
370 if (electron_arm.eq.1 .or. hadron_arm.eq.1) then
371 write(7,*) 'HMS Trials: ',hSTOP.trials
372 write(7,*) 'Slit hor/vert/corners ',hSTOP.slit_hor,hSTOP.slit_vert,hSTOP.slit_oct
373 write(7,*) 'Q1 entrance/mid/exit ',hSTOP.Q1_in,hSTOP.Q1_mid,hSTOP.Q1_out
374 write(7,*) 'Q2 entrance/mid/exit ',hSTOP.Q2_in,hSTOP.Q2_mid,hSTOP.Q2_out
375 write(7,*) 'Q3 entrance/mid/exit ',hSTOP.Q3_in,hSTOP.Q3_mid,hSTOP.Q3_out
376 write(7,*) 'Dipole entrance/exit ',hSTOP.D1_in,hSTOP.D1_out
377 write(7,*) 'Events reaching hut ',hSTOP.hut
378 write(7,*) 'DC1, DC2, Scin, Cal ',hSTOP.dc1,hSTOP.dc2,hSTOP.scin,hSTOP.cal
379 jones 1.1 write(7,*) 'Successes ',hSTOP.successes
380 write(7,*)
381 endif
382 if (electron_arm.eq.2 .or. hadron_arm.eq.2) then
383 write(7,*) 'SOS Trials: ',sSTOP.trials
384 write(7,*) 'Slit hor/vert/corners ',sSTOP.slit_vert,sSTOP.slit_hor,sSTOP.slit_oct
385 write(7,*) 'Quad entrance/mid/exit',sSTOP.quad_in,sSTOP.quad_mid,sSTOP.quad_out
386 write(7,*) 'D1 entrance/exit ',sSTOP.bm01_in,sSTOP.bm01_out
387 write(7,*) 'D2 entrance/exit ',sSTOP.bm02_in,sSTOP.bm02_out
388 write(7,*) 'Vacuum exit ',sSTOP.exit
389 write(7,*) 'Events reaching hut ',sSTOP.hut
390 write(7,*) 'DC1, DC2, Scin, Cal ',sSTOP.dc1,sSTOP.dc2,sSTOP.scin,sSTOP.cal
391 write(7,*) 'Successes ',sSTOP.successes
392 write(7,*)
393 endif
394 if (electron_arm.eq.3 .or. hadron_arm.eq.3) then
395 write(7,*) 'HRSr Trials: ',rSTOP.trials
396 write(7,*) 'Slit hor/vert ',rSTOP.slit_vert,rSTOP.slit_hor
397 write(7,*) 'Q1 entrance/mid/exit ',rSTOP.Q1_in,rSTOP.Q1_mid,rSTOP.Q1_out
398 write(7,*) 'Q2 entrance/mid/exit ',rSTOP.Q2_in,rSTOP.Q2_mid,rSTOP.Q2_out
399 write(7,*) 'Dipole entrance/exit ',rSTOP.D1_in,rSTOP.D1_out
400 jones 1.1 write(7,*) 'Q3 entrance/mid/exit ',rSTOP.Q3_in,rSTOP.Q3_mid,rSTOP.Q3_out
401 write(7,*) 'Events reaching hut ',rSTOP.hut
402 write(7,*) 'VDC1, VDC2 ',rSTOP.dc1,rSTOP.dc2
403 write(7,*) 'S1, S2, Cal ',rSTOP.s1,rSTOP.s2,rSTOP.cal
404 write(7,*)
405 endif
406 if (electron_arm.eq.4 .or. hadron_arm.eq.4) then
407 write(7,*) 'HRSl Trials: ',lSTOP.trials
408 write(7,*) 'Slit hor/vert ',lSTOP.slit_vert,lSTOP.slit_hor
409 write(7,*) 'Q1 entrance/mid/exit ',lSTOP.Q1_in,lSTOP.Q1_mid,lSTOP.Q1_out
410 write(7,*) 'Q2 entrance/mid/exit ',lSTOP.Q2_in,lSTOP.Q2_mid,lSTOP.Q2_out
411 write(7,*) 'Dipole entrance/exit ',lSTOP.D1_in,lSTOP.D1_out
412 write(7,*) 'Q3 entrance/mid/exit ',lSTOP.Q3_in,lSTOP.Q3_mid,lSTOP.Q3_out
413 write(7,*) 'Events reaching hut ',lSTOP.hut
414 write(7,*) 'VDC1, VDC2 ',lSTOP.dc1,lSTOP.dc2
415 write(7,*) 'S1, S2, Cal ',lSTOP.s1,lSTOP.s2,lSTOP.cal
416 endif
417
418 close(7)
419
420 open(unit=7, file=genfile, status='unknown')
421 jones 1.1
422 write(7,*) 'E arm Experimental Target Distributions:'
423 write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM'
424 do i=1,nHbins
425 write(7,'(3(1x,2(e11.4,1x)))')
426 > H.RECON.e.delta.min+(i-0.5)*H.RECON.e.delta.bin,
427 > H.RECON.e.delta.buf(i), H.RECON.e.yptar.min+(i-0.5)*
428 > H.RECON.e.yptar.bin, H.RECON.e.yptar.buf(i),
429 > H.RECON.e.xptar.min+(i-0.5)*H.RECON.e.xptar.bin,
430 > H.RECON.e.xptar.buf(i)
431 enddo
432
433 write(7,*) 'P arm Experimental Target Distributions:'
434 write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM'
435 do i=1,nHbins
436 write(7,'(3(1x,2(e11.4,1x)))')
437 > H.RECON.p.delta.min+(i-0.5)*H.RECON.p.delta.bin,
438 > H.RECON.p.delta.buf(i), H.RECON.p.yptar.min+(i-0.5)*
439 > H.RECON.p.yptar.bin, H.RECON.p.yptar.buf(i),
440 > H.RECON.p.xptar.min+(i-0.5)*H.RECON.p.xptar.bin,
441 > H.RECON.p.xptar.buf(i)
442 jones 1.1 enddo
443
444 write(7,*) 'Distributions of Contributing E arm Events:'
445 write(7,'(6a12)') 'delta','CONTRIB','yuptar','CONTRIB','xptar','CONTRIB'
446 do i=1,nHbins
447 write(7,'(3(1x,2(e11.4,1x)))')
448 > H.gen.e.delta.min+(i-0.5)*H.gen.e.delta.bin,
449 > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5)*
450 > H.gen.e.yptar.bin, H.gen.e.yptar.buf(i),
451 > H.gen.e.xptar.min+(i-0.5)*H.gen.e.xptar.bin, H.gen.e.xptar.buf(i)
452 enddo
453
454 write(7,*) 'Distributions of Contributing P arm Events:'
455 write(7,'(6a12)') 'delta','CONTRIB','yptar','CONTRIB','xptar','CONTRIB'
456 do i=1,nHbins
457 write(7,'(3(1x,2(e11.4,1x)))')
458 > H.gen.p.delta.min+(i-0.5)*H.gen.p.delta.bin,
459 > H.gen.p.delta.buf(i), H.gen.p.yptar.min+(i-0.5)*
460 > H.gen.p.yptar.bin, H.gen.p.yptar.buf(i),
461 > H.gen.p.xptar.min+(i-0.5)*H.gen.p.xptar.bin, H.gen.p.xptar.buf(i)
462 enddo
463 jones 1.1
464 write(7,*) 'Original E arm Events:'
465 write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN'
466 do i=1,nHbins
467 write(7,'(3(1x,2(e11.4,1x)))')
468 > H.gen.e.delta.min+(i-0.5)*H.gen.e.delta.bin,
469 > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5)*
470 > H.gen.e.yptar.bin, H.geni.e.yptar.buf(i),
471 > H.gen.e.xptar.min+(i-0.5)*H.gen.e.xptar.bin, H.geni.e.xptar.buf(i)
472 enddo
473
474 write(7,*) 'Original P arm Events:'
475 write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN'
476 do i=1,nHbins
477 write(7,'(3(1x,2(e11.4,1x)))')
478 > H.geni.p.delta.min+(i-0.5)*H.geni.p.delta.bin,
479 > H.geni.p.delta.buf(i), H.geni.p.yptar.min+(i-0.5)*
480 > H.geni.p.yptar.bin, H.geni.p.yptar.buf(i),
481 > H.geni.p.xptar.min+(i-0.5)*H.geni.p.xptar.bin, H.geni.p.xptar.buf(i)
482 enddo
483
484 jones 1.1 write(7,*) 'Original Em/Pm distributions:'
485 write(7,'(4a12)') 'Em','ORIGIN','Pm','ORIGIN'
486 do i=1,nHbins
487 write(7,'(3(1x,4(e11.4,1x)))')
488 > H.geni.Em.min+(i-0.5)*H.geni.Em.bin,
489 > H.geni.Em.buf(i), H.geni.Pm.min+(i-0.5)*
490 > H.geni.Pm.bin, H.geni.Pm.buf(i)
491 enddo
492
493 close(7)
494
495 ! Counters, etc report to the screen and to the diagnostic histogram file
496 ! call report(6,timestring1,timestring2,central,contrib,sum_sigcc)
497 open(unit=7,file=histfile,status='unknown')
498 call report(7,timestring1,timestring2,central,contrib,sum_sigcc)
499 close(unit=7)
500
501 1000 end
502
503 !-------------------------------------------------------------------
504
505 jones 1.1 subroutine inc(hist,val,weight)
506
507 implicit none
508 include 'histograms.inc'
509
510 integer*4 ibin
511 real*8 val, weight
512 record /hist_entry/ hist
513
514 ibin= nint(0.5+(val-hist.min)/hist.bin)
515 if(ibin.ge.1.and.ibin.le.nHbins)then
516 hist.buf(ibin) = hist.buf(ibin) + weight
517 endif
518
519 return
520 end
521
522 !--------------------------------------------------------------------
523
524 subroutine report(iun,timestring1,timestring2,central,contrib,sum_sigcc)
525
526 jones 1.1 implicit none
527 include 'simulate.inc'
528 include 'radc.inc'
529 include 'brem.inc'
530
531 integer*4 iun
532 real*8 sum_sigcc
533 character*23 timestring1, timestring2
534 record /contrib/ contrib
535 record /event_central/ central
536
537 ! Output of diagnostics
538
539 ! Run time
540 write(iun,'(/1x,''BEGIN Time: '',a23)') timestring1
541 write(iun,'(1x,''END Time: '',a23)') timestring2
542
543 ! Kinematics
544 write(iun,*) 'KINEMATICS:'
545 if (doing_eep) then
546 if (doing_hyd_elast) then
547 jones 1.1 write(iun,*) ' ****-------- H(e,e''p) --------****'
548 else if (doing_deuterium) then
549 write(iun,*) ' ****-------- D(e,e''p) --------****'
550 else if (doing_heavy) then
551 write(iun,*) ' ****-------- A(e,e''p) --------****'
552 else
553 stop 'I don''t have ANY idea what (e,e''p) we''re doing!!!'
554 endif
555 else if (doing_pion) then
556 if (doing_hydpi) then
557 if (targ.A .eq. 1) then
558 write(iun,*) ' ****-------- H(e,e''pi) --------****'
559 else if (targ.A .ge.3) then
560 write(iun,*) ' ****-------- A(e,e''pi) --------****'
561 endif
562 else if (doing_deutpi) then
563 write(iun,*) ' ****-------- D(e,e''pi) --------****'
564 else if (doing_hepi) then
565 write(iun,*) ' ****-------- A(e,e''pi) --------****'
566 else
567 stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!'
568 jones 1.1 endif
569 if (which_pion .eq. 0) then
570 write(iun,*) ' ****---- Default Final State ----****'
571 else if (which_pion .eq. 10) then
572 write(iun,*) ' ****---- Final State is A+pi ----****'
573 endif
574 else if (doing_dvcs) then
575 if (doing_hyddvcs) then
576 if (targ.A .eq. 1) then
577 write(iun,*) ' ****-------- H(e,e''gamma) --------****'
578 else if (targ.A .ge.3) then
579 write(iun,*) ' ****-------- A(e,e''gamma) --------****'
580 endif
581 else if (doing_deutdvcs) then
582 write(iun,*) ' ****-------- D(e,e''gamma) --------****'
583 else if (doing_hedvcs) then
584 write(iun,*) ' ****-------- A(e,e''gamma) --------****'
585 else
586 stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!'
587 endif
588 if (which_dvcs .eq. 0) then
589 jones 1.1 write(iun,*) ' ****---- Default Final State ----****'
590 else if (which_dvcs .eq. 10) then
591 write(iun,*) ' ****---- Final State is A+gamma ----****'
592 endif
593 else if (doing_kaon) then
594 if (doing_hydkaon) then
595 if (targ.A .eq. 1) then
596 write(iun,*) ' ****-------- H(e,e''K) --------****'
597 else if (targ.A .ge.3) then
598 write(iun,*) ' ****-------- A(e,e''K) --------****'
599 endif
600 else if (doing_deutkaon) then
601 write(iun,*) ' ****-------- D(e,e''K) --------****'
602 else if (doing_hekaon) then
603 write(iun,*) ' ****-------- A(e,e''K) --------****'
604 else
605 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
606 endif
607 if (which_kaon.eq.0) then
608 write(iun,*) ' ****---- producing a LAMBDA ----****'
609 else if (which_kaon.eq.1) then
610 jones 1.1 write(iun,*) ' ****---- producing a SIGMA0 ----****'
611 else if (which_kaon.eq.2) then
612 write(iun,*) ' ****---- producing a SIGMA- ----****'
613 else if (which_kaon.eq.10) then
614 write(iun,*) ' ****---- WITH BOUND LAMBDA ----****'
615 else if (which_kaon.eq.11) then
616 write(iun,*) ' ****---- WITH BOUND SIGMA0 ----****'
617 else if (which_kaon.eq.12) then
618 write(iun,*) ' ****---- WITH BOUND SIGMA- ----****'
619 else
620 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
621 endif
622 else if (doing_phsp) then
623 write(iun,*) ' ****--- PHASE SPACE - NO physics, NO radiation (may not work)---****'
624 else
625 stop 'I don''t have ANY idea what we''re doing!!!'
626 endif
627
628 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'Ebeam', Ebeam, 'MeV'
629 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)')
630 > '(dE/E)beam', dEbeam/Ebeam, '(full wid)'
631 jones 1.1 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'x-width', gen.xwid, 'cm'
632 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'y-width', gen.ywid, 'cm'
633 write(iun,'(9x,a12,'' = '',i15,2x,a16)')
634 > 'fr_pattern',targ.fr_pattern, '1=square,2=circ'
635 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr1', targ.fr1, 'cm'
636 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr2', targ.fr2, 'cm'
637
638 write(iun,*) ' '
639 write(iun,'(9x,18x,2(a15,2x))') '____E arm____','____P arm____'
640 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'angle',
641 > spec.e.theta*degrad, spec.p.theta*degrad, 'deg'
642 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'momentum',
643 > spec.e.p, spec.p.p, 'MeV/c'
644
645 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'x offset',
646 > spec.e.offset.x, spec.p.offset.x, 'cm'
647 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'y offset',
648 > spec.e.offset.y, spec.p.offset.y, 'cm'
649 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'z offset',
650 > spec.e.offset.z, spec.p.offset.z, 'cm'
651 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'xptar offset',
652 jones 1.1 > spec.e.offset.xptar, spec.p.offset.xptar, 'mr'
653 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'yptar offset',
654 > spec.e.offset.yptar, spec.p.offset.yptar, 'mr'
655
656
657 write(iun,'(6x,''Central Values: '',a5,'' = '',f15.4,2x,a9)')
658 > 'Q2', central.Q2/1.d6, '(GeV/c)^2'
659 write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'q', central.q, 'MeV/c'
660 write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'omega',
661 > central.omega, 'MeV'
662 write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'Em', central.Em, 'MeV'
663 write(iun,'(22x,a5,'' = '',f15.4,2x,a9)') 'Pm', central.Pm, 'MeV/c'
664
665 ! Target
666 write(iun,*) 'TARGET specs:'
667 9911 format(2x,2(5x,a10,' = ',e12.6,1x,a5))
668 write(iun,9911) 'A', targ.A, ' ', 'Z', targ.Z, ' '
669 write(iun,9911) 'mass', targ.mass_amu, 'amu', 'mass', targ.M,'MeV'
670 write(iun,9911) 'Mrec', targ.mrec_amu, 'amu', 'Mrec', targ.Mrec,'MeV'
671 write(iun,9911) 'Mtar_struck',targ.Mtar_struck,'MeV',
672 > 'Mrec_struck',targ.Mrec_struck,'MeV'
673 jones 1.1 write(iun,9911) 'rho', targ.rho, 'g/cm3', 'thick', targ.thick,'g/cm2'
674 write(iun,9911) 'angle', targ.angle*degrad, 'deg', 'abundancy',
675 > targ.abundancy, '%'
676 write(iun,9911) 'X0', targ.X0, 'g/cm2', 'X0_cm', targ.X0_cm,'cm'
677 write(iun,9911) 'length',targ.length,'cm','zoffset',targ.zoffset,'cm'
678 write(iun,9911) 'xoffset',targ.xoffset,'cm','yoffset',targ.yoffset,'cm'
679 write(iun,'(t12,3a15)') '__ave__', '__lo__', '__hi__'
680 9912 format(1x,a15,3f15.5,2x,a6)
681 write(iun,9912) 'Coulomb', targ.Coulomb.ave, targ.Coulomb.min,
682 > targ.Coulomb.max, 'MeV'
683 write(iun,9912) 'Eloss_beam', targ.Eloss(1).ave,
684 > targ.Eloss(1).min, targ.Eloss(1).max, 'MeV'
685 write(iun,9912) 'Eloss_e', targ.Eloss(2).ave, targ.Eloss(2).min,
686 > targ.Eloss(2).max, 'MeV'
687 write(iun,9912) 'Eloss_p', targ.Eloss(3).ave, targ.Eloss(3).min,
688 > targ.Eloss(3).max, 'MeV'
689 write(iun,9912) 'teff_beam', targ.teff(1).ave, targ.teff(1).min,
690 > targ.teff(1).max, 'radlen'
691 write(iun,9912) 'teff_e', targ.teff(2).ave, targ.teff(2).min,
692 > targ.teff(2).max, 'radlen'
693 write(iun,9912) 'teff_p', targ.teff(3).ave, targ.teff(3).min,
694 jones 1.1 > targ.teff(3).max, 'radlen'
695 9913 format(1x,a15,t25,f15.5,2x,a6)
696 write(iun,9913) 'musc_nsig_max', targ.musc_nsig_max, ' '
697 write(iun,9913) 'musc_max_beam', targ.musc_max(1)*1000., 'mr'
698 write(iun,9913) 'musc_max_e', targ.musc_max(2)*1000., 'mr'
699 write(iun,9913) 'musc_max_p', targ.musc_max(3)*1000., 'mr'
700
701 ! Flags
702 write(iun,*) 'FLAGS:'
703 write(iun,'(5x,2(2x,a19,''='',l2))') 'doing_dvcs', doing_dvcs,
704 > 'which_dvcs', which_dvcs
705 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_eep', doing_eep,
706 > 'doing_kaon', doing_kaon, 'doing_pion', doing_pion
707 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_phsp', doing_phsp,
708 > 'which_pion', which_pion, 'which_kaon', which_kaon
709 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hyd_elast', doing_hyd_elast,
710 > 'doing_deuterium', doing_deuterium, 'doing_heavy', doing_heavy
711 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hyddvcs', doing_hyddvcs,
712 > 'doing_deutdvcs', doing_deutdvcs, 'doing_hedvcs', doing_hedvcs
713 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydpi', doing_hydpi,
714 > 'doing_deutpi', doing_deutpi, 'doing_hepi', doing_hepi
715 jones 1.1 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydkaon', doing_hydkaon,
716 > 'doing_deutkaon', doing_deutkaon, 'doing_hekaon', doing_hekaon
717 write(iun,'(5x,3(2x,a19,''='',l2))') 'mc_smear', mc_smear,
718 > 'electron_arm', electron_arm, 'hadron_arm', hadron_arm
719 write(iun,'(5x,3(2x,a19,''='',l2)))') 'using_Eloss', using_Eloss,
720 > 'using_Coulomb',using_Coulomb,'deForest_flag',deForest_flag
721 write(iun,'(5x,3(2x,a19,''='',l2)))') 'correct_Eloss', correct_Eloss,
722 > 'correct_raster',correct_raster, 'doing_decay', doing_decay
723 write(iun,'(5x,2(2x,a19,''='',l2))')
724 > 'using_E_arm_montecarlo', using_E_arm_montecarlo,
725 > 'using_P_arm_montecarlo', using_P_arm_montecarlo
726 write(iun,'(7x,a11,''='',f10.3,a4)') 'ctau',ctau,'cm'
727
728 ! Counters
729 write(iun,*) 'COUNTERS:'
730 write(iun,'(12x,''Ngen (request) = '',i10)') ngen
731 write(iun,'(12x,''Ntried = '',i10)') ntried
732 write(iun,'(12x,''Ncontribute = '',i10)') ncontribute
733 write(iun,'(12x,''Nco_no_rad_prot= '',i10)') ncontribute_no_rad_proton
734 write(iun,'(12x,''-> %no_rad_prot= '',f10.3)')
735 > (100.*ncontribute_no_rad_proton/max(dfloat(ncontribute),0.1))
736 jones 1.1 write(iun,'(/1x,''INTEGRATED WEIGHTS (number of counts in delta/Em cuts!):'')')
737 write(iun,'(1x,'' MeV: wtcontr= '',e16.8)') wtcontribute/nevent
738
739 ! Radiative corrections
740 write(iun,*) 'RADIATIVE CORRECTIONS:'
741 if (.not.using_rad) then
742 write(iun,'(x,a14,''='',l3)') 'using_rad', using_rad
743 else
744 write(iun,'(4(x,a14,''='',l3))') 'use_expon',use_expon,
745 > 'include_hard',include_hard,'calc_spence',calculate_spence
746 write(iun,'(2(x,a14,''='',l3))') 'using_rad', using_rad,
747 > 'use_offshell_rad', use_offshell_rad
748 write(iun,'(4(x,a14,''='',i3))') 'rad_flag',rad_flag,
749 > 'extrad_flag', extrad_flag, 'one_tail', one_tail,
750 > 'intcor_mode', intcor_mode
751 write(iun,'(x,a14,''='',f11.3)') 'dE_edge_test',dE_edge_test
752 write(iun,'(x,a14,''='',f11.3)') 'Egamma_max', Egamma_tot_max
753
754 9914 format(1x,a18,' = ',4f11.3)
755 write(iun,*) 'Central Values:'
756 write(iun,9914) 'hardcorfac', central.rad.hardcorfac
757 jones 1.1 write(iun,9914) 'etatzai', central.rad.etatzai
758 write(iun,9914) 'frac(1:3)', central.rad.frac
759 write(iun,9914) 'lambda(1:3)', central.rad.lambda
760 write(iun,9914) 'bt(1:2)', central.rad.bt
761 write(iun,9914) 'c_int(0:3)', central.rad.c_int
762 write(iun,9914) 'c_ext(0:3)', central.rad.c_ext
763 write(iun,9914) 'c(0:3)', central.rad.c
764 write(iun,9914) 'g_int', central.rad.g_int
765 write(iun,9914) 'g_ext', central.rad.g_ext
766 write(iun,9914) 'g(0:3)', central.rad.g
767 endif
768
769 ! Miscellaneous
770 write(iun,*) 'MISCELLANEOUS:'
771 9915 format(12x,a14,' = ',e16.6,1x,a6)
772 write(iun,*) 'Note that central.sigcc is for central delta,theta,phi in both spectrometers'
773 write(iun,*) ' and may give non-physical kinematics, esp. for Hydrogen'
774 write(iun,*) 'Note also that AVE.sigcc is really AVER.weight (the two arenot exactly equal)'
775 write(iun,9915) 'CENTRAL.sigcc', central.sigcc, ' '
776 write(iun,9915) 'AVERAGE.sigcc', sum_sigcc/nevent, ' '
777 write(iun,9915) 'charge', EXPER.charge, 'mC'
778 jones 1.1 write(iun,9915) 'targetfac', targetfac, ' '
779 write(iun,9915) 'luminosity', luminosity, 'ub^-1'
780 write(iun,9915) 'luminosity', luminosity*(hbarc/100000.)**2, 'GeV^2'
781 write(iun,9915) 'genvol', genvol, ' '
782 write(iun,9915) 'normfac', normfac, ' '
783 9916 format(12x,a14,' = ',f6.1,' to ',f6.1,' in ',i4,' bins of ',f7.2,1x,a5)
784 if (doing_heavy) then
785 write(iun,'(12x,''Theory file: '',a)')
786 > theory_file(1:index(theory_file,' ')-1)
787 endif
788
789 ! Input spectrometer limits.
790
791 write(iun,*) 'Input Spectrometer Limits:'
792 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.delta.min/max',
793 > SPedge.e.delta.min,SPedge.e.delta.max,'%'
794 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.yptar.min/max',
795 > SPedge.e.yptar.min,SPedge.e.yptar.max,'rad'
796 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.xptar.min/max',
797 > SPedge.e.xptar.min,SPedge.e.xptar.max,'rad'
798 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.delta.min/max',
799 jones 1.1 > SPedge.p.delta.min,SPedge.p.delta.max,'%'
800 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.yptar.min/max',
801 > SPedge.p.yptar.min,SPedge.p.yptar.max,'rad'
802 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.xptar.min/max',
803 > SPedge.p.xptar.min,SPedge.p.xptar.max,'rad'
804
805
806 ! Edges used in generation and checking, as well as range of events found
807 ! to contribute within those edges
808
809 ! ... on GENERATED qties
810 write(iun,*) 'Limiting VERTEX values (vertex.e/p.*,Em,Pm,Trec)'
811 write(iun,*) ' USED limits are gen.e/p.*, and VERTEXedge.Em,Pm,Trec'
812 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') '______used______',
813 > '_____found______'
814 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi'
815 9917 format(1x,a18,t21,2f12.3,t50,2f10.3,2x,a5)
816 write(iun,9917) 'E arm delta', gen.e.delta.min, gen.e.delta.max,
817 > contrib.gen.e.delta.lo, contrib.gen.e.delta.hi, '%'
818 write(iun,9917) 'E arm yptar', gen.e.yptar.min*1000.,
819 > gen.e.yptar.max*1000.,
820 jones 1.1 > contrib.gen.e.yptar.lo*1000., contrib.gen.e.yptar.hi*1000.,'mr'
821 write(iun,9917) 'E arm xptar', gen.e.xptar.min*1000.,
822 > gen.e.xptar.max*1000.,
823 > contrib.gen.e.xptar.lo*1000., contrib.gen.e.xptar.hi*1000.,'mr'
824 write(iun,9917) 'P arm delta', gen.p.delta.min, gen.p.delta.max,
825 > contrib.gen.p.delta.lo, contrib.gen.p.delta.hi, '%'
826 write(iun,9917) 'P arm yptar', gen.p.yptar.min*1000.,
827 > gen.p.yptar.max*1000.,
828 > contrib.gen.p.yptar.lo*1000., contrib.gen.p.yptar.hi*1000.,'mr'
829 write(iun,9917) 'P arm xptar', gen.p.xptar.min*1000.,
830 > gen.p.xptar.max*1000.,
831 > contrib.gen.p.xptar.lo*1000., contrib.gen.p.xptar.hi*1000.,'mr'
832 write(iun,9917) 'sumEgen', gen.sumEgen.min, gen.sumEgen.max,
833 > contrib.gen.sumEgen.lo, contrib.gen.sumEgen.hi, 'MeV'
834
835 ! ... on VERTEX qties
836 ! write(iun,*) 'Limiting VERTEX values:'
837 ! write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)')
838 ! > '______used______', '_____found______'
839 ! write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi'
840 write(iun,9917) 'Trec', VERTEXedge.Trec.min, VERTEXedge.Trec.max,
841 jones 1.1 > contrib.vertex.Trec.lo, contrib.vertex.Trec.hi, 'MeV'
842 write(iun,9917) 'Em', VERTEXedge.Em.min, VERTEXedge.Em.max,
843 > contrib.vertex.Em.lo, contrib.vertex.Em.hi, 'MeV'
844 write(iun,9917) 'Pm', VERTEXedge.Pm.min, VERTEXedge.Pm.max,
845 > contrib.vertex.Pm.lo, contrib.vertex.Pm.hi, 'MeV/c'
846 if ((doing_deuterium .or. doing_pion .or. doing_kaon .or. doing_dvcs) .and. using_rad)
847 > write(iun,*) ' *** NOTE: sumEgen.min only used in GENERATE_RAD'
848
849 ! ... on TRUE qties
850 write(iun,*) 'Limiting ORIGINAL values: orig.e/p.*,Em,Pm,Trec (no edge.* limits for Pm,Trec)'
851 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)')
852 > '______used______', '_____found______'
853 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi'
854 write(iun,9917) 'E arm E', edge.e.E.min, edge.e.E.max,
855 > contrib.tru.e.E.lo, contrib.tru.e.E.hi, 'MeV'
856 write(iun,9917) 'E arm yptar', edge.e.yptar.min*1000.,
857 > edge.e.yptar.max*1000.,
858 > contrib.tru.e.yptar.lo*1000., contrib.tru.e.yptar.hi*1000.,'mr'
859 write(iun,9917) 'E arm xptar', edge.e.xptar.min*1000., edge.e.xptar.max*1000.,
860 > contrib.tru.e.xptar.lo*1000., contrib.tru.e.xptar.hi*1000., 'mr'
861 write(iun,9917) 'P arm E', edge.p.E.min, edge.p.E.max,
862 jones 1.1 > contrib.tru.p.E.lo, contrib.tru.p.E.hi, 'MeV'
863 write(iun,9917) 'P arm yptar', edge.p.yptar.min*1000.,
864 > edge.p.yptar.max*1000.,
865 > contrib.tru.p.yptar.lo*1000., contrib.tru.p.yptar.hi*1000.,'mr'
866 write(iun,9917) 'P arm xptar', edge.p.xptar.min*1000., edge.p.xptar.max*1000.,
867 > contrib.tru.p.xptar.lo*1000., contrib.tru.p.xptar.hi*1000., 'mr'
868 write(iun,9917) 'Em', max(-999999.999d0,edge.Em.min), min(999999.999d0,edge.Em.max),
869 > contrib.tru.Em.lo, contrib.tru.Em.hi, 'MeV'
870 write(iun,9917) 'Pm', 0., 0., contrib.tru.Pm.lo,
871 > contrib.tru.Pm.hi, 'MeV'
872 write(iun,9917) 'Trec', 0., 0., contrib.tru.Trec.lo,
873 > contrib.tru.Trec.hi, 'MeV'
874
875 ! ... on SPECTROMETER qties
876 write(iun,*) 'Limiting SPECTROMETER values'
877 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)')
878 > '______used______', '_____found______'
879 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi'
880 write(iun,9917) 'E arm delta', SPedge.e.delta.min,
881 > SPedge.e.delta.max, contrib.SP.e.delta.lo,
882 > contrib.SP.e.delta.hi, '%'
883 jones 1.1 write(iun,9917) 'E arm yptar', SPedge.e.yptar.min*1000.,
884 > SPedge.e.yptar.max*1000.,
885 > contrib.SP.e.yptar.lo*1000., contrib.SP.e.yptar.hi*1000., 'mr'
886 write(iun,9917) 'E arm xptar', SPedge.e.xptar.min*1000.,
887 > SPedge.e.xptar.max*1000.,
888 > contrib.SP.e.xptar.lo*1000., contrib.SP.e.xptar.hi*1000., 'mr'
889 write(iun,9917) 'P arm delta', SPedge.p.delta.min,
890 > SPedge.p.delta.max, contrib.SP.p.delta.lo,
891 > contrib.SP.p.delta.hi, '%'
892 write(iun,9917) 'P arm yptar', SPedge.p.yptar.min*1000.,
893 > SPedge.p.yptar.max*1000.,
894 > contrib.SP.p.yptar.lo*1000., contrib.SP.p.yptar.hi*1000., 'mr'
895 write(iun,9917) 'P arm xptar', SPedge.p.xptar.min*1000.,
896 > SPedge.p.xptar.max*1000.,
897 > contrib.SP.p.xptar.lo*1000., contrib.SP.p.xptar.hi*1000., 'mr'
898
899 ! ... on RADIATION qties
900 if (using_rad) then
901 write(iun,*) 'Limiting RADIATION values CONTRIBUTING to the (Em,Pm) distributions:'
902 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)')
903 > '______used______', '_____found______'
904 jones 1.1 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi'
905 write(iun,9917) 'Egamma(1)', 0., Egamma1_max,
906 > contrib.rad.Egamma(1).lo, contrib.rad.Egamma(1).hi, 'MeV'
907 write(iun,9917) 'Egamma(2)', 0., Egamma2_max, contrib.rad.Egamma(2).lo,
908 > contrib.rad.Egamma(2).hi, 'MeV'
909 write(iun,9917) 'Egamma(3)', 0., Egamma3_max, contrib.rad.Egamma(3).lo,
910 > contrib.rad.Egamma(3).hi, 'MeV'
911 write(iun,9917) 'Egamma_total', 0., Egamma_tot_max,
912 > contrib.rad.Egamma_total.lo, contrib.rad.Egamma_total.hi,'MeV'
913 endif
914
915 ! ... on slops
916 write(iun,*) 'ACTUAL and LIMITING SLOP values used/obtained:'
917 write(iun,'(t25,3a10)') '__used__', '__min__', '__max__'
918 9918 format(1x,a10,a12,t25,3f10.3,2x,a5)
919 write(iun,9918) 'slop.MC ', 'E arm delta', slop.MC.e.delta.used,
920 > slop.MC.e.delta.lo, slop.MC.e.delta.hi, '%'
921 write(iun,9918) ' ', 'E arm yptar', slop.MC.e.yptar.used*1000.,
922 > slop.MC.e.yptar.lo*1000., slop.MC.e.yptar.hi*1000., 'mr'
923 write(iun,9918) ' ', 'E arm xptar', slop.MC.e.xptar.used*1000.,
924 > slop.MC.e.xptar.lo*1000., slop.MC.e.xptar.hi*1000., 'mr'
925 jones 1.1 write(iun,9918) ' ', 'P arm delta', slop.MC.p.delta.used,
926 > slop.MC.p.delta.lo, slop.MC.p.delta.hi, '%'
927 write(iun,9918) ' ', 'P arm yptar', slop.MC.p.yptar.used*1000.,
928 > slop.MC.p.yptar.lo*1000., slop.MC.p.yptar.hi*1000., 'mr'
929 write(iun,9918) ' ', 'P arm xptar', slop.MC.p.xptar.used*1000.,
930 > slop.MC.p.xptar.lo*1000., slop.MC.p.xptar.hi*1000., 'mr'
931 write(iun,9918) 'slop.total', 'Em', slop.total.Em.used,
932 > slop.total.Em.lo, slop.total.Em.hi, 'MeV'
933 write(iun,9918) ' ', 'Pm', 0., slop.total.Pm.lo,
934 > slop.total.Pm.hi, 'MeV/c'
935
936 write(iun,'(/)')
937 return
938 end
939
940 !-----------------------------------------------------------------------
941
942 subroutine calculate_central(central)
943
944 implicit none
945 include 'simulate.inc'
946 jones 1.1 include 'radc.inc'
947
948 integer i
949 logical success
950 record /event_main/ main
951 record /event/ vertex0, recon0
952 record /event_central/ central
953
954 ! Pop in generation values for central event
955
956 if (debug(2)) write(6,*)'calc_cent: entering...'
957 main.target.x = 0.0
958 main.target.y = 0.0
959 main.target.z = 0.0
960 vertex0.Ein = Ebeam_vertex_ave
961 main.target.Coulomb = targ.Coulomb.ave
962 main.target.Eloss(1) = targ.Eloss(1).ave
963 main.target.teff(1) = targ.teff(1).ave
964 vertex0.e.delta = 0.0
965 vertex0.e.yptar = 0.0
966 vertex0.e.xptar = 0.0
967 jones 1.1 vertex0.e.theta = spec.e.theta
968 vertex0.e.phi = spec.e.phi
969 vertex0.e.P = spec.e.P*(1.+vertex0.e.delta/100.)
970 vertex0.e.E = vertex0.e.P
971 vertex0.p.delta = 0.0
972 vertex0.p.yptar = 0.0
973 vertex0.p.xptar = 0.0
974 vertex0.p.theta = spec.p.theta
975 vertex0.p.phi = spec.p.phi
976 vertex0.p.P = spec.p.P*(1.+vertex0.p.delta/100.)
977 vertex0.p.E = sqrt(vertex0.p.P**2 + Mh2)
978
979 ! Complete_recon_ev for vertex0. Note: complete_recon_ev doesn't
980 ! call radc_init_ev or calculate teff(2,3).
981
982 ! JRA: Do we want complete_ev or complete_recon_ev? Do we want to calculate
983 ! and/or dump other central values (for pion/kaon production, for example).
984
985 call complete_recon_ev(vertex0,success)
986 if (debug(2)) write(6,*)'calc_cent: done with comp_ev'
987 if (.not.success) stop 'COMPLETE_EV failed trying to complete a CENTRAL event!'
988 jones 1.1 central.Q2 = vertex0.Q2
989 central.q = vertex0.q
990 central.omega = vertex0.omega
991 central.Em = vertex0.Em
992 central.Pm = vertex0.Pm
993 main.target.teff(2) = targ.teff(2).ave
994 main.target.teff(3) = targ.teff(3).ave
995 if (debug(2)) write(6,*)'calc_cent: calling radc_init_ev'
996 if (debug(4)) write(6,*)'calc_cent: Ein =',vertex0.Ein
997 call radc_init_ev(main,vertex0)
998 if (debug(2)) write(6,*)'calc_cent: done with radc_init_ev'
999 if (using_rad) then
1000 central.rad.hardcorfac = hardcorfac
1001 central.rad.etatzai= etatzai
1002 central.rad.g_int = g_int
1003 central.rad.g_ext = g_ext
1004 do i = 0, 3
1005 if (i.gt.0) central.rad.frac(i) = frac(i)
1006 if (i.gt.0) central.rad.lambda(i) = lambda(i)
1007 if (i.gt.0.and.i.lt.3) central.rad.bt(i) = bt(i)
1008 central.rad.c_int(i) = c_int(i)
1009 jones 1.1 central.rad.c_ext(i) = c_ext(i)
1010 central.rad.c(i) = c(i)
1011 central.rad.g(i) = g(i)
1012 enddo
1013 endif
1014 if (debug(4)) write(6,*)'calc_cent: at 1'
1015
1016 do i = 1, neventfields
1017 recon0.all(i) = vertex0.all(i)
1018 enddo
1019 call complete_main(.true.,main,vertex0,recon0,success)
1020 central.sigcc = main.sigcc
1021
1022 if (debug(2)) write(6,*)'calc_cent: ending...'
1023 return
1024 end
1025
1026 !-------------------------------------------------------------------------
1027
1028 subroutine montecarlo(orig,main,recon,success)
1029
1030 jones 1.1 implicit none
1031 include 'simulate.inc'
1032
1033 * Track-coordinate and spectrometer common blocks
1034
1035 real*8 x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm,delta_E_arm
1036 real*8 x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm,delta_P_arm
1037 real*8 xfp, yfp, dxfp, dyfp
1038 real*8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2)
1039 logical success
1040 logical ok_E_arm, ok_P_arm
1041 record /event/ orig, recon
1042 record /event_main/ main
1043
1044 real*8 beta_electron
1045 parameter (beta_electron = 1.)
1046 real*8 tmpfact
1047 real*8 fry !fast raster y position.
1048 real*8 frx !fast raster x position - left as looking downstream
1049 real*8 m2 !mass for call to mc_hms(sos). Changes if decay
1050 real*8 pathlen
1051 jones 1.1 real*8 ztemp
1052 real*8 zero
1053 parameter (zero=0.0d0) !double precision zero for subroutine calls
1054
1055 ! Prepare the event for the Monte Carlo's and/or spectrometer cuts
1056
1057 success = .false.
1058 ntup.resfac = 0.0 !resfac (see simulate.inc)
1059 c if (correct_raster) then
1060 fry = -main.target.rastery
1061 frx = main.target.rasterx
1062 c else
1063 c fry = 0.0
1064 c frx = 0.0
1065 c endif
1066
1067 !BEAM MULTIPLE SCATTERING: NOT YET IMPLEMENTED, JUST GETTING IT READY!
1068
1069 ! ... multiple scattering of the beam. Generate angles for beam deflection,
1070 ! ... and then apply those angles to the electron and proton.
1071 ! ... The yptar offset goes directly to yptar of both arms (small angle approx).
1072 jones 1.1 ! ... xptar is multiplied by cos(theta) to get the xptar offsets.
1073
1074 if (mc_smear) then
1075 call target_musc(orig.Ein,beta_electron,main.target.teff(1),dang_in)
1076 else
1077 dang_in(1)=0.0 !yp offset, goes directly to yp of both arms
1078 dang_in(2)=0.0 !xp offset, mult. by cos_th for xp of both arms
1079 endif
1080
1081 if (debug(3)) write(6,*) 'mc: using p,e arm mc =',
1082 > using_P_arm_montecarlo,using_E_arm_montecarlo
1083
1084 !____ P arm ________________
1085
1086 ! Go from TRUE to SPECTROMETER quantities by computing target distortions
1087 ! ... ionization loss correction (if requested)
1088
1089 if (using_Eloss) then
1090 if (debug(3)) write(6,*)'mc: p arm stuff0 =',
1091 > orig.p.E,main.target.Eloss(3),spec.p.P
1092 main.SP.p.delta = (sqrt(abs((orig.p.E-main.target.Eloss(3))**2
1093 jones 1.1 > -Mh2))-spec.p.P) / spec.p.P*100.
1094 else
1095 if (debug(3)) write(6,*)'mc: p arm stuff1 =',orig.p.delta
1096 main.SP.p.delta = orig.p.delta
1097 endif
1098
1099 ! ... multiple scattering
1100
1101 if (mc_smear) then
1102 beta = orig.p.p/orig.p.E
1103 call target_musc(orig.p.p, beta, main.target.teff(3), dangles)
1104 else
1105 dangles(1)=0.0
1106 dangles(2)=0.0
1107 endif
1108 main.SP.p.yptar = orig.p.yptar + dangles(1) + dang_in(1)
1109 main.SP.p.xptar = orig.p.xptar + dangles(2) + dang_in(2)*spec.p.cos_th
1110
1111 ! CASE 1: Using the spectrometer Monte Carlo
1112
1113 if (using_P_arm_montecarlo) then
1114 jones 1.1
1115 ! ... change to P arm spectrometer coordinates (TRANSPORT system),
1116
1117 if (abs(cos(spec.p.phi)).gt.0.0001) then !phi not at +/- pi/2
1118 write(6,*) 'y_P_arm, z_P_arm will be incorrect if spec.p.phi <> pi/2 or 3*pi/2'
1119 write(6,*) 'spec.p.phi=',spec.p.phi,'=',spec.p.phi*180/pi,'degrees'
1120 endif
1121 delta_P_arm = main.SP.p.delta
1122 x_P_arm = -main.target.y
1123 y_P_arm = -main.target.x*spec.p.cos_th - main.target.z*spec.p.sin_th*sin(spec.p.phi)
1124 z_P_arm = main.target.z*spec.p.cos_th + main.target.x*spec.p.sin_th*sin(spec.p.phi)
1125
1126 ! ... Apply spectrometer offset (using spectrometer coordinate system).
1127 x_P_arm = x_P_arm - spec.p.offset.x
1128 y_P_arm = y_P_arm - spec.p.offset.y
1129 z_P_arm = z_P_arm - spec.p.offset.z
1130 dx_P_arm = main.SP.p.xptar - spec.p.offset.xptar
1131 dy_P_arm = main.SP.p.yptar - spec.p.offset.yptar
1132
1133 c write(6,*) 'pion BEFORE target field:'
1134 c write(6,*) 'xp, yp:', dx_P_arm, dy_P_arm
1135 jones 1.1 ! GAW -project to z=0 to compare with reconstructed target positions
1136 c in_P_arm.z = y_P_arm - z_P_arm*dy_P_arm
1137
1138 c write(6,*) 'using_tgt_field',using_tgt_field
1139 c
1140 c
1141 c write(*,*) 'sign_hms_part =' ,sign_hms_part
1142 if (using_tgt_field) then ! do target field tracking - GAW
1143
1144 c if (debug(6)) then
1145 c write(*,*) '------------------------------'
1146 c write(*,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_P_arm
1147 c endif
1148 call track_from_tgt(x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm,
1149 > sign_hms_part*spec.p.P*(1+delta_P_arm/100.),Mh,1,ok_P_arm)
1150 endif
1151 ! GAW - end 99/11/3
1152 ! ........ drift this position back to z=0, the plane through the target center
1153
1154 x_P_arm = x_P_arm - z_P_arm*dx_P_arm
1155 y_P_arm = y_P_arm - z_P_arm*dy_P_arm
1156 jones 1.1 z_P_arm = 0.0
1157
1158 c call print_coord2('virtual track z=0: ',x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm) ! GAW
1159
1160
1161 ! ... Monte Carlo through the P arm! if we succeed, continue ...
1162 ! ... Here's what's passed:
1163 ! spectrometer central momentum
1164 ! spectrometer central theta angle
1165 ! momentum delta
1166 ! x, y, z, theta, and phi in TRANSPORT coordinates
1167 ! x, y, theta, and phi at the focal plane
1168 ! radiation length of target (NOT USED!)
1169 ! particle mass (modified if the hadron decays)
1170 ! multiple scattering flag
1171 ! wcs smearing flag
1172 ! decay flag
1173 ! resmult=resolution factor for DCs (see simulate.inc)
1174 ! vertical position at target (for reconstruction w/raster MEs).
1175 ! ok flag
1176
1177 jones 1.1 m2 = Mh2
1178 pathlen = 0.0
1179 if (hadron_arm.eq.1) then
1180 call mc_hms(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm,
1181 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
1182 > m2, mc_smear, mc_smear, doing_decay,
1183 > ntup.resfac, frx, fry, ok_P_arm, pathlen, using_tgt_field
1184 > ,sign_hms_part)
1185 else if (hadron_arm.eq.2) then
1186 call mc_sos(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm,
1187 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
1188 > m2, mc_smear, mc_smear, doing_decay,
1189 > ntup.resfac, fry, ok_P_arm, pathlen)
1190 else if (hadron_arm.eq.3) then
1191 call mc_hrsr(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm,
1192 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
1193 > m2, mc_smear, mc_smear, doing_decay,
1194 > ntup.resfac, fry, ok_P_arm, pathlen)
1195 else if (hadron_arm.eq.4) then
1196 call mc_hrsl(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm,
1197 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
1198 jones 1.1 > m2, mc_smear, mc_smear, doing_decay,
1199 > ntup.resfac, fry, ok_P_arm, pathlen)
1200 else if (hadron_arm.eq.5) then
1201 call mc_calo(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm,
1202 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
1203 > m2, mc_smear, mc_smear, doing_decay,
1204 > ntup.resfac, frx, fry, ok_P_arm, pathlen, using_tgt_field,
1205 > main.target.z*cos(spec.p.theta),ztemp,drift_to_cal
1206 > )
1207 endif
1208
1209 if (.not.ok_P_arm) then
1210 if (debug(3)) write(6,*)'mc: ok_P_arm =',ok_P_arm
1211 return
1212 endif
1213
1214 c call print_coord2('recon tgt: ',x_P_arm,y_P_arm,0.,dx_P_arm,dy_P_arm) ! GAW
1215
1216 main.RECON.p.delta = delta_P_arm
1217 main.RECON.p.yptar = dy_P_arm
1218 main.RECON.p.xptar = dx_P_arm
1219 jones 1.1 main.RECON.p.z = y_P_arm
1220 main.FP.p.x = xfp
1221 main.FP.p.dx = dxfp
1222 main.FP.p.y = yfp
1223 main.FP.p.dy = dyfp
1224 main.FP.p.path = pathlen
1225
1226 ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the
1227 ! OUT record
1228
1229 else
1230 main.RECON.p.delta = main.SP.p.delta
1231 main.RECON.p.yptar = main.SP.p.yptar
1232 main.RECON.p.xptar = main.SP.p.xptar
1233 endif
1234
1235 ! Fill recon quantities.
1236
1237 recon.p.delta = main.RECON.p.delta
1238 recon.p.yptar = main.RECON.p.yptar
1239 recon.p.xptar = main.RECON.p.xptar
1240 jones 1.1 recon.p.z = main.RECON.p.z
1241 recon.p.P = spec.p.P*(1.+recon.p.delta/100.)
1242 recon.p.E = sqrt(recon.p.P**2 + Mh2)
1243 call physics_angles(spec.p.theta,spec.p.phi,recon.p.xptar,
1244 > recon.p.yptar,recon.p.theta,recon.p.phi)
1245
1246 ! ... correct for energy loss - use most probable (last flag = 4)
1247
1248 if (correct_Eloss) then
1249 call trip_thru_target (3, zero, recon.p.E,
1250 > recon.p.theta, eloss_P_arm, r,Mh,4)
1251 recon.p.E = recon.p.E + eloss_P_arm
1252 recon.p.E = max(recon.p.E,sqrt(Mh2+0.000001)) !can get P~0 when calculating hadron momentum-->P<0 after eloss
1253 recon.p.P = sqrt(recon.p.E**2-Mh2)
1254 recon.p.delta = (recon.p.P-spec.p.P)/spec.p.P*100.
1255 endif
1256
1257 !___ E arm ______________
1258
1259 ! Go from TRUE to SPECTROMETER quantities by computing target distortions
1260
1261 jones 1.1 ! ... ionization loss correction (if requested) and Coulomb deceleration
1262
1263 if (debug(3)) write(6,*)'mc: e arm stuff0 =',
1264 > orig.e.E,main.target.Eloss(2),spec.e.P
1265 main.SP.e.delta=100.* (orig.e.E - main.target.Eloss(2)
1266 > - main.target.Coulomb - spec.e.P) / spec.e.P
1267
1268 ! ... multiple scattering
1269
1270 if (mc_smear) then
1271 call target_musc(orig.e.p, beta_electron, main.target.teff(2), dangles)
1272 else
1273 dangles(1)=0.0
1274 dangles(2)=0.0
1275 endif
1276
1277 main.SP.e.yptar = orig.e.yptar + dangles(1) + dang_in(1)
1278 main.SP.e.xptar = orig.e.xptar + dangles(2) + dang_in(2)*spec.p.cos_th
1279
1280 ! CASE 1: Using the spectrometer Monte Carlo
1281
1282 jones 1.1 if (using_E_arm_montecarlo) then
1283
1284 ! ... change to E arm spectrometer coordinates (TRANSPORT system),
1285
1286 if (abs(cos(spec.e.phi)).gt.0.0001) then !phi not at +/- pi/2
1287 write(6,*) 'y_E_arm, z_E_arm will be incorrect if spec.e.phi <> pi/2 or 3*pi/2'
1288 write(6,*) 'spec.e.phi=',spec.e.phi,'=',spec.e.phi*180/pi,'degrees'
1289 endif
1290 delta_E_arm = main.SP.e.delta
1291 x_E_arm = -main.target.y
1292 y_E_arm = -main.target.x*spec.e.cos_th - main.target.z*spec.e.sin_th*sin(spec.e.phi)
1293 z_E_arm = main.target.z*spec.e.cos_th + main.target.x*spec.e.sin_th*sin(spec.e.phi)
1294
1295 ! ... Apply spectrometer offset (using spectrometer coordinate system).
1296 x_E_arm = x_E_arm - spec.e.offset.x
1297 y_E_arm = y_E_arm - spec.e.offset.y
1298 z_E_arm = z_E_arm - spec.e.offset.z
1299 dx_E_arm = main.SP.e.xptar - spec.e.offset.xptar
1300 dy_E_arm = main.SP.e.yptar - spec.e.offset.yptar
1301
1302 ! GAW -project to z=0 to compare with reconstructed target positions
1303 jones 1.1 c in_E_arm.z = y_E_arm - z_E_arm*dy_E_arm
1304
1305 if (using_tgt_field) then ! do target field tracking - GAW
1306
1307 c if (debug(6)) then
1308 c write(*,*) '------------------------------'
1309 c write(*,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_E_arm
1310 c endif
1311 c write(6,*) 'electron BEFORE target field:'
1312 c write(6,*) 'xp, yp:', dx_E_arm, dy_E_arm
1313 call track_from_tgt(x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm,
1314 > -spec.e.P*(1+delta_E_arm/100.),Me,-1,ok_E_arm)
1315 c write(6,*) 'electron AFTER target field:'
1316 c write(6,*) 'xp, yp:', dx_E_arm, dy_E_arm
1317 endif
1318 ! GAW - end 99/11/3
1319
1320 ! ........ drift this position back to z=0, the plane through the target center
1321
1322 x_E_arm = x_E_arm - z_E_arm*dx_E_arm
1323 y_E_arm = y_E_arm - z_E_arm*dy_E_arm
1324 jones 1.1 z_E_arm = 0.0
1325
1326 c call print_coord2('virtual track z=0: ',x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm) ! GAW
1327
1328 main.SP.e.z=y_E_arm
1329
1330 ! ... Monte Carlo through the E arm! if we succeed, continue ...
1331 ! ... Here's what's passed:
1332 ! spectrometer central momentum
1333 ! spectrometer central theta angle
1334 ! momentum delta
1335 ! x, y, z, theta, and phi in TRANSPORT coordinates
1336 ! x, y, theta, and phi at the focal plane
1337 ! radiation length of target (NOT USED!)
1338 ! particle mass
1339 ! multiple scattering flag
1340 ! wcs smearing flag
1341 ! resmult=resolution factor for DCs (see simulate.inc)
1342 ! decay flag (hardwired to .false. for electron arm).
1343 ! ok flag
1344
1345 jones 1.1 m2 = me2
1346 pathlen = 0.0
1347 if (electron_arm.eq.1) then
1348 call mc_hms(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm,
1349 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
1350 > me2, mc_smear, mc_smear, .false.,
1351 > tmpfact, frx, fry, ok_E_arm, pathlen,using_tgt_field)
1352 else if (electron_arm.eq.2) then
1353 call mc_sos(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm,
1354 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
1355 > me2, mc_smear, mc_smear, .false.,
1356 > tmpfact, fry, ok_E_arm, pathlen)
1357 else if (electron_arm.eq.3) then
1358 call mc_hrsr(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm,
1359 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
1360 > me2, mc_smear, mc_smear, .false.,
1361 > tmpfact, fry, ok_E_arm, pathlen)
1362 else if (electron_arm.eq.4) then
1363 call mc_hrsl(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm,
1364 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
1365 > me2, mc_smear, mc_smear, .false.,
1366 jones 1.1 > tmpfact, fry, ok_E_arm, pathlen)
1367 else if (electron_arm.eq.5) then
1368 ztemp = (recon.p.z/sin(spec.p.theta))
1369 call mc_calo(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm,
1370 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
1371 > me2, mc_smear, mc_smear, .false.,
1372 > tmpfact, frx,fry, ok_E_arm, pathlen,using_tgt_field,
1373 > main.target.z,ztemp,drift_to_cal)
1374 endif
1375 ntup.resfac=ntup.resfac+tmpfact
1376
1377 if (.not.ok_E_arm) then
1378 if (debug(3)) write(6,*)'mc: ok_E_arm =',ok_E_arm
1379 return
1380 endif
1381
1382 c call print_coord2('recon tgt: ',x_E_arm,y_E_arm,0.,dx_E_arm,dy_E_arm) ! GAW
1383
1384 main.RECON.e.delta = delta_E_arm
1385 main.RECON.e.yptar = dy_E_arm
1386 main.RECON.e.xptar = dx_E_arm
1387 jones 1.1 main.RECON.e.z = y_E_arm
1388 main.FP.e.x = xfp
1389 main.FP.e.dx = dxfp
1390 main.FP.e.y = yfp
1391 main.FP.e.dy = dyfp
1392 main.FP.e.path = pathlen
1393
1394 ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the
1395 ! OUT record
1396
1397 else
1398 main.RECON.e.delta = main.SP.e.delta
1399 main.RECON.e.yptar = main.SP.e.yptar
1400 main.RECON.e.xptar = main.SP.e.xptar
1401 endif
1402
1403 ! Fill recon quantities.
1404 recon.e.delta = main.RECON.e.delta
1405 recon.e.yptar = main.RECON.e.yptar
1406 recon.e.xptar = main.RECON.e.xptar
1407 recon.e.z = main.RECON.e.z
1408 jones 1.1 recon.e.P = spec.e.P*(1.+recon.e.delta/100.)
1409 recon.e.E = recon.e.P
1410 call physics_angles(spec.e.theta,spec.e.phi,recon.e.xptar,
1411 1 recon.e.yptar,recon.e.theta,recon.e.phi)
1412
1413 ! ... correct for energy loss and correct for Coulomb deceleration
1414
1415 if (correct_Eloss) then
1416 call trip_thru_target (2, zero, recon.e.E, recon.e.theta,
1417 & eloss_E_arm, r, Me, 4)
1418 recon.e.E = recon.e.E + eloss_E_arm
1419 endif
1420 recon.e.E = recon.e.E + targ.Coulomb.ave
1421 recon.e.P = recon.e.E
1422 recon.e.delta = (recon.e.P-spec.e.P)/spec.e.P*100.
1423
1424 ! Made it!
1425 success = .true.
1426
1427 return
1428 end
|
Return to
simc.f
CVS log
Up to [HallC] / Poltar