LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added to repository.

Browse Source
This commit is contained in:
prokscha 2004-06-17 13:05:47 +00:00
commit a8b8352190
24 changed files with 4695 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
geant3/inp/geant_lemsr.input Normal file
View File

@ -0,0 +1,114 @@
plot 1
rung 3022
trig 10
rndm 9871 59367
COMM 1 2 3 4 5 6 7 8 9
sets 5 -74 1 9 0 0 0 0 0
bfie 0. 0. 0. 0. 0.
spar ' e+' ' e-' ' gam'
cuts 0.0001 0.0001 0.01 0.01 0.00001
COMM CutGam CutEle CutNeu CutHad CutMuo
ntva 'init' 'time' 'flag' 'scie' 'scoe' 'ltrb'
geom 'mcps' 'samp' 'efla' 'cryo' 'crsh' 'saal' 'ru11'
COMM geom 'mcps' 'efla' 'ru11' 'mcp2' 'mcpa'
view 1 60 40 0 10 10 4 4
stop
This is the Format Free Input File for GEANT_LEMSR.FOR:
=======================================================
how to use the data cards:
COMM ---> comment line, no effect when reading
plot 0 ---> no graphics.
plot 1 ---> graphics on screen, and also information on detector hits.
plot 2 ---> graphics on screen and PS file.
rung 12 ---> define GEANT run number; will be appended to the NTP filename
trig 100000 ---> number of events to be processed; on PSICL1 10^5 takes
up to 3:30 min CPU time.
1 2 3 4 5 6 7 8 9
sets 5 -74 2 5 0 1 0 2 5 ---> first number : particle ID, 5 = mu+ (2 = e+)
ID = 500; user defined mu+
with anistropic decay prob.
second number: if == 0, michel-distributed at
MCP2 position
if ne 0, start position is at z=sets(2)
(here 74 cm upstream of MCP tube center)
with beam momentum
sets(3).sets(4)sets(5),
this is 2.50 MeV/c in this example
sets(6): gaussian beam divergence (in degree)
sets(7): if = 1, read initial values from file
MC$INP:GEANT_LEMSR_INPUT.DAT
sets(8): gaussian dp/p distribution,
sigma_p = sets(8).sets(9)%
spar ' e+' ' e-' ' gam' --> select secondary particles to be tracked.
BFIE 0. 0. 100. 90. 50. ---> Sample B-field (Gauss), Bx, By, Bz, and muon polarization=90%
and muon polarization=50% in sample.
ntva 'time' 'flag' ---> see comments below
view 1 40 40 0 10 10 3 3 ---> in case of graphical output define view angle,
(60)(40)(0)(10)(10)(2)(2) position and scaling of graph.
default values view(1): if == 1 then use the command lines
parameters for graphical output;
otherwise use default parameter
view(2): theta angle (degree)
view(3): phi angle
view(4): psi angle
view(5): x0, x-position in HIGZ window
view(6): y0, y-position in HIGZ window
view(7): 10*x-scaling, view(7)=3 <==> xscal=0.3
view(8): 10*y-scaling, view(8)=3 <==> yscal=0.3
geom 'mcp2' 'efla' ---> see comment below
stop ---> stops reading this file and continues the program.
c------------------------------------------------------------------------------
c
c flags for geometrical setup; it is possible to mount or remove the
c
c 'mcp2' MCP2 plates
c 'mcpa' MCP2 anode
c 'mcps' MCP2 stainless steel vacuum tube
c 'samp' sample (this is at the moment the sample holder)
c 'cryo' material of the cryostat without shielding
c 'crsh' cryo shield
c 'efla' the 100 CF flange at the end of the sample/MCP2 tube
c 'ru11' use MCP setup of Run 11 (delay line anode)
c 'ru10' use MCP setup of Run 10 and before (WSZ detector)
c has no influence for cryo geometry
c 'saal' uses Al instead of Cu for sample holder plates
c
c the volumes are activated when printing
c
c geom 'mcp2' 'efla' for example when data cards are requested
c
c---------------------------------------
c
c flags to define NTuple variables for CWN
c
c 'init' write initial parameter to NT
c 'code' initial particle ID to NT
c 'endp' end position and energy of e+
c 'time' t-mcp and t-sci
c 'scip' position when particle hits inner Sc.
c 'scie' energy loss of particles in inner sc.
c 'scop' position when particle hits outer Sc.
c 'scoe' energy loss of particles in outer sc.
c 'ltrb' energy losses in single detectors (inner and outer sc.)
c 'flag' Vol. numbers and particle codes; needed to get the single
c decay spectra
c 'nhit' Number of Hits in detectors MCP, SCI and SCO
c 'posm' Generated momentum of decay positron
c 'mcpe' energy loss in MCP detector
c 'allv' use all possible NT variables
c
c--------------------------------------------------------------------------------

35
geant3/paw/fit_prec.for Normal file
View File

@ -0,0 +1,35 @@
function fit_prec(t)
real t, t0
vector tzero(1)
common/pawpar/par(6)
c t: time [ns]
c par(1): Background [1]
c par(2): normalization [1]
c par(3): asymmetry [1]
c par(4): frequency [MHz]
c par(5): phase [degree]
c par(6): muon lifetime [microsec]
c par(7): relaxation [1/microsec]
c um richtige Fitergebnisse zu erhalten, muss der Fitroutine der Zeitnullpunkt
c des Spektrums mitgeteilt werden:
t0 = tzero(1)
t = t-t0
twopi = 6.283185
fit_prec = par(1)+ ! Background
+ par(2)* ! normalization
+ exp(-t/(1000.*par(6)))* ! decay
+ ( 1.+par(3)* ! asymmetry
c + exp(-t*(par(7)/1000.))* ! relaxation
+ cos( twopi*(par(4)/1000. *t + par(5)/360.)) ! oscillation
+ )
end

383
geant3/paw/fit_prec2_ge.kumac Normal file
View File

