#!/bin/bash

# This procedure runs Qweak Geant version 4A in "interactive mode".
# The executable is assumed to have already been created with the name
# (or a soft link to the executable) "Qweak".  It is further assumed
# that all other required files (xsneut95.dat, , etc.) are
# in the current directory (or that there are links to those files).

# Check operating system name to see if this one is supported.
OSNAME=`uname -s`
if [ $OSNAME != HP-UX -a $OSNAME != SunOS -a $OSNAME != IRIX64 ]; then
    echo The operating sytem $OSNAME is not really supported by $0.
    echo The following may or may not work.
fi

# Determine the name to use for fortran files referenced only by
# logical unit number
if [ $OSNAME = HP-UX ]; then
    FORTPRE=ftn
else
    FORTPRE=fort.
fi

# Create the command input file which is processed by Cern FFKey routines.
# Note that commands which set up physics, cuts, reactions, etc. must appear
# in this file.
echo Create the file ftn43...
cat<<-endg0 >ftn43
	C Select the scattering mode.  non-0 => forward-proton mode
	C 0 => backward electron mode.
	MODE 1
	
	C Select the number of events and the print frequency in the
	C case of a non-interactive version)
	XINT 100000 100
	
	C Set random seed for GEANT
	RNDM 1114225425 156485520
	C Random seed that returns the Sigfer subroutine error immediately:
	C RNDM 1552361693 1
	C A different random seed that returns the Sigfer subroutine error
	C immediately
	C RNDM 1286875785 1
	
	C Initialize volume statistics
	STAT 0
	
	C +----------+
	C | Geometry |
	C +----------+
	
	C Option to read the geometry from an euclid file (default is OFF)
	C (If this is ON, then some of the switches below which define
	C the Geometry are ignored.)
	KEUC ON
	
	C Z position of the target (used in UGEOM and GUKINE). This must
	C agree with the euclid file.  NOTE:  The nominal position of
	C the target center is +44 in forward proton mode and -44 in
	C backward electron mode (see euclid file for more details).
	ZTAR -650.
	
	C Z Position of the beam line window (used in GUKINE reaction 6).
	C This switch is used for a special check of the effects of
	C vacuum windows installed in the beamline.  If such a window
	C is included in the geometry file (euclid file) this value must
	C agree with its position there.
	ZWIN 250.

	C Thickness (dz) of the beam line window (dZ) (used in GUKINE
	C reaction 6). This switch is used for a special check of the
	C effects of vacuum windows installed in the beamline.  If such
	C a window is included in the geometry file (euclid file) this
	C value must agree with its position there.
	TWIN 0.3
	
	C The following options are used only if the geometry is NOT
	C defined using the euclid file (KEUC OFF) (i.e. the geometry
	C is defined in UGEOM)
	
	C HTAR defines the target half length (in cm) (default = 10.)
	C (Used only if KEUC OFF) NOT USED
	HTAR 17.5
	
	C Number of octants (default = 8) (Used only if KEUC OFF)
	NUMB_OCT 8
	
	C FPD determines the controls of the focal plane detector
	C (Used only if KEUC OFF)
	C     0		No detector
	C     1		Single (unsegmented) focal plane detector per sector.
	C     2		Segmented focal plane detector (default).
	FPD  2
	
	C CED determines the controls the cryostadt exit window detector
	C type (Used only if KEUC OFF)
	C     0		No detector
	C     1		Segmented exit window detector (default).
	CED 0

	C +----------------+
	C | Magnetic Field |
	C +----------------+
	
	C Magnetic field type
	C     0		No field.
	C     1		Field is toroidal in the volume covered by the
	C		coils
	C     2		Field interpolation with the Q_weak map
	C     3		Field interpolation with a map
	C     4		Field interpolation with the Q_weak and minitorus maps 
	MAGF 4

	C Magnetic field value.  This is a constant used to scale
	C the magnetic field.  If a map is used (MAGF = 3), then it
	C is a constant that multiplies all fields in the map.
	C for the current maps (November 10, 1998) a value of
	C 0.9804 should be used for the proton case.  Otherwise, this
	C is a field in Tesla (?)
	C BFIL -16.95
	C BFIL -18.95
	C BFIL -0.9804

         
C	   BFIL 0  
C          BFIL 0.7232
C          BFIL 0.9017857
C          BFIL 0.95
C	   BFIL 1.0804
C	   BFIL 1.0
C	   BFIL 0.0

