Added to repository.
This commit is contained in:
parent
e38ca8310e
commit
c7720cc66b
110
mcv3k/inp/lem2004_2750_040_ag_bc200.in
Normal file
110
mcv3k/inp/lem2004_2750_040_ag_bc200.in
Normal file
@ -0,0 +1,110 @@
|
|||||||
|
SEEDS 15130905.
|
||||||
|
STRAG_ON
|
||||||
|
SCATT_ON
|
||||||
|
PARTICLE 105.6593 0.0001 5.
|
||||||
|
BEAM_P 27.50 0.0400 0.10 32.00
|
||||||
|
BEAM_XY 2.6 2.6 0. 3.
|
||||||
|
|
||||||
|
UFORMATTED
|
||||||
|
|
||||||
|
LAYER 3. 1.39 0.0003 1.0 0.0
|
||||||
|
Mylar wind
|
||||||
|
COMPONENT 1. 1.00794 376.
|
||||||
|
COMPONENT 6. 12.011 470.
|
||||||
|
COMPONENT 8. 15.9994 188.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 0. 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
LAYER 2. 1.032 0.020 4. 0.0
|
||||||
|
Beamcounter (PVTolu)
|
||||||
|
COMPONENT 1. 1.00794 138.
|
||||||
|
COMPONENT 6. 12.011 125.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 0.00031 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
LAYER 1. 4.510 0.003 3.0 0.0
|
||||||
|
Ti Vacuum window
|
||||||
|
COMPONENT 22. 47.880 1.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 3.1 1. 3.
|
||||||
|
|
||||||
|
LAYER 1. 2.702 0.0025 1.0 0.0
|
||||||
|
LN2 Shield, Al
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 5.3 2. -2.6 2.6 -2.6 2.6
|
||||||
|
|
||||||
|
LAYER 1. 2.702 0.0025 1.0 0.0
|
||||||
|
LHe Shield, Al
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 6.0 2. -2.1 2.1 -2.1 2.1
|
||||||
|
|
||||||
|
LAYER 1. 10.49 0.0125 59.0 0.0
|
||||||
|
Target,Ag subs.
|
||||||
|
COMPONENT 47. 107.8682 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 10.2 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
STOPPED 41.
|
||||||
|
|
||||||
|
LAYER 1. 1.771 0.001 1.0 0.0
|
||||||
|
10mu Ar layer
|
||||||
|
COMPONENT 18. 39.948 1.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 10.213 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
STOPPED 41.
|
||||||
|
|
||||||
|
|
||||||
|
LAYER 1. 2.702 .001 1.0 0.0
|
||||||
|
Beam Stop (MCP)
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 115.5 1. 2.0
|
||||||
|
|
||||||
|
OUTPUT 40. 2. 1000000. 100.
|
||||||
|
|
||||||
|
COMPUTE 10000. 1100000.
|
||||||
|
|
||||||
|
EXIT
|
||||||
|
|
||||||
|
|
||||||
|
************************************************************************
|
||||||
|
|
||||||
|
3mu Mylar window
|
||||||
|
250mu beam counter
|
||||||
|
|
||||||
|
LEM muE4: S1 und Ti-Fenster naeher an Moderator geschoben
|
||||||
|
125um Ag Moderator
|
||||||
|
|
||||||
|
neue Abstaende:
|
||||||
|
|
||||||
|
S1-Ti: 3.1cm (statt 4.3cm)
|
||||||
|
S1-AlN2: 5.3cm
|
||||||
|
S1-AlHe: 6.0cm
|
||||||
|
S1-Mode: 10.2cm
|
||||||
|
S1-MCP1: 115.5cm (statt 82.1cm)
|
||||||
|
|
||||||
|
Titanfenster 30um
|
||||||
|
Ag Substrat 125um, 59 Schichten a 2.119mu
|
||||||
|
zusaetzlicher Flansch 22mm zwischen KL und Spiegel;
|
||||||
|
zusätzliches CF160 Rohr (33.4cm) zwischen Spiegel
|
||||||
|
und MCP1
|
||||||
|
|
||||||
|
Target Position 10.2cm von Mylarfolie statt 11.9cm
|
||||||
|
|
||||||
|
Gesamtstrecke = 115.5 cm statt 82.1 cm
|
||||||
|
|
||||||
|
minimale Energie = 100 eV
|
||||||
|
|
||||||
|
Mylar folie 2um statt 3um
|
||||||
|
Al Fenster 25 um statt 30 um
|
||||||
|
|
||||||
|
Dicke MCP1 auf 0.1 cm gesetzt, damit alles stoppt in MCP1
|
||||||
|
|
||||||
|
Max. Impuls auf 32 MeV/c gesetzt
|
||||||
|
|
||||||
|
rechteckige Strahlflaeche 2.6*2.6cm^2; paralleler Strahl
|
||||||
|
|
||||||
|
|
99
mcv3k/inp/run11_2750_040_ag.in
Normal file
99
mcv3k/inp/run11_2750_040_ag.in
Normal file
@ -0,0 +1,99 @@
|
|||||||
|
SEEDS 15130905.
|
||||||
|
STRAG_ON
|
||||||
|
SCATT_ON
|
||||||
|
PARTICLE 105.6593 0.0001 5.
|
||||||
|
BEAM_P 27.50 0.0400 0.10 32.00
|
||||||
|
BEAM_XY 2.6 2.6 0. 3.
|
||||||
|
|
||||||
|
UFORMATTED
|
||||||
|
|
||||||
|
LAYER 3. 1.39 0.0002 1.0 0.0
|
||||||
|
Mylar wind
|
||||||
|
COMPONENT 1. 1.00794 376.
|
||||||
|
COMPONENT 6. 12.011 470.
|
||||||
|
COMPONENT 8. 15.9994 188.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 0. 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
LAYER 2. 1.032 0.02 5. 0.0
|
||||||
|
Beamcounter (PVTolu)
|
||||||
|
COMPONENT 1. 1.00794 138.
|
||||||
|
COMPONENT 6. 12.011 125.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 0.00031 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
LAYER 1. 4.510 0.003 3.0 0.0
|
||||||
|
Ti Vacuum window
|
||||||
|
COMPONENT 22. 47.880 1.
|
||||||
|
BUILD_TAB
|
||||||
|
GEOMETRY 4.3 1. 3.
|
||||||
|
|
||||||
|
LAYER 1. 2.702 0.0025 1.0 0.0
|
||||||
|
LN2 Shield, Al
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 5.3 2. -2.6 2.6 -2.6 2.6
|
||||||
|
|
||||||
|
LAYER 1. 2.702 0.0025 1.0 0.0
|
||||||
|
LHe Shield, Al
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 6.0 2. -2.1 2.1 -2.1 2.1
|
||||||
|
|
||||||
|
LAYER 1. 10.49 0.0125 59.0 0.0
|
||||||
|
Target,Ag subs.
|
||||||
|
COMPONENT 47. 107.8682 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 10.2 2. -1.5 1.5 -1.5 1.5
|
||||||
|
|
||||||
|
stopped 41.
|
||||||
|
|
||||||
|
LAYER 1. 2.702 .001 1.0 0.0
|
||||||
|
Beam Stop (MCP)
|
||||||
|
COMPONENT 13. 26.981539 1.
|
||||||
|
BUILD_TAB 1.E-5
|
||||||
|
GEOMETRY 82.1 1. 2.0
|
||||||
|
|
||||||
|
OUTPUT 40. 2. 1000000. 100.
|
||||||
|
|
||||||
|
COMPUTE 100000. 1100000.
|
||||||
|
|
||||||
|
EXIT
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
************************************************************************
|
||||||
|
|
||||||
|
RUN 11 piE3: S1 und Ti-Fenster naeher an Moderator geschoben
|
||||||
|
125um Ag Moderator
|
||||||
|
|
||||||
|
neue Abstaende:
|
||||||
|
|
||||||
|
S1-Ti: 4.3cm
|
||||||
|
S1-AlN2: 5.3cm
|
||||||
|
S1-AlHe: 6.0cm
|
||||||
|
S1-Mode: 10.2cm
|
||||||
|
S1-MCP1: 82.1cm
|
||||||
|
|
||||||
|
Titanfenster 30um
|
||||||
|
Ag Substrat 125um, 59 Schichten a 2.119mu
|
||||||
|
zusaetzlicher Flansch 22mm zwischen KL und Spiegel; d.h.
|
||||||
|
MCP1 jetzt 2.2 cm weiter von Target entfernt;
|
||||||
|
|
||||||
|
Target Position 10.2cm von Mylarfolie statt 11.9cm
|
||||||
|
|
||||||
|
Gesamtstrecke = 82.1 cm statt 83.8 cm
|
||||||
|
|
||||||
|
minimale Energie = 100 eV
|
||||||
|
|
||||||
|
Mylar folie 2um statt 3um
|
||||||
|
Al Fenster 25 um statt 30 um
|
||||||
|
|
||||||
|
|
||||||
|
Dicke MCP1 auf 0.1 cm gesetzt, damit alles stoppt in MCP1
|
||||||
|
|
||||||
|
Max. Impuls auf 32 MeV/c gesetzt
|
||||||
|
|
||||||
|
rechteckige Strahlflaeche 2.6*2.6cm^2; paralleler Strahl
|
||||||
|
(wie Run9)
|
||||||
|
|
9
mcv3k/paw/gauss_ener.for
Normal file
9
mcv3k/paw/gauss_ener.for
Normal file
@ -0,0 +1,9 @@
|
|||||||
|
function gauss_ener(x)
|
||||||
|
common /pawpar/ par(3)
|
||||||
|
sqpihalf=1.25331414
|
||||||
|
S2=1.414213562
|
||||||
|
fact1=10.*par(1)/sqpihalf/par(3)/erfc(-par(2)/s2/par(3))
|
||||||
|
c faktor 10 wegen 10keV binning im histogramm
|
||||||
|
gauss_ener=fact1*exp(-(x-par(2))**2/2/par(3)**2)
|
||||||
|
end
|
||||||
|
|
128
mcv3k/paw/mcv3k_ana.kumac
Normal file
128
mcv3k/paw/mcv3k_ana.kumac
Normal file
@ -0,0 +1,128 @@
|
|||||||
|
* ~mc/mcv3k/paw/mcv3k_ana.kumac
|
||||||
|
*
|
||||||
|
* kumac to create some histograms from NTUPLE files
|
||||||
|
* of the MCV3K simulations of the muonium experiment
|
||||||
|
* at PSI performed at piE1 in May 1994.
|
||||||
|
* one input parameter - the filename of the NTUPLE file which
|
||||||
|
* must have the extension .NT - is required without extension.
|
||||||
|
* the relevant variables are E0 (energy in MeV) and T0 (time
|
||||||
|
* of flight between S1 and MCP2 ins ns), the NTUPLE id is 111.
|
||||||
|
* 1.6 ns are substracted from the TOF's to correct for the
|
||||||
|
* TOF S1 - Target foil. This is necessary since the experimental
|
||||||
|
* TOF's are also corrected for this time.
|
||||||
|
*
|
||||||
|
* the new histograms are
|
||||||
|
*
|
||||||
|
* - 3 TOF's with different binnings (0.5, 1.0 and 2.0 ns)
|
||||||
|
* in the intervall [0,200]ns
|
||||||
|
* - 3 energy histograms in the intervals [0,1500] keV with 10keV
|
||||||
|
* bins
|
||||||
|
* [0,100]keV with 1keV bins and
|
||||||
|
* [0,30]keV with 1keV bins
|
||||||
|
* the purpose of the last two ones is to get information
|
||||||
|
* on the fraction of muons having energies in the 10 keV region
|
||||||
|
* to be able to make some predictions on muonium formation
|
||||||
|
*
|
||||||
|
*
|
||||||
|
* no cuts are necessary since the NTUPLE data contain only those
|
||||||
|
* muons hitting the MCP2. The distances and dimensions of the
|
||||||
|
* various elementes from the Mylar foil and S1 to the MCP2
|
||||||
|
* were rechecked fro the MCV3K simulation and they should be ok now.
|
||||||
|
*
|
||||||
|
* at the end the new histograms are written to a .RZ file with the
|
||||||
|
* same name as the NTUPLE file had.
|
||||||
|
*
|
||||||
|
* TP, 20-July-1995
|
||||||
|
*
|
||||||
|
*
|
||||||
|
* added a second parameter to indicate if a Postscript output should
|
||||||
|
* be done ( if [2]=PS then post script)
|
||||||
|
*
|
||||||
|
*
|
||||||
|
* TP, 18-dec-1995 renamed from RUN7_ANA.KUMAC
|
||||||
|
* TP, 09-feb-1998 changed bin size of ID1 from 257ps to 258ps
|
||||||
|
* TP, 14-sep-2000 Unix version
|
||||||
|
*-----------------------------------------------------------------------------
|
||||||
|
*
|
||||||
|
* open the NTUPLE file
|
||||||
|
*
|
||||||
|
close 0 | first close all open units
|
||||||
|
file = $mcv3k_out/[1].nt
|
||||||
|
hi/file 20 [file] 4096
|
||||||
|
message NTUPLE file = [file]
|
||||||
|
*
|
||||||
|
* book now the histograms needed
|
||||||
|
*
|
||||||
|
do i =1,6 | delete existing histograms
|
||||||
|
hi/del [i]
|
||||||
|
enddo
|
||||||
|
1d 1 'TOF - 1.7ns, 0.258ns binning' 501 -.129 129.129
|
||||||
|
1d 2 'TOF - 1.7ns, 0.5ns binning' 401 -.25 200.25
|
||||||
|
1d 3 'TOF - 1.7ns, 1.0ns binning' 201 -.5 200.5
|
||||||
|
1d 4 'energy in keV' 301 -5. 3005.
|
||||||
|
1d 5 'energy below 100 keV' 101 -.5 100.5
|
||||||
|
1d 6 'energy below 30 keV' 31 -.5 30.5
|
||||||
|
*
|
||||||
|
* fill the histograms
|
||||||
|
*
|
||||||
|
opt liny
|
||||||
|
set ysiz 28
|
||||||
|
zone 2 3
|
||||||
|
opt nfil
|
||||||
|
title [1]
|
||||||
|
*
|
||||||
|
* open PostScript file if desired
|
||||||
|
*
|
||||||
|
if [2] .eq. 'PS' .or. [2] .eq. 'ps' then
|
||||||
|
for/file 66 $mcv3k_out/[1].ps
|
||||||
|
meta 66 -111
|
||||||
|
*
|
||||||
|
endif
|
||||||
|
nt/pl 111.(t0-1.7) id>1 ! -1 | id>1 necessary because of a little
|
||||||
|
nt/pl 111.(t0-1.7) id>1 ! -2 | mistake in the routine which creates
|
||||||
|
nt/pl 111.(t0-1.7) id>1 ! -3 | the ntuple files
|
||||||
|
nt/pl 111.1000*e0 id>1 ! -4 | energy in keV
|
||||||
|
nt/pl 111.1000*e0 id>1 ! -5
|
||||||
|
nt/pl 111.1000*e0 id>1 ! -6
|
||||||
|
*
|
||||||
|
* write histograms to RZ files
|
||||||
|
*
|
||||||
|
n = 0
|
||||||
|
do i = 1, 6
|
||||||
|
temp = $HINFO([i],'ENTRIES')
|
||||||
|
if ([temp].eq.0) then
|
||||||
|
n = [n] + 1
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
if ( [n] .eq. 6) then
|
||||||
|
goto ret
|
||||||
|
endif
|
||||||
|
*
|
||||||
|
set ysiz
|
||||||
|
opt fit
|
||||||
|
set fit 111111
|
||||||
|
v/cre par(3) r 200 450 200
|
||||||
|
* hi/fit 4 gauss_ener.for l 3 par
|
||||||
|
hi/file 21 $mcv3k_out/[1].rz ! n | "! n" means new file
|
||||||
|
*
|
||||||
|
if [2] .eq. 'PS' .or. [2] .eq. 'ps' then
|
||||||
|
close 66
|
||||||
|
message PostScript file $mcv3k_out/[1].ps
|
||||||
|
endif
|
||||||
|
*
|
||||||
|
do i = 1,6
|
||||||
|
hrout [i]
|
||||||
|
enddo
|
||||||
|
message Histograms written to $mcv3k_out/[1].rz
|
||||||
|
close 20
|
||||||
|
close 21
|
||||||
|
set ysiz
|
||||||
|
opt file
|
||||||
|
title ' '
|
||||||
|
ret:
|
||||||
|
return
|
||||||
|
*
|
||||||
|
*------------------------------------------------------------------------------
|
||||||
|
*
|
||||||
|
* EOF ~/mc/mcv3k/paw/mcv3k_ana.kumac
|
||||||
|
*
|
130
mcv3k/paw/plot_range.kumac
Normal file
130
mcv3k/paw/plot_range.kumac
Normal file
@ -0,0 +1,130 @@
|
|||||||
|
*
|
||||||
|
* $mcv3k_root\paw\plot_range.kumac
|
||||||
|
*
|
||||||
|
* Plot mu+ ranges in Al up to p_in = 29.8 MeV/c and fit
|
||||||
|
* a1xp**3.5 to the data.
|
||||||
|
* Ranges calculated with mcv3k.
|
||||||
|
* Histograms saved in $mcv3_root\data\al_range.rz.
|
||||||
|
*
|
||||||
|
* requires one input parameter:
|
||||||
|
*
|
||||||
|
* [1] material, for example al
|
||||||
|
*
|
||||||
|
* TP, 28-Nov-2000, PSI LEM group
|
||||||
|
*
|
||||||
|
*-----------------------------------------------------------------------------
|
||||||
|
*
|
||||||
|
if ( [#] = 0 ) then
|
||||||
|
mess
|
||||||
|
mess Need input parameter to build input file name.
|
||||||
|
mess For example:
|
||||||
|
mess
|
||||||
|
mess plot_range al
|
||||||
|
mess
|
||||||
|
mess Input file name will be $mcv3k_root/data/al_range.rz
|
||||||
|
mess
|
||||||
|
exitm
|
||||||
|
endif
|
||||||
|
*
|
||||||
|
mat = [1]
|
||||||
|
*
|
||||||
|
*------------------------------------------------------------------------------
|
||||||
|
*
|
||||||
|
close 0 | close all open units
|
||||||
|
hi/de 0 | delete all histograms in //pawc
|
||||||
|
v/de * | delete all vectors
|
||||||
|
*
|
||||||
|
v/crea hid(100) | vector for histogram ID's
|
||||||
|
v/crea x0(2) | for defining graphic window
|
||||||
|
v/crea y0(2) | same for y
|
||||||
|
*
|
||||||
|
*-------------------------------
|
||||||
|
*
|
||||||
|
file = $mcv3k_root/data/[mat]_range.rz | file with range histograms
|
||||||
|
if ( $fexist([file]) = 0) then
|
||||||
|
mess
|
||||||
|
mess File [file] does not exist ?!
|
||||||
|
mess
|
||||||
|
exitm
|
||||||
|
endif
|
||||||
|
*
|
||||||
|
* open file
|
||||||
|
*
|
||||||
|
hi/file 20 [file]
|
||||||
|
hrin 0
|
||||||
|
close 20
|
||||||
|
nhist = 0
|
||||||
|
do i = 1, 300 | look for histograms; ID = 10 x momentum [MeV/c]
|
||||||
|
if ( $hexist([i]) ) then
|
||||||
|
mess Found histogram [i]
|
||||||
|
nhist = [nhist] + 1
|
||||||
|
v/inp hid([nhist]) [i]
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
*
|
||||||
|
* now we are ready to extract information from histograms
|
||||||
|
*
|
||||||
|
v/crea r([nhist]) | vector for mean ranges
|
||||||
|
v/crea er([nhist]) | vector for RMS range
|
||||||
|
v/crea p([nhist]) | vector for momenta
|
||||||
|
v/crea ep([nhist]) r 8*0. | error in momentum = 0
|
||||||
|
*
|
||||||
|
do i = 1, [nhist]
|
||||||
|
id = hid([i])
|
||||||
|
mom = [id]/10.
|
||||||
|
range = $hinfo([id],'MEAN')
|
||||||
|
rms = $hinfo([id],'RMS')
|
||||||
|
v/inp r([i]) [range]
|
||||||
|
v/inp er([i]) [rms]
|
||||||
|
v/inp p([i]) [mom]
|
||||||
|
enddo
|
||||||
|
xmax = [mom]*1.15
|
||||||
|
ymax = [range]*1.07
|
||||||
|
v/inp x0(1) 0.
|
||||||
|
v/inp y0(1) 0.
|
||||||
|
v/inp x0(2) [xmax]
|
||||||
|
v/inp y0(2) [ymax]
|
||||||
|
*
|
||||||
|
* plot
|
||||||
|
*
|
||||||
|
set ysiz
|
||||||
|
zone
|
||||||
|
opt grid
|
||||||
|
mat = $upper([mat])
|
||||||
|
title [mat]//', [m]^+! mean range, [D]p=0'
|
||||||
|
*
|
||||||
|
*
|
||||||
|
*
|
||||||
|
graph 2 x0 y0 aw
|
||||||
|
hpl/err p r ep er [nhist] 20 0.3
|
||||||
|
*
|
||||||
|
* fit to data par(1)xp**3.5
|
||||||
|
*
|
||||||
|
v/crea par(1) r 1.
|
||||||
|
opt nstat
|
||||||
|
opt fit
|
||||||
|
set fit 111
|
||||||
|
v/fit p r er $mcv3k_root/paw/range.f s 1 par
|
||||||
|
scale = par(1)
|
||||||
|
max = p(1)
|
||||||
|
func/pl [scale]*x**3.5 0. [max] s
|
||||||
|
atitle 'p (MeV/c)' 'mean range ([m]m)'
|
||||||
|
*
|
||||||
|
xpos = 0.05*x0(2)
|
||||||
|
ypos = 0.95*r([nhist])
|
||||||
|
str = $mcv3k_root/paw/plot"-#range.kumac
|
||||||
|
text [xpos] [ypos] [str] 0.2
|
||||||
|
ypos = 0.90*r([nhist])
|
||||||
|
text [xpos] [ypos] 'Fit function "J# P1 x p^3.5!' .25
|
||||||
|
*
|
||||||
|
*
|
||||||
|
*
|
||||||
|
title ' '
|
||||||
|
opt ngri
|
||||||
|
opt stat
|
||||||
|
*
|
||||||
|
*-------------------------------------------------------------------------------
|
||||||
|
*
|
||||||
|
|
||||||
|
|
||||||
|
|
50
mcv3k/src/Makefile
Normal file
50
mcv3k/src/Makefile
Normal file
@ -0,0 +1,50 @@
|
|||||||
|
#
|
||||||
|
# ~/mc/mcv3k/src/Makefile
|
||||||
|
#
|
||||||
|
# makefile to compile and link the MCV3K Monte Carlo program
|
||||||
|
#
|
||||||
|
# TP, 12-Se-2000, PSI
|
||||||
|
#
|
||||||
|
# TP, 21-Nov-2001: extended to Linux
|
||||||
|
#
|
||||||
|
#--------------------------------------------------------------------------
|
||||||
|
#
|
||||||
|
# for Digital Unix
|
||||||
|
#
|
||||||
|
#FF = f77
|
||||||
|
#CFLAGS = -c
|
||||||
|
#------------------------------
|
||||||
|
# for Linux
|
||||||
|
#
|
||||||
|
FF = g77
|
||||||
|
CFLAGS = -c -Wno-globals -finit-local-zero -fno-automatic
|
||||||
|
|
||||||
|
cern_lib = $(CERN_ROOT)/lib
|
||||||
|
|
||||||
|
BIN_DIR = $(HOME)/bin
|
||||||
|
OBJ = mcv3k.o mcv3k_sub.o
|
||||||
|
.f.o:
|
||||||
|
$(FF) $(CFLAGS) $*.f
|
||||||
|
|
||||||
|
all: mcv3k readmcv3k
|
||||||
|
|
||||||
|
mcv3k: $(OBJ)
|
||||||
|
$(FF) -o $(BIN_DIR)/mcv3k $(OBJ) $(cern_lib)/libmathlib.a \
|
||||||
|
$(cern_lib)/libkernlib.a
|
||||||
|
|
||||||
|
readmcv3k: read_mcv3k_out_nt.o
|
||||||
|
$(FF) -o $(BIN_DIR)/read_mcv3k_out_nt read_mcv3k_out_nt.o $(cern_lib)/libpacklib.a \
|
||||||
|
$(cern_lib)/libkernlib.a
|
||||||
|
|
||||||
|
|
||||||
|
clean:
|
||||||
|
rm -f *.o *~
|
||||||
|
|
||||||
|
install:
|
||||||
|
cp -pv define_mcv3k run_mcv3k $(BIN_DIR)
|
||||||
|
|
||||||
|
#
|
||||||
|
#--------------------------------------------------------------------------
|
||||||
|
#
|
||||||
|
# EOF
|
||||||
|
#
|
31
mcv3k/src/define_mcv3k
Executable file
31
mcv3k/src/define_mcv3k
Executable file
@ -0,0 +1,31 @@
|
|||||||
|
#!/bin/bash
|
||||||
|
#-----------------------------------------------
|
||||||
|
# ~/mc/mcv3k/define_mcv3k
|
||||||
|
#-----------------------------------------------
|
||||||
|
#
|
||||||
|
# script to define environmental variables and
|
||||||
|
# aliases to run the MCV3K Monte Carlo program.
|
||||||
|
#
|
||||||
|
# TP, 13-Sep-2000, PSI
|
||||||
|
#
|
||||||
|
#-----------------------------------------------
|
||||||
|
#
|
||||||
|
# first define directory variables
|
||||||
|
#
|
||||||
|
export mcv3k_root=$HOME/simulation/mcv3k
|
||||||
|
export mcv3k_bin=$HOME/bin
|
||||||
|
export mcv3k_src=$HOME/simulation/mcv3k/src
|
||||||
|
export mcv3k_inp=$HOME/simulation/mcv3k/inp
|
||||||
|
export mcv3k_out=$HOME/simulation/mcv3k/data
|
||||||
|
export mcv3k_log=$HOME/simulation/mcv3k/log
|
||||||
|
export mcv3k_paw=$HOME/simulation/mcv3k/paw
|
||||||
|
#
|
||||||
|
# define aliases
|
||||||
|
#
|
||||||
|
alias mcv3k='$mcv3k_bin/mcv3k'
|
||||||
|
alias mcv3knt='$mcv3k_bin/read_mcv3knt'
|
||||||
|
#
|
||||||
|
echo " "
|
||||||
|
echo "Variables and aliases defined for MCV3K . . ."
|
||||||
|
echo " "
|
||||||
|
#
|
2264
mcv3k/src/mcv3k.f
Normal file
2264
mcv3k/src/mcv3k.f
Normal file
File diff suppressed because it is too large
Load Diff
1386
mcv3k/src/mcv3k_sub.f
Normal file
1386
mcv3k/src/mcv3k_sub.f
Normal file
File diff suppressed because it is too large
Load Diff
51
mcv3k/src/read_mcv3k_out.f
Normal file
51
mcv3k/src/read_mcv3k_out.f
Normal file
@ -0,0 +1,51 @@
|
|||||||
|
c read_mcv3k_out.for
|
||||||
|
c
|
||||||
|
program read_mcv3k_out
|
||||||
|
c
|
||||||
|
c program to read the unformatted output file of the MCV3K_TEST program
|
||||||
|
c at the moment the "format" is
|
||||||
|
c 12 entries per record
|
||||||
|
c
|
||||||
|
c 1.entry number of particle I*4
|
||||||
|
c 2.entry number of layer I*4
|
||||||
|
c 3.-5.entry x,y,z R*4
|
||||||
|
c 6.-8.entry px,py,pz R*4
|
||||||
|
c 9.entry kinetic energy R*4
|
||||||
|
c 10.entry theta R*4
|
||||||
|
c 11.entry phi R*4
|
||||||
|
c 12.entry total time R*4
|
||||||
|
c
|
||||||
|
c later these data should be written to a NTUPLE to be analysed using PAW
|
||||||
|
c
|
||||||
|
c TP, PSI, 2-jun-1995
|
||||||
|
c
|
||||||
|
c-----------------------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
implicit none ! inquire variable definition
|
||||||
|
c
|
||||||
|
integer*4 lun ! logical unit for data file
|
||||||
|
integer*4 id, il ! particle and layer number
|
||||||
|
real*4 x,y,z,px,py,pz,ekin ! location, momentum and energy of particle
|
||||||
|
real*4 theta, phi,time ! direction angulars and total time of
|
||||||
|
c ! flight
|
||||||
|
c
|
||||||
|
c-----------------------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
write(6,*) 'Enter the LUN of the MCV3K output file : '
|
||||||
|
read(5,*) lun
|
||||||
|
do while ( .true. )
|
||||||
|
read(lun,err=11,end=10) id,il,x,y,z,px,py,pz,ekin,theta,phi,time
|
||||||
|
write(6,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
|
||||||
|
1 ekin,theta,phi,time
|
||||||
|
|
||||||
|
c
|
||||||
|
enddo
|
||||||
|
10 write(6,'('' Succesfully read data from LUN '',i)') lun
|
||||||
|
call exit
|
||||||
|
11 write(6,'('' Error during reading...'')')
|
||||||
|
call exit
|
||||||
|
end
|
||||||
|
c
|
||||||
|
c------------------------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
c EOF READ_MCV3K_OUT.FOR
|
277
mcv3k/src/read_mcv3k_out_nt.f
Normal file
277
mcv3k/src/read_mcv3k_out_nt.f
Normal file
@ -0,0 +1,277 @@
|
|||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c ~/mc/mcv3k/src/read_mcv3k_out_nt.f
|
||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
program read_mcv3k_out_nt
|
||||||
|
c
|
||||||
|
c program to read the unformatted output file of the MCV3K_TEST program
|
||||||
|
c data are written into a ntuple to use PAW to display data
|
||||||
|
c
|
||||||
|
c at the moment the data "format" is
|
||||||
|
c 12 entries per record
|
||||||
|
c
|
||||||
|
c 1.entry number of particle I*4
|
||||||
|
c 2.entry number of layer I*4
|
||||||
|
c 3.-5.entry x,y,z R*4
|
||||||
|
c 6.-8.entry px,py,pz R*4
|
||||||
|
c 9.entry kinetic energy R*4
|
||||||
|
c 10.entry theta R*4
|
||||||
|
c 11.entry phi R*4
|
||||||
|
c 12.entry total time R*4
|
||||||
|
c
|
||||||
|
c at the moment particle data at the exit of seven layers are written
|
||||||
|
c to the data file (the name must be FOR0##.DAT);
|
||||||
|
c
|
||||||
|
c TP, PSI, 13-jun-1995
|
||||||
|
c
|
||||||
|
c TP, 26-jun-1995: modified number of NTUPLE entries for
|
||||||
|
c use in analysing of the piE5 data
|
||||||
|
c
|
||||||
|
c TP, 27-jun-1995: modified again number of ntuple entries
|
||||||
|
c only the data in front of the last layer are
|
||||||
|
c selected to minimize file space
|
||||||
|
c there is a little error now: the first entry
|
||||||
|
c of the NTUPLE will be ID=1 and all other variables
|
||||||
|
c equal to 0, if the first ID read from file is not
|
||||||
|
c 1
|
||||||
|
c prompt now for the input file name without extension
|
||||||
|
c ( the default extension is .DAT !). the root file
|
||||||
|
c name will be used to create the output filename
|
||||||
|
c with the extension .NT
|
||||||
|
c
|
||||||
|
c TP, 10-Aug-1995: read filename without extension from command
|
||||||
|
c line using the CLI$GET_VALUE function
|
||||||
|
c
|
||||||
|
c TP, 30-Jan-1996: read initial and final data from file, so, add data
|
||||||
|
c to NTUPLE
|
||||||
|
c
|
||||||
|
c TP, 30-Jan-1996: changed back
|
||||||
|
c TP, 13-Sep-2000: changed to run on Windows NT 4.0 and Unix
|
||||||
|
c
|
||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
implicit none ! inquire variable definition
|
||||||
|
c
|
||||||
|
character*80 fname ! string containing filename without
|
||||||
|
integer*4 length ! extension and length
|
||||||
|
c
|
||||||
|
integer*4 id, il ! particle and layer number
|
||||||
|
real*4 x,y,z,px,py,pz,ekin ! location, momentum and energy of
|
||||||
|
c ! particle
|
||||||
|
real*4 theta, phi,time ! direction angulars and total time of
|
||||||
|
c ! flight
|
||||||
|
c
|
||||||
|
c for the NTUPLE
|
||||||
|
c
|
||||||
|
c first declarations for the initialisation of the HBOOK package
|
||||||
|
c
|
||||||
|
integer*4 MAXPAGES
|
||||||
|
parameter (MAXPAGES = 1000)
|
||||||
|
integer*4 ipawc(128 * MAXPAGES)
|
||||||
|
c
|
||||||
|
common /pawc/ ipawc
|
||||||
|
c
|
||||||
|
integer*4 iquest(100)
|
||||||
|
common /quest/iquest !needed to get enough space for NTUPLE
|
||||||
|
c
|
||||||
|
c for handling the NTUPLE file
|
||||||
|
c
|
||||||
|
integer*4 lunin /20/, lunnt /21/ ,istat
|
||||||
|
character*80 infile, ntfile
|
||||||
|
c
|
||||||
|
c for booking and filling of the NTUPLE
|
||||||
|
c
|
||||||
|
integer*4 nvar, nprime
|
||||||
|
c parameter nvar = 71 !number of variables in NTUPLE
|
||||||
|
parameter (nvar = 11) !number of variables in NTUPLE
|
||||||
|
parameter (nprime = 1000) !NPRIME words are allocated
|
||||||
|
character*8 tags(nvar) !short name of each variable
|
||||||
|
real*4 xtup(nvar) !array with the current data
|
||||||
|
c !points
|
||||||
|
data tags
|
||||||
|
1 / 'ID',
|
||||||
|
2 'X0', 'Y0', 'Z0', 'PX0', 'PY0', 'PZ0', 'E0', 'TH0', 'PH0', 'T0'/
|
||||||
|
c 3 'X1', 'Y1', 'Z1', 'PX1', 'PY1', 'PZ1', 'E1', 'TH1', 'PH1', 'T1'/
|
||||||
|
c 5 'X2', 'Y2', 'Z2', 'PX2', 'PY2', 'PZ2', 'E2', 'TH2', 'PH2', 'T2'/
|
||||||
|
c 6 'X3', 'Y3', 'Z3', 'PX3', 'PY3', 'PZ3', 'E3', 'TH3', 'PH3', 'T3'/
|
||||||
|
c 7 'X4', 'Y4', 'Z4', 'PX4', 'PY4', 'PZ4', 'E4', 'TH4', 'PH4', 'T4',
|
||||||
|
c 8 'X5', 'Y5', 'Z5', 'PX5', 'PY5', 'PZ5', 'E5', 'TH5', 'PH5', 'T5',
|
||||||
|
c 9 'X6', 'Y6', 'Z6', 'PX6', 'PY6', 'PZ6', 'E6', 'TH6', 'PH6', 'T6'
|
||||||
|
c 1 /
|
||||||
|
c
|
||||||
|
c other variables
|
||||||
|
c
|
||||||
|
integer*4 temp, n, i
|
||||||
|
logical*1 cont
|
||||||
|
c
|
||||||
|
c-----------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
c read first the filename for data input from an MCV3K run and make
|
||||||
|
c the NTUPLE filename
|
||||||
|
c
|
||||||
|
write(6,'(/,$,a)')
|
||||||
|
1 'Enter the filename with the MCV3K data without extension : '
|
||||||
|
read(5,'(a)') fname
|
||||||
|
call str_trim (fname, length)
|
||||||
|
c
|
||||||
|
c write(6,'('' length of string = '',i)') length
|
||||||
|
c
|
||||||
|
infile = fname(1:length)//'.dat'
|
||||||
|
ntfile = fname(1:length)//'.nt'
|
||||||
|
c
|
||||||
|
c open input file
|
||||||
|
c
|
||||||
|
open (lunin, file=infile, status = 'OLD', form = 'UNFORMATTED',
|
||||||
|
1 err= 1)
|
||||||
|
goto 2
|
||||||
|
1 continue
|
||||||
|
write(6,'('' %OPEN_MCV3K-E- error opening input file '',a)')
|
||||||
|
1 infile
|
||||||
|
call exit
|
||||||
|
2 continue
|
||||||
|
write(6,'('' %OPEN_MCV3K-I- successfully opened input file '',a)')
|
||||||
|
1 infile
|
||||||
|
c
|
||||||
|
c initialize HBOOK
|
||||||
|
c
|
||||||
|
call hlimit ( 128 * MAXPAGES )
|
||||||
|
c
|
||||||
|
c open NTUPLE file and book NTUPLE
|
||||||
|
c
|
||||||
|
iquest(10) = 65000 ! to allocate enough space for the NTUPLE
|
||||||
|
c
|
||||||
|
call hropen(lunnt,'MCV3K',ntfile,'NQ',4096,istat)
|
||||||
|
c
|
||||||
|
if (istat.eq.0) then ! everything is all right
|
||||||
|
write(6,
|
||||||
|
1 '('' %READ_MCV3K-I- succesfully opened NTUPLE file '',a)') ntfile
|
||||||
|
else
|
||||||
|
write(6,'('' %READ_MCV3K-E- error opening NTUPLE file '',a)')
|
||||||
|
1 ntfile
|
||||||
|
call exit
|
||||||
|
endif
|
||||||
|
c
|
||||||
|
call hbookn(111,'MCV3KDAT',nvar,'//MCV3K',nprime,tags)
|
||||||
|
c
|
||||||
|
temp = 1 ! initialize the temporary storage with the
|
||||||
|
c ! the first particle number
|
||||||
|
id = 1
|
||||||
|
n = 0
|
||||||
|
c
|
||||||
|
c endless loop on reading the file and filling the ntuple
|
||||||
|
c
|
||||||
|
do while ( .true. )
|
||||||
|
c
|
||||||
|
xtup(1) = temp
|
||||||
|
cont = .true.
|
||||||
|
c
|
||||||
|
do while ( cont ) ! read until a particle is lost
|
||||||
|
c ! then, the array XTUP will be written
|
||||||
|
c ! into the NTUPLE, and then XTUP will
|
||||||
|
c ! be cleared and the data of the next
|
||||||
|
c ! particle will be read
|
||||||
|
c
|
||||||
|
read(lunin,err=11,end=10) id,il,x,y,z,px,py,pz,ekin,theta,
|
||||||
|
1 phi,time
|
||||||
|
c write(6,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
|
||||||
|
c 1 ekin,theta,phi,time
|
||||||
|
c
|
||||||
|
c copy data into the ntuple array xtup
|
||||||
|
c
|
||||||
|
if ( id .eq. temp) then ! still data from one particle
|
||||||
|
c
|
||||||
|
xtup( 2 + n*10 ) = x
|
||||||
|
xtup( 3 + n*10 ) = y
|
||||||
|
xtup( 4 + n*10 ) = z
|
||||||
|
xtup( 5 + n*10 ) = px
|
||||||
|
xtup( 6 + n*10 ) = py
|
||||||
|
xtup( 7 + n*10 ) = pz
|
||||||
|
xtup( 8 + n*10 ) = ekin
|
||||||
|
xtup( 9 + n*10 ) = theta
|
||||||
|
xtup( 10 + n*10 ) = phi
|
||||||
|
xtup( 11 + n*10 ) = time
|
||||||
|
n = n + 1
|
||||||
|
else
|
||||||
|
call hfn (111, xtup)
|
||||||
|
c
|
||||||
|
c reset the xtup array to be ready for the next particle
|
||||||
|
c
|
||||||
|
do i = 1, nvar
|
||||||
|
xtup (i) = 0.
|
||||||
|
enddo
|
||||||
|
c
|
||||||
|
c we just read the initial paramter of the following
|
||||||
|
c particle, so, write them to the XTUP array
|
||||||
|
c
|
||||||
|
xtup( 2 ) = x
|
||||||
|
xtup( 3 ) = y
|
||||||
|
xtup( 4 ) = z
|
||||||
|
xtup( 5 ) = px
|
||||||
|
xtup( 6 ) = py
|
||||||
|
xtup( 7 ) = pz
|
||||||
|
xtup( 8 ) = ekin
|
||||||
|
xtup( 9 ) = theta
|
||||||
|
xtup( 10) = phi
|
||||||
|
xtup( 11) = time
|
||||||
|
temp = id
|
||||||
|
cont = .false.
|
||||||
|
n = 1
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
c
|
||||||
|
enddo
|
||||||
|
10 write(6,'('' Succesfully read data from file '',a)') infile
|
||||||
|
close (lunin)
|
||||||
|
call hfn (111, xtup) ! save also data from the last particle
|
||||||
|
call end_ntuple( lunnt, ntfile )
|
||||||
|
call exit
|
||||||
|
11 write(6,'('' Error during reading...'')')
|
||||||
|
close (lunin)
|
||||||
|
call end_ntuple( lunnt, ntfile )
|
||||||
|
call exit
|
||||||
|
c
|
||||||
|
end
|
||||||
|
c
|
||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
subroutine end_ntuple(lunnt, ntfile)
|
||||||
|
c
|
||||||
|
c routine to close the ntuple file
|
||||||
|
c
|
||||||
|
implicit none
|
||||||
|
c
|
||||||
|
integer*4 lunnt ! LUN of Ntuple file
|
||||||
|
character*(*) ntfile
|
||||||
|
integer*4 i
|
||||||
|
c
|
||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
call hrout (111, i, ' ')
|
||||||
|
call hrend ('MCV3K')
|
||||||
|
close (lunnt)
|
||||||
|
write(6,'('' %READ_MCV3K-I- closed NTUPLE file '',a,/)') ntfile
|
||||||
|
return
|
||||||
|
end
|
||||||
|
c
|
||||||
|
c------------------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
subroutine str_trim(string, length)
|
||||||
|
c
|
||||||
|
c routine to get the length of string without blanks
|
||||||
|
c
|
||||||
|
implicit none
|
||||||
|
c
|
||||||
|
character*(*) string
|
||||||
|
integer*4 i, length
|
||||||
|
c
|
||||||
|
i = 1
|
||||||
|
do while (string(i:i) .ne. ' ')
|
||||||
|
i = i + 1
|
||||||
|
enddo
|
||||||
|
length = i-1
|
||||||
|
return
|
||||||
|
end
|
||||||
|
c
|
||||||
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||||
|
c
|
||||||
|
c EOF READ_MCV3K_OUT.FOR
|
78
mcv3k/src/run_mcv3k
Executable file
78
mcv3k/src/run_mcv3k
Executable file
@ -0,0 +1,78 @@
|
|||||||
|
#!/bin/csh
|
||||||
|
#--------------------------------------
|
||||||
|
# $(HOME)/bin/run_mcv3k
|
||||||
|
#--------------------------------------
|
||||||
|
#
|
||||||
|
# shell script to run the MCV3K Monte Carlo program by Werner
|
||||||
|
# Jacobs for the simulation of energy loss, energy straggling
|
||||||
|
# and multiple scattering. Via xlsbatch it can be run in batch
|
||||||
|
# mode on the Digital Unix farm at PSI. After calculation data
|
||||||
|
# can be written to Ntuple file to be analyzed using PAW.
|
||||||
|
#
|
||||||
|
# Environment variables must be defined with
|
||||||
|
# source ~/mc/mcv3k/define_mcv3k or by putting this line into the
|
||||||
|
# .login script.
|
||||||
|
#
|
||||||
|
# Needs up to three input variables:
|
||||||
|
#
|
||||||
|
# $1 : the input filename without extension, expected to be located
|
||||||
|
# in mc_inp, with extension .in.
|
||||||
|
# $2 : the number of fortan file with unformatted output, for example,
|
||||||
|
# 40, this will assume an existing file fort.40 in mc_out.
|
||||||
|
# $3 : the number of fortran file with the formatted stop distribution.
|
||||||
|
#
|
||||||
|
#------------------------------------------------------------------------------
|
||||||
|
#
|
||||||
|
# TP, 13-Sep-2000, PSI
|
||||||
|
#
|
||||||
|
#-------------------------------------------------------------------------------
|
||||||
|
#
|
||||||
|
|
||||||
|
define_mcv3k
|
||||||
|
|
||||||
|
if ( $1 == "") then
|
||||||
|
echo " "
|
||||||
|
echo " MCV3K-E- You have to enter at least the input filename !"
|
||||||
|
echo " "
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
#
|
||||||
|
set infile = $mcv3k_inp/$1.in
|
||||||
|
set outfile = $mcv3k_out/$1.out
|
||||||
|
set datfile = $mcv3k_out/$1.dat
|
||||||
|
set stofile = $mcv3k_out/$1.stop
|
||||||
|
#
|
||||||
|
# check if file exists
|
||||||
|
#
|
||||||
|
if ( ! -e $infile ) then
|
||||||
|
echo " "
|
||||||
|
echo " MCV3K-E- file $infile does not exist ! "
|
||||||
|
echo " "
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
#
|
||||||
|
# now we are ready to run MCV3K; standard input will be re-directed
|
||||||
|
# to infile, standard output will be re-directed to outfile
|
||||||
|
#
|
||||||
|
cd $mcv3k_out
|
||||||
|
$mcv3k_bin/mcv3k < $infile > $outfile
|
||||||
|
#
|
||||||
|
# the simulation finished, now look, if we shall create a NT file
|
||||||
|
#
|
||||||
|
if ( $2 != "") then
|
||||||
|
mv fort.$2 $datfile
|
||||||
|
$mcv3k_bin/read_mcv3knt <<end
|
||||||
|
$datfile:r
|
||||||
|
end
|
||||||
|
#
|
||||||
|
# <<end means: read input from script file until end
|
||||||
|
# $datfile:r means: filename until period; the read_mcv3knt programs
|
||||||
|
# requires filenames without extension, it expects .dat
|
||||||
|
#
|
||||||
|
endif
|
||||||
|
#
|
||||||
|
if ( $3 != "") then
|
||||||
|
mv fort.$3 $stofile
|
||||||
|
endif
|
||||||
|
#
|
||||||
|
exit
|
Loading…
x
Reference in New Issue
Block a user