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
|
401 gaskelld 1.1 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
422 gaskelld 1.1 ! ... 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 > one_tail.eq.-3 .or. one_tail.eq.-1
443 gaskelld 1.1 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 1 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 endif
464 gaskelld 1.1
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
485 gaskelld 1.1 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 ! ... Read in (and normalize) momentum distribution. Normalizing to one
506 gaskelld 1.1 ! ... 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 enddo
527 gaskelld 1.1 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 tmpfile='benharsf_56.dat'
548 gaskelld 1.1 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 ! ... Read in input file for kaon cross section model (saghai model).
569 gaskelld 1.1 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 do iread=1,10
590 gaskelld 1.1 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 read(3,*) dum1,dum2
611 gaskelld 1.1 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 ! ... some announcements
632 gaskelld 1.1
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 elseif (doing_deutsemi) then
653 gaskelld 1.1 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 endif
674 gaskelld 1.1 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 endif
695 gaskelld 1.1 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 write(6,*) ' ****-------- H(e,e''K) --------****'
716 gaskelld 1.1 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 else if (which_kaon.eq.12) then
737 gaskelld 1.1 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 stop
758 gaskelld 1.1 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 else if (electron_arm.eq.6) then
779 gaskelld 1.1 write(6,*) ' SHMS is detecting electrons (LSA TUNE)'
|
816 gaskelld 1.1 else
817 write(6,*) '**WARNING: SIMC works best when both spectrometers are in the same hall!!'
818 endif
819
820 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 gaskelld 1.1 subroutine regallvars()
838
839 implicit none
840 include 'simulate.inc'
841 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).eq.1) 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 gaskelld 1.1 * EXPERIMENT
859
860 ierr = regparmdouble('Ebeam',Ebeam,0)
861 ierr = regparmdouble('dEbeam',dEbeam,0)
862 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
|
898 gaskelld 1.1
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 * 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 gaskelld 1.1 * 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 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 gaskelld 1.1 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 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 gaskelld 1.1 ierr = regparmdouble('sigc_kin_ind',sigc_kin_ind,0.0)
|