C	   BFIL 1.048
C	   BFIL 1.025
C	   BFIL 0.980

C	   BFIL 1.040
	   BFIL 1.000
C	   BFIL 1.126
C	   BFIL 1.161

	C Minitorus field value.  This is a constant used to scale
	C the magnetic field.  If minitor map is used (MAGF = 4), then it
	C is a constant that multiplies all fields in the minitorus map.
	C for the current maps (September 2003) a value of
	C ??? should be used.
C        MFIL = .1325= 50/(26*3.81^2)
C        MFIL = -.091875
C        MFIL = -.059625

C   factor     MFIL
C    .45       .059625
C    .50       .06625
C    .55       .072875
C    .60       .0795
C    .65       .086125

C    .44       .0583
C    .47       .062275
C    .50       .06625
C    .53       .070225 
C    .56       .0742 
C    .59       .078175
C    .62       .08215 
C    .68       .0901
C    .71       .094075
C    .74       .09805
C    .77       .102025

C        MFIL = -.0583
C        MFIL = -.062275
C        MFIL = -.06625
C        MFIL = -.070225
C        MFIL = -.0742
        MFIL = -.078175
C        MFIL = -.08215
C        MFIL = -.0901
C        MFIL = -.094075
C        MFIL = -.09805
C        MFIL = -.102025
C        MFIL = -.1325


C        MFIL = 0.

C  Optimum current density of ~470 A/cm^2 for down
C        MFIL = 0.0795
C        MFIL = -0.0795

	C Minitorus z position - must be less than -330 or minitorus 
        C field will overlap main torus field.
