1 gaskelld 1.1 subroutine dbase_read(H)
2
3 ! dbase_read reads the input file and sets the run-time flags for the run.
4 ! Note that there are four INDEPENDENT ways to run SIMC.
5 !
6 ! 1. doing_eep: (e,e'p) subcases:doing_hyd_elast, doing_deuterium, doing_heavy
7 !
8 ! 2. doing_kaon:(e,e'K) subcases:doing_hydkaon, doing_deutkaon,doing_hekaon.
9 ! which_kaon= 0/ 1/ 2 for Lambda/Sigam0/Sigma- quasifree.
10 ! which_kaon=10/11/12 for Lambda/Sigam0/Sigma- coherent production.
11 ! For bound hypernuclear states, set doing_hydkaon true for ALL
12 ! targets (i.e. treat as a heavy proton) but with targ.Mtar_struck
13 ! and targ.Mrec_struck set appropriatly.
14 !
15 ! 3. doing_pion:(e,e'p) subcases:doing_hydpi, doing_deutpi,doing_hepi.
16 ! which_pion = 0( 1) gives pi+ (pi-) quasifree production.
17 ! which_pion = 10(11) gives pi+ (pi-) coherent production.
18 ! This sets doing_hydpi true for ALL targets (i.e. treat as
19 ! a heavy proton) but with targ.Mtar_struck and targ.Mrec_struck
20 ! set appropriatly.
21 !
22 gaskelld 1.1 ! 4. doing_phsp:Generate acceptance with radiation and cross section disabled,
23 ! use doing_kaon or doing_pion to set hadron mass, then
24 ! set doing_eep,doing_kaon and doing_pion to FALSE.
25 ! May or may not be properly implemented.
26 !
27 ! 5. doing_delta:(e,e'p)pi subcases:doing_hyddelta,doing_deutdelta,doing_hedelta.
28 ! Identical to doing_pion, except for changes in mass parameters and
29 ! cross section model. Initial implementation will be for proton target
30 ! only (generation well be general, but cross section model will be
31 ! a 'shortcut' version to start with.
32 ! 6. doing_semi: H(e,e'pi)X (doing_semipi) and H(e,e'k)X (doing_semika)
33 ! 7. doing_rho: H(e,e'rho)p
34
35 implicit none
36 include 'radc.inc'
37 include 'histograms.inc'
38 include 'simulate.inc'
39
40 real*8 dum1,dum2,dum3,dum4,dum5,dum6,dum7
41 real*8 mrec_guess
42
43 gaskelld 1.1 integer*4 iread, iq2, iang
44 integer*4 i, j, k, ii
45 integer*4 ierr, thload, thbook
46 logical success
47 character filename*80,tmpfile*80
48 character dbase_file*60 !needs to be shorter than filename
49
50 type (histograms):: H
51
52 ! ... Register CTP vars
53
54 call regallvars
55
56 ! ... read the dbase file.
57
58 dbase_file=' '
59 extra_dbase_file=' '
60 write(6,*) 'Enter the input filename (assumed to be in infiles directory)'
61 read(5,'(a)') dbase_file
62
63 j=index(dbase_file,'/')
64 gaskelld 1.1 if (j.eq.0) then !no path given
65 write(filename,'(a)') 'infiles/'//dbase_file
66 else
67 write(filename,'(a)') dbase_file
68 endif
69 i=index(filename,'.')
70 j=index(filename,'/')
71 if(i.eq.0) then !add .inp if not included in filename
72 i=index(filename,' ')
73 if(i+2.le.len(filename)) write(filename(i:),'(''.inp'')')
74 endif
75 write(6,'(a10,a69)')'filename=',filename
76 if (i.gt.1) base=filename(j+1:i-1)
77
78 ! ... load and book input file
79
80 ierr = thload(filename)
81 if (ierr.ne.0) stop 'Loading problem! Not going to try again...wouldnt be prudent.'
82 ierr = thbook()
83 if (ierr.ne.0) stop ' Booking problem! Not going to try again...wouldnt be prudent'
84
85 gaskelld 1.1 ! ... read the secondary dbase file.
86
87 if (extra_dbase_file.ne.' ') then !new filename from dbase file
88 write(filename,'(a)') 'infiles/'//extra_dbase_file
89 i=index(filename,'.')
90 j=index(filename,'/')
91 if(i.eq.0) then !add .inp if not included in filename
92 i=index(filename,' ')
93 i f(i+2.le.len(filename)) write(filename(i:),'(''.inp'')')
94 endif
95 write(6,'(a10,a69)')'filename=',filename
96 ierr = thload(filename)
97 if (ierr.ne.0) stop ' Loading problem! Not going to try again...wouldnt be prudent.'
98 ierr = thbook()
99 if (ierr.ne.0) stop ' Booking problem! Not going to try again...wouldnt be prudent'
100 endif !extra dbase input file
101
102 C DJG: Ugly hack! This must come before the test on doing_pion
103 if(doing_pion .and. doing_semi) then
104 doing_semipi=.true.
105 doing_semika=.false.
106 gaskelld 1.1 doing_pion=.false.
107 endif
108 if(doing_kaon .and. doing_semi) then
109 doing_semika=.true.
110 doing_semipi=.false.
111 doing_kaon=.false.
112 endif
113 C DJG:
114
115 ! ... dbase field experiment.
116 if (doing_pion) then
117 Mh=Mpi
118 if (nint(targ%A).eq.1 .and. which_pion.eq.1)
119 > stop 'Pi- production from Hydrogen not allowed!'
120 if (nint(targ%A).le.2 .and. which_pion.ge.10)
121 > stop 'Coherent production from Hydrogen/Deuterium not allowed!'
122 if (nint(targ%A).eq.3 .and. which_pion.eq.11)
123 > stop 'Coherent Pi- production from 3He not allowed!'
124 if (nint(targ%A).eq.4 .and. which_pion.ge.10)
125 > stop 'Coherent production from 4He not allowed!'
126 doing_hydpi = (nint(targ%A).eq.1)
127 gaskelld 1.1 doing_deutpi = (nint(targ%A).eq.2)
128 doing_hepi = (nint(targ%A).ge.3)
129 doing_eep = .false.
130
131 * quasifree production is default (which_pion=0,1 for pi+/pi0)
132 * Treat (A+pi) final state as production from heavy proton (which_pion=10,11)
133 if (which_pion.ge.10) then
134 doing_hydpi = .true.
135 doing_deutpi = .false.
136 doing_hepi = .false.
137 endif
138
139 else if (doing_kaon) then
140 Mh=Mk
141 if (nint(targ%A).eq.1 .and. which_kaon.eq.2)
142 > stop 'Sigma- production from Hydrogen not allowed!'
143 if (nint(targ%A).le.2 .and. which_kaon.ge.10)
144 > stop 'Coherent production from Hydrogen/Deuterium not allowed!'
145 doing_hydkaon = (nint(targ%A).eq.1)
146 doing_deutkaon = (nint(targ%A).eq.2)
147 doing_hekaon = (nint(targ%A).ge.3)
148 gaskelld 1.1 doing_eep = .false.
149
150 * quasifree production is default (which_kaon=0,1,2 for lambda/sigma/sigma-).
151 * Treat (A+pi) final state as production from heavy proton (+10 to which_kaon)
152 if (which_kaon.ge.10) then
153 doing_hydkaon = .true.
154 doing_deutkaon = .false.
155 doing_hekaon = .false.
156 endif
157
158 else if (doing_delta) then
159 Mh=Mp
160 if (nint(targ%A).ge.2)
161 > write(6,*) 'WARNING: Delta cross section model only set up for proton target!'
162 doing_hyddelta = (nint(targ%A).eq.1)
163 doing_deutdelta = (nint(targ%A).eq.2)
164 doing_hedelta = (nint(targ%A).ge.3)
165 doing_eep = .false.
166
167 else if (doing_semi) then
168 if(doing_semipi) then
169 gaskelld 1.1 Mh=Mpi
170 elseif(doing_semika) then
171 Mh=Mk
172 endif
173 doing_hydsemi = (nint(targ%A).eq.1)
174 doing_deutsemi = (nint(targ%A).eq.2)
175 if(doing_hydsemi.and.do_fermi) then
176 write(6,*) 'WARNING: Cannot do Fermi motion for Hydrogen!'
177 write(6,*) 'Do you mean to be running deuterium?'
178 do_fermi = .false.
179 endif
180 doing_eep=.false.
181
182 else if (doing_rho) then
183 Mh = Mrho
184 doing_hydrho = (nint(targ%A).eq.1)
185 doing_deutrho = (nint(targ%A).eq.2)
186 doing_herho = (nint(targ%A).eq.3)
187 doing_eep=.false.
188
189 else !doing_eep if nothing else set.
190 gaskelld 1.1 Mh=Mp
191 doing_eep = .true.
192 doing_hyd_elast = (nint(targ%A).eq.1)
193 doing_deuterium = (nint(targ%A).eq.2)
194 doing_heavy = (nint(targ%A).ge.3)
195 endif
196
197 Mh2 = Mh*Mh
198 Mh2_final = Mh2
199
200 if (doing_phsp) then
201 rad_flag=0
202 doing_eep=.false. !need to set one of these to determine Mh,
203 doing_pion=.false. !but doing_phsp is independent of others.
204 doing_kaon=.false.
205 doing_delta=.false.
206 doing_rho=.false.
207 endif
208
209 ! ... dbase field kinematics_main
210
211 gaskelld 1.1 dEbeam=Ebeam*dEbeam/100.
212 spec%e%theta=abs(spec%e%theta)/degrad
213 spec%e%cos_th=cos(spec%e%theta)
214 spec%e%sin_th=sin(spec%e%theta)
215 spec%p%theta=abs(spec%p%theta)/degrad
216 spec%p%cos_th=cos(spec%p%theta)
217 spec%p%sin_th=sin(spec%p%theta)
218
219 ! ... Get phi for spectrometers. HMS/HRSr are on right, phi=3*pi/2
220 if (electron_arm.eq.1 .or. electron_arm.eq.3 .or. electron_arm.eq.7) then
221 spec%e%phi = 3*pi/2.
222 else if (electron_arm.eq.2 .or. electron_arm.eq.4 .or.
223 > electron_arm.eq.5 .or. electron_arm.eq.6 .or. electron_arm.eq.8) then
224 spec%e%phi = pi/2.
225 else
226 stop 'I dont know what phi should be for the electron arm'
227 endif
228
229 if (hadron_arm.eq.1 .or. hadron_arm.eq.3) then
230 spec%p%phi = 3*pi/2.
231 else if (hadron_arm.eq.2 .or. hadron_arm.eq.4 .or.
232 gaskelld 1.1 > hadron_arm.eq.5 .or. hadron_arm.eq.6) then
233 spec%p%phi = pi/2.
234 else
235 stop 'I dont know what phi should be for the hadron arm'
236 endif
237
238 ! ... dbase field target
239
240 targ%N = targ%A - targ%Z
241 targ%M = targ%mass_amu*amu
242 targ%Mrec = targ%mrec_amu*amu
243
244 ! ... improve precision for hydrogen and deuterium targets. For deuterium,
245 ! ... need to wait to fill targ.Mrec until we know which nucleon was struck.
246
247 if (nint(targ%A).eq.1) then !proton target
248 if (abs(targ%M-Mp).gt.0.1) write(6,*) 'WARNING: forcing target mass to Mp!!!'
249 if (abs(targ%Mrec).gt.0.1) write(6,*) 'WARNING: forcing recoil mass to zero!!!'
250 targ%M = Mp
251 targ%Mrec = 0.
252 else if (nint(targ%A).eq.2) then !deuterium target
253 gaskelld 1.1 if (abs(targ%M-Md).gt.0.1) write(6,*) 'WARNING: forcing target mass to Md!!!'
254 targ%M = Md
255 else !heavy target
256 Mrec_guess = targ%M - Mp !approx A-1 mass, don't know mtar_struck yet
257 if (abs(targ%Mrec-Mrec_guess).gt.100.) then !assume OK if within 100MeV.
258 write(6,*) 'targ.Mrec=',targ%Mrec,' I was expecting closer to',Mrec_guess
259 write(6,*) 'I AM REPLACING YOUR TARG.MREC WITH MY EXPECTATION!!!!!'
260 targ%Mrec = Mrec_guess
261 endif
262 endif
263
264
265 ! ... set target/recoil masses.
266 if (doing_eep) then !Strike proton, no 'recoil' particle
267 targ%Mtar_struck = Mp
268 targ%Mrec_struck = 0.0
269 sign_hadron=1.0
270 else if (doing_delta) then !Strike (and detect) proton, pion 'recoil'
271 targ%Mtar_struck = Mp
272 targ%Mrec_struck = Mpi
273 sign_hadron=1.0
274 gaskelld 1.1 else if(doing_semi) then ! For now, just assuming proton mass
275 targ%Mtar_struck = Mp
276 targ%Mrec_struck = Mp !must have at LEAST a recoiling proton
277 !Probably more...
278 if(doing_deutsemi) then
279 targ%Mtar_struck = (Mp+Mn)/2.0
280 targ%Mrec_struck = (Mp+Mn)/2.0
281 endif
282 if(doing_hplus) then
283 sign_hadron=1.0
284 else
285 sign_hadron=-1.0
286 endif
287
288 else if(doing_rho) then
289 targ%Mtar_struck = Mp
290 targ%Mrec_struck = Mp
291 if(doing_hplus) then
292 sign_hadron=1.0
293 else
294 sign_hadron=-1.0
295 gaskelld 1.1 endif
296
297
298 ! ... for normal production, Strike p (n), recoil partile is n(p).
299 ! ... for bound final state, use targ.Mrec if it appears to be OK (same A
300 ! ... A as target, with one n->p or p->n.
301
302 else if (doing_pion) then
303 if (which_pion .eq. 0) then
304 targ%Mtar_struck = Mp ! H(e,e'pi+)n or D(e,e'pi+)nn
305 targ%Mrec_struck = Mn
306 sign_hadron=1.0
307 else if (which_pion .eq. 1) then
308 targ%Mtar_struck = Mn ! D(e,e'pi-) pp
309 targ%Mrec_struck = Mp
310 sign_hadron=-1.0
311 else if (which_pion .eq. 10) then
312 targ%Mtar_struck = targ%M ! A(e,e'pi+)A'
313 targ%Mrec_struck = targ%Mrec
314 Mrec_guess = targ%M - Mp + Mn
315 sign_hadron=1.0
316 gaskelld 1.1 else if (which_pion .eq. 11) then
317 targ%Mtar_struck = targ%M ! A(e,e'pi-)A'
318 targ%Mrec_struck = targ%Mrec
319 Mrec_guess = targ%M - Mn + Mp
320 sign_hadron=-1.0
321 else
322 stop 'Bad value for which_pion'
323 endif
324
325 if (which_pion .ge. 10) then !Coherent pion production.
326 if (abs(targ%Mrec_struck-Mrec_guess).gt.100.) then
327 targ%Mrec_struck = Mrec_guess
328 write(6,*) 'For coherent pion production, targ.mrec_amu should'
329 write(6,*) 'be the mass of the final nucleus.'
330 write(6,*) 'I AM REPLACING THE MASS FROM THE INPUT FILE WITH MY GUESS:'
331 write(6,*) 'M_final = M_A + M_N(struck) - M_N(recoil) =',Mrec_guess
332 endif
333 targ%Mrec = 0. !no 'A-1' recoil system for coherent production.
334 endif
335
336
337 gaskelld 1.1 ! ... for normal production, Strike p(n), recoil is Lambda/Sigma/Sigma-.
338 ! ... for bound final state, use targ.Mrec if it appears to be OK (same A
339 ! ... A as target, with one n->p or p->n.
340
341 else if (doing_kaon) then !Strike p(n), recoil is L0/S0(S-)
342 if (which_kaon .eq. 0) then ! p(e,eK+)Lambda0
343 targ%Mtar_struck = Mp
344 targ%Mrec_struck = Mlambda
345 sign_hadron=1.0
346 else if (which_kaon .eq. 1) then ! p(e,eK+)Sigma0
347 targ%Mtar_struck = Mp
348 targ%Mrec_struck = Msigma0
349 sign_hadron=1.0
350 else if (which_kaon .eq. 2) then ! n(e,eK+)Sigma-
351 targ%Mtar_struck = Mn
352 targ%Mrec_struck = Msigma_minus
353 sign_hadron=1.0
354 else if (which_kaon .eq. 10) then !A(e,eK+)A_Lambda0 (bound hypernucleon)
355 targ%Mtar_struck = targ%M
356 targ%Mrec_struck = targ%Mrec
357 Mrec_guess = targ%M - Mp + Mlambda
358 gaskelld 1.1 sign_hadron=1.0
359 else if (which_kaon .eq. 11) then !A(e,eK+)A_Sigma0 (bound hypernucleon)
360 targ%Mtar_struck = targ%M
361 targ%Mrec_struck = targ%Mrec
362 Mrec_guess = targ%M - Mp + Msigma0
363 sign_hadron=1.0
364 else if (which_kaon .eq. 12) then !A(e,eK+)A_Sigma- (bound hypernucleon)
365 targ%Mtar_struck = targ%M
366 targ%Mrec_struck = targ%Mrec
367 Mrec_guess = targ%M - Mn + Msigma_minus
368 sign_hadron=1.0
369 else
370 stop 'Bad value for which_kaon'
371 endif
372
373 if (which_kaon .ge. 10) then !Coherent hypernucleus final state.
374 if (abs(targ%Mrec_struck-Mrec_guess).gt.100.) then
375 targ%Mrec_struck = Mrec_guess
376 write(6,*) 'For coherent kaon production, targ.mrec_amu should'
377 write(6,*) 'be the mass of the final hypernucleus.'
378 write(6,*) 'I AM REPLACING THE MASS FROM THE INPUT FILE WITH MY GUESS:'
379 gaskelld 1.1 write(6,*) 'M_hypernucleus = M_A - M_N + M_Y =',Mrec_guess
380 endif
381 targ%Mrec = 0. !no 'A-1' recoil system for coherent production.
382 endif
383 endif
384
385 ! ... Now that we know targ.Mtar_struck, we can improve the precision of
386 ! ... targ.Mrec for deuterium target.
387
388 if (nint(targ%A).eq.2) then !force Mrec.
389 Mrec_guess = Mp + Mn - targ%Mtar_struck
390 if (abs(targ%Mrec-Mrec_guess).gt.0.1) write(6,*)
391 > 'WARNING: forcing recoil mass to M_Nucleon'
392 targ%Mrec = Mrec_guess
393 endif
394
395
396 targ%thick=targ%thick/1000.
397 targ%length=targ%thick/targ%rho
398 targ%angle=targ%angle/degrad
399 targ_Bangle = targ_Bangle/degrad
400 gaskelld 1.1 targ_Bphi = targ_Bphi/degrad
401 if(targ%Z.lt.2.4) then !liquid target
402 if (abs(targ%angle) .gt. 0.0001) then
403 targ%angle = 0.0
404 write(6,*) 'Forcing target angle to zero for cryotarget.'
405 endif
406 if (targ%can.ne.1 .and. targ%can.ne.2) stop 'bad targ.can value'
407 endif
408 if(sin(targ%angle) .gt. 0.85) then
409 write(6,*) 'BAD targ.angle (0 is perp. to beam, +ve is rotated towards SOS)'
410 write(6,*) 'If the target really is >60 degrees from normal to the beam, then'
411 write(6,*) 'comment out this error checking code.'
412 stop
413 endif
414
415 ! ... dbase field spect_offset
416
417 spec%e%offset%xptar=spec%e%offset%xptar/1000.
418 spec%e%offset%yptar=spec%e%offset%yptar/1000.
419 spec%p%offset%xptar=spec%p%offset%xptar/1000.
420 spec%p%offset%yptar=spec%p%offset%yptar/1000.
421 gaskelld 1.1
422 ! ... dbase fields e_arm_accep, p_arm_accep. Convert angles to radians.
423
424 if (SPedge%e%delta%min.le.-100.0) SPedge%e%delta%min=-99.99
425 if (SPedge%p%delta%min.le.-100.0) SPedge%p%delta%min=-99.99
426 SPedge%e%yptar%min = SPedge%e%yptar%min/1000.
427 SPedge%e%yptar%max = SPedge%e%yptar%max/1000.
428 SPedge%e%xptar%min = SPedge%e%xptar%min/1000.
429 SPedge%e%xptar%max = SPedge%e%xptar%max/1000.
430 SPedge%p%yptar%min = SPedge%p%yptar%min/1000.
431 SPedge%p%yptar%max = SPedge%p%yptar%max/1000.
432 SPedge%p%xptar%min = SPedge%p%xptar%min/1000.
433 SPedge%p%xptar%max = SPedge%p%xptar%max/1000.
434
435 ! ... dbase field simulate
436
437 ! ... Set flags for e,e',hadron radiation tails.
438
439 doing_tail(1) = one_tail.eq.0 .or. one_tail.eq.1 .or.
440 > one_tail.eq.-2 .or. one_tail.eq.-3
441 doing_tail(2) = one_tail.eq.0 .or. one_tail.eq.2 .or.
442 gaskelld 1.1 > one_tail.eq.-3 .or. one_tail.eq.-1
443 doing_tail(3) = one_tail.eq.0 .or. one_tail.eq.3 .or.
444 > one_tail.eq.-1 .or. one_tail.eq.-2
445
446 if (.not.doing_tail(1) .or. .not.doing_tail(2)) then
447 write(6,*) '**WARNING: Shutting off tails 1 or 2 not well defined'
448 endif
449
450 if (abs(one_tail).gt.3 .and. using_rad)
451 > stop 'Moron! one_tail>3 turns radiation off, but using_rad wants it on.'
452
453 if (.not.using_rad) then
454 do i = 1, 3
455 doing_tail(i) = .false.
456 enddo
457 endif
458
459 if (Egamma_gen_max.gt.0.01) then !use hardwired limits if gen_max > 0
460 hardwired_rad = .true.
461 else
462 hardwired_rad = .false.
463 gaskelld 1.1 endif
464
465 ! ... change angles to radians
466 gen%e%yptar%min=gen%e%yptar%min/1000.
467 gen%e%yptar%max=gen%e%yptar%max/1000.
468 gen%p%yptar%min=gen%p%yptar%min/1000.
469 gen%p%yptar%max=gen%p%yptar%max/1000.
470 gen%e%xptar%min=gen%e%xptar%min/1000.
471 gen%e%xptar%max=gen%e%xptar%max/1000.
472 gen%p%xptar%min=gen%p%xptar%min/1000.
473 gen%p%xptar%max=gen%p%xptar%max/1000.
474
475 ! ... set some flags
476 using_E_arm_montecarlo = spect_mode.ne.1 .and. spect_mode.ne.-1
477 using_P_arm_montecarlo = spect_mode.ne.1 .and. spect_mode.ne.-2
478
479 if (doing_pion .or. doing_kaon .or. doing_delta .or.
480 > (cuts%Em%min.eq.cuts%Em%max) ) then
481 cuts%Em%min = -1.d6
482 cuts%Em%max = 1.d6
483 endif
484 gaskelld 1.1
485 if (abs(deForest_flag).gt.1) stop 'Idiot! check setting of deForest_flag'
486
487 if (correct_Eloss .and. .not.using_Eloss) then
488 write(6,*) 'using_Eloss is false, SO WE ARE FORCING CORRECT_ELOSS TO BE FALSE!!!'
489 correct_Eloss = .false.
490 endif
491
492 if (nint(targ%Z).eq.1) using_Coulomb=.false. !no coulomb corr. for Z=1
493
494 ! ... target
495
496 call target_init(using_Eloss, using_Coulomb, spec%e%theta, spec%p%theta,
497 > spec%p%P, Mh2, Ebeam, spec%e%P)
498
499 ! ... Read in the theory file (for A(e,e'p))
500 if (doing_deuterium .or. doing_heavy) then
501 call theory_init(success)
502 if (.not.success) stop 'THEORY_INIT failed!'
503 endif
504
505 gaskelld 1.1 ! ... Read in (and normalize) momentum distribution. Normalizing to one
506 ! ... may not be correct if there is strength beyond p=p_max (MeV/c)
507
508 if(doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.doing_deutsemi) then
509 if(doing_deutpi .or. doing_deutkaon .or. doing_deutsemi) open(1,file='deut.dat',status='old',form='formatted')
510 if(doing_hepi .or. doing_hekaon) then
511 if (nint(targ%A).eq.3) then
512 open(1,file='he3.dat',status='old',form='formatted')
513 else if (nint(targ%A).eq.4) then
514 open(1,file='he4.dat',status='old',form='formatted')
515 else if (nint(targ%A).eq.12) then
516 open(1,file='c12.dat',status='old',form='formatted')
517 else
518 write(6,*) 'No Momentum Distribution for A = ',targ%A
519 write(6,*) 'Defaulting to carbon momentum distribution'
520 open(1,file='c12.dat',status='old',form='formatted')
521 endif
522 endif
523 do ii=1,2000
524 read(1,*,end=999) pval(ii),mprob(ii)
525 nump=ii
526 gaskelld 1.1 enddo
527 999 continue
528 do ii=1,nump
529 mprob(ii)=mprob(ii)/mprob(nump)
530 enddo
531 close(1)
532 endif
533
534 ! ... Read in (and normalize) spectral function for A>2 electroproduction.
535 ! ... Normalizing to one may not be correct if there is strength beyond
536 ! ... Em_max,Pm_max.
537
538 if(doing_hepi.or.doing_hekaon) then
539 if (nint(targ%A).eq.3) then
540 write(6,*) 'Using the mod version of 3He S.F. rather than Paris.'
541 tmpfile='benharsf_3mod.dat'
542 else if (nint(targ%A).eq.4) then
543 tmpfile='benharsf_4.dat'
544 else if (nint(targ%A).eq.12) then
545 tmpfile='benharsf_12.dat'
546 else if (nint(targ%A).eq.56) then
547 gaskelld 1.1 tmpfile='benharsf_56.dat'
548 else if (nint(targ%A).eq.197) then
549 tmpfile='benharsf_197.dat'
550 else
551 write(6,*) 'No spectral function for A = ',targ%A
552 write(6,*) 'Defaulting to Carbon S.F.'
553 tmpfile='benharsf_12.dat'
554 endif
555 ! Choos proton or neutron spectral function based on targ.Mtar_struck
556 if (abs(targ%Mtar_struck-Mp).le.1.d-6) then
557 call sf_lookup_init(tmpfile,.true.) !proton S.F.
558 else if (abs(targ%Mtar_struck-Mn).le.1.d-6) then
559 call sf_lookup_init(tmpfile,.false.) !neutron S.F.
560 else
561 write(6,*) 'targ%Mtar_struck = ',targ%Mtar_struck
562 write(6,*) 'targ%Mtar_struck not equal to Mp or Mn'
563 write(6,*) 'DEFAULTING TO PROTON SPECTRAL FUNCTION!!!'
564 call sf_lookup_init(tmpfile,.true.) !proton S.F.
565 endif
566 endif
567
568 gaskelld 1.1 ! ... Read in input file for kaon cross section model (saghai model).
569 if (doing_kaon) then !read in model xsec info.
570 open(1,file='weightc.dat',status='old',form='formatted')
571 do i=1,50
572 do j=1,20
573 read(1,*) weightc(j,i)
574 enddo
575 enddo
576 close(1)
577
578 open(1,file='weightd.dat',status='old',form='formatted')
579 do i=1,30
580 do j=1,40
581 do k=1,8
582 read(1,*) weightd(k,j,i)
583 enddo
584 enddo
585 enddo
586 close(1)
587
588 open(3,file='saghai_proton.dat',status='old',form='formatted')
589 gaskelld 1.1 do iread=1,10
590 do iq2=1,11
591 read(3,*) dum1,dum2
592 do iang=1,19
593 read(3,4005) dum3,dum4,dum5,dum6,dum7
594 read(3,4005) zrff1(iread,iq2,iang),ziff1(iread,iq2,iang),
595 1 zrff2(iread,iq2,iang),ziff2(iread,iq2,iang),
596 2 zrff3(iread,iq2,iang),ziff3(iread,iq2,iang)
597 read(3,4005) zrff4(iread,iq2,iang),ziff4(iread,iq2,iang),
598 1 zrff5(iread,iq2,iang),ziff5(iread,iq2,iang),
599 2 zrff6(iread,iq2,iang),ziff6(iread,iq2,iang)
600 enddo
601 enddo
602 enddo
603 4005 format(6e12.4)
604 close(3)
605
606 open(3,file='saghai_sigma0.dat',status='old',form='formatted')
607 do iread=1,20
608 do iq2=1,10
609 do iang=1,19
610 gaskelld 1.1 read(3,*) dum1,dum2
611 read(3,4005) dum3,dum4,dum5,dum6,dum7
612 read(3,4005) zsrff1(iread,iq2,iang),zsiff1(iread,iq2,iang),
613 1 zsrff2(iread,iq2,iang),zsiff2(iread,iq2,iang),
614 2 zsrff3(iread,iq2,iang),zsiff3(iread,iq2,iang)
615 read(3,4005) zsrff4(iread,iq2,iang),zsiff4(iread,iq2,iang),
616 1 zsrff5(iread,iq2,iang),zsiff5(iread,iq2,iang),
617 2 zsrff6(iread,iq2,iang),zsiff6(iread,iq2,iang)
618 enddo
619 enddo
620 enddo
621 close(3)
622
623 endif !kaon input files.
624
625
626
627 ! ... initialize limits on generation, cuts, edges
628
629 call limits_init(H)
630
631 gaskelld 1.1 ! ... some announcements
632
633 if (doing_eep) then
634 if (doing_hyd_elast) then
635 write(6,*) ' ****-------- H(e,e''p) --------****'
636 else if (doing_deuterium) then
637 write(6,*) ' ****-------- D(e,e''p) --------****'
638 else if (doing_heavy) then
639 write(6,*) ' ****-------- A(e,e''p) --------****'
640 else
641 stop 'I don''t have ANY idea what (e,e''p) we''re doing!!!'
642 endif
643
644 else if (doing_semi) then
645 if (doing_semipi) then
646 if(doing_hydsemi) then
647 if(doing_hplus) then
648 write(6,*) ' ****-------- H(e,e''pi+)X --------****'
649 else
650 write(6,*) ' ****-------- H(e,e''pi-)X --------****'
651 endif
652 gaskelld 1.1 elseif (doing_deutsemi) then
653 if(doing_hplus) then
654 write(6,*) ' ****-------- D(e,e''pi+)X --------****'
655 else
656 write(6,*) ' ****-------- D(e,e''pi-)X --------****'
657 endif
658 endif
659
660 else if (doing_semika) then
661 if (doing_hydsemi) then
662 if(doing_hplus) then
663 write(6,*) ' ****-------- H(e,e''K+)X --------****'
664 else
665 write(6,*) ' ****-------- H(e,e''K-)X --------****'
666 endif
667 elseif(doing_deutsemi) then
668 if(doing_hplus) then
669 write(6,*) ' ****-------- D(e,e''K+)X --------****'
670 else
671 write(6,*) ' ****-------- D(e,e''K-)X --------****'
672 endif
673 gaskelld 1.1 endif
674 else
675 stop 'I don''t have ANY idea what (e,e''m)X we''re doing!!!'
676 endif
677
678 else if (doing_delta) then
679 if (doing_hyddelta) then
680 write(6,*) ' ****-------- H(e,e''p)pi --------****'
681 else if (doing_deutpi) then
682 write(6,*) ' ****-------- D(e,e''p)pi --------****'
683 else if (doing_hepi) then
684 write(6,*) ' ****-------- A(e,e''p)pi --------****'
685 else
686 stop 'I don''t have ANY idea what (e,e''p)pi we''re doing!!!'
687 endif
688 else if (doing_pion) then
689 if (doing_hydpi) then
690 if (targ%A .eq. 1) then
691 write(6,*) ' ****-------- H(e,e''pi) --------****'
692 else if (targ%A .ge. 3) then
693 write(6,*) ' ****-------- A(e,e''pi) --------****'
694 gaskelld 1.1 endif
695 else if (doing_deutpi) then
696 write(6,*) ' ****-------- D(e,e''pi) --------****'
697 else if (doing_hepi) then
698 write(6,*) ' ****-------- A(e,e''pi) --------****'
699 else
700 stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!'
701 endif
702 if (which_pion.eq.0 .or. which_pion.eq.10) then
703 write(6,*) ' ****------- pi+ production -------****'
704 else if (which_pion.eq.1 .or. which_pion.eq.11) then
705 write(6,*) ' ****------- pi- production -------****'
706 endif
707 if (which_pion.eq.0 .or. which_pion.eq.1) then
708 write(6,*) ' ****---- Quasifree Production ----****'
709 else if (which_pion.eq.10 .or. which_pion.eq.11) then
710 write(6,*) ' ****---- Coherent Production ----****'
711 endif
712 else if (doing_kaon) then
713 if (doing_hydkaon) then
714 if (targ%A .eq. 1) then
715 gaskelld 1.1 write(6,*) ' ****-------- H(e,e''K) --------****'
716 else if (targ%A .ge. 3) then
717 write(6,*) ' ****-------- A(e,e''K) --------****'
718 endif
719 else if (doing_deutkaon) then
720 write(6,*) ' ****-------- D(e,e''K) --------****'
721 else if (doing_hekaon) then
722 write(6,*) ' ****-------- A(e,e''K) --------****'
723 else
724 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
725 endif
726 if (which_kaon.eq.0) then
727 write(6,*) ' ****----- LAMBDA Quasifree -----****'
728 else if (which_kaon.eq.1) then
729 write(6,*) ' ****----- SIGMA0 Quasifree -----****'
730 else if (which_kaon.eq.2) then
731 write(6,*) ' ****----- SIGMA- Quasifree -----****'
732 else if (which_kaon.eq.10) then
733 write(6,*) ' ****----- LAMBDA Coherent -----****'
734 else if (which_kaon.eq.11) then
735 write(6,*) ' ****----- SIGMA0 Coherent -----****'
736 gaskelld 1.1 else if (which_kaon.eq.12) then
737 write(6,*) ' ****----- SIGMA- Coherent -----****'
738 else
739 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
740 endif
741 else if (doing_rho) then
742 if (doing_hydrho) then
743 if (targ%A .eq. 1) then
744 write(6,*) ' ****-------- H(e,e''rho) --------****'
745 else if (targ%A .ge. 3) then
746 write(6,*) ' ****-------- A(e,e''rho) --------****'
747 write(6,*) ' **** ---- Not yet implemented -----****'
748 stop
749 endif
750 else if (doing_deutrho) then
751 write(6,*) ' ****-------- D(e,e''rho) --------****'
752 write(6,*) ' **** ---- Not yet implemented -----****'
753 stop
754 else if (doing_herho) then
755 write(6,*) ' ****-------- A(e,e''rho) --------****'
756 write(6,*) ' **** ---- Not yet implemented -----****'
757 gaskelld 1.1 stop
758 else
759 stop 'I don''t have ANY idea what (e,e''rho) we''re doing!!!'
760 endif
761
762 else if (doing_phsp) then
763 write(6,*) ' ****-------- PHASE SPACE - NO physics, NO radiation --------****'
764 else
765 stop 'I don''t have ANY idea what we''re doing!!!'
766 endif
767
768 if (electron_arm.eq.1) then
769 write(6,*) ' HMS is detecting electrons'
770 else if (electron_arm.eq.2) then
771 write(6,*) ' SOS is detecting electrons'
772 else if (electron_arm.eq.3) then
773 write(6,*) ' HRSR is detecting electrons'
774 else if (electron_arm.eq.4) then
775 write(6,*) ' HRSL is detecting electrons'
776 else if (electron_arm.eq.5) then
777 write(6,*) ' SHMS is detecting electrons (SSA TUNE)'
778 gaskelld 1.1 else if (electron_arm.eq.6) then
779 write(6,*) ' SHMS is detecting electrons (LSA TUNE)'
780 else if (electron_arm.eq.7) then
781 write(6,*) ' BIGCAL is detecting electrons (HMS side of beam)'
782 else if (electron_arm.eq.8) then
783 write(6,*) ' BIGCAL is detecting electrons (SOS side of beam)'
784 endif
785
786 if (hadron_arm.eq.1) then
787 write(6,*) ' HMS is detecting hadrons'
788 else if (hadron_arm.eq.2) then
789 write(6,*) ' SOS is detecting hadrons'
790 else if (hadron_arm.eq.3) then
791 write(6,*) ' HRSR is detecting hadrons'
792 else if (hadron_arm.eq.4) then
793 write(6,*) ' HRSL is detecting hadrons'
794 else if (hadron_arm.eq.5) then
795 write(6,*) ' SHMS is detecting hadrons (SSA TUNE)'
796 else if (hadron_arm.eq.6) then
797 write(6,*) ' SHMS is detecting hadrons (LSA TUNE)'
798 else if (hadron_arm.eq.7 .or. hadron_arm.eq.8) then
799 gaskelld 1.1 write(6,*) ' Cannot use Bigcal for hadrons'
800 stop
801 endif
802
803 if (hadron_arm.eq.electron_arm) then
804 write(6,*) '**WARNING: SIMC works best with two DIFFERENT spectrometers!!!'
805 else if ((hadron_arm+electron_arm).eq.3) then
806 write(6,*) 'Rolf welcomes you to Hall C at Jefferson Lab - Get to work!'
807 else if (electron_arm.eq.5 .or. hadron_arm.eq.5 .or.
808 > electron_arm.eq.6 .or. hadron_arm.eq.6) then
809 write(6,*) 'Grumpy Dutch Giant welcomes you to Hall C++ at JLab12 - Get to work!'
810 else if ((hadron_arm+electron_arm).eq.7) then
811 write(6,*) 'Kees welcomes you to Hall A at Jefferson Lab - Enjoy your stay!'
812 else if ( electron_arm.eq.7 .and. hadron_arm .eq. 5) then
813 write(6,*) ' Bigcal and SHMS'
814 else if ( electron_arm.eq.8 .and. hadron_arm .eq. 1) then
815 write(6,*) ' Bigcal and HMS'
816 else
817 write(6,*) '**WARNING: SIMC works best when both spectrometers are in the same hall!!'
818 endif
819
820 gaskelld 1.1 if (doing_phsp) write(6,*) '**WARNING: phsp option is not really implemented - buyer beware!!'
821 if (doing_kaon .and. .not. doing_decay) write(6,*) 'NOTE: not doing decay, so a decay weight will be applied to WEIGHT'
822 if (doing_semika .and. .not. doing_decay) write(6,*) 'NOTE: not doing decay, so a decay weight will be applied to WEIGHT'
823 if(doing_deutsemi.and.do_fermi) write(6,*) 'NOTE: Fermi motion enabled for semi-incusive production from deuterium'
824 if (.not.using_rad) write(6,*) 'NOTE: Will NOT be applying radiative corrections'
825 if (.not.using_E_arm_montecarlo) write(6,*) 'NOTE: Will NOT be running events through the E arm Monte Carlo'
826 if (.not.using_P_arm_montecarlo) write(6,*) 'NOTE: Will NOT be running events through the P arm Monte Carlo'
827 if (.not.using_Eloss) write(6,*) 'NOTE: Will NOT be calculating energy loss in the target'
828 if (.not.using_Coulomb) write(6,*) 'NOTE: Will NOT be calculating Coulomb correction (default for Hydrogen target)'
829 if (.not.correct_Eloss) write(6,*) 'NOTE: Will NOT correct reconstructed data for energy loss'
830 if (.not.correct_raster) write(6,*) 'NOTE: Will NOT use raster terms in reconstruction'
831
832 return
833 end
834
835 !---------------------------------------------------------------
836
837 subroutine regallvars()
838
839 implicit none
840 include 'simulate.inc'
841 gaskelld 1.1 include 'radc.inc'
842
843 integer*4 ierr
844 integer*4 regparmint, regparmdouble
845 integer*4 regparmintarray, regparmdoublearray
846 integer*4 regparmstring
847 ! integer*4 regparmreal
848
849 * RESTSW
850
851 if (debug(2)) write(6,*)'regallvars: entering'
852 ierr = regparmint('mc_smear',mc_smear,0)
853 ierr = regparmint('electron_arm',electron_arm,0)
854 ierr = regparmint('hadron_arm',hadron_arm,0)
855 ierr = regparmint('use_first_cer',use_first_cer,0)
856 ierr = regparmstring('extra_dbase_file',extra_dbase_file,0)
857
858 * EXPERIMENT
859
860 ierr = regparmdouble('Ebeam',Ebeam,0)
861 ierr = regparmdouble('dEbeam',dEbeam,0)
862 gaskelld 1.1 ierr = regparmdouble('EXPER%charge',EXPER%charge,0)
863 ierr = regparmint('doing_kaon',doing_kaon,0)
864 ierr = regparmint('which_kaon',which_kaon,0)
865 ierr = regparmint('doing_pion',doing_pion,0)
866 ierr = regparmint('which_pion',which_pion,0)
867 ierr = regparmint('doing_delta',doing_delta,0)
868 ierr = regparmint('doing_semi', doing_semi,0)
869 ierr = regparmint('doing_hplus', doing_hplus,1)
870 ierr = regparmint('doing_rho',doing_rho,0)
871 ierr = regparmint('doing_decay',doing_decay,0)
872 ierr = regparmdouble('ctau',ctau,0)
873
874 * DEBUG
875
876 ierr = regparmintarray('debug',debug,6,0)
877
878 * TARGET
879
880 ierr = regparmdouble('targ%A',targ%A,0)
881 ierr = regparmdouble('targ%Z',targ%Z,0)
882 ierr = regparmdouble('targ%mass_amu',targ%mass_amu,0)
883 gaskelld 1.1 ierr = regparmdouble('targ%mrec_amu',targ%mrec_amu,0)
884 ierr = regparmdouble('targ%rho',targ%rho,0)
885 ierr = regparmdouble('targ%thick',targ%thick,0)
886 ierr = regparmdouble('targ%xoffset',targ%xoffset,0)
887 ierr = regparmdouble('targ%yoffset',targ%yoffset,0)
888 ierr = regparmint('targ%fr_pattern',targ%fr_pattern,0)
889 ierr = regparmdouble('targ%fr1',targ%fr1,0)
890 ierr = regparmdouble('targ%fr2',targ%fr2,0)
891 ierr = regparmdouble('targ%zoffset',targ%zoffset,0)
892 ierr = regparmdouble('targ%angle',targ%angle,0)
893 ierr = regparmdouble('targ%abundancy',targ%abundancy,0)
894 ierr = regparmint('targ%can',targ%can,0)
895 ierr = regparmdouble('targ_Bangle',targ_Bangle,0)
896 ierr = regparmdouble('targ_Bphi',targ_Bphi,0)
897 ierr = regparmdouble('targ_pol',targ_pol,0)
898
899 * E_ARM_MAIN
900
901 ierr = regparmdouble('spec%e%P',spec%e%P,0)
902 ierr = regparmdouble('spec%e%theta',spec%e%theta,0)
903
904 gaskelld 1.1 * P_ARM_MAIN
905
906 ierr = regparmdouble('spec%p%P',spec%p%P,0)
907 ierr = regparmdouble('spec%p%theta',spec%p%theta,0)
908
909 * E_ARM_OFFSET
910
911 ierr = regparmdouble('gen%ywid',gen%ywid,0)
912 ierr = regparmdouble('gen%xwid',gen%xwid,0)
913 ierr = regparmdouble('spec%e%offset%x',spec%e%offset%x,0)
914 ierr = regparmdouble('spec%e%offset%y',spec%e%offset%y,0)
915 ierr = regparmdouble('spec%e%offset%z',spec%e%offset%z,0)
916 ierr = regparmdouble('spec%e%offset%xptar',spec%e%offset%xptar,0)
917 ierr = regparmdouble('spec%e%offset%yptar',spec%e%offset%yptar,0)
918
919 * P_ARM_OFFSET
920
921 ierr = regparmdouble('spec%p%offset%x',spec%p%offset%x,0)
922 ierr = regparmdouble('spec%p%offset%y',spec%p%offset%y,0)
923 ierr = regparmdouble('spec%p%offset%z',spec%p%offset%z,0)
924 ierr = regparmdouble('spec%p%offset%xptar',spec%p%offset%xptar,0)
925 gaskelld 1.1 ierr = regparmdouble('spec%p%offset%yptar',spec%p%offset%yptar,0)
926
927 * MISC2INT
928
929 ierr = regparmint('use_expon',use_expon,0)
930 ierr = regparmint('one_tail',one_tail,0)
931 ierr = regparmint('intcor_mode',intcor_mode,0)
932
933 * SIMULATE
934
935 ierr = regparmint('ngen',ngen,0)
936 ierr = regparmint('hard_cuts',hard_cuts,0)
937 ierr = regparmint('using_rad',using_rad,0)
938 ierr = regparmint('spect_mode',spect_mode,0)
939 ierr = regparmint('doing_phsp',doing_phsp,0)
940 ierr = regparmdouble('cuts%Em%min',cuts%Em%min,0)
941 ierr = regparmdouble('cuts%Em%max',cuts%Em%max,0)
942 ierr = regparmint('using_Eloss',using_Eloss,0)
943 ierr = regparmint('correct_Eloss',correct_eloss,0)
944 ierr = regparmint('correct_raster',correct_raster,0)
945 ierr = regparmint('deForest_flag',deForest_flag,0)
946 gaskelld 1.1 ierr = regparmint('rad_flag',rad_flag,0)
947 ierr = regparmint('extrad_flag',extrad_flag,0)
948 ierr = regparmdoublearray('lambda',lambda,3,0)
949 ierr = regparmint('using_cit_generation',using_cit_generation,0)
950 ierr = regparmint('Nntu',Nntu,0)
951 ierr = regparmint('using_Coulomb',using_Coulomb,0)
952 ierr = regparmdouble('dE_edge_test',dE_edge_test,0)
953 ierr = regparmint('use_offshell_rad',use_offshell_rad,0)
954 ierr = regparmdouble('Egamma_gen_max',Egamma_gen_max,0)
955 ierr = regparmint('do_fermi',do_fermi,0)
956 ierr = regparmdouble('pt_b_param',pt_b_param,0)
957 ierr = regparmint('sigc_flag',sigc_flag,0)
958 ierr = regparmint('sigc_nbin',sigc_nbin,10)
959 ierr = regparmdouble('sigc_kin_min',sigc_kin_min,0)
960 ierr = regparmdouble('sigc_kin_max',sigc_kin_max,1.0)
961 ierr = regparmdouble('sigc_kin_ind',sigc_kin_ind,0.0)
962 ierr = regparmint('using_tgt_field',using_tgt_field,0)
963 ierr = regparmstring('tgt_field_file',tgt_field_file,0)
964 ierr = regparmdouble('drift_to_cal',drift_to_cal,0)
965
966 * E_ARM_ACCEPT
967 gaskelld 1.1
968 ierr = regparmdouble('SPedge%e%delta%min',SPedge%e%delta%min,0)
969 ierr = regparmdouble('SPedge%e%delta%max',SPedge%e%delta%max,0)
970 ierr = regparmdouble('SPedge%e%yptar%min',SPedge%e%yptar%min,0)
971 ierr = regparmdouble('SPedge%e%yptar%max',SPedge%e%yptar%max,0)
972 ierr = regparmdouble('SPedge%e%xptar%min',SPedge%e%xptar%min,0)
973 ierr = regparmdouble('SPedge%e%xptar%max',SPedge%e%xptar%max,0)
974
975 * P_ARM_ACCEPT
976
977 ierr = regparmdouble('SPedge%p%delta%min',SPedge%p%delta%min,0)
978 ierr = regparmdouble('SPedge%p%delta%max',SPedge%p%delta%max,0)
979 ierr = regparmdouble('SPedge%p%yptar%min',SPedge%p%yptar%min,0)
980 ierr = regparmdouble('SPedge%p%yptar%max',SPedge%p%yptar%max,0)
981 ierr = regparmdouble('SPedge%p%xptar%min',SPedge%p%xptar%min,0)
982 ierr = regparmdouble('SPedge%p%xptar%max',SPedge%p%xptar%max,0)
983
984 if (debug(2)) write(6,*)'regallvars: ending'
985 return
986 end
|
Return to
dbase.f
CVS log
Up to [HallC] / simc_gfortran