@ -0,0 +1,383 @@
if [1].EQ.'?' .OR. [1].EQ.' ' then
mess '==============================================================================='
mess 'FIT_PREC IDLeft B fromFit toFit t0'
mess
mess 'Kumac to fit the GEANT simulated decay spectra of L and R.'
mess 'If IDleft is given, then it is assumed that'
mess ' IDright = IDleft+3,'
mess
mess 'B : magnetic field in gauss'
mess 'fromFit : lower boundary for fit'
mess 'toFit : upper boundary for fit'
mess 't0 : time zero'
mess
mess '==============================================================================='
exitm
endif
*
*===========================================================================
current = $hcdir()
* define Histo ID's
idl = [1]
idt = [idl]+2
idr = [idl]+3
idb = [idl]+4
*
* check if histograms exist in PAW memory
*
cd //pawc
if ( $hexist([idl]) .eq. 0 ) then
mess
mess Did not found Left histogram ID [idl]
mess
exitm
endif
if ( $hexist([idt]) .eq. 0 ) then
mess
mess Did not found Top histogram ID [idt]
mess
exitm
endif
if ( $hexist([idr]) .eq. 0 ) then
mess
mess Did not found Right histogram ID [idr]
mess
exitm
endif
if ( $hexist([idb]) .eq. 0 ) then
mess
mess Did not found Bottom histogram ID [idb]
mess
exitm
endif
shift
* transfer magnetic field strength:
B = [1]
* get region of histo to be used for fit:
if ($index([2],'.').NE.0) then
von = $rsigma([2])
else
von = [2]
endif
if ($index([3],'.').NE.0) then
bis = $rsigma([3])
else
bis = [3]
endif
* get rebinning factor and fit options:
t0 = [4]
if ( [t0] .eq. ' ' .or. [t0] .eq. '!' ) then
mess
mess Missing time zero !
mess
exitm
endif
if ( $vexist(tzero) ) then
v/de tzero
endif
v/cre tzero(1)
v/inp tzero(1) [t0]
* check if added histos are to be used and read in histos:
* do the fits:
freq = [B] * 0.013554 | precession frequency in MHz
freqmin = 0.8*[freq]
freqmax = 1.2*[freq]
freqstep = 0.005*[freq]
*
if ( $vexist(par) ) then
v/de par
endif
if ( $vexist(pmin) ) then
v/de pmin
endif
if ( $vexist(pmax) ) then
v/de pmax
endif
if ( $vexist(step) ) then
v/de step
endif
if ( $vexist(errpar) ) then
v/de errpar
endif
*
*
*
* BKG N0 Asym mu+ freq. Phase tau_mu lambda
* =======================================================
*
* v/cre par(7) r 0. 500. 0.25 0. 0. 2.19703 0.
* v/cre pmin(7) r -10. 0. 0. 0. -400. 0. 0.
* v/cre pmax(7) r 100. 10000. 0.60 0. +400. 3. 100.
* v/cre step(7) r 0. 10. 0.01 0. 5. 0. 0.
* v/cre errpar(7) r 1. 1. 1. 1. 1. 1. 1.
v/cre par(6) r 0. 500. 0.25 [freq] 0. 2.19703
v/cre pmin(6) r -10. 0. 0. [freqmin] -400. 0.
v/cre pmax(6) r 100. 10000. 0.40 [freqmax] +400. 3.
v/cre step(6) r 0. 10. 0.01 [freqstep] 5. 0.
v/cre errpar(6) r 0. 1. 1. 1. 1. 0.
*
tdec0 = tzero(1) | vector tzero filled in FIT_PREC.KUMAC
ttt = Fit from [von] to [bis], t0 = [tdec0]
title [ttt]
*
*
hi/fit [idl]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar
*
* store fit result and error for A_Left, frequency and relaxation rate
*
a_l = par(3)
a_lerr = errpar(3)
a_l_weight = [a_l]/([a_lerr]*[a_lerr])
inv_a_lerr = 1./([a_lerr]*[a_lerr])
f_l = par(4)
f_lerr = errpar(4)
f_l_weight = [f_l]/([f_lerr]*[f_lerr])
inv_f_lerr = 1./([f_lerr]*[f_lerr])
* r_l = par(7)
* r_lerr = errpar(7)
* r_l_weight = [r_l]/([r_lerr]*[r_lerr])
* inv_r_lerr = 1./([r_lerr]*[r_lerr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_l = par(8)
* beta_lerr = errpar(8)
* beta_l_weight = [beta_l]/([beta_lerr]*[beta_lerr])
* inv_beta_lerr = 1./([beta_lerr]*[beta_lerr])
* endif
*
*
* initialize the phase for the top detector
*
goto right
*
phase = par(5)
phase = [phase]+90.
v/inp par(5) [phase]
hi/fit [idt]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar
*
* store fit result and error for A_top
*
a_t = par(3)
a_terr = errpar(3)
a_t_weight = [a_t]/([a_terr]*[a_terr])
inv_a_terr = 1./([a_terr]*[a_terr])
f_t = par(4)
f_terr = errpar(4)
f_t_weight = [f_t]/([f_terr]*[f_terr])
inv_f_terr = 1./([f_terr]*[f_terr])
* r_t = par(7)
* r_terr = errpar(7)
* r_t_weight = [r_t]/([r_terr]*[r_terr])
* inv_r_terr = 1./([r_terr]*[r_terr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_t = par(8)
* beta_terr = errpar(8)
* beta_t_weight = [beta_t]/([beta_terr]*[beta_terr])
* inv_beta_terr = 1./([beta_terr]*[beta_terr])
* endif
*
*
*
* initialize the phase for the right detector
*
right:
phase = par(5)
phase = [phase]+180.
v/inp par(5) [phase]
hi/fit [idr]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar
* store fit result and error for A_right
*
a_r = par(3)
a_rerr = errpar(3)
a_r_weight = [a_r]/([a_rerr]*[a_rerr])
inv_a_rerr = 1./([a_rerr]*[a_rerr])
f_r = par(4)
f_rerr = errpar(4)
f_r_weight = [f_r]/([f_rerr]*[f_rerr])
inv_f_rerr = 1./([f_rerr]*[f_rerr])
* r_r = par(7)
* r_rerr = errpar(7)
* r_r_weight = [r_r]/([r_rerr]*[r_rerr])
* inv_r_rerr = 1./([r_rerr]*[r_rerr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_r = par(8)
* beta_rerr = errpar(8)
* beta_r_weight = [beta_r]/([beta_rerr]*[beta_rerr])
* inv_beta_rerr = 1./([beta_rerr]*[beta_rerr])
* endif
*
*
*
* initialize the phase for the bottom detector
*
goto end
*
phase = par(5)
phase = [phase]+90.
v/inp par(5) [phase]
hi/fit [idb]([von]:[bis]) paw$dir:fit_prec.for bl0 6 par step pmin pmax errpar
*
* store fit result and error for A_bottom
*
a_b = par(3)
a_berr = errpar(3)
a_b_weight = [a_b]/([a_berr]*[a_berr])
inv_a_berr = 1./([a_berr]*[a_berr])
f_b = par(4)
f_berr = errpar(4)
f_b_weight = [f_b]/([f_berr]*[f_berr])
inv_f_berr = 1./([f_berr]*[f_berr])
* r_b = par(7)
* r_berr = errpar(7)
* r_b_weight = [r_b]/([r_berr]*[r_berr])
* inv_r_berr = 1./([r_berr]*[r_berr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_b = par(8)
* beta_berr = errpar(8)
* beta_b_weight = [beta_b]/([beta_berr]*[beta_berr])
* inv_beta_berr = 1./([beta_berr]*[beta_berr])
* endif
*
*
* plot result:
end:
zone 2 2
hi/pl [idl] e
hi/pl [idt] e
hi/pl [idr] e
hi/pl [idb] e
*
*--------------------------------
*
* calculate the weighted mean and the "inner" and "outer" standard deviation
* fo asymmetry
*
sum1 = [a_l_weight] + [a_r_weight]
sum2 = [inv_a_lerr] + [inv_a_rerr]
a_mean = [sum1]/[sum2]
in_a_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* now calculate the outer standard deviation
*
dev_l = ([a_mean]-[a_l])*([a_mean]-[a_l])*[inv_a_lerr]
dev_r = ([a_mean]-[a_r])*([a_mean]-[a_r])*[inv_a_rerr]
sum3 = [dev_l] + [dev_r]
out_a_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*----------------------------------
*
* mean of frequency
*
sum1 = [f_l_weight] + [f_r_weight]
sum2 = [inv_f_lerr] + [inv_f_rerr]
f_mean = [sum1]/[sum2]
in_f_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
dev_l = ([f_mean]-[f_l])*([f_mean]-[f_l])*[inv_f_lerr]
dev_r = ([f_mean]-[f_r])*([f_mean]-[f_r])*[inv_f_rerr]
sum3 = [dev_l] + [dev_r]
out_f_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*--------------------------------------
*
* mean of relaxation
*
* sum1 = [r_l_weight] + [r_t_weight] + [r_r_weight] + [r_b_weight]
* sum2 = [inv_r_lerr] + [inv_r_terr] + [inv_r_rerr] + [inv_r_berr]
* r_mean = [sum1]/[sum2]
* in_r_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
* dev_l = ([r_mean]-[r_l])*([r_mean]-[r_l])*[inv_r_lerr]
* dev_t = ([r_mean]-[r_t])*([r_mean]-[r_t])*[inv_r_terr]
* dev_r = ([r_mean]-[r_r])*([r_mean]-[r_r])*[inv_r_rerr]
* dev_b = ([r_mean]-[r_b])*([r_mean]-[r_b])*[inv_r_berr]
* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
* out_r_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
* mean of stretched exponent
*
*
* if ( [opt] .eq. 'strexp' ) then
* sum1 = [beta_l_weight] + [beta_t_weight] + [beta_r_weight] + [beta_b_weight]
* sum2 = [inv_beta_lerr] + [inv_beta_terr] + [inv_beta_rerr] + [inv_beta_berr]
* beta_mean = [sum1]/[sum2]
* in_beta_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
* dev_l = ([beta_mean]-[beta_l])*([beta_mean]-[beta_l])*[inv_beta_lerr]
* dev_t = ([beta_mean]-[beta_t])*([beta_mean]-[beta_t])*[inv_beta_terr]
* dev_r = ([beta_mean]-[beta_r])*([beta_mean]-[beta_r])*[inv_beta_rerr]
* dev_b = ([beta_mean]-[beta_b])*([beta_mean]-[beta_b])*[inv_beta_berr]
* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
* out_beta_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
* endif
*
*----------------------------------------
*
* output fit parameter definitions to text window:
mess
mess 'par(1): Background [1]'
mess 'par(2): normalization [1]'
mess 'par(3): asymmetry [1]'
mess 'par(4): frequency [MHz]'
mess 'par(5): phase [degree]'
mess 'par(6): muon lifetime [microSec]'
mess 'par(7): relaxation [1/microSec]'
mess
mess
mess A = [a_mean]
mess A_err_in = [in_a_mean]
mess A_err_out = [out_a_mean]
mess
mess Freq. = [f_mean]
mess F_err_in = [in_f_mean]
mess F_err_out = [out_f_mean]
mess
* mess Relax. = [r_mean]
* mess r_err_in = [in_r_mean]
* mess r_err_out = [out_r_mean]
* mess
* if ( [opt] .eq. 'strexp' ) then
* mess
* mess Beta = [beta_mean]
* mess Beta_err_in = [in_beta_mean]
* mess Beta_err_out = [out_beta_mean]
* mess
* endif
* restore previous settings:
cd [current]
zone
title ' '
set ysiz
set xsiz

415
geant3/paw/fit_prec_ge.kumac Normal file
View File

@ -0,0 +1,415 @@
if [1].EQ.'?' .OR. [1].EQ.' ' then
mess '==============================================================================='
mess 'FIT_PREC IDLeft B fromFit toFit t0'
mess
mess 'Kumac to fit the GEANT simulated decay spectra of L,T,R,B.'
mess 'If IDleft is given, then it is assumed that'
mess ' IDtop = IDleft+2,'
mess ' IDright = IDleft+3,'
mess ' IDbottom = IDleft+4.'
mess
mess 'B : magnetic field in gauss'
mess 'fromFit : lower boundary for fit'
mess 'toFit : upper boundary for fit'
mess 't0 : time zero'
mess
mess '==============================================================================='
exitm
endif
*
*===========================================================================
current = $hcdir()
fitopt = 'bl0' | fit option b=parameters with boundaries
| l=log likelihood fit (Poisson stat.)
| 0=do not plot fit result
* define Histo ID's
idl = [1]
idt = [idl]+2
idr = [idl]+3
idb = [idl]+4
*
* check if histograms exist in PAW memory
*
cd //pawc
if ( $hexist([idl]) .eq. 0 ) then
mess
mess Did not found Left histogram ID [idl]
mess
exitm
endif
if ( $hexist([idt]) .eq. 0 ) then
mess
mess Did not found Top histogram ID [idt]
mess
exitm
endif
if ( $hexist([idr]) .eq. 0 ) then
mess
mess Did not found Right histogram ID [idr]
mess
exitm
endif
if ( $hexist([idb]) .eq. 0 ) then
mess
mess Did not found Bottom histogram ID [idb]
mess
exitm
endif
shift
* transfer magnetic field strength:
B = [1]
* get region of histo to be used for fit:
if ($index([2],'.').NE.0) then
von = $rsigma([2])
else
von = [2]
endif
if ($index([3],'.').NE.0) then
bis = $rsigma([3])
else
bis = [3]
endif
* get rebinning factor and fit options:
t0 = [4]
if ( [t0] .eq. ' ' .or. [t0] .eq. '!' ) then
mess
mess Missing time zero !
mess
exitm
endif
if ( $vexist(tzero) ) then
v/de tzero
endif
v/cre tzero(1)
v/inp tzero(1) [t0]
* check if added histos are to be used and read in histos:
* do the fits:
freq = [B] * 0.013554 | precession frequency in MHz
freqmin = 0.8*[freq]
freqmax = 1.2*[freq]
freqstep = 0.005*[freq]
*
if ( $vexist(par) ) then
v/de par
endif
if ( $vexist(pmin) ) then
v/de pmin
endif
if ( $vexist(pmax) ) then
v/de pmax
endif
if ( $vexist(step) ) then
v/de step
endif
if ( $vexist(errpar) ) then
v/de errpar
endif
*
*
*
* BKG N0 Asym mu+ freq. Phase tau_mu lambda
* =======================================================
*
* v/cre par(7) r 0. 500. 0.25 0. 0. 2.19703 0.
* v/cre pmin(7) r -10. 0. 0. 0. -400. 0. 0.
* v/cre pmax(7) r 100. 10000. 0.60 0. +400. 3. 100.
* v/cre step(7) r 0. 10. 0.01 0. 5. 0. 0.
* v/cre errpar(7) r 1. 1. 1. 1. 1. 1. 1.
v/cre par(6) r 0. 500. 0.25 [freq] 90. 2.19703
v/cre pmin(6) r -10. 0. 0. [freqmin] -400. 0.
v/cre pmax(6) r 100. 10000. 0.40 [freqmax] +400. 3.
v/cre step(6) r 0. 10. 0.01 [freqstep] 5. 0.
v/cre errpar(6) r 0. 1. 1. 1. 1. 0.
*
tdec0 = tzero(1) | vector tzero filled in FIT_PREC.KUMAC
ttt = Fit from [von] to [bis], t0 = [tdec0]
* title [ttt]
*
*
hi/fit [idl]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
* store fit result and error for A_Left, frequency and relaxation rate
*
a_l = par(3)
a_lerr = errpar(3)
a_l_weight = [a_l]/([a_lerr]*[a_lerr])
inv_a_lerr = 1./([a_lerr]*[a_lerr])
f_l = par(4)
f_lerr = errpar(4)
f_l_weight = [f_l]/([f_lerr]*[f_lerr])
inv_f_lerr = 1./([f_lerr]*[f_lerr])
* r_l = par(7)
* r_lerr = errpar(7)
* r_l_weight = [r_l]/([r_lerr]*[r_lerr])
* inv_r_lerr = 1./([r_lerr]*[r_lerr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_l = par(8)
* beta_lerr = errpar(8)
* beta_l_weight = [beta_l]/([beta_lerr]*[beta_lerr])
* inv_beta_lerr = 1./([beta_lerr]*[beta_lerr])
* endif
*
*
* initialize the phase for the top detector
*
phase = par(5)
phase = [phase]+90.
v/inp par(5) [phase]
hi/fit [idt]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
* store fit result and error for A_top
*
a_t = par(3)
a_terr = errpar(3)
a_t_weight = [a_t]/([a_terr]*[a_terr])
inv_a_terr = 1./([a_terr]*[a_terr])
f_t = par(4)
f_terr = errpar(4)
f_t_weight = [f_t]/([f_terr]*[f_terr])
inv_f_terr = 1./([f_terr]*[f_terr])
* r_t = par(7)
* r_terr = errpar(7)
* r_t_weight = [r_t]/([r_terr]*[r_terr])
* inv_r_terr = 1./([r_terr]*[r_terr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_t = par(8)
* beta_terr = errpar(8)
* beta_t_weight = [beta_t]/([beta_terr]*[beta_terr])
* inv_beta_terr = 1./([beta_terr]*[beta_terr])
* endif
*
*
*
* initialize the phase for the right detector
*
phase = par(5)
phase = [phase]+90.
v/inp par(5) [phase]
hi/fit [idr]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
* store fit result and error for A_right
*
a_r = par(3)
a_rerr = errpar(3)
a_r_weight = [a_r]/([a_rerr]*[a_rerr])
inv_a_rerr = 1./([a_rerr]*[a_rerr])
f_r = par(4)
f_rerr = errpar(4)
f_r_weight = [f_r]/([f_rerr]*[f_rerr])
inv_f_rerr = 1./([f_rerr]*[f_rerr])
* r_r = par(7)
* r_rerr = errpar(7)
* r_r_weight = [r_r]/([r_rerr]*[r_rerr])
* inv_r_rerr = 1./([r_rerr]*[r_rerr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_r = par(8)
* beta_rerr = errpar(8)
* beta_r_weight = [beta_r]/([beta_rerr]*[beta_rerr])
* inv_beta_rerr = 1./([beta_rerr]*[beta_rerr])
* endif
*
*
*
* initialize the phase for the bottom detector
*
phase = par(5)
phase = [phase]+90.
v/inp par(5) [phase]
hi/fit [idb]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
* store fit result and error for A_bottom
*
a_b = par(3)
a_berr = errpar(3)
a_b_weight = [a_b]/([a_berr]*[a_berr])
inv_a_berr = 1./([a_berr]*[a_berr])
f_b = par(4)
f_berr = errpar(4)
f_b_weight = [f_b]/([f_berr]*[f_berr])
inv_f_berr = 1./([f_berr]*[f_berr])
* r_b = par(7)
* r_berr = errpar(7)
* r_b_weight = [r_b]/([r_berr]*[r_berr])
* inv_r_berr = 1./([r_berr]*[r_berr])
*
* if ( [opt] .eq. 'strexp' ) then
* beta_b = par(8)
* beta_berr = errpar(8)
* beta_b_weight = [beta_b]/([beta_berr]*[beta_berr])
* inv_beta_berr = 1./([beta_berr]*[beta_berr])
* endif
*
*
* plot result:
pics kum
set ysiz 25
zone 2 3
hi/pl [idl] e
hi/pl [idt] e
hi/pl [idr] e
hi/pl [idb] e
*
*--------------------------------
*
* calculate the weighted mean and the "inner" and "outer" standard deviation
* fo asymmetry
*
sum1 = [a_l_weight] + [a_t_weight] + [a_r_weight] + [a_b_weight]
sum2 = [inv_a_lerr] + [inv_a_terr] + [inv_a_rerr] + [inv_a_berr]
a_mean = [sum1]/[sum2]
in_a_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* now calculate the outer standard deviation
*
dev_l = ([a_mean]-[a_l])*([a_mean]-[a_l])*[inv_a_lerr]
dev_t = ([a_mean]-[a_t])*([a_mean]-[a_t])*[inv_a_terr]
dev_r = ([a_mean]-[a_r])*([a_mean]-[a_r])*[inv_a_rerr]
dev_b = ([a_mean]-[a_b])*([a_mean]-[a_b])*[inv_a_berr]
sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
out_a_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*----------------------------------
*
* mean of frequency
*
sum1 = [f_l_weight] + [f_t_weight] + [f_r_weight] + [f_b_weight]
sum2 = [inv_f_lerr] + [inv_f_terr] + [inv_f_rerr] + [inv_f_berr]
f_mean = [sum1]/[sum2]
in_f_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
dev_l = ([f_mean]-[f_l])*([f_mean]-[f_l])*[inv_f_lerr]
dev_t = ([f_mean]-[f_t])*([f_mean]-[f_t])*[inv_f_terr]
dev_r = ([f_mean]-[f_r])*([f_mean]-[f_r])*[inv_f_rerr]
dev_b = ([f_mean]-[f_b])*([f_mean]-[f_b])*[inv_f_berr]
sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
out_f_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*--------------------------------------
*
* mean of relaxation
*
* sum1 = [r_l_weight] + [r_t_weight] + [r_r_weight] + [r_b_weight]
* sum2 = [inv_r_lerr] + [inv_r_terr] + [inv_r_rerr] + [inv_r_berr]
* r_mean = [sum1]/[sum2]
* in_r_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
* dev_l = ([r_mean]-[r_l])*([r_mean]-[r_l])*[inv_r_lerr]
* dev_t = ([r_mean]-[r_t])*([r_mean]-[r_t])*[inv_r_terr]
* dev_r = ([r_mean]-[r_r])*([r_mean]-[r_r])*[inv_r_rerr]
* dev_b = ([r_mean]-[r_b])*([r_mean]-[r_b])*[inv_r_berr]
* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
* out_r_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
* mean of stretched exponent
*
*
* if ( [opt] .eq. 'strexp' ) then
* sum1 = [beta_l_weight] + [beta_t_weight] + [beta_r_weight] + [beta_b_weight]
* sum2 = [inv_beta_lerr] + [inv_beta_terr] + [inv_beta_rerr] + [inv_beta_berr]
* beta_mean = [sum1]/[sum2]
* in_beta_mean = $sigma(sqrt(1./[sum2])) | inner standard deviation
*
* no calculate the outer standard deviation
*
* dev_l = ([beta_mean]-[beta_l])*([beta_mean]-[beta_l])*[inv_beta_lerr]
* dev_t = ([beta_mean]-[beta_t])*([beta_mean]-[beta_t])*[inv_beta_terr]
* dev_r = ([beta_mean]-[beta_r])*([beta_mean]-[beta_r])*[inv_beta_rerr]
* dev_b = ([beta_mean]-[beta_b])*([beta_mean]-[beta_b])*[inv_beta_berr]
* sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
* out_beta_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
* endif
*
*----------------------------------------
*
* output fit parameter definitions to text window:
mess
mess 'par(1): Background [1]'
mess 'par(2): normalization [1]'
mess 'par(3): asymmetry [1]'
mess 'par(4): frequency [MHz]'
mess 'par(5): phase [degree]'
mess 'par(6): muon lifetime [microSec]'
mess 'par(7): relaxation [1/microSec]'
mess
mess
mess A = [a_mean]
mess A_err_in = [in_a_mean]
mess A_err_out = [out_a_mean]
mess
mess Freq. = [f_mean]
mess F_err_in = [in_f_mean]
mess F_err_out = [out_f_mean]
mess
* mess Relax. = [r_mean]
* mess r_err_in = [in_r_mean]
* mess r_err_out = [out_r_mean]
* mess
* if ( [opt] .eq. 'strexp' ) then
* mess
* mess Beta = [beta_mean]
* mess Beta_err_in = [in_beta_mean]
* mess Beta_err_out = [out_beta_mean]
* mess
* endif
* output to graphics window (quick and dirty)
v/cre x(2) r 0 10
v/cre y(2) r 0 10
set xtic 0.0001
set ytic 0.0001
set vsiz 0.0001
zone 1 3 3 s
graph 2 x y w
*
asym = $format([a_mean],F8.4)
in_err = $format([in_a_mean],F7.4)
out_err = $format([out_a_mean],F7.4)
*
text 1 8 'A = '//[asym]//'+-'//[in_err]//'+-'//[out_err] 0.4
*
freq = $format([f_mean],F8.5)
in_err = $format([in_f_mean],F7.5)
out_err = $format([out_f_mean],F7.5)
*
app = ' MHz'
text 1 5 '[n] = '//[freq]//'+-'//[in_err]//'+-'//[out_err]//[app] 0.4
* restore previous settings:
cd [current]
zone
* title ' '
set xtic
set ytic
set vsiz
set ysiz
set xsiz

158
geant3/paw/geant_dohist.kumac Normal file
View File

@ -0,0 +1,158 @@
*
* quicky KUMAC to get decay spectra from GEANT simulations
* (GEANT_LEMSR.KUMAC)
*
* TP, 03-Feb-1999, PSI
*
* [1] Geant Runnumber; for histogram titles
*
*----------------------------------------------------------------------
*
if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then
mess
mess Missing GEANT RunNumber as parameter.
mess Needed for Histogram titles.
mess
mess Syntax:
mess
mess GEANT_DOHIST RunNo
mess
exitm
endif
runno = [1]
*
* cuts for e+ and e- and gammas
*
cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec.
cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec.
cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec.
*
cut $11 p0>30.
cut $10 desci>0.and.desco>0
cut $1 de_left_i>0.and.de_left_o>0
cut $2 de_top_i>0.and.de_top_o>0
cut $3 de_rite_i>0.and.de_rite_o>0
cut $4 de_bot_i>0.and.de_bot_o>0
*
cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out
cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out
cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out
cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out
*
1d 1 'GEANT '//[runno]//' volno, dE gt 0' 256 -.5 255.5
1d 2 'GEANT '//[runno]//' volno, dE gt 0, e+ bit or e- bit required' 256 -.5 255.5
1d 3 'GEANT '//[runno]//' volno, dE ge 0, only gammas' 256 -.5 255.5
1d 4 'GEANT '//[runno]//' volno, dE gt 0, only gammas' 256 -.5 255.5
1d 5 'GEANT '//[runno]//' volno, volno gt 0' 256 -.5 255.5
*goto fine
1d 404 'GEANT '//[runno]//' only left, de gt 0, 50ns bins' 220 0. 11000.
1d 406 'GEANT '//[runno]//' only top, de gt 0, 50ns bins' 220 0. 11000.
1d 407 'GEANT '//[runno]//' only right, de gt 0, 50ns bins' 220 0. 11000.
1d 408 'GEANT '//[runno]//' only bottom, de gt 0, 50ns bins' 220 0. 11000.
1d 504 'GEANT '//[runno]//' left , de gt 0, 50ns bins' 220 0. 11000.
1d 506 'GEANT '//[runno]//' top , de gt 0, 50ns bins' 220 0. 11000.
1d 507 'GEANT '//[runno]//' right, de gt 0, 50ns bins' 220 0. 11000.
1d 508 'GEANT '//[runno]//' bottom, de gt 0, 50ns bins' 220 0. 11000.
1d 604 'GEANT '//[runno]//' left , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 606 'GEANT '//[runno]//' top , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 607 'GEANT '//[runno]//' right , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 608 'GEANT '//[runno]//' bottom, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 704 'GEANT '//[runno]//' left , only gamma, de ge 0, 50ns bins' 220 0. 11000.
1d 706 'GEANT '//[runno]//' top , only gamma, de ge 0, 50ns bins' 220 0. 11000.
1d 707 'GEANT '//[runno]//' right , only gamma, de ge 0, 50ns bins' 220 0. 11000.
1d 708 'GEANT '//[runno]//' bottom, only gamma, de ge 0, 50ns bins' 220 0. 11000.
* goto cont
*
fine:
1d 1404 'GEANT '//[runno]//' only left, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1406 'GEANT '//[runno]//' only top, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1407 'GEANT '//[runno]//' only right, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1408 'GEANT '//[runno]//' only bottom, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1504 'GEANT '//[runno]//' left , de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1506 'GEANT '//[runno]//' top , de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1507 'GEANT '//[runno]//' right, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1508 'GEANT '//[runno]//' bottom, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1604 'GEANT '//[runno]//' left , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1606 'GEANT '//[runno]//' top , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1607 'GEANT '//[runno]//' right , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1608 'GEANT '//[runno]//' bottom, e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
*
mess Plot 1=Volno dE>0
nt/pl 111.volno volno>0.and.$10 -1
mess Plot 2=Volno dE>0 and e+ or e- bit reuqired
nt/pl 111.volno volno>0.and.$10.and.($21.or.$22) -2
mess Plot 3=Volno dE>=0, only gammas
nt/pl 111.volno volno>0.and.partcode=17 -3
mess Plot 4=Volno dE>0, only gammas
nt/pl 111.volno volno>0.and.partcode=17.and.$10 -4
mess Plot 5=Volno volno>0
nt/pl 111.volno volno>0 -5
mess NTPlot left=404
nt/pl 111.tsci volno=17.and.$10 -404
mess NTPlot top =406
nt/pl 111.tsci volno=34.and.$10 -406
mess NTPlot right=407
nt/pl 111.tsci volno=68.and.$10 -407
mess NTPlot bottom=408
nt/pl 111.tsci volno=136.and.$10 -408
mess NTPlot left=1404
nt/pl 111.tsci volno=17.and.$10 -1404
mess NTPlot top =1406
nt/pl 111.tsci volno=34.and.$10 -1406
mess NTPlot right=1407
nt/pl 111.tsci volno=68.and.$10 -1407
mess NTPlot bottom=1408
nt/pl 111.tsci volno=136.and.$10 -1408
mess NTPlot left id=504
nt/pl 111.tsci $1.and.$5 -504
mess NTPlot top id=506
nt/pl 111.tsci $2.and.$6 -506
mess NTPlot rite id=507
nt/pl 111.tsci $3.and.$7 -507
mess NTPlot bot id=508
nt/pl 111.tsci $4.and.$8 -508
mess NTPlot left id=604
nt/pl 111.tsci $1.and.$5.and.($21.or.$22) -604
mess NTPlot top id=606
nt/pl 111.tsci $2.and.$6.and.($21.or.$22) -606
mess NTPlot rite id=607
nt/pl 111.tsci $3.and.$7.and.($21.or.$22) -607
mess NTPlot bot id=608
nt/pl 111.tsci $4.and.$8.and.($21.or.$22) -608
mess NTPlot left id=704
nt/pl 111.tsci partcode=17.and.$5 -704
mess NTPlot top id=706
nt/pl 111.tsci partcode=17.and.$6 -706
mess NTPlot rite id=707
nt/pl 111.tsci partcode=17.and.$7 -707
mess NTPlot bot id=708
nt/pl 111.tsci partcode=17.and.$8 -708
mess NTPlot left id=1504
nt/pl 111.tsci $1.and.$5 -1504
mess NTPlot top id=1506
nt/pl 111.tsci $2.and.$6 -1506
mess NTPlot rite id=1507
nt/pl 111.tsci $3.and.$7 -1507
mess NTPlot bot id=1508
nt/pl 111.tsci $4.and.$8 -1508
mess NTPlot left id=1604
nt/pl 111.tsci $1.and.$5.and.($21.or.$22) -1604
mess NTPlot top id=1606
nt/pl 111.tsci $2.and.$6.and.($21.or.$22) -1606
mess NTPlot rite id=1607
nt/pl 111.tsci $3.and.$7.and.($21.or.$22) -1607
mess NTPlot bot id=1608
nt/pl 111.tsci $4.and.$8.and.($21.or.$22) -1608

58
geant3/paw/geant_dohist1.kumac Normal file
View File

@ -0,0 +1,58 @@
*
* quicky KUMAC to get decay spectra from GEANT simulations
* (GEANT_LEMSR.KUMAC)
*
* TP, 03-Feb-1999, PSI
*
* [1] Geant Runnumber; for histogram titles
*
*
* decay spectra, only 1 detector has hit, e+ e- bit required
*----------------------------------------------------------------------
*
if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then
mess
mess Missing GEANT RunNumber as parameter.
mess Needed for Histogram titles.
mess
mess Syntax:
mess
mess GEANT_DOHIST RunNo
mess
exitm
endif
runno = [1]
*
* cuts for e+ and e- and gammas
*
cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec.
cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec.
cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec.
*
cut $11 p0>30.
cut $10 desci>0.and.desco>0
cut $1 de_left_i>0.and.de_left_o>0
cut $2 de_top_i>0.and.de_top_o>0
cut $3 de_rite_i>0.and.de_rite_o>0
cut $4 de_bot_i>0.and.de_bot_o>0
*
cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out
cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out
cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out
cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out
*
1d 614 'GEANT '//[runno]//' only left , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 616 'GEANT '//[runno]//' only top , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 617 'GEANT '//[runno]//' only right , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 618 'GEANT '//[runno]//' only bottom, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
mess NTPlot left id=614
nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and.volno.eq.17 -614
mess NTPlot top id=616
nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and.volno.eq.34 -616
mess NTPlot rite id=617
nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and.volno.eq.68 -617
mess NTPlot bot id=618
nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and.volno.eq.136 -618

72
geant3/paw/geant_dohist2.kumac Normal file
View File

@ -0,0 +1,72 @@
*
* quicky KUMAC to get decay spectra from GEANT simulations
* (GEANT_LEMSR.KUMAC)
*
* TP, 03-Feb-1999, PSI
*
* [1] Geant Runnumber; for histogram titles
*
*
* decay spectra, only 1 detector has hit, e+ e- bit required
*----------------------------------------------------------------------
*
if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then
mess
mess Missing GEANT RunNumber as parameter.
mess Needed for Histogram titles.
mess
mess Syntax:
mess
mess GEANT_DOHIST RunNo
mess
exitm
endif
runno = [1]
*
* cuts for e+ and e- and gammas
*
cut $20 btest(int(partcode),0).and.btest(int(partcode),4) |require gamma in both detec.
cut $21 btest(int(partcode),1).and.btest(int(partcode),5) |require e+ in both detec.
cut $22 btest(int(partcode),2).and.btest(int(partcode),6) |require e- in both detec.
*
cut $11 p0>30.
cut $10 desci>0.and.desco>0
cut $1 de_left_i>0.and.de_left_o>0
cut $2 de_top_i>0.and.de_top_o>0
cut $3 de_rite_i>0.and.de_rite_o>0
cut $4 de_bot_i>0.and.de_bot_o>0
*
cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out
cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out
cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out
cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out
*
1d 624 'GEANT '//[runno]//' left, dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 626 'GEANT '//[runno]//' top , dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 627 'GEANT '//[runno]//' right, dE others=0 , e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 628 'GEANT '//[runno]//' bottom, dE others=0, e+ or e- bit, de gt 0, 50ns bins' 220 0. 11000.
1d 1624 'GEANT '//[runno]//' left, dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1626 'GEANT '//[runno]//' top , dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1627 'GEANT '//[runno]//' right, dE others=0 , e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
1d 1628 'GEANT '//[runno]//' bottom, dE others=0, e+ or e- bit, de gt 0, 2ns bins' 5500 0.5 11000.5
mess NTPlot left id=624
nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and..not.($2.or.$3.or.$4) -624
mess NTPlot top id=626
nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and..not.($1.or.$3.or.$4) -626
mess NTPlot rite id=627
nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and..not.($1.or.$2.or.$4) -627
mess NTPlot bot id=628
nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and..not.($1.or.$2.or.$3) -628
mess NTPlot left id=1624
nt/pl 111.tsci $1.and.$5.and.($21.or.$22).and..not.($2.or.$3.or.$4) -1624
mess NTPlot top id=1626
nt/pl 111.tsci $2.and.$6.and.($21.or.$22).and..not.($1.or.$3.or.$4) -1626
mess NTPlot rite id=1627
nt/pl 111.tsci $3.and.$7.and.($21.or.$22).and..not.($1.or.$2.or.$4) -1627
mess NTPlot bot id=1628
nt/pl 111.tsci $4.and.$8.and.($21.or.$22).and..not.($1.or.$2.or.$3) -1628

92
geant3/paw/geant_eloss.kumac Normal file
View File

@ -0,0 +1,92 @@
*
* quicky KUMAC to get energy loss spectra from GEANT simulations
* (GEANT_ELOSS.KUMAC)
*
* TP, 03-Feb-1999, PSI
*
* [1] Geant Runnumber; for histogram titles
*
*----------------------------------------------------------------------
*
if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!') then
mess
mess Missing GEANT RunNumber as parameter.
mess Needed for Histogram titles.
mess
mess Syntax:
mess
mess GEANT_DOHIST RunNo
mess
exitm
endif
runno = [1]
*
* cuts for e+ and e- and gammas
*
cut $20 btest(int(partcode),0).and.btest(int(partcode),4) | bits for gamma's
cut $21 btest(int(partcode),1).and.btest(int(partcode),5) | bits for e+
cut $22 btest(int(partcode),2).and.btest(int(partcode),6) | bits for e-
* cut $30 partcode=17 | only gammas
*
* cut $11 p0>30.
* cut $10 desci>0.and.desco>0
* cut $1 de_left_i>0.and.de_left_o>0
* cut $2 de_top_i>0.and.de_top_o>0
* cut $3 de_rite_i>0.and.de_rite_o>0
* cut $4 de_bot_i>0.and.de_bot_o>0
*
* cut $5 btest(int(volno),0).and.btest(int(volno),4) | bits for left in and out
* cut $6 btest(int(volno),1).and.btest(int(volno),5) | bits for top in and out
* cut $7 btest(int(volno),2).and.btest(int(volno),6) | bits for right in and out
* cut $8 btest(int(volno),3).and.btest(int(volno),7) | bits for bottom in and out
*
1d 100 'GEANT '//[runno]//' desci gt 0' 200 0. 10.
1d 101 'GEANT '//[runno]//' descie+ gt 0' 200 0. 10.
1d 102 'GEANT '//[runno]//' descie- gt 0' 200 0. 10.
1d 103 'GEANT '//[runno]//' desciga gt 0' 200 0. 1.
1d 110 'GEANT '//[runno]//' desci gt 0, require e+ or e- bits' 200 0. 10.
1d 200 'GEANT '//[runno]//' desco gt 0' 200 0. 10.
1d 201 'GEANT '//[runno]//' descoe+ gt 0' 200 0. 10.
1d 202 'GEANT '//[runno]//' descoe- gt 0' 200 0. 10.
1d 203 'GEANT '//[runno]//' descoga gt 0' 200 0. 1.
1d 210 'GEANT '//[runno]//' desco gt 0, require e+ or e- bits' 200 0. 10.
*
* inner detectors
*
mess Plot 100=Desci gt 0
nt/pl 111.desci desci>0 -100
mess Plot 101=Descie+ gt 0
nt/pl 111.descipos descipos>0 -101
mess Plot 102=Descie- gt 0
nt/pl 111.desciele desciele>0 -102
mess Plot 103=Descigamma gt 0
nt/pl 111.descigam descigam>0 -103
mess Plot 110=Desci gt 0, e+ or e- bits required
nt/pl 111.desci desci>0.and.($21.or.$22) -110
*
* outer detectors
*
mess Plot 200=Desco gt 0
nt/pl 111.desco desco>0 -200
mess Plot 201=Descoe+ gt 0
nt/pl 111.descopos descopos>0 -201
mess Plot 202=Descoe- gt 0
nt/pl 111.descoele descoele>0 -202
mess Plot 203=Descogamma gt 0
nt/pl 111.descogam descogam>0 -203
mess Plot 210=Desco gt 0, e+ or e- bits required
nt/pl 111.desco desco>0.and.($21.or.$22) -210

196
geant3/paw/geant_part.kumac Normal file
View File

@ -0,0 +1,196 @@
*
* paw$dir:geant_part.kumac
*
* KUMAC to copy first 400 ns of GEANT decay spectra to
* histogram.
*
* [1] ID of left histogram; other ID's will be reconstructed
*
* TP, 23-feb-1999
*
*----------------------------------------------------------------------------
*
if ( [1] .eq. ' ' .or. [1] .eq. '?' .or. [1] .eq. '!' ) then
mess
mess GEANT_PART IDleft
mess
mess Copies first 400ns of existing GEANT decay spectra
mess to new IDs; NewID = OldID+1000
mess
exitm
endif
idl = [1]
idt = [idl]+2
idr = [idt]+1
idb = [idr]+1
idln = [idl]+1000
idtn = [idt]+1000
idrn = [idr]+1000
idbn = [idb]+1000
*
upp = 400. | upper bin center of new histograms
*
if ( $vexist(data) ) then
v/de data
endif
*
* get histo information
*
if ( $hexist([idl]) .eq. 0 ) then
mess
mess Left Histogram ID = [idl] does not exist !
mess
exitm
endif
if ( $hexist([idt]) .eq. 0 ) then
mess
mess Top Histogram ID = [idt] does not exist !
mess
exitm
endif
if ( $hexist([idr]) .eq. 0 ) then
mess
mess Right Histogram ID = [idr] does not exist !
mess
exitm
endif
if ( $hexist([idb]) .eq. 0 ) then
mess
mess Bottom Histogram ID = [idb] does not exist !
mess
exitm
endif
*
*----------------------------------------------------------------------------
*
* L E F T
*
xmin = $hinfo([idl], 'XMIN')
xmax = $hinfo([idl], 'XMAX')
nbin = $hinfo([idl], 'XBINS')
bin = ([xmax] - [xmin])/[nbin]
*
mess Left ID = [idl] with bin size [bin]
*
* copy histogram to vector
*
v/cre data([nbin]
get_vec/con [idl] data
*
* construct new upper limit
*
nbin = [upp]/[bin] + 1
xmax = [xmin] + [nbin]*[bin]
mess New LeftId = [idln] with xmax = [xmax] and nbin = [nbin]
*
* book histogram
*
if ( $hexist([idln]) ) then
hi/de [idln]
endif
1d [idln] 'part of left, id = '//[idl] [nbin] [xmin] [xmax]
put_vec/con [idln] data(1:[nbin])
v/de data
*
*---------------------------------------------------------------------------
*
* T O P
*
xmin = $hinfo([idt], 'XMIN')
xmax = $hinfo([idt], 'XMAX')
nbin = $hinfo([idt], 'XBINS')
bin = ([xmax] - [xmin])/[nbin]
*
mess Top ID = [idt] with bin size [bin]
*
* copy histogram to vector
*
v/cre data([nbin]
get_vec/con [idt] data
*
* construct new upper limit
*
nbin = [upp]/[bin] + 1
xmax = [xmin] + [nbin]*[bin]
mess New TopId = [idtn] with xmax = [xmax] and nbin = [nbin]
*
* book histogram
*
if ( $hexist([idtn]) ) then
hi/de [idtn]
endif
1d [idtn] 'part of top, id = '//[idt] [nbin] [xmin] [xmax]
put_vec/con [idtn] data(1:[nbin])
v/de data
*
*---------------------------------------------------------------------------
*
* Right
*
xmin = $hinfo([idr], 'XMIN')
xmax = $hinfo([idr], 'XMAX')
nbin = $hinfo([idr], 'XBINS')
bin = ([xmax] - [xmin])/[nbin]
*
mess Right ID = [idr] with bin size [bin]
*
* copy histogram to vector
*
v/cre data([nbin]
get_vec/con [idr] data
*
* construct new upper limit
*
nbin = [upp]/[bin] + 1
xmax = [xmin] + [nbin]*[bin]
mess New RightId = [idrn] with xmax = [xmax] and nbin = [nbin]
*
* book histogram
*
if ( $hexist([idrn]) ) then
hi/de [idrn]
endif
1d [idrn] 'part of Right, id = '//[idr] [nbin] [xmin] [xmax]
put_vec/con [idrn] data(1:[nbin])
v/de data
*
*------------------------------------------------------------------------------
*
* Bottom
*
xmin = $hinfo([idb], 'XMIN')
xmax = $hinfo([idb], 'XMAX')
nbin = $hinfo([idb], 'XBINS')
bin = ([xmax] - [xmin])/[nbin]
*
mess Bottom ID = [idb] with bin size [bin]
*
* copy histogram to vector
*
v/cre data([nbin]
get_vec/con [idb] data
*
* construct new upper limit
*
nbin = [upp]/[bin] + 1
xmax = [xmin] + [nbin]*[bin]
mess New BottomId = [idbn] with xmax = [xmax] and nbin = [nbin]
*
* book histogram
*
if ( $hexist([idbn]) ) then
hi/de [idbn]
endif
1d [idbn] 'part of Bottom, id = '//[idb] [nbin] [xmin] [xmax]
put_vec/con [idbn] data(1:[nbin])
v/de data
*
*------------------------------------------------------------------------------
*
zone 2 2
hi/pl [idln]
hi/pl [idtn]
hi/pl [idrn]
hi/pl [idbn]
zone
*

62
geant3/paw/readge.kumac Normal file
View File

@ -0,0 +1,62 @@
*
* READGE.KUMAC
*
* KUMAC to open NT file from GEANT_LEMSR simulation on next free
* LUN.
*
* Input parameter:
*
* [1]: GEant Runnumber
*
* TP, 10-feb-1999, PSI
*
* TP, 06-Apr-2001, PSI: Unix, Linux and Windows version
*
*---------------------------------------------------------------------------
*
if ([1] .eq. '?' .or. [1] .eq. ' ' .or. [1] .eq. '!') then
mess
mess READGE kumac to open GEANT_LEMSR NTUPLE file.
mess
mess Syntax:
mess
mess READGE runNr
mess
exitm
endif
runnr = [1]
*
* determine operating system to build the filename;
* only Windows or Unix/Linux
*
os = $OS
*
if ( [os] .eq. 'Windows') then
file = 'c:\users\thomas\mc\geant\data\geant_lemsr_'//[runnr]//'.nt'
goto check_file
endif
if ( [os] .eq. 'UNIX') then
* if ( $OSTYPE .eq. 'Linux' ) then
file = /scratch/$USER/geant/geant_lemsr_[runnr].nt
goto check_file
* else
* file = $HOME//'/mc/geant/data/geant_lemsr_'//[runnr]//'.nt'
* goto check_file
* endif
endif
*
mess Operating system [os] not supported
exitm
*
check_file:
if ( $fexist([file]) ) then
hi/file 0 [file] 4096
else
mess
mess File [file] does not exist...
mess
endif
*
*

50
geant3/src/lemsr/Makefile Normal file
View File

@ -0,0 +1,50 @@
#
# $HOME/mc/geant/src/lemsr/Makefile
#
# makefile to compile and link the program geant_lemsr to
# simulate decaying muons at the LE-muSR spectrometer at PSI
# NTuple file to be analysed by PAW.
#
# TP, 15-Sep-2000, PSI
# 09-Apr-2001
#
#--------------------------------------------------------------------------
#
# The following lines contain specific switches for different
# Unix systems. Outcomment the unwanted.
#
#--------------------------------------------------------------------------
#
# For Digital Unix
#
#FF = f77
#
#--------------------------------------------------------------------------
#
# For Linux
#
FF = g77
CFLAGS = -c -I$(CERN_ROOT)/include/ -Wno-globals -finit-local-zero
#CFLAGS += -fno-automatic # not really necessary
#
#--------------------------------------------------------------------------
#
OBJ = geant_lemsr_main.o get_direction.o gudcay.o gukine.o \
guout.o gustep.o gutrack.o gutreve.o michel.o rotate.o ugeom.o
.F.o:
$(FF) $(CFLAGS) $*.F
geant_lemsr: $(OBJ)
$(FF) -o $(HOME)/mc/geant/bin/geant_lemsr $(OBJ) \
`cernlib geant321 pawlib graflib packlib mathlib`
# `cernlib geant321 pawlib graflib/Motif packlib mathlib`
#
# the program cernlib searches the cernlib libraries as well as
# needed system libraries
#
clean:
rm -f *.o *~ \#* .#*
#--------------------------------------------------------------------------
#
# EOF
#

124
geant3/src/lemsr/common.inc Normal file
View File

@ -0,0 +1,124 @@
c
c ~/mc/geant/lemsr/src/common.inc
c
c GEANT_LEMSR program,
c include common blocks other than NTuple variables
c
c TP, 05-Apr-2001, PSI
c
c 17-Dec-2002 TP extended bfield array from 4 to 5,
c bfield(5)=polarisation in sample
c
c-----------------------------------------------------------------------
c
c check OS
c
#if defined( _WIN32 )
#define OS_WIN
#else
#define OS_UNIX
#endif
c
c parameter
c
character*(*) IN_DIR, OUT_DIR
#if defined(OS_WIN)
parameter IN_DIR = 'c:\users\thomas\mc\geant\inp\'
parameter OUT_DIR = 'c:\users\thomas\mc\geant\data\'
#else
parameter (IN_DIR = '/afs/psi.ch/user/p/prokscha/mc/geant/inp/')
c parameter (OUT_DIR = '/afs/psi.ch/user/p/prokscha/mc/geant/data/')
parameter (OUT_DIR = '/scratch/prokscha/geant/')
#endif
c
character*(*) INFILE, PARFILE
parameter (INFILE = IN_DIR//'geant_lemsr_input.dat')
parameter (PARFILE = IN_DIR//'geant_lemsr.input')
c
c
c INFILE: read input data from file instead of throwing them
c in the program
c
c PARFILE: read input parameter from this file
c
real*4 MU_MASS, TWO_PI, POS_MASS, SQR_POS_MASS,
1 MAX_TOTAL_ENERGY, TAU_MUON, MU_GYRO, R_MCP
parameter (MU_MASS = .105658389 ) !muon mass, GeV/c^2
parameter (POS_MASS = 0.51099906 ) ! e+ mass, MeV/c^2
parameter (SQR_POS_MASS = 0.26112 ) ! squared positron mass
parameter (MAX_TOTAL_ENERGY = 52.8304302) ! MeV
c
parameter (TAU_MUON = 2.19703E-6 ) ! muon life time in s
parameter (MU_GYRO = 13554.) ! Hz/G
c
parameter (R_MCP = 1.00 ) ! radius of MCP in cm
parameter (TWO_PI = 6.283185307)
c
c-------------------------------------------------------------------
c
real*4 bfield(5) ! new data card with B-field components
c ! bfield(4) = polarization (%)
c ! bfield(5) = polarization (%) in sample
integer*4 ix1 !random generator seed
c
c count successfull and rejected MICHEL-energy throws
c
integer*4 cou, rejcou
c
c------------------------------------------------------------------------------
c
c flags for geometrical setup; it is possible to mount or remove the
c
c 'mcp2' MCP2 plates
c 'mcpa' MCP2 anode
c 'mcps' MCP2 stainless steel vacuum tube
c 'samp' sample (this is at the moment the sample holder)
c 'cryo' material of the cryostat without shielding
c 'crsh' cryo shield
c 'efla' the 100 CF flange at the end of the sample/MCP2 tube
c
c the volumes are activated when printing
c
c geom 'mcp2' 'efla' for example when data cards are requested
c
logical*1 l_mcp2
logical*1 l_mcpa
logical*1 l_mcps
logical*1 l_samp
logical*1 l_cryo
logical*1 l_crsh
logical*1 l_efla
logical*1 l_run11
character*4 samp_mat
c
c------------------------------------------------------------------------
c
c flags for secondary particles
c
integer*4 lspar(10) ! new data card defined by CALL FFKEY to
c ! select seoncdary particles
c ! ' e+', ' e-' and ' gam' possible.
c
logical*1 l_pos
logical*1 l_ele
logical*1 l_gam
c
c-------------------------------------------------------------------------
c
c flags for output of percentage of processed events
c
integer*4 percent(9)
c
c-------------------------------------------------------------------------
c
c--- the common blocks
c
common /magnetic_field/ bfield
common /seeds/ ix1
common /counter/ cou, rejcou
common /geo_flags/ l_mcp2, l_mcpa, l_mcps, l_samp, l_cryo,
1 l_crsh, l_efla, l_run11, samp_mat
common /part_flag/ l_pos, l_ele, l_gam
common /proc_event/ percent
c-----------------------------------------------------------------------

45
geant3/src/lemsr/cwn.inc Normal file
View File

@ -0,0 +1,45 @@
c
c cwn.inc
c
c include variable declaration for cwn Ntuples
c
c TP, 15-Sep-2000, PSI
c
c-----------------------------------------------------------------------------
c
c CWN variables
c
real*4 nt_x0, nt_y0, nt_z0, nt_px0, nt_py0, nt_pz0, nt_p0
integer*4 nt_ipart
real*4 nt_xe, nt_ye, nt_ze, nt_eposend
real*4 nt_tmcp, nt_tsci ! times are integer, use 1ns bins
real*4 nt_xsci, nt_ysci, nt_zsci
real*4 nt_desci, nt_descipos, nt_descigam, nt_desciele
real*4 nt_xsco, nt_ysco, nt_zsco
real*4 nt_desco, nt_descopos, nt_descogam, nt_descoele
integer*4 nt_volno, nt_partcode
integer*4 nt_mcpnhits, nt_scinhits, nt_sconhits
real*4 nt_pxpos, nt_pypos, nt_pzpos
real*4 nt_demcp
real*4 nt_descil, nt_descit, nt_descir, nt_descib,
1 nt_descol, nt_descot, nt_descor, nt_descob
c
c the common blocks
c
common /de_ltrb/ nt_descil, nt_descit, nt_descir, nt_descib,
1 nt_descol, nt_descot, nt_descor, nt_descob
common /init_par/ nt_x0, nt_y0, nt_z0, nt_px0, nt_py0, nt_pz0,
1 nt_p0
common /partcode/ nt_ipart
common /end_pos/ nt_xe, nt_ye, nt_ze, nt_eposend
common /times/ nt_tmcp, nt_tsci
common /sci_pos/ nt_xsci, nt_ysci, nt_zsci
common /sci_elos/ nt_desci, nt_descipos, nt_descigam, nt_desciele
common /sco_pos/ nt_xsco, nt_ysco, nt_zsco
common /sco_elos/ nt_desco , nt_descopos, nt_descogam, nt_descoele
common /flags/ nt_volno, nt_partcode
common /nhits/ nt_mcpnhits, nt_scinhits, nt_sconhits
common /pos_mom/ nt_pxpos, nt_pypos, nt_pzpos
common /de_mcp/ nt_demcp
c
c-----------------------------------------------------------------------------

816
geant3/src/lemsr/geant_lemsr_main.F Normal file
View File

@ -0,0 +1,816 @@
#include "geant321/pilot.h"
c-----------------------------------------------
c ~/mc/geant/src/lemsr/geant_lemsr_main.F
c-----------------------------------------------
c
program geant_lemsr
c
c Program very much "grown", so many things could be programmed in
c a better way...
c
c program to simulate the response of the low energy muon muSR
c detector setup at PSI. Some parts are taken from Nai_simul.for.
c It is intended to steer the simulation by format free inputs
c (routine gffgo); this can be done interactively or data sets
c can be read from file.
c
c Uses Run11 setup geometry.
c
c For precession: at the moment use the total TOF for the calculation
c of the precession angle. When the muon is started as particle code=500
c and upstream of the target, the precession is calculated as if the
c muon is flying all the time through the B-field.
c
c TP 18-Jan-1999 PSI
c
c TP, 4-Feb-1999: CW Ntuples
c TP, 9-Feb-1999: store energy loss per detector in NTuple
c TP,16-Feb-1999: time in NTUPLE as REAL*4 instead of INTEGER
c TP,26-Feb-1999: If particles start upstream gaussian divergence of
c beam now possible.
c
c*****************************************************************************
c
c TP, 11-May-2000 PSI
c Windows NT 4.0 version; uses include files which
c I created manually from the VMS GEANT321.TLB text library.
c
c TP, 15-Sep-2000 PSI
c Unix version Used fsplit to separate program into
c subroutines and functions.
c
c
c TP, 05-Apr-2001 PSI
c Unix, NT, and Linux (to be checked)
c version; use the Compaq Fortran fpp
c precompiler on NT
c
c 17-Dec-2002 TP extend bfield array to account for different
c polarization in sample (=bfield(5))
c gaussian muon momentum distribution possible
c via sets(8) and sets(9):
c sigma_p = p0 * (sets(8)/100.+sets(9)/1000.)
c
c------------------------------------------------------------------
c
implicit none
c
c include some files from $CERN_ROOT/include/
c
c to avoid explicitly including the gt***.inc file
c define CERNLIB_TYPE; the gt***.inc files are then
c automatically included in the gc***.inc files.
c
#define CERNLIB_TYPE
#include "geant321/gclist.inc"
#include "geant321/gcflag.inc"
#include "geant321/gcphys.inc"
#include "geant321/gccuts.inc"
#include "geant321/gckine.inc"
c
#include "cwn.inc"
#include "common.inc"
c
c
c the particle definition (variable IPART): see GEANT manual CONS300-1
c
c particle IPART
c
c gamma 1
c positron 2
c electron 3
c neutrino 4
c muon+ 5
c muon- 6
c pion0 7
c pion+ 8
c pion- 9
c ...
c
integer*4 NG, NH
parameter (NG = 2000000)
parameter (NH = 1000000)
integer*4 ipawc,q
integer*4 iwktyp
integer*4 iquest(100)
c
common /quest/iquest !needed to get enough space for NTUPLE
common /pawc/ipawc(NH)
common /gcbank/q(NG)
c
c parameter for drawing (used in GEANT routine GDRAW)
c
real*4 gdtheta,gdphi,gdpsi
real*4 gdx0,gdy0,gdsx,gdsy
c
real*4 xn,yn,scal,r(2)
integer*4 i
c
c for opening NTUPLE file and NTUPLE
c
integer*4 istat
character*72 ntfile
c
c the run number copied from GEAN variable IDRUN
c
character*4 runno
c
c for graphic output to PostScript file
c
character*72 psfile
c
integer*4 lunin /4/
integer*4 ipres
real*4 ub(1)
c
c-----------------------------------------------------------------
c
c flags for NTUPLE variable definitions
c
c 'init' write initial parameter to NT
c 'code' initial particle ID to NT
c 'endp' end position and energy of e+
c 'time' t-mcp and t-sci
c 'scip' position when particle hits inner Sc.
c 'scie' energy loss of particles in inner sc.
c 'scop' position when particle hits outer Sc.
c 'scoe' energy loss of particles in outer sc.
c 'flag' Vol. numbers and particle codes; needed to get the single
c decay spectra
c 'nhit' Number of Hits in detectors MCP, SCI and SCO
c 'posm' Generated momentum of decay positron
c 'mcpe' energy loss in MCP detector
c 'allv' use all possible NT variables
c 'ltrb' energy loss in single detectors (inner and outer separate)
c
integer*4 ntvar(20) ! new data card defined by CALL FFKEY to
c ! select NTuple variables
c
logical*1 l_init_par /.false./
logical*1 l_code /.false./
logical*1 l_endp /.false./
logical*1 l_time /.false./
logical*1 l_scip /.false./
logical*1 l_scie /.false./
logical*1 l_scop /.false./
logical*1 l_scoe /.false./
logical*1 l_flag /.false./
logical*1 l_nhit /.false./
logical*1 l_posm /.false./
logical*1 l_mcpe /.false./
logical*1 l_allv /.false./
logical*1 l_ltrb /.false./
c
c------------------------------------------------------------------
c
c start of program code
c
c------------------------------------------------------------------
c
c general initializations
c
c initial random generator seeds
c
ix1 = 343156821
c
do i = 1, 9
percent(i) = 0
enddo
c
l_mcp2 = .false.
l_mcpa = .false.
l_mcps = .false.
l_samp = .false.
l_cryo = .false.
l_crsh = .false.
l_efla = .false.
l_run11 = .false.
samp_mat = ' '
c
l_pos = .false.
l_ele = .false.
l_gam = .false.
c
c-----------------------------------------
c
call gzebra(NG)
call hlimit(-NH)
call ginit ! initialize GEANT variables
call gzinit ! initialize GEANT data structure
call gpart ! GEANT standard particles
call gspart(500, 'muon+', 5, 0.1056584, 1., 2.19703E-6, ub, 1)
c
c------------------------------------------------------------------
c
c open parfile in any case, the user can decide if data should
c be read from parfile or entered by terminal
c
open(lunin, file=PARFILE, status='OLD',err=1003)
write(6,*) ' '
write(6,*) '%GEANT-LEMSR-I- Successfully opened data card file '
write(6,*)
write(6,*) PARFILE
c
c read data cards
c
write(6,*) ' '
write(6,*) 'Type in the data cards now... finish with ''STOP'''
write(6,*) '--------------------------------------------------'
write(6,*) ' '
write(6,*) 'READ 4 : reads input data from '//
1 'LUN=4--> file '
write(6,*) ' '
write(6,*) PARFILE
write(6,*) ' '
c write(6,*) 'RUNG 24 : user run number'
c write(6,*) ' '
c write(6,*) 'PLOT 0 : no graphical output'
c write(6,*) 'PLOT 1 : output on workstation only'
c write(6,*) 'PLOT 2 : output on workstation and metafile'
c write(6,*) ' '
c write(6,*) 'SETS 2 0 : particle is e+(=2), michel '//
c 1 'distributed'
c write(6,*) ' start MCP2 detector'
c write(6,*) 'SETS 2 -74 27 5 1 : beam positron with 27.51 MeV/c '//
c 1 'at z=-74cm'
c write(6,*) 'SETS 5 -74 2 0 0 : beam muon(=5) with 2.00 MeV/c'
c write(6,*) ' particle started 74cm upstream of'
c write(6,*) ' of center of MCP2 tube.'
c write(6,*) ' '
c write(6,*) 'TRIG 100 : 100 events to be generated'
c write(6,*) ' '
c write(6,*) 'GEOM ''mcp2'' ''samp'' : mcp2 and cryo are mounted.'
c write(6,*) ' '
c write(6,*) 'SPAR '' e+'' '' e-'' '' gam'' : select sec. '//
c 1 'particles.'
c write(6,*) ' '
c write(6,*) 'BFIE 0. 0. 100. 90. : Sample B-field (Gauss), '//
c 1 'Bx, By, Bz,'
c write(6,*) ' and Polarization = 90 % here.'
c write(6,*)
c write(6,*) 'NTVA ''init'' ''time'' ''flag'' : define NT '//
c 1 'variables.'
c write(6,*) 'See geant_lemsr.input for more information.'
c write(6,*)
c write(6,*) 'View 1 Thet Phi Psi u0 v0 10*scal_x 10*scal_y : '//
c 1 'for graphical output'
c write(6,*) 'Def. 60 40 0 10 10 2 2'
c write(6,*) ' '
c
c---------------------------------------------------------------
c
c define new data card for particles
c
call ffkey('SPAR', lspar, 10, 'M')
c
c define new data card for magnetic field components
c
call ffkey('BFIE', bfield, 5, 'R')
c
c define new data card for NTUPLE variables
c
call ffkey('NTVA', ntvar, 20, 'M')
c
c read data card
c
call gffgo
c
close(lunin)
c
write(6,*) ' '
write(6,*) '%GEANT-LEMSR-I- Closed file ', PARFILE
write(6,*)
c
write(6,*) 'B-field components Bx, By, Bz = ', bfield(1),
1 bfield(2), bfield(3)
write(6,*) ' '
write(6,*) 'Initial muon polarization = ', bfield(4),' %.'
write(6,*)
c
if ( bfield(1) .ne. 0. ) then
write(6,*) '******************************************'
write(6,*) ' '
write(6,*) 'GEANT-LEMSR-E- B_x not yet implemented...'
write(6,*) 'Sorry, - - - - > E X I T'
write(6,*) ' '
write(6,*) '******************************************'
call exit
endif
if ( bfield(2) .ne. 0. .and. bfield(3) .ne. 0. ) then
write(6,*) '******************************************'
write(6,*) ' '
write(6,*) 'GEANT-LEMSR-E- B_y >0 and B_z >0 not yet '//
1 'implemented...'
write(6,*) 'Sorry, - - - - > E X I T'
write(6,*) ' '
write(6,*) '******************************************'
call exit
endif
c
if ( lsets(2) .ne. 0 .and. lsets(7) .ne. 1) then
write(6,*) ' '
write(6,*) 'particles of type ',lsets(1),' start at z0 = ',
1 lsets(2)
write(6,*) ' '
else if ( lsets(7) .eq. 1 ) then
write(6,*) ' '
write(6,*) 'particles of type ',lsets(1),' start at z0 '//
1 'read from file'
write(6,*) ' '
else
write(6,*) ' '
write(6,*) 'particles of type ',lsets(1),' start at z0 = '//
1 '1.40002 cm.'
write(6,*) ' '
endif
c
c-------------------------------------
c
c look for particles to be created during tracking
c
ipres = 0
call glook(' e+',lspar,10,ipres)
if ( ipres .ne. 0 ) then
l_pos = .true.
write(6,*) ' '
write(6,*) 'secondary particle: e+...'
endif
call glook(' e-',lspar,10,ipres)
if ( ipres .ne. 0 ) then
l_ele = .true.
write(6,*) 'secondary particle: e-...'
endif
call glook(' gam',lspar,10,ipres)
if ( ipres .ne. 0 ) then
l_gam = .true.
write(6,*) 'secondary particle: gamma...'
endif
c
c------------------------------------------------------------
c
c now look for geometrical setup variables
c
call glook('mcp2',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_mcp2 = .true.
write(6,*) 'Setup with: MCP2'
endif
c
call glook('mcpa',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_mcpa = .true.
write(6,*) 'Setup with: MCP2 anode and support ring'
endif
c
call glook('ru11',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_run11 = .true.
if ( l_mcp2 .or. l_mcpa ) then
write(6,*) 'Setup with: MCP2 Run 11 geometry, delay '//
1 'line anode'
endif
else
if ( l_mcp2 .or. l_mcpa ) then
write(6,*) 'Setup with: MCP2 Run 10 geometry, WSZ anode'
endif
endif
c
call glook('mcps',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_mcps = .true.
write(6,*) 'Setup with: MCP2 stainless steel vacuum tube'
endif
c
call glook('samp',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_samp = .true.
write(6,*) 'Setup with: sample holder'
endif
c
call glook('saal',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
samp_mat = 'saal'
if ( l_samp ) then
write(6,*) 'Setup with: Al - sample holder'
endif
else
if ( l_samp ) then
write(6,*) 'Setup with: Cu - sample holder'
endif
endif
call glook('cryo',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_cryo = .true.
write(6,*) 'Setup with: cryostat'
endif
c
call glook('crsh',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_crsh = .true.
write(6,*) 'Setup with: He shield'
endif
c
call glook('efla',lgeom,20,ipres)
if ( ipres .ne. 0 ) then
l_efla = .true.
write(6,*) 'Setup with: 100 CF flange at the end'
endif
c
c----------------------------------------------------
c
c look for NT variables
c
call glook('init', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_init_par = .true.
write(6,*) 'NT variable: initial position and momenta.'
endif
c
call glook('code', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_code = .true.
write(6,*) 'NT variable: Code of primary particle.'
endif
c
call glook('endp', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_endp = .true.
write(6,*) 'NT variable: End position and end energy of e+.'
endif
c
call glook('time', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_time = .true.
write(6,*) 'NT variable: t-MCP2 and t-SCI.'
endif
call glook('scip', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_scip = .true.
write(6,*) 'NT variable: Position of particle when '//
1 'hitting inner Sc..'
endif
c
call glook('scie', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_scie = .true.
write(6,*) 'NT variable: Energy deposited in inner Sc.'
endif
c
call glook('scop', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_scop = .true.
write(6,*) 'NT variable: Position of particle when '//
1 'hitting outer Sc..'
endif
c
call glook('scoe', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_scoe = .true.
write(6,*) 'NT variable: Energy deposited in outer Sc.'
endif
c
call glook('flag', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_flag = .true.
write(6,*) 'NT variable: Volume number and particle codes.'
write(6,*) 'Note: the Vol. No. is necessary to get '//
1 'single decay spectra.'
endif
c
call glook('nhit', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_nhit = .true.
write(6,*) 'NT variable: Hits per Track in MCP, Sci Sco.'
endif
c
call glook('posm', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_posm = .true.
write(6,*) 'NT variable: inital e+ momentum.'
endif
c
call glook('mcpe', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_mcpe = .true.
write(6,*) 'NT variable: Energy deposited in MCP2.'
endif
c
call glook('allv', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_allv = .true.
write(6,*) 'NT variable: all variables.'
endif
c
call glook('ltrb', ntvar,20,ipres)
if ( ipres .ne. 0 ) then
l_ltrb = .true.
write(6,*) 'NT variable: energies deposited in single '//
1 'scintill.'
endif
c
c----------------------------------------------------
c
write(runno,'(I4)',err=1002) idrun ! copy run number to string
c
if ( lsets(2) .ne. 0 .and. lsets(7) .ne. 1) then
if ( lsets(3) .gt. 0 ) then
pkine(1) = float(lsets(3)) + float(lsets(4))/10. +
1 float(lsets(5))/100.
else
pkine(1) = float(lsets(3)) - float(lsets(4))/10. -
1 float(lsets(5))/100.
endif
write(6,'(a,1x,f10.4)') ' %GEANT-LEMSR-I- beam momentum '//
1 'set to ', pkine(1)
endif
if ( lsets(6) .ne. 0 ) then
write(6,*) ' '
write(6,*) 'Gaussian beam divergence set to (degree) ',
1 lsets(6)
write(6,*) ' '
endif
c
c define volumes and detectors
c
call ugeom
c
c compute the cross-sections for the different tracking media
c
call gphysi
c
c make MARS invisible
c
call gsatt('MARS','SEEN',0)
c
call ggclos ! end of the geometry, must be called
c
c print some information about volumes and tracking media on screen
c
c call gprint('VOLU',0)
c call gptmed(0)
c call gpmate(17)
c call gpmate(18)
c
c-----------------------------------------------------------------------------
c
c from GEXAM.FOR example :
c
c the card PLOT is used in the following way:
c no PLOT card --> LPLOT(1)=0 --> no graphic output
c PLOT 1 --> LPLOT(1)=1 --> plotting onto workstation
c PLOT 2 --> LPLOT(1)=2 --> plottinig onto workstation
c and PostScript file
if (lplot(1) .gt. 0) then
c
if ( lview(1) .eq. 1 ) then
gdtheta = lview(2)
gdphi = lview(3)
gdpsi = lview(4)
gdx0 = lview(5)
gdy0 = lview(6)
gdsx = lview(7)/10.
gdsy = lview(8)/10.
else
gdtheta = 60.
gdphi = 40.
gdpsi = 0.
gdx0 = 10.
gdy0 = 10.
gdsx = 0.2
gdsy = 0.2
endif
c
c
call iginit(0) ! init HIGZ
c
c ask for worktsation type
c this defines the geometry of the X11 window, according to
c HIGZ_WINDOWS.DAT
c
write(6,'(''0Enter the workstation type : '')')
read(5,*) iwktyp
call hplint(iwktyp)
call igset('2BUF',11.)
c
if (lplot(1) .gt. 1) then
c
c open graphics output file
c
psfile = OUT_DIR//'geant_lemsr_'//runno//'.ps'
open(unit=10, file=psfile, status='NEW',
1 form='FORMATTED',
2 access='SEQUENTIAL', err=1001)
write(6,*) '%GEANT-LEMSR-I- succesfully opened '//
1 'metafile '
write(6,*)
write(6,*) psfile
c
c open and activate workstation for metafile output
c workstation ID = 5
c goes to file LUN = 10
c type = -4112 = PostScript, Landscape, A4 (--> HIGZ manual)
c
call iopwk(5, 10, -4112)
call iacwk(5)
c
c find out and select the appropriate transformation and viewport
c
call igqwk(5,'MXDS',r)
xn = min( 1., r(1)/r(2) )
yn = min( 1., r(2)/r(1) )
call isvp( 5, 0., xn, 0., yn)
scal = 20. / min( xn, yn )
call iswn( 5, 0., scal, 0., scal)
call iselnt(5)
endif ! if LPLOT(1) .GT. 1
c
c initialize the GEANT's drawing package before the first plot
c
call gdinit
c
c draw the chamber
c
call gdraw('MARS',gdtheta,gdphi,gdpsi,gdx0,gdy0,gdsx,gdsy)
call gdaxis(-10.,15.,-20.,8.)
call gdscal(gdx0, gdy0)
endif ! IF LPLOT(1) .GT. 0
c
c------------------------------------------------------------------------------
c
c open NTUPLE file and book NTUPLE
c
ntfile = OUT_DIR//'geant_lemsr_'//runno//'.nt'
c
iquest(10) = 65000 ! to allocate enough space for NTUPLE
call hropen(80,'GEANT',ntfile,'NQ',4096,istat)
if (istat.eq.0) then
write(6,*) '%GEANT-LEMSR-I- succesfully opened NTUPLE file '
write(6,*) ' '
write(6,*) ntfile
else
write(6,'('' %GEANT-LEMSR-E- error opening NTUPLE file '',/,a)')
1 ntfile
call exit
endif
c
c book CWN
c
call hbnt(111,'GEANT_LEMSR', 'D') ! 'D' disk resident NTUPLE
c
c now define the variable blocks for CWN
c
if ( l_init_par .or. l_allv ) then
call hbname(111,'init_par', nt_x0,
1 'x0:r*4, y0:r*4, z0:r*4,'//
2 'px0:r*4, py0:r*4, pz0:r*4, p0:r*4')
endif
c
if ( l_code .or. l_allv ) then
call hbname(111,'partcode', nt_ipart, 'Ipart:u*4:9')
endif
if ( l_endp .or. l_allv ) then
call hbname(111,'endpos' , nt_xe,
1 'xe:r*4, ye:r*4, ze:r*4, epos:r*4')
endif
if ( l_time .or. l_allv ) then
call hbname(111,'times' , nt_tmcp, 'tmcp:r*4, tsci:r*4')
endif
if ( l_scip .or. l_allv ) then
call hbname(111,'sci_pos' , nt_xsci,
1 'xsci:r*4, ysci:r*4, zsci:r*4')
endif
if ( l_scie .or. l_allv ) then
call hbname(111,'sci_elos', nt_desci,
1 'desci:r*4, descipos:r*4, descigam:r*4, desciele:r*4')
endif
c
if ( l_scop .or. l_allv ) then
call hbname(111,'sco_pos' , nt_xsco,
1 'xsco:r*4, ysco:r*4, zsco:r*4')
endif
if ( l_scoe .or. l_allv ) then
call hbname(111,'sco_elos', nt_desco,
1 'desco:r*4, descopos:r*4, descogam:r*4, descoele:r*4')
endif
c
if ( l_flag .or. l_allv ) then
call hbname(111,'flags' , nt_volno, 'volno:u*4:8,'//
1 'partcode:u*4:7')
endif
if ( l_nhit .or. l_allv ) then
call hbname(111,'nhits' , nt_mcpnhits,
1 'mcpnhits:u*4:8, scinhits:u*4:8, sconhits:u*4:8')
endif
if ( l_posm .or. l_allv ) then
call hbname(111,'pos_mom' , nt_pxpos,
1 'pxpos:r*4, pypos:r*4, pzpos:r*4')
endif
if ( l_mcpe .or. l_allv ) then
call hbname(111,'de_mcp' , nt_demcp, 'demcp:r*4')
endif
c
if ( l_ltrb .or. l_allv ) then
call hbname(111,'de_ltrb' , nt_descil, 'de_left_i:r*4,'//
1 'de_top_i:r*4, de_rite_i:r*4, de_bot_i:r*4,'//
2 'de_left_o:r*4,'//
3 'de_top_o:r*4, de_rite_o:r*4, de_bot_o:r*4')
endif
c
if ( lsets(7) .eq. 1 ) then ! read initial data from file
c
write(6,*) ' '
write(6,*) ' lsets(7) = 1 ---> I will read data from file ',
1 INFILE
open(70, file=INFILE, form='UNFORMATTED', status='OLD',
1 err=2000)
write(6,*) ' '
write(6,*) '%GEANT_LEMSR-I- Successfully opened file ', INFILE
read(70) nevent
write(6,*) ' '
write(6,*) ' Number of events in file = ', nevent
write(6,*) ' '
endif
c
c----------------------------------------------------------------- c
c
c E N D OF I N I T I A L I S A T I O N
c c
c-----------------------------------------------------------------
c
c now do the tracking
c
call grun
c
c update all workstations
c
if ( lplot(1) .gt. 0 ) call iuwk(0,1)
c
c accept some input to have time to look at the picture
c
if ( lplot(1) .gt. 0 ) then
write(6,*) ' '
write(6,*) ' type anything to continue '
read(5,*)
endif
if ( lplot(1) .gt. 1 ) then
call iclwk(5) ! advance paper
close(10)
endif
c
c end GEANT and HPLOT
c
call glast
c
if ( lplot(1) .gt. 0 ) call igend
c
if ( lsets(7) .eq. 1) then
close(70)
write(6,*) ' '
write(6,*) '%GEANT-LEMSR-I- closed input file ', INFILE
write(6,*) ' '
endif
c
c write data and close NTUPLE file
c
call hrout(111,i,' ')
call hrend('GEANT')
close(80)
write(6,*) '%GEANT-LEMSR-I- closed NTUPLE file ', ntfile
write(6,*) '%GEANT-LEMSR-I- results of NEUMANN rejection for '//
1 'MICHEL spectrum : '
write(6,'(t10,'' rejected trials : '',i10)') rejcou
write(6,'(t10,'' succes. trials : '',i10)') cou
write(6,*) ' '
write(6,'(t5,'' RNDM generator seed at the end : '')')
write(6,'(t10,'' ix1 = '',i15)') ix1
write(6,*) ' '
call exit
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c error messages
c
1000 continue
write(6,'(''0%GEANT-LEMSR-E- error reading GDRAW parameter...'')')
call exit
c
1001 continue
write(6,'(''0%GEANT-LEMSR-E- error opening '',a)') psfile
call exit
c
1002 continue
write(6,'(''0%GEANT-LEMSR-E- error write RunNo to string...'')')
write(6,*) ' may be, you entered a run number with more than '
write(6,*) ' 4 digits...'
call exit
c
1003 continue
write(6,'(''0%GEANT-LEMSR-E- error opening input file '',a)')
1 PARFILE
call exit
c
2000 continue
write(6,*) '%GEANT-LEMSR-E- error opening input file ', INFILE
call exit
end

77
geant3/src/lemsr/get_direction.F Normal file
View File

@ -0,0 +1,77 @@
c
c----------------------------------------------------------------------------
c
subroutine get_direction(costheta, sintheta, cosphi, sinphi,
1 pos_tot)
c
c
c T. Prokscha, 15-Sep-2000 PSI
c
c Unix Version
c
c------------------------------------------------------------------------------
c
c generates the initial direction of the decay-e+
c
implicit none
c
#define CERNLIB_TYPE
#include "geant321/gctrak.inc"
#include "common.inc"
c
real*4 pos_tot ! total e+ energy including mass
real*4 p ! polarization == bfield(4)
real*4 costheta,phi ! decay direction in rotated system
real*4 sintheta
real*4 cosphi, sinphi
real*4 eps ! pos_tot/MAX_TOTAL_ENERGY
real*4 Dw ! asymmetry factor
real*4 random(3)
real*4 x,y,z,r
c
c-----------------------------------------------------------------------------
c
x = vect(1)
y = vect(2)
z = vect(3)
r = sqrt(x**2 + y**2)
if (1.4 .le. z .and. z .le. 1.41 .and. r .le. 3.5) then
p = bfield(5) / 100.
else
p = bfield(4) / 100.
endif
c
#if defined (OS_UNIX)
call ranlux(random, 3) ! random generator from Mathlib
#else
random(1) = ran(ix1)
random(2) = ran(ix1)
random(3) = ran(ix1)
#endif
if (p.eq.0.0) then
c
costheta = 1. - 2.*random(1)
c
else
c
c see TP Logbook 7, p. 167 for getting the equations
c
eps = pos_tot / MAX_TOTAL_ENERGY
Dw = (2. * eps - 1.) / (3. - 2. * eps) * p
costheta = 1/Dw*(-1. +
1 sqrt( 1. - 2*Dw*( 2.*random(2) - 1. ) + Dw**2))
endif
c
sintheta = sqrt(1. - costheta**2)
c
phi=random(3) * TWO_PI
c
cosphi = cos(phi)
sinphi = sin(phi)
c
c now we have the decay direction in the rest system of the muon
c return.
c
return
end

110
geant3/src/lemsr/gudcay.F Normal file
View File

@ -0,0 +1,110 @@
c
#include "geant321/pilot.h"
c---------------------------------------------------------------------------
c
subroutine gudcay
c
c in case of muon started with partcode = 500, throw michel spectrum
c and anisotropic angular distribution
c
c T. Prokscha, 15-Sep-2000 PSI
c
c Unix version
c
c TP, 05-Apr-2001 PSI Unix, NT, Linux
c
c---------------------------------------------------------------------------
c
implicit none
c
#define CERNLIB_TYPE
#include "geant321/gctrak.inc"
#include "geant321/gckine.inc"
#include "geant321/gcking.inc"
c
#include "common.inc"
c
c Michel spectrum
c
real*4 pos_tot ! total positron energy (including rest mass)
real*4 p0 ! total positron momentum (MeV)
real*4 px0, py0, pz0 ! components of p0
c
c decay direction in muon rest system (Spin parallel to x-axis !!!)
c
real*4 costheta, sintheta, cosphi, sinphi
c
c 4-momentum of muon particle (if decay in flight); if decay in flight
c do a Lorentz-transformation to get momenta in real space
c
real*4 decay_rest(4), beta(4), decay_lab(4), gamma, beta_tot
real*4 decay_rot(3) ! copy of momentum comp. of decay_lab(4)
c
c----------------------------------------------------------------------------
c
c get the Michel spectrum
c
call michel(pos_tot, p0)
c
call get_direction(costheta, sintheta, cosphi, sinphi, pos_tot)
c
c get the single momentum components in the laboratory system
c the muon spin is directed along the x-axis which points to the
c LEFT detector.)
c
px0 = p0 * costheta
py0 = p0 * sintheta * cosphi
pz0 = p0 * sintheta * sinphi
c
decay_rest(1) = px0/1000. ! GEANT uses GeV/c instead of MeV/c
decay_rest(2) = py0/1000.
decay_rest(3) = pz0/1000.
decay_rest(4) = pos_tot/1000.
c
c for the Lorentz transformation, maybe not important, because the
c mu+ are slow
c
gamma = getot / MU_MASS
beta_tot = vect(7) / getot ! this is p/E, vect(7) = p
beta(1) = -vect(4) * beta_tot ! vect(4) = px/p
beta(2) = -vect(5) * beta_tot ! vect(5) = py/p
beta(3) = -vect(6) * beta_tot ! vect(6) = pz/p
c
c negative sign because we need the velocity of
c the lab-system with respect to the rest-system
c
beta(4) = gamma
c
call gloren(beta, decay_rest, decay_lab)
c
c if there is a B-field present we have to rotate the momentum vector
c at the moment only Bz implemented; exit program if Bx or By are
c given.
c
c
if ( bfield(1) .ne. 0. .or. bfield(2) .ne. 0. or. bfield(3) .ne.
1 0.) then
decay_rot(1) = decay_lab(1)
decay_rot(2) = decay_lab(2)
decay_rot(3) = decay_lab(3)
c
call rotate(decay_rot, tofg)
c
decay_lab(1) = decay_rot(1)
decay_lab(2) = decay_rot(2)
decay_lab(3) = decay_rot(3)
endif
c
ngkine = 1
tofd(1) = 0.
gkin(1,1) = decay_lab(1)
gkin(2,1) = decay_lab(2)
gkin(3,1) = decay_lab(3)
gkin(4,1) = decay_lab(4)
gkin(5,1) = 2
gpos(1,1) = vect(1)
gpos(2,1) = vect(2)
gpos(3,1) = vect(3)
c
return
end

229
geant3/src/lemsr/gukine.F Normal file
View File

@ -0,0 +1,229 @@
*
#include "geant321/pilot.h"
*
c------------------------------------------------------------------------------
c
subroutine gukine
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version
c
c TP, 05-Apr-2001 PSI: Unix, NT, Linux
c
c------------------------------------------------------------------------------
c
implicit none
c
#define CERNLIB_TYPE
#include "geant321/gckine.inc"
#include "geant321/gclist.inc"
#include "geant321/gctrak.inc"
c
#include "cwn.inc"
#include "common.inc"
c
c-----------------------------------------------------------------------------
c
c local variables
c
real*4 decay_rot(3) ! rotation of muon spin,
c ! rotated e+-momenta
real*4 x0, y0, z0, px0, py0, pz0, p0, pos_tot
real*4 radius, phi0, cosphi, sinphi, costheta, sintheta
real*4 theta2d, theta_beam, sintheta_beam, costheta_beam
real*4 proj, phi_beam, sinphi_beam, cosphi_beam, sigma_beam
integer*4 ntbeam,nttarg,nubuf,i
parameter (nubuf = 1)
real*4 ubuf(nubuf)
real*4 random_1(3), random_2(7)
real*4 sigma_p
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ipart = lsets(1)
c
c-----------------------------------------------------------------------------
c
c throw decay point at MCP, isotropic distribution and MICHEL-spectrum
c
if ( lsets(2) .eq. 0 .and. lsets(7) .ne. 1) then
c
#if defined (OS_UNIX)
call ranlux(random_1,3)
#else
do i = 1, 3
random_1(i) = ran(ix1)
enddo
#endif
c
c------------------
c
c the decay point at MCP or sample
c
c z0 = +1.6 ! distance of MCP center to center of MARS
c z0 = +1.48002 ! the Hamamatsu MCP's end is at z=1.48cm
c
c set MCP2 to 1.52cm; doing this we have the same initial z-position
c for the e+ with MCP mounting and sample mounting:
c
c Cu ring of sample (z/2 = 0.2cm) at z = 1.6 cm
c MCP Hamamatsu Z-stack (z/2 = 0.12cm) at z = 1.52 cm
c
z0 = +1.40002 ! the Hamamatsu, Galileo MCP's end and sample end is at z=1.40cm now
c ! z=1.40002 corresponds to an implatation depth of the
c ! mu+ of 200 nm.
c
radius = R_MCP * sqrt(random_1(1)) ! uniform distribution on MCP
c radius = 4.7
phi0 = TWO_PI * random_1(2)
x0 = radius * cos(phi0)
y0 = radius * sin(phi0)
c
c---------------------------
c
c the decay time
c
tofg = -TAU_MUON*alog(1.-random_1(3))
c
c the michel spectrum
c
call michel(pos_tot, p0)
c
c anisotropic angular distribution, P=100%
c
call get_direction(costheta, sintheta, cosphi, sinphi,
1 pos_tot)
c
c get the single momentum components in the laboratory system
c (the muon spin is directed along the x-axis which points to the
c LEFT detector.)
c
if ( bfield(1) .ne. 0. or. bfield(2) .ne. 0. or. bfield(3)
1 .ne. 0.) then
c
c rotate when there is a B-field present; the subroutine
c ROTATE(...) will return without any action if bfield(4)=pol.=0.
c
decay_rot(1) = p0 * costheta
decay_rot(2) = p0 * sintheta * cosphi
decay_rot(3) = p0 * sintheta * sinphi
c
call rotate(decay_rot, tofg)
c
px0 = decay_rot(1)
py0 = decay_rot(2)
pz0 = decay_rot(3)
c
else
c
px0 = p0 * costheta
py0 = p0 * sintheta * cosphi
pz0 = p0 * sintheta * sinphi
c
endif
c
c-------------------------
c
c read input data from file
c
else if ( lsets(7) .eq. 1 ) then
c
read(70) x0, y0, z0, px0, py0, pz0
z0 = z0-0.3
c
c-------------------------
c
c for beam particles hitting MCP with momentum defined by
c lsets(3), lsets(4), lsets(5)
c
else
#if defined (OS_UNIX)
call ranlux(random_2,7)
#else
do i = 1, 5
random_2(i) = ran(ix1)
enddo
#endif
c
c z0 = -74. ! start in vacuum upstream from MCP
z0 = float(lsets(2)) ! start in vacuum at z=z0;
radius = R_MCP * sqrt(random_2(1)) ! uniform distribution on MCP
phi0 = TWO_PI * random_2(2)
x0 = radius * cos(phi0)
y0 = radius * sin(phi0)+0.25
c
c if sets(8) > 0 throw gaussian momentum distribution
c
p0 = pkine(1)
if ( lsets(8) .gt. 0 ) then
sigma_p = p0 * (lsets(8)/100. + lsets(9)/1000.)
theta2d = sigma_p *
1 sqrt( -2. * alog(1.-random_2(6)))
proj = sin( TWO_PI * random_2(7))
p0 = p0 + theta2d*proj
endif
c
c throw THETA gaussian distributed to get some divergence of
c the beam
c
sigma_beam = float(lsets(6)) ! in degree
c
if ( sigma_beam .eq. 0. ) then
px0 = 0.
py0 = 0.
pz0 = p0
else
sigma_beam = sigma_beam / 57.29577951 ! in rad
theta2d = sigma_beam *
1 sqrt( -2. * alog(1.-random_2(3)))
proj = sin( TWO_PI * random_2(4))
theta_beam = proj * theta2d
sintheta_beam= sin(theta_beam)
costheta_beam= cos(theta_beam)
c
phi_beam = TWO_PI * random_2(5)
sinphi_beam = sin(phi_beam)
cosphi_beam = cos(phi_beam)
c
px0 = p0 * sintheta_beam * cosphi_beam
py0 = p0 * sintheta_beam * sinphi_beam
pz0 = p0 * costheta_beam
c
endif
c
endif
c
c-----------------------------------------------------------------------------
c
c type*,' ivert = ',ivert
c
vert(1) = x0
vert(2) = y0
vert(3) = z0
pvert(1) = px0/1000. ! GEANT uses GeV/c
pvert(2) = py0/1000.
pvert(3) = pz0/1000.
nt_x0 = x0
nt_y0 = y0
nt_z0 = z0
nt_px0 = px0
nt_py0 = py0
nt_pz0 = pz0
nt_p0 = p0
nt_ipart = ipart
c
c write(6,*) x0,y0,z0,px0,py0,pz0
c
c store "vertex" (this means in this case the initial position and
c momentum of the paricle) data to the GEANT data structures
c JVERTX and JKINE. It has to be done !
c
call gsvert(vert,ntbeam,nttarg,ubuf,nubuf,ivert)
call gskine(pvert,ipart,ivert,ubuf,nubuf,itra)
c call gpvert(ivert)
c call gpkine(itra)
c
return
end

439
geant3/src/lemsr/guout.F Normal file
View File

@ -0,0 +1,439 @@
#include "geant321/pilot.h"
c
c-----------------------------------------------------------------------------
c
subroutine guout
c
c this routine writes the data at the end of one event to NTUPLE.
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix Version, part of the lemsr_geant program
c
c TP, 05-Apr-2001, PSI Unix, NT, Linux version
c
c-----------------------------------------------------------------------------
c
implicit none
c
#define CERNLIB_TYPE
#include "geant321/gcsets.inc"
#include "geant321/gcflag.inc"
#include "geant321/gclist.inc"
c
c ntuple variables and other common variablea
c
#include "cwn.inc"
#include "common.inc"
c
c-------------------------------------------------------------------------
c
integer*4 i, j, NVDIM, NHMAX, NHDIM
c
parameter (NVDIM = 1)
parameter (NHMAX = 100)
parameter (NHDIM = 7)
integer*4 numvs(NVDIM)
integer*4 itra(NHMAX)
integer*4 nbv(NVDIM,NHMAX)
real*4 hits(NHDIM,NHMAX)
integer*4 nhits, maxind
c
real*4 detot, descipos, descigam, desciele ! energy loss total,
c ! and in Sc. due
c ! to e+, gammas, e-
c
real*4 descil, descit, descir, descib ! energy loss in one detec.
c
integer*4 bit
integer*4 volno ! volume number; this is to decide which
c ! scintillator was hit:
c ! volno = 1 (LEFT) ----> write 1 to NTuple variable
c ! = 2 (TOP) ----> write 2
c ! = 3 (RIGHT) ----> write 4
c ! = 4 (BOTTOM)----> write 8
c ---> sum variable entry if more than one detector
c was hit.
c ! for outer scintillators, start LEFT ---> 16
c and so on.
c
integer*4 partcode ! store in NTuple which particles had contributed
c ! partcode = 1 (Gamma) ---> write 1 to NT variable
c ! partcode = 2 (e+) ---> write 2 to NT
c ! partcode = 3 (e-) ---> write 4 to NT
c ---> sum NT variable values if two or three
c different particles were generated.
c for outer scintillators start with
c GAMMA ---> 8
c
real*4 ubuf(1)
integer*4 nubuf /1/
integer*4 nvert, part
real*4 vert(3), pvert(3) ! e+ vertex values
c
c-----------------------------------------------------------------------------
c
if ( lsets(1) .eq. 5 .or. lsets(1) .eq. 500) then
c ! if initial particle was a muon,
c ! then, look for e+ momentum which
c ! is created as 2. vertex
c ! "2" is therefore the 1. argument
c ! of gfkine
c
call gfkine(2, vert, pvert, part, nvert, ubuf, nubuf)
c
nt_pxpos = 1000. * pvert(1) ! MeV/c
nt_pypos = 1000. * pvert(2) ! MeV/c
nt_pzpos = 1000. * pvert(3) ! MeV/c
c
endif
c
c now look for detector hits to get information on energy loss,
c particle and which detector has a hit, and fill NTUPLE variuables
c
c------------------------------
c
c the MCP2 detector
c
do i = 1, NHDIM ! reset first HITS(,) array
do j = 1, NHMAX
hits(i, j) = 0.
enddo
enddo
detot = 0.
c
call gfhits('MDET','DMCP',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,
1 hits, nhits)
c
if ( nhits .gt. 0 ) then ! ~MCP2 Nhits>0
if ( nhits .le. NHMAX ) then
maxind = nhits
else
maxind = NHMAX
endif
do j = 1, maxind
c
detot = detot + hits(6,j) ! sum the energy loss per step
c
c write(6,'(1x,7E10.3)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
c + hits(5,j),hits(6,j),hits(7,j)
enddo
c
c now fill MCP relevant data to NTUPLE variables
c
nt_tmcp = int(hits(4,1)*1.e9) ! decay time at MCP in ns
nt_demcp = 1000. * detot ! total energy loss in MCP in MeV
c
endif ! ~ endif MCP2 hits >0
c
nt_mcpnhits = nhits
c
c-------------------------------------
c
c the inner scintillators
c
do i = 1, NHDIM ! reset first HITS(,) array
do j = 1, NHMAX
hits(i, j) = 0.
enddo
enddo
c
nhits = 0
detot = 0.
descipos = 0.
descigam = 0.
desciele = 0.
volno = 0
partcode = 0
descil = 0.
descit = 0.
descir = 0.
descib = 0.
c
call gfhits('ISSC','SCIS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits,
1 nhits)
c
c write(6,*) 'Number of Hits in inner Sc.: ',nhits
c
if ( nhits .gt. 0 ) then ! ~ hit in inner Sc.
c
if ( nhits .le. NHMAX ) then
maxind = nhits
else
maxind = NHMAX
endif
do j = 1, maxind
detot = detot + hits(6,j) ! sum the energy loss per step
if ( hits(7,j) .eq. 1. ) then ! particle is gamma
descigam = descigam + hits(6,j)
endif
if ( hits(7,j) .eq. 2. ) then ! particle is positron
descipos = descipos + hits(6,j)
endif
if ( hits(7,j) .eq. 3. ) then ! particle is electron
desciele = desciele + hits(6,j)
endif
c
c now look for volume
c
bit = nbv(1,j) - 1
volno = ibset(volno, bit)
c
c sum up the energy deposited in one detector
c
if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector
descil = descil + hits(6,j)
endif
if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector
descit = descit + hits(6,j)
endif
if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector
descir = descir + hits(6,j)
endif
if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector
descib = descib + hits(6,j)
endif
c
c look for particles
c
bit = int(hits(7,j)) - 1
partcode = ibset(partcode, bit)
c
c write(6,'(1x,7E10.3,2I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno
c
enddo
c
nt_tsci = hits(4,1)*1.e9 ! decay time at Sci. in ns
nt_xsci = hits(1,1)
nt_ysci = hits(2,1)
nt_zsci = hits(3,1)
nt_desci = detot * 1000.
nt_descipos = descipos * 1000.
nt_descigam = descigam * 1000.
nt_desciele = desciele * 1000.
nt_descil = descil * 1000.
nt_descit = descit * 1000.
nt_descir = descir * 1000.
nt_descib = descib * 1000.
nt_volno = volno
nt_partcode = partcode
c
endif ! ~ hit in inner Sc.
c
nt_scinhits = nhits
c
c end of inner sc.
c
c-------------------------------
c
c the outer scintillators
c
do i = 1, NHDIM ! reset first HITS(,) array
do j = 1, NHMAX
hits(i, j) = 0.
enddo
enddo
c
nhits = 0
detot = 0.
descipos = 0.
descigam = 0.
desciele = 0.
descil = 0.
descit = 0.
descir = 0.
descib = 0.
c
c
call gfhits('OSSC','SCOS',NVDIM,NHDIM,NHMAX,0,numvs,itra,nbv,hits,
1 nhits)
c
c write(6,*) 'Number of Hits in outer Sc.: ',nhits
c
if ( nhits .gt. 0 ) then ! ~ hit in outer Sc.
c
if ( nhits .le. NHMAX ) then
maxind = nhits
else
maxind = NHMAX
endif
do j = 1, maxind
detot = detot + hits(6,j) ! sum the energy loss per step
if ( hits(7,j) .eq. 1. ) then ! particle is gamma
descigam = descigam + hits(6,j)
endif
if ( hits(7,j) .eq. 2. ) then ! particle is positron
descipos = descipos + hits(6,j)
endif
if ( hits(7,j) .eq. 3. ) then ! particle is electron
desciele = desciele + hits(6,j)
endif
c
c now look for volume
c
bit = nbv(1,j) - 1 + 4
volno = ibset(volno, bit)
c
c sum up the energy deposited in one detector
c
if ( nbv(1, j) .eq. 1 ) then ! hit is in left detector
descil = descil + hits(6,j)
endif
if ( nbv(1, j) .eq. 2 ) then ! hit is in top detector
descit = descit + hits(6,j)
endif
if ( nbv(1, j) .eq. 3 ) then ! hit is in right detector
descir = descir + hits(6,j)
endif
if ( nbv(1, j) .eq. 4 ) then ! hit is in bottom detector
descib = descib + hits(6,j)
endif
c
c look for particles
c
bit = int(hits(7,j)) - 1 + 4
partcode = ibset(partcode, bit)
c
c write(6,'(1x,7E10.3,3I5)') hits(1,j),hits(2,j),hits(3,j),hits(4,j),
c + hits(5,j),hits(6,j),hits(7,j),nbv(1,j),volno,idtype
enddo
c
c
c nt_tsco = hits(4,1)*1.e9 ! decay time at Sci. in ns
nt_xsco = hits(1,1)
nt_ysco = hits(2,1)
nt_zsco = hits(3,1)
nt_desco = detot * 1000.
nt_descopos = descipos * 1000.
nt_descogam = descigam * 1000.
nt_descoele = desciele * 1000.
nt_descol = descil * 1000.
nt_descot = descit * 1000.
nt_descor = descir * 1000.
nt_descob = descib * 1000.
nt_volno = volno
nt_partcode = partcode
c
endif ! ~ hit in inner Sc.
c
nt_sconhits = nhits
c
c end out. scintillators
c
c---------------------------------
c
c fill NTUPLE
c
call hfnt(111)
c
c reset NTUPLE variables
c
nt_x0 = 0.
nt_y0 = 0.
nt_z0 = 0.
nt_px0 = 0.
nt_py0 = 0.
nt_pz0 = 0.
nt_p0 = 0.
nt_ipart = 0
nt_xe = 0.
nt_ye = 0.
nt_ze = 0.
nt_eposend = 0.
nt_tmcp = 0
nt_tsci = 0
nt_xsci = 0.
nt_ysci = 0.
nt_zsci = 0.
nt_desci = 0.
nt_descipos = 0.
nt_descigam = 0.
nt_desciele = 0.
c
nt_xsco = 0.
nt_ysco = 0.
nt_zsco = 0.
nt_desco = 0.
nt_descopos = 0.
nt_descogam = 0.
nt_descoele = 0.
c
nt_volno = 0
nt_partcode = 0
nt_mcpnhits = 0
nt_scinhits = 0
nt_sconhits = 0
c
nt_pxpos = 0.
nt_pypos = 0.
nt_pzpos = 0.
nt_demcp = 0.
c
nt_descil = 0.
nt_descit = 0.
nt_descir = 0.
nt_descib = 0.
c
nt_descol = 0.
nt_descot = 0.
nt_descor = 0.
nt_descob = 0.
c
c-----------------------------------
c
c some information on progress of simulation
c
if ( float(idevt)/float(nevent) .ge. 0.1 .and. percent(1) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 10%'
percent(1) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.2 .and. percent(2) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 20%'
percent(2) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.3 .and. percent(3) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 30%'
percent(3) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.4 .and. percent(4) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 40%'
percent(4) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.5 .and. percent(5) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 50%'
percent(5) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.6 .and. percent(6) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 60%'
percent(6) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.7 .and. percent(7) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 70%'
percent(7) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.8 .and. percent(8) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 80%'
percent(8) = 1
endif
if ( float(idevt)/float(nevent) .ge. 0.9 .and. percent(9) .eq. 0)
1 then
write(6,*) ' fraction of processed events: 90%'
percent(9) = 1
endif
if ( idevt .eq. nevent) then
write(6,*) ' fraction of processed events: 100%'
endif
c
c-----------------------------
c
return
end

161
geant3/src/lemsr/gustep.F Normal file
View File

@ -0,0 +1,161 @@
#include "geant321/pilot.h"
c
c----------------------------------------------------------------------
c
subroutine gustep
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix Version, part of the geant_lemsr program
c
c TP, 05-Apr-2001, PSI Unix, NT, Linux
c
c----------------------------------------------------------------------
c
implicit none
c
c include again some GEANT variables
c
#define CERNLIB_TYPE
#include "geant321/gckine.inc"
#include "geant321/gctrak.inc"
#include "geant321/gctmed.inc"
#include "geant321/gcking.inc"
#include "geant321/gcsets.inc"
#include "geant321/gcvolu.inc"
#include "geant321/gclist.inc"
#include "geant321/gcflag.inc"
c
c NTuple variable and other common variables
c
#include "cwn.inc"
#include "common.inc"
c
integer*4 j, part
integer*4 ihit
c
real*4 hits(7)
c
c------------------------------------------------------------------
c
c store secondary particle data to GEANT structure JSTAK
c
do j = 1,ngkine
part = gkin(5,j)
if ( part .eq. 4 ) then ! neutrinos, discard them
iflgk(j) = -1
elseif ( part .eq. 2) then ! save e+ in JKINE and transport
if ( l_pos) then
iflgk(j) = 1
else
iflgk(j) = -1
endif
elseif ( part .eq. 1) then
if ( l_gam ) then
iflgk(j) = 0 ! simply transport the rest
else
iflgk(j) = -1
endif
elseif ( part .eq. 3) then
if ( l_ele ) then
iflgk(j) = 0 ! simply transport the rest
else
iflgk(j) = -1
endif
else ! discard other part. in any
iflgk(j) = -1 ! case
endif
enddo
call gsking(0) ! store secondary particle data to JSTAK
c
c
c print*,' particle type = ',ipart
c print*,' tracking medium number = ',numed
c print*,' isvol = ',isvol
c print*,' current track position = ',vect(1),vect(2),vect(3)
c print*,' current momentum = ',vect(4),vect(5),vect(6)
c print*,' total momentum = ',vect(7)
c print*,' current kinetic energy = ',gekin
c print*,' number of steps for = ',nstep
c print*,' total energy lost = ',destep
c print*,' current TOF = ',tofg
c print*,' istop = ',istop
c print*,' igauto = ',idecad
c print*,' inwvol = ',inwvol
c print*,'---------------------'
c
c
if ( lplot(1) .gt. 0) call gdcxyz
c
c store spatial,momentum and energy loss data to GEANT data structure
c JXYZ
c
c call gsxyz
c call gpcxyz
c
c look where the e+ stops
c
if (ipart.eq.2.and.(inwvol.eq.3.or.istop.ge.1)) then
c
nt_xe = vect(1)
nt_ye = vect(2)
nt_ze = vect(3)
nt_eposend = 1000. * gekin
c
endif
c
c look now for detector hits
c
if ( idtype .eq. 1) then
c
hits(1) = vect(1)
hits(2) = vect(2)
hits(3) = vect(3)
hits(4) = tofg
hits(5) = gekin
hits(6) = destep
hits(7) = float(ipart)
call gsahit(iset,idet,itra,numbv,hits,ihit) ! MCP
if ( lplot(1) .gt. 0 ) then ! draw and plot hit
call gphits('MDET','DMCP')
call gdhits('MDET','DMCP',0,854,.3)
endif
c
endif
c
if ( idtype .eq. 2) then
c
hits(1) = vect(1)
hits(2) = vect(2)
hits(3) = vect(3)
hits(4) = tofg
hits(5) = gekin
hits(6) = destep
hits(7) = float(ipart)
call gsahit(iset,idet,itra,numbv,hits,ihit) ! inner Scint.
if ( lplot(1) .gt. 0 ) then ! plot hit
call gphits('ISSC','SCIS')
call gdhits('ISSC','SCIS',0,852,.3)
endif
c
endif
c
if ( idtype .eq. 3) then
c
hits(1) = vect(1)
hits(2) = vect(2)
hits(3) = vect(3)
hits(4) = tofg
hits(5) = gekin
hits(6) = destep
hits(7) = float(ipart)
call gsahit(iset,idet,itra,numbv,hits,ihit) ! outer Scint.
c
if ( lplot(1) .gt. 0 ) then ! plot hit
call gphits('OSSC','SCOS')
call gdhits('OSSC','SCOS',0,842,.3)
endif
c
endif
return
end

19
geant3/src/lemsr/gutrack.F Normal file
View File

@ -0,0 +1,19 @@
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine gutrack
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version, part of the geant_lemsr program
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
implicit none
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c type*,' routine GUTRACK'
call gtrack
return
end

29
geant3/src/lemsr/gutreve.F Normal file
View File

@ -0,0 +1,29 @@
c#include "geant321/pilot.h"
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
subroutine gutreve
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version, part of the geant_lemsr program
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
implicit none
c
c include again some GEANT variables
c
c INTEGER IKINE,ITRA,ISTAK,IVERT,IPART,ITRTYP,NAPART,IPAOLD
c REAL PKINE,AMASS,CHARGE,TLIFE,VERT,PVERT
c
c#include "geant321/gtkine.inc"
c#include "geant321/gckine.inc"
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c type*,' routine GUTREVE'
call gtreve
c
return
end

56
geant3/src/lemsr/michel.F Normal file
View File

@ -0,0 +1,56 @@
c
c----------------------------------------------------------------------------
c
subroutine michel(pos_tot, p0)
c
c routine to throw the e+ energy spectrum in muon decay.
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version
c
c TP, 05-Apr-2001, PSI Unix, NT, Linux
c
c----------------------------------------------------------------------------
c
implicit none
c
#include "common.inc"
c
real*4 pos_tot, p0
real*4 eta, temp, rho, ratio
real*4 random(2)
c
c---------------------------
c
c MICHEL distributed e+ energy by using NEUMANN rejection
c "rejection function" is 2.25 * eta
c
10 continue
#if defined (OS_UNIX)
call ranlux(random, 2)
#else
random(1) = ran(ix1)
random(2) = ran(ix2)
#endif
eta = sqrt(random(1))
temp = 2. * eta**2 * (3. - 2. * eta)
ratio = temp / (2.25 * eta)
rho = random(2)
if ( rho .gt. ratio ) then
rejcou = rejcou + 1
goto 10
else
cou = cou + 1
endif
c
pos_tot = MAX_TOTAL_ENERGY * eta
c
if ( pos_tot .ge. POS_MASS ) then
p0 = sqrt(pos_tot**2 - SQR_POS_MASS)
else
p0 = 0.
endif
c
return
end

52
geant3/src/lemsr/rotate.F Normal file
View File

@ -0,0 +1,52 @@
c
c----------------------------------------------------------------------------
c
subroutine rotate(decay_rot, time)
c
c routine to rotate the e+-momentum vector in case of an applied
c bfield; BFIELD is a 4-component array with bfield(4) = polarization.
c If Polarization = 0 return.
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version, part of the geant_lemsr program
c
c-----------------------------------------------------------------------------
c
implicit none
c
#include "common.inc"
c
real*4 decay_rot(3), time
real*4 omega, prec_angle, costh, sinth, cosph, sinph
c
c------------------------------------------------
c
if ( bfield(4) .eq. 0. ) return
c
if ( bfield(2) .ne. 0. ) then
omega = bfield(2) * TWO_PI * MU_GYRO ! cycle/nsec
prec_angle = omega * time
costh = cos(prec_angle)
sinth = -sin(prec_angle) ! clockwise rotation
cosph = 1.
sinph = 0.
c
call gdrot(decay_rot, costh, sinth, cosph, sinph)
c
endif
c
if ( bfield(3) .ne. 0. ) then
omega = bfield(3) * TWO_PI * MU_GYRO ! cycle/nsec
prec_angle = omega * time
cosph = cos(prec_angle)
sinph = -sin(prec_angle) ! clockwise rotation
costh = 1.
sinth = 0.
c
call gdrot(decay_rot, costh, sinth, cosph, sinph)
c
endif
c
return
end

903
geant3/src/lemsr/ugeom.F Normal file
View File

@ -0,0 +1,903 @@
#include "geant321/pilot.h"
c
c-----------------------------------------------------------------------
c
subroutine ugeom
c
c this subroutine defines the geometrical set-up.
c The z-axis is defines as the low-energy mu+ beam axis.
c This means x-axis --> horizontal, pointing to the left
c y-axis --> vertical, upwards.
c
c T. Prokscha, 15-Sep-2000, PSI
c
c Unix version, part of the geant_lemsr program
c
c TP, 05-Apr-2001, PSI Unix, NT, Linux
c
c------------------------------------------------------------------------
c
implicit none
c
#include "common.inc"
c
#define CERNLIB_TYPE
#include "geant321/gctmed.inc"
#include "geant321/gcsets.inc"
#include "geant321/gclist.inc"
#include "geant321/gcflag.inc"
#include "geant321/gcphys.inc"
#include "geant321/gccuts.inc"
#include "geant321/gckine.inc"
c
c-----------------------------------------------------------------------------
c
real*4 par_box(3)
real*4 par_tube(3)
real*4 par_tubs(5) ! segment of a tube
real*4 par_cone(5)
c
real*4 theta1,theta2,theta3,phi1,phi2,phi3 !rotation angles
integer*4 ivol !for the side arm
c
character*4 namesv(4) !for detector descr.
character*4 namesh(7) ! to define hit
integer*4 nbitsh(7) ! paramters on MCP
real*4 orig(7),fact(7)
c
real*4 atomw(9),zatom(9),wmat(9),dens ! to define MCP material
c
integer*4 mat_type ! material type for sample plate
real*4 off
real*4 he_shield_thickness
integer*4 Air, Vacuum, Al, Fe, Cu, Steel, Brass, Sapphire,
1 Macor, MCP, NaI, Scint ! medium numbers
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
data namesh /'x ','y ','z ','t ','e ','elos', 'part'/
data nbitsh /7*32/
data orig /3*100.,0.,0.,0.,0./
data fact /3*100.,1.e+9,1000000.,1000000.,1./
c
c the values for x,y,z,t,e and eloss will be internally converted
c by GEANT to integers;
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c MCP material data are taken from J.L.WIZA,Microchannel Plate Detectors,
c Nucl.Instr.Meth. 162 (1979), 587-601
c
c data dens /4./ ! density of MCP glass
data dens /2./ ! about half of the MCP is "empty" due to
c ! the MCP channels
data zatom /82.,8.,14.,19.,37.,56.,33.,55.,11./
data wmat /.479,.258,.182,.042,.018,.013,.004,.002,.001/
data atomw /207.19,15.999,28.086,39.10,85.47,137.33,74.92,
1 132.91,22.99/
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define materials
c
call gmate ! standard GEANT materials
c
c some additional materials
c
c define LEAD GLASS (MCP), NaI
c
call gsmixt(17,'MCPGLASS',atomw,zatom,dens,9,wmat)
c
c define NaI
c
zatom(1) = 11.
atomw(1) = 22.99
zatom(2) = 53
atomw(2) = 126.9
wmat(1) = 1.
wmat(2) = 1.
dens = 3.67
call gsmixt(18,'NaICALOR',atomw,zatom,dens,-2,wmat)
c
c define Scintillator NE102A
c
zatom(1) = 1.
atomw(1) = 1.0079
zatom(2) = 6.
atomw(2) = 12.011
wmat(1) = 1.104
wmat(2) = 1.
dens = 1.032
call gsmixt(19,'SCINT',atomw,zatom,dens,-2,wmat)
c
c stainless steel
c
dens = 7.930
zatom(1) = 26.
zatom(2) = 24.
zatom(3) = 28.
atomw(1) = 55.85
atomw(2) = 52.00
atomw(3) = 58.70
wmat(1) = 0.71
wmat(2) = 0.18
wmat(3) = 0.11
call gsmixt(20,'STEEL',atomw,zatom,dens,3,wmat)
c
c brass
c
dens = 8.67
zatom(1) = 29.
zatom(2) = 30.
atomw(1) = 63.55
atomw(2) = 65.38
wmat(1) = .80
wmat(2) = .19
c
call gsmixt(21,'BRASS',atomw,zatom,dens,2,wmat)
c
c define Sapphire == AL2O3
c
zatom(1) = 13.
atomw(1) = 26.98
zatom(2) = 8.
atomw(2) = 16.00
wmat(1) = 2.
wmat(2) = 3.
dens = 3.985
c
call gsmixt(22,'SAPPHIRE',atomw,zatom,dens,-2,wmat)
c
c define MACOR as SiO2 with higher density; the correct compound is
c
c (SiO2)46 (Al2O3)16 (MgO)17 (K2O)10 (B2O3)7
c
zatom(1) = 14.
atomw(1) = 28.0855
zatom(2) = 8.
atomw(2) = 16.00
wmat(1) = 1.
wmat(2) = 2.
dens = 2.52
c
call gsmixt(23,'MACOR',atomw,zatom,dens,-2,wmat)
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define tracking medium numbers corresponding to the material numbers
c stored in JMATE by calling GMATE
c
c 9 = Aluminum
c 10 = Iron
c 11 = Copper
c 15 = Air
c 16 = Vacuum
c 17 = MCP
c 18 = NaI <--- not used for present setup
c 19 = Scintillator
c 20 = Stainless steel
c 21 = brass <--- not used
c 22 = sapphire
c 23 = macor
c
isvol = 0
ifield = 0
c epsil = 0.0001 ! boundary crossing precision set to 1 micron
epsil = 0.001 ! boundary crossing precision set to 10 micron
c
call gstmed( 1,'Air $', 15, isvol, ifield, fieldm,
1 tmaxfd, STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 2,'Vacuum $', 16, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 3,'Aluminum $', 9, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 4,'Iron $', 10, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 5,'Copper $', 11, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 6,'Steel $', 20, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 7,'Brass $', 21, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 8,'Sapphire $', 22, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed( 9,'Macor $', 23, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
c
isvol = 1 ! sensitive volume; this means this volume is a detector
call gstmed(10,'MCPDetector$', 17, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed(11,'NaIDetector$', 18, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
call gstmed(12,'Scintillato$', 19, isvol, ifield, fieldm, tmaxfd,
1 STEMAX, deemax, epsil, stmin, 0, 0)
c
Air = 1
Vacuum = 2
Al = 3
Fe = 4
Cu = 5
Steel = 6
Brass = 7
Sapphire = 8
Macor = 9
MCP = 10
NaI = 11
Scint = 12
c
c redefine kinetic cut energies for muons in all media from its
c default value = 10 MeV to 10keV =0.00001 GeV
c
c ---> define cut energies via CUTS input cards
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define MOTHER VOLUME MARS
c
par_box(1) = 100.
par_box(2) = 100.
par_box(3) = 100.
c
call gsvolu('MARS','BOX ', Air, par_box, 3, ivol)
if (ivol.lt.0) goto 900 !this means an error during
c !call to GSVOLU
c------------------------------------
c
c vacuum tube of sample or MCP2, data from Balzers offer
c
par_tube(1) = 0.
par_tube(2) = 7.65
par_tube(3) = 16.2
c
call gsvolu('MCPV','TUBE', Vacuum, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c now stainless steel part of tube
c
par_tube(1) = 7.65
par_tube(2) = 7.95
par_tube(3) = 16.2
c
call gsvolu('MCPS','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c 160 CF flange upstream of MCP2 tube
c
par_tube(1) = 7.95
par_tube(2) = 10.125
par_tube(3) = 1.1
c
call gsvolu('F160','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c 100 CF flange at the end
c
par_tube(1) = 0.
par_tube(2) = 7.7
par_tube(3) = 1.0
c
call gsvolu('F100','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c 200 CF flange upstream of MCP2 tube to connect to gate valve chamber
c
par_tube(1) = 7.65
par_tube(2) = 10.325
par_tube(3) = 1.2
c
call gsvolu('F200','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c par_tube(1) = 0.
c par_tube(2) = 7.95
c par_tube(3) = 1.2
c
c call gsvolu('T200','TUBE', Vacuum, par_tube, 3, ivol)
c if (ivol.lt.0) goto 900
c
c the gate valve chamber
c
c
par_tube(1) = 0.
par_tube(2) = 10.325
par_tube(3) = 9.25
c
call gsvolu('GATV','TUBE', Vacuum, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c for simplicity choose outer diameter of the Gate valve chamber equal
c to the outer diameter of 200 CF flange
c
par_tube(1) = 10.325
par_tube(2) = 12.65
par_tube(3) = 9.25
c
call gsvolu('GATS','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c the L3 tube, vacuum
c
par_tube(1) = 0.
par_tube(2) = 10.0
par_tube(3) = 23.2
c
call gsvolu('L3VA','TUBE', Vacuum, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c the L3 tube, stainless steel
c
par_tube(1) = 10.0
par_tube(2) = 10.3
par_tube(3) = 23.2
c
call gsvolu('L3ST','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c conical anode
c
c the downstream part of HV piece
c
par_cone(1) = 2.5 ! half length dz
par_cone(2) = 4.8 ! inside radius at -dz
par_cone(3) = 6.2 ! outside radius at -dz
par_cone(4) = 3.3 ! inside radius at dz
par_cone(5) = 3.6 ! outside radius at dz
c
call gsvolu('RA-E','CONE', Steel, par_cone, 5, ivol)
if (ivol.lt.0) goto 900
c
c the middle part, still on HV
c
par_cone(1) = 2.0 ! half length dz
par_cone(2) = 5.8 ! inside radius at -dz
par_cone(3) = 6.2 ! outside radius at -dz
par_cone(4) = 4.8 ! inside radius at dz
par_cone(5) = 6.2 ! outside radius at dz
c
call gsvolu('RA-M','CONE', Steel, par_cone, 5, ivol)
if (ivol.lt.0) goto 900
c
c ground cylinder, stainless steel
c
par_tube(1) = 5.8
par_tube(2) = 6.2
par_tube(3) = 5.8
c
call gsvolu('RA-G','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c------------------------------------------------------
c
c sample cryo things
c
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
c
c sample holder
c
c 1. Cu plate (sample holder) on Cold finger, 0.5cm
c
par_tube(1) = 0.
par_tube(2) = 3.5
par_tube(3) = 0.25
c
if ( samp_mat .eq. 'saal' ) then
mat_type = Al ! Al
else
mat_type = Cu ! CU
endif
c
call gsvolu('SAH1','TUBE',mat_type,par_tube,3,ivol)
c print*,'SAH1 ivol = ',ivol
if (ivol.lt.0) goto 900
c
c 2. Cu plate (0.4cm), sample is mounted on this plate
c
par_tube(1) = 0.
par_tube(2) = 3.5
par_tube(3) = 0.2
c
call gsvolu('SAH2','TUBE',mat_type,par_tube,3,ivol)
c print*,'SAH2 ivol = ',ivol
if (ivol.lt.0) goto 900
c
c now select new color if sample plates are Al
c
if ( mat_type .eq. Al ) then
call gsatt('SAH1', 'COLO', 4)
call gsatt('SAH2', 'COLO', 4)
endif
c
c Sample mounting ring (0.1cm), 4cm inner diameter
c
par_tube(1) = 2.
par_tube(2) = 3.5
par_tube(3) = 0.05
c
call gsvolu('SAH3','TUBE',mat_type,par_tube,3,ivol)
c print*,'SAH3 ivol = ',ivol
if ( mat_type .eq. Al ) then
call gsatt('SAH3', 'COLO', 4) ! set color for aluminum
endif
if (ivol.lt.0) goto 900
if ( l_samp) then
write(6,*) ' '
if ( mat_type .eq. Cu ) then
write(6,*) ' Sample holder ring: Cu.'
else if ( mat_type .eq. Al ) then
write(6,*) ' Sample holder ring: Al.'
else
write(6,*) ' Sample holder ring: ', mat_type
endif
write(6,*) ' '
endif
c
c
c Sapphire mounted between 1. and 2. Cu plate, 0.6cm thick, 6cm diameter
c
par_tube(1) = 0.
par_tube(2) = 3.0
par_tube(3) = 0.3
c
call gsvolu('SAPH','TUBE', Sapphire, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
call gsatt('SAPH', 'COLO', 3) ! green color for saphhire
c
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
c
c cryo
c
c 1 cm plate, "cold finger"
c
par_tube(1) = 0.
par_tube(2) = 2.75
par_tube(3) = 0.5
c
call gsvolu('COFI','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c end plate of cryo (3cm diameter, 0.7cm thick)
c
par_tube(1) = 0.
par_tube(2) = 1.5
par_tube(3) = 0.35
c
call gsvolu('CRY1','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c "Waermetauscher", assume that it has 1cm open cylinder ?
c 5 cm long
c
par_tube(1) = 0.5
par_tube(2) = 1.5
par_tube(3) = 2.5
c
call gsvolu('CRY2','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c the mounting ring for He-shield
c
par_tube(1) = 3.8
par_tube(2) = 4.7
par_tube(3) = 0.55
c
call gsvolu('CRY3','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c add 2mm thick plate for mounting ring
c this is just to close the downstream side, I don't know if
c dimensions and position are right.
c
par_tube(1) = 1.5
par_tube(2) = 3.8
par_tube(3) = 0.1
c
call gsvolu('CRY4','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c call gsatt('CRY4', 'COLO', 3)
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
c
c He shield
c
par_tube(1) = 4.7
par_tube(2) = 4.8
par_tube(3) = 4.5
c
he_shield_thickness = par_tube(2) - par_tube(1)
if ( l_crsh ) then
write(6,*) ' '
write(6,*) 'Thickness of Cu-He-shield: ',he_shield_thickness,
1 ' cm.'
write(6,*) ' '
endif
c
call gsvolu('CRSH','TUBE', Cu, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c------------------------------------------------------
c
c MCP detector
c
par_tube(1) = 0.
par_tube(2) = 2.5
c
if ( l_run11 ) then
par_tube(3) = 0.15 ! Galileo Z-stack, 2.4mm totallength
! Galileo would have 3.0mm totallength
else
par_tube(3) = 0.12 ! hamamatsu Z-stack, 2.4mm totallength
endif
c
if ( l_mcp2 ) then
call gsvolu('DMCP','TUBE', MCP, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
endif
c
c MCP2 mounting and anode up to Run11
c
if ( .not. l_run11) then
c
c macor ring
c
par_tube(1) = 2.5
par_tube(2) = 3.5
par_tube(3) = 0.55
call gsvolu('MCPM', 'TUBE', Macor, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c anode, use MACOR
c
par_tube(1) = 0.
par_tube(2) = 3.5
par_tube(3) = 0.05
c
call gsvolu('MCPA', 'TUBE', Macor, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c stainless steel ring
c
par_tube(1) = 2.5
par_tube(2) = 3.5
par_tube(3) = 0.05
call gsvolu('MCSR', 'TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c stainless steel support ring
c
par_tube(1) = 2.75
par_tube(2) = 4.25
par_tube(3) = 0.2
call gsvolu('MCSS', 'TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c--------------------------------------------
c
else
c
c Run 11 MCP detector
c
c ceramic rings, use MACOR as material
c
par_tube(1) = 2.4
par_tube(2) = 3.25
par_tube(3) = 0.075
call gsvolu('MCPM', 'TUBE', Macor, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c anode in Run11,
c
par_box(1) = 3.65
par_box(2) = 3.65
par_box(3) = 0.4
call gsvolu('MCPA','BOX ', Steel, par_box, 3, ivol)
if (ivol.lt.0) goto 900
c
c now the anode is not completely made of steel but there
c the are "vacuum cylinders".
c
par_tube(1) = 0.
par_tube(2) = 2.75
par_tube(3) = 0.15
c
call gsvolu('ANVA','TUBE', Vacuum, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c stainless steel ring for MCP2 mounting
c
par_box(1) = 3.65
par_box(2) = 3.65
par_box(3) = 0.1
call gsvolu('MCSR','BOX ', Steel, par_box, 3, ivol)
if (ivol.lt.0) goto 900
c
par_tube(1) = 0
par_tube(2) = 2.75
par_tube(3) = 0.1
call gsvolu('MCVR','TUBE', Vacuum, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
c stainless steel support ring
c
par_tube(1) = 4.0
par_tube(2) = 4.8
par_tube(3) = 0.25
call gsvolu('MCSS','TUBE', Steel, par_tube, 3, ivol)
if (ivol.lt.0) goto 900
c
endif ! if .not. l_run11
c
c------------------------------------------------------
c
c inner scintillator
c
par_tubs(1) = 9.0 ! inner radius
par_tubs(2) = 9.5 ! outer radius
par_tubs(3) = 13.0 ! half length in z
par_tubs(4) = -45. ! starting angle
par_tubs(5) = +45. ! ending angle
call gsvolu('SCIS','TUBS', Scint, par_tubs, 5, ivol)
if (ivol.lt.0) goto 900
c
c outer scintillator
c
par_tubs(1) = 9.6 ! inner radius
par_tubs(2) = 10.1 ! outer radius
par_tubs(3) = 13.0 ! half length in z
par_tubs(4) = -45. ! starting angle
par_tubs(5) = +45. ! ending angle
call gsvolu('SCOS','TUBS', Scint, par_tubs, 5, ivol)
if (ivol.lt.0) goto 900
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define now the positions of vacuum chambers and detectors
c within the MARS
c
c z defines the beam direction
c
if ( l_mcps) then
call gspos('MCPS',1,'MARS',0.,0.,0.,0,'ONLY')
endif
if ( l_efla ) then
call gspos('F100',1,'MARS',0.,0.,17.2,0,'ONLY')
endif
c
call gspos('MCPV',1,'MARS',0.,0.,0.,0,'ONLY')
call gspos('F160',1,'MARS',0.,0.,-15.10,0,'ONLY')
c call gspos('T200',1,'MARS',0.,0.,-17.40,0,'ONLY')
call gspos('GATV',1,'MARS',0.,0.,-25.45,0,'ONLY')
call gspos('GATS',1,'MARS',0.,0.,-25.45,0,'ONLY')
call gspos('F200',1,'GATV',0.,0., 8.05,0,'ONLY')
call gspos('L3VA',1,'MARS',0.,0.,-57.90,0,'ONLY')
call gspos('L3ST',1,'MARS',0.,0.,-57.90,0,'ONLY')
c
c Ring anode
c
if ( l_run11) then
call gspos('RA-E',1,'MCPV',0.,0., -8.6,0,'ONLY')
call gspos('RA-M',1,'MCPV',0.,0.,-13.1,0,'ONLY')
call gspos('RA-G',1,'GATV',0.,0.,+3.35,0,'ONLY')
endif
c
c position of sample cryo things
c
c off = 0.25
off = 0.0
if ( l_samp ) then
call gspos('SAH3',1,'MCPV',0.,0.,+1.35,0,'ONLY')
call gspos('SAH2',1,'MCPV',0.,0.,+1.6,0,'ONLY')
call gspos('SAPH',1,'MCPV',0.,0.,+2.1,0,'ONLY')
call gspos('SAH1',1,'MCPV',0.,0.,+2.65+off,0,'ONLY')
endif
c
c off = 0.5
off = 0.0
if ( l_cryo ) then ! parts of the CryoVac cryostat
call gspos('COFI',1,'MCPV',0.,0.,+3.40+off,0,'ONLY')
call gspos('CRY1',1,'MCPV',0.,0.,+4.25+off,0,'ONLY')
call gspos('CRY2',1,'MCPV',0.,0.,+6.75+off,0,'ONLY')
call gspos('CRY3',1,'MCPV',0.,0.,+5.10+off,0,'ONLY')
call gspos('CRY4',1,'MCPV',0.,0.,+5.10+off,0,'ONLY')
endif
c
if ( l_crsh ) then
call gspos('CRSH',1,'MCPV',0.,0.,+1.6,0,'ONLY')
endif
c
c position of MCP
c
if ( l_mcp2 ) then
if ( l_run11) then
call gspos('DMCP',1,'MCPV',0.,0.,+1.55,0,'ONLY') ! Galileo MCP
else
call gspos('DMCP',1,'MCPV',0.,0.,+1.52,0,'ONLY') ! Hamamatsu MCP
endif
endif
c
c anode and stainless steel ring
c
if ( l_mcpa ) then
c
if ( .not. l_run11 ) then
c
c Run 10
c
call gspos('MCPM', 1, 'MCPV', 0., 0., +1.85, 0, 'ONLY') ! MACOR ring
call gspos('MCPA', 1, 'MCPV', 0., 0., +3.45, 0, 'ONLY') ! MACOR anodee
call gspos('MCSS', 1, 'MCPV', 0., 0., +3.90, 0, 'ONLY') ! steel support ring
call gspos('MCSR', 1, 'MCPV', 0., 0., +2.45, 0, 'ONLY') ! guard ring
call gspos('MCSR', 2, 'MCPV', 0., 0., +2.95, 0, 'ONLY') ! guard ring
else
c
c Run 11
c
call gspos('MCPM',1,'MCPV',0.,0.,+1.325,0,'ONLY') ! ceramic ring
call gspos('MCPM',2,'MCPV',0.,0.,+1.775,0,'ONLY') ! ceramic ring
call gspos('MCSR',1,'MCPV',0.,0.,+2.00,0,'ONLY') ! steel plate
call gspos('MCVR',1,'MCSR',0.,0.,+0.00,0,'ONLY') ! hole in plate
call gspos('MCPA',1,'MCPV',0.,0.,+3.10,0,'ONLY') ! the anode
call gspos('ANVA',1,'MCPA',0.,0.,-0.25,0,'ONLY') ! part of anode
call gspos('ANVA',2,'MCPA',0.,0.,+0.25,0,'ONLY') ! part of anode
call gspos('MCSS',1,'MCPV',0.,0.,+5.75,0,'ONLY') ! support ring
c
endif
c
endif
c
c position of the four scintillators
c
c--------------------------
c
c the LEFT detector
c
call gspos('SCIS',1,'MARS',0.,0.,0.,0,'ONLY')
call gspos('SCOS',1,'MARS',0.,0.,0.,0,'ONLY')
c
c--------------------------
c
c the TOP detector; for comments in rotation matrix, see RIGHT detector
c below
c
theta1 = 90.
phi1 = 90.
theta2 = 90.
phi2 = 180.
theta3 = 0.
phi3 = 0.
c
c GSROTM(97...) defines rotation matrix 97
c
call gsrotm(97,theta1,phi1,theta2,phi2,theta3,phi3)
c
call gspos('SCIS',2,'MARS',0.,0.,0.,97,'ONLY')
call gspos('SCOS',2,'MARS',0.,0.,0.,97,'ONLY')
c
c--------------------------------
c
c the RIGHT detector
c
c Now, I can get 3 copies of the SCIS and SCOS by defining a rotation
c matrix to place the other detectors in the MARS
c
c THETA1 = polar angle of new x axis
c PHI1 = azimuth of new x axis
c THETA2 = polar angle of new y axis
c PHI2 = azimuth of new y axis
c THETA3 = polar angle of new z axis
c PHI3 = azimuth of new z axis
c
c SCIS and SCOS above correspond to the LEFT detector of the LEM
c setup; now I make a copy to get the RIGHT detector
c
c z remains unchanged ==> THETA3 = 0, PHI3 = 0
c the x axis becomes -x ==> THETA1 = 90, PHI1 = 180
c the y axis becomes -y ==> THETA2 = 90, PHI2 = 270
c
c the unity matrix is
c
c theta1 = 90, phi1 = 0
c theta2 = 90, phi2 = 90
c theta3 = 0 , phi3 = 0
c
theta1 = 90.
phi1 = 180.
theta2 = 90.
phi2 = 270.
theta3 = 0.
phi3 = 0.
c GSROTM(98...) defines rotation matrix 98
c
call gsrotm(98,theta1,phi1,theta2,phi2,theta3,phi3)
c
call gspos('SCIS',3,'MARS',0.,0.,0.,98,'ONLY')
call gspos('SCOS',3,'MARS',0.,0.,0.,98,'ONLY')
c
c--------------------------
c
c the BOTTOM detector; for comments in rotation matrix, see RIGHT detector
c above
c
theta1 = 90.
phi1 = 270.
theta2 = 90.
phi2 = 0. ! == 360
theta3 = 0.
phi3 = 0.
c
c GSROTM(99...) defines rotation matrix 99
c
call gsrotm(99,theta1,phi1,theta2,phi2,theta3,phi3)
c
call gspos('SCIS',4,'MARS',0.,0.,0.,99,'ONLY')
call gspos('SCOS',4,'MARS',0.,0.,0.,99,'ONLY')
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define basic detector parameter and allocation of the detector
c volume to a SET using GSDET; a SET becomes useful in large
c detectors where it is possible to group together several detetectors
c like the different components of a calorimeter;
c in this simple case it isn't necessary, but GEANT aquires it
c
namesv(1) = 'DMCP'
c
c the following variables have been "included"
c
nvname = 1
numbv(1) = 1
idtype = 1 ! user defined detector ID
if ( l_mcp2 ) then
call gsdet('MDET','DMCP',nvname,namesv,numbv,idtype,100,
1 100,iset,idet)
c
c define now the hit parameters using GSDETH
c MCP : the interesting variables I want to store are
c
c - position x,y,z
c - time t
c - energy e
c - energyloss eloss
c - particle type ipart
c
call gsdeth('MDET','DMCP',7,namesh,nbitsh,orig,fact)
endif
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c define Scintillator detector
c
nvname = 1
namesv(1) = 'SCIS'
numbv(1) = 3
idtype = 2 ! user defined detector ID
call gsdet('ISSC','SCIS',nvname,namesv,numbv,idtype,100,100,
1 iset,idet)
call gsdeth('ISSC','SCIS',7,namesh,nbitsh,orig,fact)
nvname = 1
namesv(1) = 'SCOS'
numbv(1) = 3
idtype = 3 ! user defined detector ID
c
call gsdet('OSSC','SCOS',nvname,namesv,numbv,idtype,100,100,
1 iset, idet)
call gsdeth('OSSC','SCOS',7,namesh,nbitsh,orig,fact)
c
c call gpsets('MDET','DMCP')
c call gpsets('ISSC','SCIS')
c call gpsets('OSSC','SCOS')
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
return
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c error handling
c
900 write(6,'(''0GEANT_LEMSR-E- Error in routine GSVOLU...'',/)')
call exit
c
end