C        ZPOS = -425.0
C        ZPOS = -465.0
C  For neven's collimator
C       ZPOS = -505.0
       ZPOS = -495.31

	C Maximum field value (Used only when a field map is selected
	C i.e. MAGF = 3)
	BMAX 5000.

	C +--------------------------+
	C | Reactions and kinematics |
	C +--------------------------+

	C REAC selects the initial reaction type
	C     0		Event generated from a file (In this the value used
	C		in the TRIG command MUST agree with the number of
	C		events in the file.  There is no check!)
	C		The structure of the event file is (for each event):
	C		  1)  The number of particles generated in the event
	C		  2)  The Q**2 for the event.
	C		Then for each particle in the event:
	C		  1)  The Geant particle type number
	C		  2)  The 3 momentum components at the event vertex,
	C		      Px Py Pz
	C     1		Electromagnetic shower/ Cerenkov Studies
	C     2		electron-proton and electron-neutron elastic       
	C		scattering
	C     3		electron-proton and electron-neutron inelastic
	C		scattering: e+p ==> e+n+(pi+)
	C     4		electron-proton and electron-neutron inelastic
	C		scattering: e+p ==> e+p+(pi+)+(pi-)
	C     5		electron-proton and electron-neutron inelastic
	C		scattering: e+p ==> e+p+(pi0)
	C     8		Moller electrons
	C     9		elastic e+Al ==> e+Al
	C     >10	REAC-10 is the particle type, X in an inclusive
	C		(e,X) reaction
	REAC 2
	
	C ISOT selects the weighting used for an elastic reaction (REAC=2)
	C     0		GENBOD is used to select the initial kinematics
	C		and weight.
	C     1		An "isotropic" distribution is used with weights
	C		based on the cross section.
	ISOT 1

	C IPTY defines the interaction point type.
	C     0		Z position of IP randomly chosen along the
	C		length of the target
	C     1		Z position fixed at end of target away from
	C		detectors.
	C     2		Z position fixed at center of target
	C     3		Z position fixed at detector end of target.
	C     4		Z position from dowstream endcap of target.
	C     5		Z position from upstream endcap of target.
	IPTY 0
	
	C Tracking of primary particles.  Logical values, TRAP(i) are defined
	C for each Geant particle type (up #14 - proton).  If a particle
	C is created in the primary reaction and the corresponding TRAP logical
	C is ON, the particle will be tracked (all are tracked by default).
	C Track only the proton--#14:
	TRAP 1=OFF
	TRAP 2=OFF
	TRAP 3=ON
	TRAP 4=OFF
	TRAP 5=OFF
	TRAP 6=OFF
	TRAP 7=OFF
	TRAP 8=OFF
	TRAP 9=OFF
	TRAP 10=OFF
	TRAP 11=OFF
	TRAP 12=OFF
	TRAP 13=OFF
	TRAP 14=OFF

	C Tracking of secondary particles  Logical values, TRAS(i) are defined
	C for each Geant particle type (up to #14 - proton).  If a particle
	C is created as a secondary particle (during the shower following the
	C initial reaction) and the corresponding TRAS value is ON, the
	C particle will be tracked (all are tracked by default).
	C Do not Track any secondaries:
	 TRAS 1=OFF
	 TRAS 2=OFF
	 TRAS 3=ON
	  TRAS 4=OFF
	  TRAS 5=OFF
	  TRAS 6=OFF
	  TRAS 7=OFF
	  TRAS 8=OFF
	  TRAS 9=OFF
	  TRAS 10=OFF
	  TRAS 11=OFF
	  TRAS 12=OFF
	  TRAS 13=OFF
	  TRAS 14=OFF

	C Incident electron beam momentum (GeV/c).  For forward-proton
	C mode this should be 3.0 GeV.  For backward electron mode, the
	C value depends on Q**2 as follows:
	C    Q**2	PLAB
	C    0.2	0.335
	C    0.3	0.428
	C    0.4	0.512
	C    0.5	0.590
	C    1.0	0.935
           PLAB 1.165
	CPROT        PLAB 5.0

	C Electron beam position (cm)
	BPOX 0.
	BPOY 0.
	
	C Electron beam x and y width (SIGMA) (cm)
	BEAX -0.2
	BEAY -0.2
	
	C Beam divergence.  If > 0, then a uniform distribution is used
	C up to the value given (assumed to be in degrees).  If < 0,
	C a Gaussian distribution is used with the value given as the FWHM.
	C (This is used only for the electromagnetic shower reaction,
	C REAC = 1)
	BDIV 0.
	
	C ROTR allows one to rotate all trajectories produced in
	C a reaction to place them at a given PHI angle
	C specified in degrees.  Angles greater than 360 or
	C less than -360 mean no rotation (the default).
	ROTR 361.
		
	C BDIR defines the direction of the beam.  With the introduction
	C of the geometry definition in a euclid+ file, this command
	C is essentially obsolete.
	C     0		Beam in the +Z direction (default, and the usual
	C		direction for proton mode or when using a euclid+
	C		file in either mode)
	C     1		Beam in the -Z direction
	BDIR 0

	C Cut on THETA angle for proton elastic (REAC = 2) or on X THETA
	C for (e,X) inclusive reaction (REAC > 10) (in rad).  This must
	C be different from 0 (>0.03) for the inclusive reaction (for
	C the weight calculation). Note that setting limits on THETA
	C of 70+-7.5 gives a limit on Q**2 of 0.1 to 0.5 for elastic.
	C Congruent with q**2 of .1 to .5 (as stated above, 62.5-77.5 deg):
	TLIM 0. 3.14159
	
	C Cut on THETA angles for electron for the elastic and inelastic case (REAC = 2 and 3)
	C (in rad)
	C Theta min and max for photon rates from ep - mollers
	C	 EROT  .005555 1.238

	C7-28 What was downloaded:
C        EROT  0.1047198 0.2094395
C        EROT  0.0872665 0.2094395
C 3 to 16 degrees for extreme rays thru collimator
C        EROT  0.0523599 0.2792526

C 3 to 19 degrees for extreme rays thru first collimator
C        EROT  0.0523599 0.3665191

C 0.01 to 50 degrees to test
        EROT  0.0001745 0.8726646

C Big theta, phi
C        EROT  0.01745 .3491

C       Theta = 6 degrees
C        EROT  0.104 0.105 

C       Theta = 10 degrees
C        EROT  0.174 0.175 


	CPROT        EROT  1.5 3.14
	C	 EROT  0.0872665 0.1919862

	C Theta max for photon rates from ep - mollers
	C	 EROT  1.1 1.4
	C Theta min for photon rates from ep - mollers
	C	 EROT  .004 .006

	C Limit on Q**2 for the elastic case (REAC = 2)
	QLIM 0.0 0.15
	CPROT	QLIM 0.5 10.15

	C Limit for proton PHI (GENBOD generation) or electron PHI
	C (isotropic generation) for elastic and inelastic reaction (REAC = 2 and 3) or
	C on X PHI in (e,x) inclusive reaction (REAC > 10) (in rad.)

C Standard phi range         
	PLIM -0.270526 0.270526

C 3 octants phi
C	PLIM -0.9774 0.9774

C Big theta, phi
C	PLIM -0.87266 0.87266

C	PLIM -1.18 1.18


C	PLIM -0.01 0.01

	C7-11        PLIM -0.2530727 0.2530727
	C PROT         PLIM 2.89 3.39
	C PLIM -0.2617994 0.2617994
	C PLIM -0.2443461 0.2443461
	C PLIM -0.2094395 0.2094395
	C PLIM -0. 0.
	C PLIM 0.3839724 0.4014257
	C PLIM -0.1919862 0.1919862

	C CUTS allows one to change kinetic energy cuts.
	C The order of cuts and defaults is:
	C     CUTGAM   Cut for gammas (0.001 GeV)
	C     CUTELE   Cut for electrons (0.001 GeV)
	C     CUTNEU   Cut for neutral hadrons (0.01 GeV)
	C     CUTHAD   Cut for charged hadrons (0.01 GeV)
	C     CUTMUO   Cut for muons (0.01 GeV)
	C     BCUTE    Cut for electron brems. (CUTGAM)
	C     BCUTM    Cut for muon brems. (CUTGAM)
	C     DCUTE    Cut for electron delta-rays. (10 TeV)
	C     DCUTM    Cut for muon delta-rays. (10 TeV)
	C     PPCUTM   Cut for e+e- pairs by muons. (0.01 GeV)
	C     TOFMAX   Time of flight cut (1.E+10 sec)
	C     GCUTS    5 user words (0.)
	C
	C CUTS  0.00050 0.00050 0.00050 0.00050 0.00050
	C CUTS  0.0050 0.0050 0.0050 0.0050 0.0050
	CCUTS  0.00001 .001 0.01 0.01 0.01
	CUTS  0.001 .001 0.01 0.01 0.01
	
	C +------------------------------+
	C | Reactions Included in Shower |
	C +------------------------------+
	C Include Rayleigh scattering of photons
	C     0		No Rayleigh scattering
	C     1		Rayleigh scattering (default)
	RAYL 0
	
	C DCAY controls decay process
	C     0		No decay effect
	C     1		Decay with generation of secondaries
	C     2		Decay without generation of secondaries
        DCAY 1

	C ANNI controls positron annihilation process
	C     0		No positron annihilation effect
	C     1		Positron annihilation with generation of
	C		secondaries
	C     2		same without generation of secondaries
        ANNI 1

	C PFIS controls photofission process
	C     0		No photofission
	C     1		Photofission with generation of secondaries
	C     2		Same without generation of secondaries
        PFIS 1

	C BREM controls Bremsstrahlung process
	C     0		No Bremsstrahlung effect
	C     1		Bremsstrahlung with generation of secondaries
	C     2		Same without generation of secondaries
	BREM 1

	C COMP controls Compton scattering process
	C     0		No Compton scattering
	C     1		Compton scattering with generation of secondaries
	C     2		Same without generation of secondaries
	COMP 1

	C DRAY controls delta rays process
	C     0		No delta rays effect
	C     1		Delta rays with generation of secondaries
	C     2		Same without generation of secondaries
	C (Note that this is automatically set to 0 if LOSS=2)
	DRAY 1

	C STRA controls straggling for thin layers
	C     0		No straggling in effect
	C     1		Straggling for thin layers in effect
	STRA 0
	
	C SYNC controls synchrotron radiation generation
	C     0		No synchrotron radiation
	C     1		Synchrotron radiation included
	SYNC 0
	
	C HADR controls hadronic interactions
	C     0		No hadronic interactions.
	C     1		Hadronic interactions. Secondaries processed.
	C     2		Hadronic interactions. No secondaries stored.
	HADR 1

	C LOSS controls energy loss
	C     0		No energy loss
	C     1		Unrestricted energy loss fluctuations
	C     2		Complete energy loss fluctuations
	C     3		Same as 1
	C     4		No energy loss fluctuations
	LOSS 2
	
	C MULS controls multiple scattering
	C     0		No multiple scattering
	C     1		Moliere or Coulomb scattering
	C     2		Moliere or Coulomb scattering
	C     3		Gaussian scattering
        MULS 1

	C MUNU controls muon nuclear interactions
	C     0		No muon-nuclear interactions
	C     1		Nuclear interactions. Secondaries processed.
	C     2		Nuclear interactions. Secondaries not processed.
	MUNU 1

	C PAIR controls the pair production mechanism
	C     0		No pair production.
	C     1		Pair production. Secondaries processed.
	C     2		Pair production. No secondaries stored.
	PAIR 1

	C PHOT controls the photo electric effect
	C     0		No photo electric effect
	C     1		Photo effect. Electron processed.
	C     2		Photo effect. No electron stored.
	PHOT 0
	
	C +-----------------+
	C | Control of Hits |
	C +-----------------+
	C GEOM (This is used only if the geometry is NOT defined using a
	C "euclid file" (KEUC = OFF))
	GEOM 'GSCI' 'GDET' 'TARG' 'GLOB' 'BEAM'
	
	C SETS (This is used only if the geometry is NOT defined using a
	C "euclid file" (KEUC = OFF))
	C SETS 'GSCI' 'GDET' 'GLOB' 'TARG' 'DECA'
	SETS 'DECA'
	
	C PRIN is used to selectively print the data structures.  Additional
	C keys may be supplied to print the following
	C     MATE	Materials (standard and user defined)
	C     PART	Particle constants
	C     ROTM	Rotation matrices
	C     TMED	Tracking media
	C     VOLU	Volumes
	C     SETS	Detectors defined for hits and digitizations
	C     JXYZ	Track parameters at each step
	C     VERT	Primary vertices at each event
	C     KINE	Primary tracks kinematics at each event
	C     HITS	Hits at each event
	C     DIGI	Digitisations at each event
	C PRIN 'VOLU' 'MATE' 'TMED' 'VERT' 'HITS' 'KINE' 'ROTM' 'JXYZ' 'DIGI'
	C PRIN 'VERT' 'HITS' 'KINE' 'JXYZ' 'DIGI'
	
	C PLOT 'GLOB' 'GDET'
	
	C AUTO controls automatic computation of STMIN, STEMAX, DEEMAX
	C and TMAXFD for tracking material definition.
	C     0		Tracking media parameters taken from the argument
	C		list of GSTMED (or the euclid TMED command)
	C     1		Tracking media parameters calculated by GEANT
	AUTO 1
	
	C DEBU is used to turn on debugging output.  Three arguments can
	C be provided (defaults are 0 for all three).
	C     IDEMIN	First event to debug
	C     IDEMAX	Last event to debug
	C     ITEST	Print control frequency
	 DEBU 1 100 
	
	C SWIT sets switches to control debug output. The following values
	C are recognized:
	C SWIT	Value	Meaning
	C  1	2	Call GPGKIN to list particles generated
	C		during step
	C  2	1	Call GSXYZ to store point in JXYZ
	C	2	Call GPCXYZ to list step by step printed
	C		debug output
	C	3	Call GDCXYZ to do interactive drawing of
	C		trajectories (see switch 4, value = 3)
	C	4	The particle trajectory is drawn when the particle
	C		stops via GDTRAK
	C  4	3	Do not call GDCXYZ to draw neutral particles
	C  9	12345	Print debug info from GGCLOS (memory management?)
	C 	123456789 Print debug info from GTNEXT
	SWIT 1 2
	SWIT 2 4
	SWIT 4 1
	C SWIT 9 123456789

	STOP
	endg0

# Create the file containing "string" input (typically file names).
echo Create the file ftn48...
cat<<-endg0 >ftn48
	ep_m59_collu4_b1000_he_r5.nt, Name of the n-tuple
	sn-dir.dat,  	Name of the field directory file
	celeg_event.dat,	Name of the event file
	opt_collu4.euclid,   Name of the euclid file
	sm_gap_mintor_r5.dat,   Name of the minitorus map file
	endg0

# The files xsneut95.dat and flukaaf.dat must be available.  The
# program will look first in the current directory, and then in the
# directory $CERN_ROOT/lib/ (where CERN_ROOT is an environment variable)
# for these files.  If the files are actually stored elsewhere, you
# might want to create soft links to them here.
#ln -s /cern/95a/lib/xsneut95.dat xsneut95.dat

# Now run Qweak.  Note that we disable the automatic reading of
# the user's ~/.pawlogon.kumac because some paw commands are misinterpreted
# by interactive Geant.
./Qweak -n

# Clean-up -- remove old (temporary) input and useless output files
#rm ftn43
#rm ftn48
rm paw.metafile
rm last.kumac
rm last.kumacold
#rm flukaerr.dat

# Remove links (if any) that were created to get to the data files
#rm xsneut95.dat

# Save the random seed file
mv ${FORTPRE}47 Qweak.ranout