Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JGedanken.sh
Go to the documentation of this file.
1#!/usr/bin/env zsh
2#
3# \author mdejong
4version=1.0
5script=${0##*/}
6
7# ------------------------------------------------------------------------------------------
8#
9# Auxiliary script for Gedanken experiment.
10#
11# ------------------------------------------------------------------------------------------
12
13if [ -z $JPP_DIR ]; then
14 echo "Variable JPP_DIR undefined."
15 exit
16fi
17
18source $JPP_DIR/setenv.sh $JPP_DIR
19
20zmodload zsh/mathfunc
21
22# graphics
23
24set_variable: FORMAT GRAPHICS_FORMAT gif
25set_variable+ BATCH GRAPHICS_BATCH -B
26
27if do_usage $*; then
28 usage "$script <input file>"\
29 "\nAn example input file \"gedanken.txt\" is created with option \"?\"."
30fi
31
32if (( $# == 1 )); then
33 set_variable INPUT_FILE $1
34else
35 fatal "Wrong number of options."
36fi
37
38if [[ $INPUT_FILE == "?" ]]; then
39
40 cat>gedanken.txt<<EOF
41
42# general
43
44set_variable DEBUG 2
45set_variable WORKDIR ${TMPDIR:-/tmp}/gedanken
46set_variable NUMBER_OF_EVENTS 1000000
47set_variable NUMBER_OF_STEPS 20 # shower elongation
48
49# location of optical module
50
51set_variable R_M 25.0
52set_variable Z_M 50.0
53
54# particle
55
56set_variable E_GEV 20 # energy
57set_variable POS 0.0 0.0 0.0 # position
58set_variable DIR 0.0 0.0 1.0 # direction
59set_variable TYPE 13 # PDG
60EOF
61
62 return
63fi
64
65if [[ -f $INPUT_FILE ]]; then
66 source $INPUT_FILE
67else
68 fatal "Invalid input file <$INPUT_FILE>."
69fi
70
71if [[ ! -d $WORKDIR ]]; then
72 mkdir -p $WORKDIR
73fi
74
75set_variable DETECTOR $WORKDIR/detector.detx
76
77let "D_M = sqrt($R_M*$R_M + $Z_M*$Z_M)"
78let "CD = $Z_M / $D_M"
79
80if [[ ! -f $DETECTOR ]]; then
81
82 JModule \
83 -D 1001 \
84 -M 1 \
85 -P "$R_M 0.0 $Z_M" \
86 -o $DETECTOR \
87 -d $DEBUG
88fi
89
90if [[ ! -f $WORKDIR/gedanken.$TYPE.root ]]; then
91
92 $JPP_DIR/examples/JSirene/JGedanken \
93 -o $WORKDIR/gedanken.$TYPE.root \
94 -P "$POS" \
95 -D "$DIR" \
96 -E $E_GEV \
97 -T $TYPE \
98 -n $NUMBER_OF_EVENTS \
99 -a $DETECTOR \
100 -d 0 --!
101fi
102
103if [[ ! -f $WORKDIR/sirene.$TYPE.root ]]; then
104
105 JSirene \
106 -a $DETECTOR \
107 -f $WORKDIR/gedanken.$TYPE.root \
108 -o $WORKDIR/sirene.$TYPE.root \
109 -F $JPP_DATA/I%p.dat \
110 -N 0 \
111 -d $DEBUG --!
112fi
113
114printf "Generating histograms..."
115
116$JPP_DIR/examples/JSirene/JDomino \
117 -a $DETECTOR \
118 -f $WORKDIR/sirene.$TYPE.root \
119 -o $WORKDIR/domino.root \
120 -d 0
121
122printf "\n"
123
124
125set_variable PDF $JPP_DATA/J%p.dat
126set_variable EPS 1.0E-6
127set_variable PI $((acos(-1)))
128
129# constrain angle between [epsilon, pi - epsilon]
130#
131function constrain()
132{
133 if (( $1 > $PI - $EPS )); then
134 echo $(($PI - $EPS))
135 elif (( $1 < $EPS )); then
136 echo $EPS
137 else
138 echo $1
139 fi
140}
141
142typeset -A TITLE
143
144for id in {1..31}; do
145
146 printf "Generating PDFs PMT %2d\r" $id
147
148 getPMT -a $DETECTOR -p $id | read ID X Y Z DX DY DZ T0 STATUS
149
150 let "THETA = acos($DZ)"
151
152 if (( $DX > +$EPS )); then
153 let "PHI = atan(fabs($DY/$DX))"
154 elif (( $DX < -$EPS )); then
155 let "PHI = $PI - atan(fabs($DY/$DX))"
156 else
157 let "PHI = 0.5*$PI"
158 fi
159
160 let "THETA = $(constrain $THETA)"
161 let "PHI = $(constrain $PHI)"
162
163 TITLE[$id]="$(printf "PMT[%d] (%5.3f,%5.3f)" $id $THETA $PHI)"
164
165 if (( $TYPE == -13 || $TYPE == +13 || $TYPE == -15 || $TYPE == +15 )); then
166
167 JPlotPDF \
168 -f ${PDF/\%/1} \
169 -f ${PDF/\%/2} \
170 -f ${PDF/\%/3} \
171 -f ${PDF/\%/4} \
172 -f ${PDF/\%/5} \
173 -f ${PDF/\%/6} \
174 -D "$THETA $PHI" \
175 -R $R_M -E $E_GEV -z $Z_M \
176 -H "860 -10 850" \
177 -o $WORKDIR/f1\[${id}\].root
178
179 JDrawPDF \
180 -F 1 \
181 -F 2 \
182 -F 3 \
183 -F 4 \
184 -F 5 \
185 -F 6 \
186 -D "$THETA $PHI" \
187 -R $R_M -E $E_GEV -z $Z_M \
188 -H "860 -10 850" \
189 -o $WORKDIR/g1\[${id}\].root
190
191 else
192
193 # equivalent electro-magnetic energy
194
195 let "E = $(getPythia -P $TYPE -E $E_GEV)"
196
197 JPlotPDG \
198 -f ${PDF/\%/13} \
199 -f ${PDF/\%/14} \
200 -D "$THETA $PHI" \
201 -R $D_M -E $E -c $CD \
202 -N $NUMBER_OF_STEPS \
203 -H "860 -10 850" \
204 -o $WORKDIR/f1\[${id}\].root
205
206 JDrawPDG \
207 -F 13 \
208 -F 14 \
209 -D "$THETA $PHI" \
210 -R $D_M -E $E -c $CD \
211 -H "860 -10 850" \
212 -o $WORKDIR/g1\[${id}\].root
213 fi
214done
215
216printf "\n"
217
218
219if (( $TYPE == -13 || $TYPE == +13 || $TYPE == -15 || $TYPE == +15 )); then
220
221
222 JLatex \
223 -p "0.1 0.8 E" \
224 -p "0.3 0.8 =" \
225 -p "0.4 0.8 $E_GEV [GeV]" \
226 -p "0.1 0.7 R" \
227 -p "0.3 0.7 =" \
228 -p "0.4 0.7 $(printf "%5.1f [m]" $R_M)" \
229 -p "0.1 0.6 z" \
230 -p "0.3 0.6 =" \
231 -p "0.4 0.6 $(printf "%5.1f" $Z_M)" \
232 -@ "align = 11" \
233 -@ "size = 30" \
234 -o $WORKDIR/t1.root
235
236else
237
238 JLatex \
239 -p "0.1 0.8 E" \
240 -p "0.3 0.8 =" \
241 -p "0.4 0.8 $E_GEV [GeV]" \
242 -p "0.1 0.7 D" \
243 -p "0.3 0.7 =" \
244 -p "0.4 0.7 $(printf "%5.1f [m]" $D_M)" \
245 -p "0.1 0.6 cos(#theta)" \
246 -p "0.3 0.6 =" \
247 -p "0.4 0.6 $(printf "%5.2f" $CD)" \
248 -@ "align = 11" \
249 -@ "size = 30" \
250 -o $WORKDIR/t1.root
251fi
252
253JDraw \
254 -f $WORKDIR/t1.root:\.\* \
255 -o $WORKDIR/t1.$FORMAT $BATCH
256
257let "YMAX = 0.0"
258
259for id in {1..31}; do
260
261 let "Y = $(JPrintResult -f $WORKDIR/f1\[${id}\].root:h0 -F GetMaximum)"
262
263 if (( $Y > $YMAX )); then
264 let "YMAX = $Y"
265 fi
266done
267
268let "YMAX = 10**($(printf "%1.0f" $((log10($YMAX) + 0.5))))"
269let "YMIN = $YMAX * 1.0e-4"
270
271for id in {1..31}; do
272
273 printf "Generating graphics PMT %2d\r" $id
274
275 JPlot1D \
276 -f $WORKDIR/"domino.root:h\[${id}\]" \
277 -f $WORKDIR/f1\[${id}\].root:h0 \
278 -f $WORKDIR/g1\[${id}\].root:h0 \
279 -y "$YMIN $YMAX" -Y \
280 -> "#Deltat [ns]" \
281 -\^ "dP/dt npe/ns" \
282 -T "$TITLE[$id]" \
283 -O "HIST][" \
284 -o $WORKDIR/pdf_${id}.$FORMAT $BATCH
285done
286
287printf "\n"
288
289montage \
290 -tile 8x4 \
291 -geometry +0+0 \
292 $WORKDIR/pdf_{1..31}.$FORMAT \
293 $WORKDIR/t1.$FORMAT \
294 pdf.$TYPE.$FORMAT >& /dev/null
295
296
297for id in {1..31}; do
298
299 printf "Generating graphics PMT %2d\r" $id
300
301 JSum1D \
302 -f $WORKDIR/"domino.root:h\[${id}\]" \
303 -o $WORKDIR/"Domino.root" \
304 -B
305
306 JSum1D \
307 -f $WORKDIR/f1\[${id}\].root:h0 \
308 -o $WORKDIR/F1.root \
309 -B
310
311 JSum1D \
312 -f $WORKDIR/g1\[${id}\].root:h0 \
313 -o $WORKDIR/G1.root \
314 -B
315
316 let "Y = $(JPrintResult -f $WORKDIR/F1.root:h0 -F GetMaximum)"
317 let "YMAX = 10**($(printf "%1.1f" $((log10($Y)))) + 0.2)"
318 let "YMIN = $YMAX * 1.0e-1"
319
320 JPlot1D \
321 -f $WORKDIR/"Domino.root:h\[${id}\]" \
322 -f $WORKDIR/F1.root:h0 \
323 -f $WORKDIR/G1.root:h0 \
324 -y "$YMIN $YMAX" -Y \
325 -> "#Deltat [ns]" \
326 -\^ "#int P" \
327 -T "$TITLE[$id]" \
328 -o $WORKDIR/PDF_${id}.$FORMAT $BATCH
329
330 rm -f $WORKDIR/Domino.root
331 rm -f $WORKDIR/F1.root
332 rm -f $WORKDIR/G1.root
333done
334
335printf "\n"
336
337montage \
338 -tile 8x4 \
339 -geometry +0+0 \
340 $WORKDIR/PDF_{1..31}.$FORMAT \
341 $WORKDIR/t1.$FORMAT \
342 PDF.$TYPE.$FORMAT >& /dev/null
Generated by doxygen 1.12.0