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 record /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 1 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) 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) 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
270 else if (doing_delta) then !Strike (and detect) proton, pion 'recoil'
271 targ.Mtar_struck = Mp
272 targ.Mrec_struck = Mpi
273
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
283 else if(doing_rho) then
284 targ.Mtar_struck = Mp
285 targ.Mrec_struck = Mp
286
287
288 ! ... for normal production, Strike p (n), recoil partile is n(p).
289 ! ... for bound final state, use targ.Mrec if it appears to be OK (same A
290 ! ... A as target, with one n->p or p->n.
291
292 else if (doing_pion) then
293 if (which_pion .eq. 0) then
294 targ.Mtar_struck = Mp
295 gaskelld 1.1 targ.Mrec_struck = Mn
296 else if (which_pion .eq. 1) then
297 targ.Mtar_struck = Mn
298 targ.Mrec_struck = Mp
299 else if (which_pion .eq. 10) then
300 targ.Mtar_struck = targ.M
301 targ.Mrec_struck = targ.Mrec
302 Mrec_guess = targ.M - Mp + Mn
303 else if (which_pion .eq. 11) then
304 targ.Mtar_struck = targ.M
305 targ.Mrec_struck = targ.Mrec
306 Mrec_guess = targ.M - Mn + Mp
307 else
308 stop 'Bad value for which_pion'
309 endif
310
311 if (which_pion .ge. 10) then !Coherent pion production.
312 if (abs(targ.Mrec_struck-Mrec_guess).gt.100.) then
313 targ.Mrec_struck = Mrec_guess
314 write(6,*) 'For coherent pion production, targ.mrec_amu should'
315 write(6,*) 'be the mass of the final nucleus.'
316 gaskelld 1.1 write(6,*) 'I AM REPLACING THE MASS FROM THE INPUT FILE WITH MY GUESS:'
317 write(6,*) 'M_final = M_A + M_N(struck) - M_N(recoil) =',Mrec_guess
318 endif
319 targ.Mrec = 0. !no 'A-1' recoil system for coherent production.
320 endif
321
322
323 ! ... for normal production, Strike p(n), recoil is Lambda/Sigma/Sigma-.
324 ! ... for bound final state, use targ.Mrec if it appears to be OK (same A
325 ! ... A as target, with one n->p or p->n.
326
327 else if (doing_kaon) then !Strike p(n), recoil is L0/S0(S-)
328 if (which_kaon .eq. 0) then ! p(e,eK+)Lambda0
329 targ.Mtar_struck = Mp
330 targ.Mrec_struck = Mlambda
331 else if (which_kaon .eq. 1) then ! p(e,eK+)Sigma0
332 targ.Mtar_struck = Mp
333 targ.Mrec_struck = Msigma0
334 else if (which_kaon .eq. 2) then ! n(e,eK+)Sigma-
335 targ.Mtar_struck = Mn
336 targ.Mrec_struck = Msigma_minus
337 gaskelld 1.1 else if (which_kaon .eq. 10) then !A(e,eK+)A_Lambda0 (bound hypernucleon)
338 targ.Mtar_struck = targ.M
339 targ.Mrec_struck = targ.Mrec
340 Mrec_guess = targ.M - Mp + Mlambda
341 else if (which_kaon .eq. 11) then !A(e,eK+)A_Sigma0 (bound hypernucleon)
342 targ.Mtar_struck = targ.M
343 targ.Mrec_struck = targ.Mrec
344 Mrec_guess = targ.M - Mp + Msigma0
345 else if (which_kaon .eq. 12) then !A(e,eK+)A_Sigma- (bound hypernucleon)
346 targ.Mtar_struck = targ.M
347 targ.Mrec_struck = targ.Mrec
348 Mrec_guess = targ.M - Mn + Msigma_minus
349 else
350 stop 'Bad value for which_kaon'
351 endif
352
353 if (which_kaon .ge. 10) then !Coherent hypernucleus final state.
354 if (abs(targ.Mrec_struck-Mrec_guess).gt.100.) then
355 targ.Mrec_struck = Mrec_guess
356 write(6,*) 'For coherent kaon production, targ.mrec_amu should'
357 write(6,*) 'be the mass of the final hypernucleus.'
358 gaskelld 1.1 write(6,*) 'I AM REPLACING THE MASS FROM THE INPUT FILE WITH MY GUESS:'
359 write(6,*) 'M_hypernucleus = M_A - M_N + M_Y =',Mrec_guess
360 endif
361 targ.Mrec = 0. !no 'A-1' recoil system for coherent production.
362 endif
363 endif
364
365 ! ... Now that we know targ.Mtar_struck, we can improve the precision of
366 ! ... targ.Mrec for deuterium target.
367
368 if (nint(targ.A).eq.2) then !force Mrec.
369 Mrec_guess = Mp + Mn - targ.Mtar_struck
370 if (abs(targ.Mrec-Mrec_guess).gt.0.1) write(6,*)
371 > 'WARNING: forcing recoil mass to M_Nucleon'
372 targ.Mrec = Mrec_guess
373 endif
374
375
376 targ.thick=targ.thick/1000.
377 targ.length=targ.thick/targ.rho
378 targ.angle=targ.angle/degrad
379 gaskelld 1.1 if(targ.Z.lt.2.4) then !liquid target
380 if (abs(targ.angle) .gt. 0.0001) then
381 targ.angle = 0.0
382 write(6,*) 'Forcing target angle to zero for cryotarget.'
383 endif
384 if (targ.can.ne.1 .and. targ.can.ne.2) stop 'bad targ.can value'
385 endif
386 if(sin(targ.angle) .gt. 0.85) then
387 write(6,*) 'BAD targ.angle (0 is perp. to beam, +ve is rotated towards SOS)'
388 write(6,*) 'If the target really is >60 degrees from normal to the beam, then'
389 write(6,*) 'comment out this error checking code.'
390 stop
391 endif
392
393 ! ... dbase field spect_offset
394
395 spec.e.offset.xptar=spec.e.offset.xptar/1000.
396 spec.e.offset.yptar=spec.e.offset.yptar/1000.
397 spec.p.offset.xptar=spec.p.offset.xptar/1000.
398 spec.p.offset.yptar=spec.p.offset.yptar/1000.
399
400 gaskelld 1.1 ! ... dbase fields e_arm_accep, p_arm_accep. Convert angles to radians.
401
402 if (SPedge.e.delta.min.le.-100.0) SPedge.e.delta.min=-99.99
403 if (SPedge.p.delta.min.le.-100.0) SPedge.p.delta.min=-99.99
404 SPedge.e.yptar.min = SPedge.e.yptar.min/1000.
405 SPedge.e.yptar.max = SPedge.e.yptar.max/1000.
406 SPedge.e.xptar.min = SPedge.e.xptar.min/1000.
407 SPedge.e.xptar.max = SPedge.e.xptar.max/1000.
408 SPedge.p.yptar.min = SPedge.p.yptar.min/1000.
409 SPedge.p.yptar.max = SPedge.p.yptar.max/1000.
410 SPedge.p.xptar.min = SPedge.p.xptar.min/1000.
411 SPedge.p.xptar.max = SPedge.p.xptar.max/1000.
412
413 ! ... dbase field simulate
414
415 ! ... Set flags for e,e',hadron radiation tails.
416
417 doing_tail(1) = one_tail.eq.0 .or. one_tail.eq.1 .or.
418 > one_tail.eq.-2 .or. one_tail.eq.-3
419 doing_tail(2) = one_tail.eq.0 .or. one_tail.eq.2 .or.
420 > one_tail.eq.-3 .or. one_tail.eq.-1
421 gaskelld 1.1 doing_tail(3) = one_tail.eq.0 .or. one_tail.eq.3 .or.
422 > one_tail.eq.-1 .or. one_tail.eq.-2
423
424 if (.not.doing_tail(1) .or. .not.doing_tail(2)) then
425 write(6,*) '**WARNING: Shutting off tails 1 or 2 not well defined'
426 endif
427
428 if (abs(one_tail).gt.3 .and. using_rad)
429 1 stop 'Moron! one_tail>3 turns radiation off, but using_rad wants it on.'
430
431 if (.not.using_rad) then
432 do i = 1, 3
433 doing_tail(i) = .false.
434 enddo
435 endif
436
437 if (Egamma_gen_max.gt.0.01) then !use hardwired limits if gen_max > 0
438 hardwired_rad = .true.
439 else
440 hardwired_rad = .false.
441 endif
442 gaskelld 1.1
443 ! ... change angles to radians
444 gen.e.yptar.min=gen.e.yptar.min/1000.
445 gen.e.yptar.max=gen.e.yptar.max/1000.
446 gen.p.yptar.min=gen.p.yptar.min/1000.
447 gen.p.yptar.max=gen.p.yptar.max/1000.
448 gen.e.xptar.min=gen.e.xptar.min/1000.
449 gen.e.xptar.max=gen.e.xptar.max/1000.
450 gen.p.xptar.min=gen.p.xptar.min/1000.
451 gen.p.xptar.max=gen.p.xptar.max/1000.
452
453 ! ... set some flags
454 using_E_arm_montecarlo = spect_mode.ne.1 .and. spect_mode.ne.-1
455 using_P_arm_montecarlo = spect_mode.ne.1 .and. spect_mode.ne.-2
456
457 if (doing_pion .or. doing_kaon .or. doing_delta .or.
458 > (cuts.Em.min.eq.cuts.Em.max) ) then
459 cuts.Em.min = -1.d6
460 cuts.Em.max = 1.d6
461 endif
462
463 gaskelld 1.1 if (abs(deForest_flag).gt.1) stop 'Idiot! check setting of deForest_flag'
464
465 if (correct_Eloss .and. .not.using_Eloss) then
466 write(6,*) 'using_Eloss is false, SO WE ARE FORCING CORRECT_ELOSS TO BE FALSE!!!'
467 correct_Eloss = .false.
468 endif
469
470 if (nint(targ.Z).eq.1) using_Coulomb=.false. !no coulomb corr. for Z=1
471
472 ! ... target
473
474 call target_init(using_Eloss, using_Coulomb, spec.e.theta, spec.p.theta,
475 > spec.p.P, Mh2, Ebeam, spec.e.P)
476
477 ! ... Read in the theory file (for A(e,e'p))
478 if (doing_deuterium .or. doing_heavy) then
479 call theory_init(success)
480 if (.not.success) stop 'THEORY_INIT failed!'
481 endif
482
483 ! ... Read in (and normalize) momentum distribution. Normalizing to one
484 gaskelld 1.1 ! ... may not be correct if there is strength beyond p=p_max (MeV/c)
485
486 if(doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.doing_deutsemi) then
487 if(doing_deutpi .or. doing_deutkaon .or. doing_deutsemi) open(1,file='deut.dat',status='old',form='formatted')
488 if(doing_hepi .or. doing_hekaon) then
489 if (nint(targ.A).eq.3) then
490 open(1,file='he3.dat',status='old',form='formatted')
491 else if (nint(targ.A).eq.4) then
492 open(1,file='he4.dat',status='old',form='formatted')
493 else if (nint(targ.A).eq.12) then
494 open(1,file='c12.dat',status='old',form='formatted')
495 else
496 write(6,*) 'No Momentum Distribution for A = ',targ.A
497 write(6,*) 'Defaulting to carbon momentum distribution'
498 open(1,file='c12.dat',status='old',form='formatted')
499 endif
500 endif
501 do ii=1,2000
502 read(1,*,end=999) pval(ii),mprob(ii)
503 nump=ii
504 enddo
505 gaskelld 1.1 999 continue
506 do ii=1,nump
507 mprob(ii)=mprob(ii)/mprob(nump)
508 enddo
509 close(1)
510 endif
511
512 ! ... Read in (and normalize) spectral function for A>2 electroproduction.
513 ! ... Normalizing to one may not be correct if there is strength beyond
514 ! ... Em_max,Pm_max.
515
516 if(doing_hepi.or.doing_hekaon) then
517 if (nint(targ.A).eq.3) then
518 write(6,*) 'Using the mod version of 3He S.F. rather than Paris.'
519 tmpfile='benharsf_3mod.dat'
520 else if (nint(targ.A).eq.4) then
521 tmpfile='benharsf_4.dat'
522 else if (nint(targ.A).eq.12) then
523 tmpfile='benharsf_12.dat'
524 else if (nint(targ.A).eq.56) then
525 tmpfile='benharsf_56.dat'
526 gaskelld 1.1 else if (nint(targ.A).eq.197) then
527 tmpfile='benharsf_197.dat'
528 else
529 write(6,*) 'No spectral function for A = ',targ.A
530 write(6,*) 'Defaulting to Carbon S.F.'
531 tmpfile='benharsf_12.dat'
532 endif
533 ! Choos proton or neutron spectral function based on targ.Mtar_struck
534 if (abs(targ.Mtar_struck-Mp).le.1.d-6) then
535 call sf_lookup_init(tmpfile,.true.) !proton S.F.
536 else if (abs(targ.Mtar_struck-Mn).le.1.d-6) then
537 call sf_lookup_init(tmpfile,.false.) !neutron S.F.
538 else
539 write(6,*) 'targ.Mtar_struck = ',targ.Mtar_struck
540 write(6,*) 'targ.Mtar_struck not equal to Mp or Mn'
541 write(6,*) 'DEFAULTING TO PROTON SPECTRAL FUNCTION!!!'
542 call sf_lookup_init(tmpfile,.true.) !proton S.F.
543 endif
544 endif
545
546 ! ... Read in input file for kaon cross section model (saghai model).
547 gaskelld 1.1 if (doing_kaon) then !read in model xsec info.
548 open(1,file='weightc.dat',status='old',form='formatted')
549 do i=1,50
550 do j=1,20
551 read(1,*) weightc(j,i)
552 enddo
553 enddo
554 close(1)
555
556 open(1,file='weightd.dat',status='old',form='formatted')
557 do i=1,30
558 do j=1,40
559 do k=1,8
560 read(1,*) weightd(k,j,i)
561 enddo
562 enddo
563 enddo
564 close(1)
565
566 open(3,file='saghai_proton.dat',status='old',form='formatted')
567 do iread=1,10
568 gaskelld 1.1 do iq2=1,11
569 read(3,*) dum1,dum2
570 do iang=1,19
571 read(3,4005) dum3,dum4,dum5,dum6,dum7
572 read(3,4005) zrff1(iread,iq2,iang),ziff1(iread,iq2,iang),
573 1 zrff2(iread,iq2,iang),ziff2(iread,iq2,iang),
574 2 zrff3(iread,iq2,iang),ziff3(iread,iq2,iang)
575 read(3,4005) zrff4(iread,iq2,iang),ziff4(iread,iq2,iang),
576 1 zrff5(iread,iq2,iang),ziff5(iread,iq2,iang),
577 2 zrff6(iread,iq2,iang),ziff6(iread,iq2,iang)
578 enddo
579 enddo
580 enddo
581 4005 format(6e12.4)
582 close(3)
583
584 open(3,file='saghai_sigma0.dat',status='old',form='formatted')
585 do iread=1,20
586 do iq2=1,10
587 do iang=1,19
588 read(3,*) dum1,dum2
589 gaskelld 1.1 read(3,4005) dum3,dum4,dum5,dum6,dum7
590 read(3,4005) zsrff1(iread,iq2,iang),zsiff1(iread,iq2,iang),
591 1 zsrff2(iread,iq2,iang),zsiff2(iread,iq2,iang),
592 2 zsrff3(iread,iq2,iang),zsiff3(iread,iq2,iang)
593 read(3,4005) zsrff4(iread,iq2,iang),zsiff4(iread,iq2,iang),
594 1 zsrff5(iread,iq2,iang),zsiff5(iread,iq2,iang),
595 2 zsrff6(iread,iq2,iang),zsiff6(iread,iq2,iang)
596 enddo
597 enddo
598 enddo
599 close(3)
600
601 endif !kaon input files.
602
603
604
605 ! ... initialize limits on generation, cuts, edges
606
607 call limits_init(H)
608
609 ! ... some announcements
610 gaskelld 1.1
611 if (doing_eep) then
612 if (doing_hyd_elast) then
613 write(6,*) ' ****-------- H(e,e''p) --------****'
614 else if (doing_deuterium) then
615 write(6,*) ' ****-------- D(e,e''p) --------****'
616 else if (doing_heavy) then
617 write(6,*) ' ****-------- A(e,e''p) --------****'
618 else
619 stop 'I don''t have ANY idea what (e,e''p) we''re doing!!!'
620 endif
621
622 else if (doing_semi) then
623 if (doing_semipi) then
624 if(doing_hydsemi) then
625 if(doing_hplus) then
626 write(6,*) ' ****-------- H(e,e''pi+)X --------****'
627 else
628 write(6,*) ' ****-------- H(e,e''pi-)X --------****'
629 endif
630 elseif (doing_deutsemi) then
631 gaskelld 1.1 if(doing_hplus) then
632 write(6,*) ' ****-------- D(e,e''pi+)X --------****'
633 else
634 write(6,*) ' ****-------- D(e,e''pi-)X --------****'
635 endif
636 endif
637
638 else if (doing_semika) then
639 if (doing_hydsemi) then
640 if(doing_hplus) then
641 write(6,*) ' ****-------- H(e,e''K+)X --------****'
642 else
643 write(6,*) ' ****-------- H(e,e''K-)X --------****'
644 endif
645 elseif(doing_deutsemi) then
646 if(doing_hplus) then
647 write(6,*) ' ****-------- D(e,e''K+)X --------****'
648 else
649 write(6,*) ' ****-------- D(e,e''K-)X --------****'
650 endif
651 endif
652 gaskelld 1.1 else
653 stop 'I don''t have ANY idea what (e,e''m)X we''re doing!!!'
654 endif
655
656 else if (doing_delta) then
657 if (doing_hyddelta) then
658 write(6,*) ' ****-------- H(e,e''p)pi --------****'
659 else if (doing_deutpi) then
660 write(6,*) ' ****-------- D(e,e''p)pi --------****'
661 else if (doing_hepi) then
662 write(6,*) ' ****-------- A(e,e''p)pi --------****'
663 else
664 stop 'I don''t have ANY idea what (e,e''p)pi we''re doing!!!'
665 endif
666 else if (doing_pion) then
667 if (doing_hydpi) then
668 if (targ.A .eq. 1) then
669 write(6,*) ' ****-------- H(e,e''pi) --------****'
670 else if (targ.A .ge. 3) then
671 write(6,*) ' ****-------- A(e,e''pi) --------****'
672 endif
673 gaskelld 1.1 else if (doing_deutpi) then
674 write(6,*) ' ****-------- D(e,e''pi) --------****'
675 else if (doing_hepi) then
676 write(6,*) ' ****-------- A(e,e''pi) --------****'
677 else
678 stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!'
679 endif
680 if (which_pion.eq.0 .or. which_pion.eq.10) then
681 write(6,*) ' ****------- pi+ production -------****'
682 else if (which_pion.eq.1 .or. which_pion.eq.11) then
683 write(6,*) ' ****------- pi- production -------****'
684 endif
685 if (which_pion.eq.0 .or. which_pion.eq.1) then
686 write(6,*) ' ****---- Quasifree Production ----****'
687 else if (which_pion.eq.10 .or. which_pion.eq.11) then
688 write(6,*) ' ****---- Coherent Production ----****'
689 endif
690 else if (doing_kaon) then
691 if (doing_hydkaon) then
692 if (targ.A .eq. 1) then
693 write(6,*) ' ****-------- H(e,e''K) --------****'
694 gaskelld 1.1 else if (targ.A .ge. 3) then
695 write(6,*) ' ****-------- A(e,e''K) --------****'
696 endif
697 else if (doing_deutkaon) then
698 write(6,*) ' ****-------- D(e,e''K) --------****'
699 else if (doing_hekaon) then
700 write(6,*) ' ****-------- A(e,e''K) --------****'
701 else
702 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
703 endif
704 if (which_kaon.eq.0) then
705 write(6,*) ' ****----- LAMBDA Quasifree -----****'
706 else if (which_kaon.eq.1) then
707 write(6,*) ' ****----- SIGMA0 Quasifree -----****'
708 else if (which_kaon.eq.2) then
709 write(6,*) ' ****----- SIGMA- Quasifree -----****'
710 else if (which_kaon.eq.10) then
711 write(6,*) ' ****----- LAMBDA Coherent -----****'
712 else if (which_kaon.eq.11) then
713 write(6,*) ' ****----- SIGMA0 Coherent -----****'
714 else if (which_kaon.eq.12) then
715 gaskelld 1.1 write(6,*) ' ****----- SIGMA- Coherent -----****'
716 else
717 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!'
718 endif
719 else if (doing_rho) then
720 if (doing_hydrho) then
721 if (targ.A .eq. 1) then
722 write(6,*) ' ****-------- H(e,e''rho) --------****'
723 else if (targ.A .ge. 3) then
724 write(6,*) ' ****-------- A(e,e''rho) --------****'
725 write(6,*) ' **** ---- Not yet implemented -----****'
726 stop
727 endif
728 else if (doing_deutrho) then
729 write(6,*) ' ****-------- D(e,e''rho) --------****'
730 write(6,*) ' **** ---- Not yet implemented -----****'
731 stop
732 else if (doing_herho) then
733 write(6,*) ' ****-------- A(e,e''rho) --------****'
734 write(6,*) ' **** ---- Not yet implemented -----****'
735 stop
736 gaskelld 1.1 else
737 stop 'I don''t have ANY idea what (e,e''rho) we''re doing!!!'
738 endif
739
740 else if (doing_phsp) then
741 write(6,*) ' ****-------- PHASE SPACE - NO physics, NO radiation --------****'
742 else
743 stop 'I don''t have ANY idea what we''re doing!!!'
744 endif
745
746 if (electron_arm.eq.1) then
747 write(6,*) ' HMS is detecting electrons'
748 else if (electron_arm.eq.2) then
749 write(6,*) ' SOS is detecting electrons'
750 else if (electron_arm.eq.3) then
751 write(6,*) ' HRSR is detecting electrons'
752 else if (electron_arm.eq.4) then
753 write(6,*) ' HRSL is detecting electrons'
754 else if (electron_arm.eq.5) then
755 write(6,*) ' SHMS is detecting electrons (SSA TUNE)'
756 else if (electron_arm.eq.6) then
757 gaskelld 1.1 write(6,*) ' SHMS is detecting electrons (LSA TUNE)'
758 endif
759
760 if (hadron_arm.eq.1) then
761 write(6,*) ' HMS is detecting hadrons'
762 else if (hadron_arm.eq.2) then
763 write(6,*) ' SOS is detecting hadrons'
764 else if (hadron_arm.eq.3) then
765 write(6,*) ' HRSR is detecting hadrons'
766 else if (hadron_arm.eq.4) then
767 write(6,*) ' HRSL is detecting hadrons'
768 else if (hadron_arm.eq.5) then
769 write(6,*) ' SHMS is detecting hadrons (SSA TUNE)'
770 else if (hadron_arm.eq.6) then
771 write(6,*) ' SHMS is detecting hadrons (LSA TUNE)'
772 endif
773
774 if (hadron_arm.eq.electron_arm) then
775 write(6,*) '**WARNING: SIMC works best with two DIFFERENT spectrometers!!!'
776 else if ((hadron_arm+electron_arm).eq.3) then
777 write(6,*) 'Rolf welcomes you to Hall C at Jefferson Lab - Get to work!'
778 gaskelld 1.1 else if (electron_arm.eq.5 .or. hadron_arm.eq.5 .or.
779 > electron_arm.eq.6 .or. hadron_arm.eq.6) then
780 write(6,*) 'Grumpy Dutch Giant welcomes you to Hall C++ at JLab12 - Get to work!'
781 else if ((hadron_arm+electron_arm).eq.7) then
782 write(6,*) 'Kees welcomes you to Hall A at Jefferson Lab - Enjoy your stay!'
783 else
784 write(6,*) '**WARNING: SIMC works best when both spectrometers are in the same hall!!'
785 endif
786
787 if (doing_phsp) write(6,*) '**WARNING: phsp option is not really implemented - buyer beware!!'
788 if (doing_kaon .and. .not. doing_decay) write(6,*) 'NOTE: not doing decay, so a decay weight will be applied to WEIGHT'
789 if (doing_semika .and. .not. doing_decay) write(6,*) 'NOTE: not doing decay, so a decay weight will be applied to WEIGHT'
790 if(doing_deutsemi.and.do_fermi) write(6,*) 'NOTE: Fermi motion enabled for semi-incusive production from deuterium'
791 if (.not.using_rad) write(6,*) 'NOTE: Will NOT be applying radiative corrections'
792 if (.not.using_E_arm_montecarlo) write(6,*) 'NOTE: Will NOT be running events through the E arm Monte Carlo'
793 if (.not.using_P_arm_montecarlo) write(6,*) 'NOTE: Will NOT be running events through the P arm Monte Carlo'
794 if (.not.using_Eloss) write(6,*) 'NOTE: Will NOT be calculating energy loss in the target'
795 if (.not.using_Coulomb) write(6,*) 'NOTE: Will NOT be calculating Coulomb correction (default for Hydrogen target)'
796 if (.not.correct_Eloss) write(6,*) 'NOTE: Will NOT correct reconstructed data for energy loss'
797 if (.not.correct_raster) write(6,*) 'NOTE: Will NOT use raster terms in reconstruction'
798
799 gaskelld 1.1 return
800 end
801
802 !---------------------------------------------------------------
803
804 subroutine regallvars()
805
806 implicit none
807 include 'simulate.inc'
808 include 'radc.inc'
809
810 integer*4 ierr
811 integer*4 regparmint, regparmdouble
812 integer*4 regparmintarray, regparmdoublearray
813 integer*4 regparmstring
814 ! integer*4 regparmreal
815
816 * RESTSW
817
818 if (debug(2).eq.1) write(6,*)'regallvars: entering'
819 ierr = regparmint('mc_smear',mc_smear,0)
820 gaskelld 1.1 ierr = regparmint('electron_arm',electron_arm,0)
821 ierr = regparmint('hadron_arm',hadron_arm,0)
822 ierr = regparmint('use_first_cer',use_first_cer,0)
823 ierr = regparmstring('extra_dbase_file',extra_dbase_file,0)
824
825 * EXPERIMENT
826
827 ierr = regparmdouble('Ebeam',Ebeam,0)
828 ierr = regparmdouble('dEbeam',dEbeam,0)
829 ierr = regparmdouble('EXPER.charge',EXPER.charge,0)
830 ierr = regparmint('doing_kaon',doing_kaon,0)
831 ierr = regparmint('which_kaon',which_kaon,0)
832 ierr = regparmint('doing_pion',doing_pion,0)
833 ierr = regparmint('which_pion',which_pion,0)
834 ierr = regparmint('doing_delta',doing_delta,0)
835 ierr = regparmint('doing_semi', doing_semi,0)
836 ierr = regparmint('doing_hplus', doing_hplus,1)
837 ierr = regparmint('doing_rho',doing_rho,0)
838 ierr = regparmint('doing_decay',doing_decay,0)
839 ierr = regparmdouble('ctau',ctau,0)
840
841 gaskelld 1.1 * DEBUG
842
843 ierr = regparmintarray('debug',debug,5,0)
844
845 * TARGET
846
847 ierr = regparmdouble('targ.A',targ.A,0)
848 ierr = regparmdouble('targ.Z',targ.Z,0)
849 ierr = regparmdouble('targ.mass_amu',targ.mass_amu,0)
850 ierr = regparmdouble('targ.mrec_amu',targ.mrec_amu,0)
851 ierr = regparmdouble('targ.rho',targ.rho,0)
852 ierr = regparmdouble('targ.thick',targ.thick,0)
853 ierr = regparmdouble('targ.xoffset',targ.xoffset,0)
854 ierr = regparmdouble('targ.yoffset',targ.yoffset,0)
855 ierr = regparmint('targ.fr_pattern',targ.fr_pattern,0)
856 ierr = regparmdouble('targ.fr1',targ.fr1,0)
857 ierr = regparmdouble('targ.fr2',targ.fr2,0)
858 ierr = regparmdouble('targ.zoffset',targ.zoffset,0)
859 ierr = regparmdouble('targ.angle',targ.angle,0)
860 ierr = regparmdouble('targ.abundancy',targ.abundancy,0)
861 ierr = regparmint('targ.can',targ.can,0)
862 gaskelld 1.1
863 * E_ARM_MAIN
864
865 ierr = regparmdouble('spec.e.P',spec.e.P,0)
866 ierr = regparmdouble('spec.e.theta',spec.e.theta,0)
867
868 * P_ARM_MAIN
869
870 ierr = regparmdouble('spec.p.P',spec.p.P,0)
871 ierr = regparmdouble('spec.p.theta',spec.p.theta,0)
872
873 * E_ARM_OFFSET
874
875 ierr = regparmdouble('gen.ywid',gen.ywid,0)
876 ierr = regparmdouble('gen.xwid',gen.xwid,0)
877 ierr = regparmdouble('spec.e.offset.x',spec.e.offset.x,0)
878 ierr = regparmdouble('spec.e.offset.y',spec.e.offset.y,0)
879 ierr = regparmdouble('spec.e.offset.z',spec.e.offset.z,0)
880 ierr = regparmdouble('spec.e.offset.xptar',spec.e.offset.xptar,0)
881 ierr = regparmdouble('spec.e.offset.yptar',spec.e.offset.yptar,0)
882
883 gaskelld 1.1 * P_ARM_OFFSET
884
885 ierr = regparmdouble('spec.p.offset.x',spec.p.offset.x,0)
886 ierr = regparmdouble('spec.p.offset.y',spec.p.offset.y,0)
887 ierr = regparmdouble('spec.p.offset.z',spec.p.offset.z,0)
888 ierr = regparmdouble('spec.p.offset.xptar',spec.p.offset.xptar,0)
889 ierr = regparmdouble('spec.p.offset.yptar',spec.p.offset.yptar,0)
890
891 * MISC2INT
892
893 ierr = regparmint('use_expon',use_expon,0)
894 ierr = regparmint('one_tail',one_tail,0)
895 ierr = regparmint('intcor_mode',intcor_mode,0)
896
897 * SIMULATE
898
899 ierr = regparmint('ngen',ngen,0)
900 ierr = regparmint('hard_cuts',hard_cuts,0)
901 ierr = regparmint('using_rad',using_rad,0)
902 ierr = regparmint('spect_mode',spect_mode,0)
903 ierr = regparmint('doing_phsp',doing_phsp,0)
904 gaskelld 1.1 ierr = regparmdouble('cuts.Em.min',cuts.Em.min,0)
905 ierr = regparmdouble('cuts.Em.max',cuts.Em.max,0)
906 ierr = regparmint('using_Eloss',using_Eloss,0)
907 ierr = regparmint('correct_Eloss',correct_eloss,0)
908 ierr = regparmint('correct_raster',correct_raster,0)
909 ierr = regparmint('deForest_flag',deForest_flag,0)
910 ierr = regparmint('rad_flag',rad_flag,0)
911 ierr = regparmint('extrad_flag',extrad_flag,0)
912 ierr = regparmdoublearray('lambda',lambda,3,0)
913 ierr = regparmint('using_cit_generation',using_cit_generation,0)
914 ierr = regparmint('Nntu',Nntu,0)
915 ierr = regparmint('using_Coulomb',using_Coulomb,0)
916 ierr = regparmdouble('dE_edge_test',dE_edge_test,0)
917 ierr = regparmint('use_offshell_rad',use_offshell_rad,0)
918 ierr = regparmdouble('Egamma_gen_max',Egamma_gen_max,0)
919 ierr = regparmint('do_fermi',do_fermi,0)
920 ierr = regparmdouble('pt_b_param',pt_b_param,0)
921 ierr = regparmint('sigc_flag',sigc_flag,0)
922 ierr = regparmint('sigc_nbin',sigc_nbin,10)
923 ierr = regparmdouble('sigc_kin_min',sigc_kin_min,0)
924 ierr = regparmdouble('sigc_kin_max',sigc_kin_max,1.0)
925 gaskelld 1.1 ierr = regparmdouble('sigc_kin_ind',sigc_kin_ind,0.0)
926
927 * E_ARM_ACCEPT
928
929 ierr = regparmdouble('SPedge.e.delta.min',SPedge.e.delta.min,0)
930 ierr = regparmdouble('SPedge.e.delta.max',SPedge.e.delta.max,0)
931 ierr = regparmdouble('SPedge.e.yptar.min',SPedge.e.yptar.min,0)
932 ierr = regparmdouble('SPedge.e.yptar.max',SPedge.e.yptar.max,0)
933 ierr = regparmdouble('SPedge.e.xptar.min',SPedge.e.xptar.min,0)
934 ierr = regparmdouble('SPedge.e.xptar.max',SPedge.e.xptar.max,0)
935
936 * P_ARM_ACCEPT
937
938 ierr = regparmdouble('SPedge.p.delta.min',SPedge.p.delta.min,0)
939 ierr = regparmdouble('SPedge.p.delta.max',SPedge.p.delta.max,0)
940 ierr = regparmdouble('SPedge.p.yptar.min',SPedge.p.yptar.min,0)
941 ierr = regparmdouble('SPedge.p.yptar.max',SPedge.p.yptar.max,0)
942 ierr = regparmdouble('SPedge.p.xptar.min',SPedge.p.xptar.min,0)
943 ierr = regparmdouble('SPedge.p.xptar.max',SPedge.p.xptar.max,0)
944
945 if (debug(2).eq.1) write(6,*)'regallvars: ending'
946 gaskelld 1.1 return
947 end
|
Return to
dbase.f
CVS log
Up to [HallC] / simc_semi