Jpp 20.0.0-27-g39925593c-D
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
103set_variable CDF $JPP_DATA/I%p.dat.gz
104
105if [[ ! -f $WORKDIR/sirene.$TYPE.root ]]; then
106
107 JSirene \
108 -a $DETECTOR \
109 -f $WORKDIR/gedanken.$TYPE.root \
110 -o $WORKDIR/sirene.$TYPE.root \
111 -F $CDF \
112 -N 0 \
113 -d $DEBUG --!
114fi
115
116printf "Generating histograms..."
117
118$JPP_DIR/examples/JSirene/JDomino \
119 -a $DETECTOR \
120 -f $WORKDIR/sirene.$TYPE.root \
121 -o $WORKDIR/domino.root \
122 -d 0
123
124printf "\n"
125
126
127set_variable PDF $JPP_DATA/J%p.dat.gz
128set_variable EPS 1.0E-6
129set_variable PI $((acos(-1)))
130
131# constrain angle between [epsilon, pi - epsilon]
132#
133function constrain()
134{
135 if (( $1 > $PI - $EPS )); then
136 echo $(($PI - $EPS))
137 elif (( $1 < $EPS )); then
138 echo $EPS
139 else
140 echo $1
141 fi
142}
143
144typeset -A TITLE
145
146for id in {1..31}; do
147
148 printf "Generating PDFs PMT %2d\r" $id
149
150 getPMT -a $DETECTOR -p $id | read ID X Y Z DX DY DZ T0 STATUS
151
152 let "THETA = acos($DZ)"
153
154 if (( $DX > +$EPS )); then
155 let "PHI = atan(fabs($DY/$DX))"
156 elif (( $DX < -$EPS )); then
157 let "PHI = $PI - atan(fabs($DY/$DX))"
158 else
159 let "PHI = 0.5*$PI"
160 fi
161
162 let "THETA = $(constrain $THETA)"
163 let "PHI = $(constrain $PHI)"
164
165 TITLE[$id]="$(printf "PMT[%d] (%5.3f,%5.3f)" $id $THETA $PHI)"
166
167 if (( $TYPE == -13 || $TYPE == +13 || $TYPE == -15 || $TYPE == +15 )); then
168
169 JPlotPDF \
170 -f ${PDF/\%/1} \
171 -f ${PDF/\%/2} \
172 -f ${PDF/\%/3} \
173 -f ${PDF/\%/4} \
174 -f ${PDF/\%/5} \
175 -f ${PDF/\%/6} \
176 -D "$THETA $PHI" \
177 -R $R_M -E $E_GEV -z $Z_M \
178 -H "860 -10 850" \
179 -o $WORKDIR/f1\[${id}\].root
180
181 JDrawPDF \
182 -F 1 \
183 -F 2 \
184 -F 3 \
185 -F 4 \
186 -F 5 \
187 -F 6 \
188 -D "$THETA $PHI" \
189 -R $R_M -E $E_GEV -z $Z_M \
190 -H "860 -10 850" \
191 -o $WORKDIR/g1\[${id}\].root
192
193 else
194
195 # equivalent electro-magnetic energy
196
197 let "E = $(getPythia -P $TYPE -E $E_GEV)"
198
199 JPlotPDG \
200 -f ${PDF/\%/13} \
201 -f ${PDF/\%/14} \
202 -D "$THETA $PHI" \
203 -R $D_M -E $E -c $CD \
204 -N $NUMBER_OF_STEPS \
205 -H "860 -10 850" \
206 -o $WORKDIR/f1\[${id}\].root
207
208 JDrawPDG \
209 -F 13 \
210 -F 14 \
211 -D "$THETA $PHI" \
212 -R $D_M -E $E -c $CD \
213 -H "860 -10 850" \
214 -o $WORKDIR/g1\[${id}\].root
215 fi
216done
217
218printf "\n"
219
220
221if (( $TYPE == -13 || $TYPE == +13 || $TYPE == -15 || $TYPE == +15 )); then
222
223
224 JLatex \
225 -p "0.1 0.8 E" \
226 -p "0.3 0.8 =" \
227 -p "0.4 0.8 $E_GEV [GeV]" \
228 -p "0.1 0.7 R" \
229 -p "0.3 0.7 =" \
230 -p "0.4 0.7 $(printf "%5.1f [m]" $R_M)" \
231 -p "0.1 0.6 z" \
232 -p "0.3 0.6 =" \
233 -p "0.4 0.6 $(printf "%5.1f" $Z_M)" \
234 -@ "align = 11" \
235 -@ "size = 30" \
236 -o $WORKDIR/t1.root
237
238else
239
240 JLatex \
241 -p "0.1 0.8 E" \
242 -p "0.3 0.8 =" \
243 -p "0.4 0.8 $E_GEV [GeV]" \
244 -p "0.1 0.7 D" \
245 -p "0.3 0.7 =" \
246 -p "0.4 0.7 $(printf "%5.1f [m]" $D_M)" \
247 -p "0.1 0.6 cos(#theta)" \
248 -p "0.3 0.6 =" \
249 -p "0.4 0.6 $(printf "%5.2f" $CD)" \
250 -@ "align = 11" \
251 -@ "size = 30" \
252 -o $WORKDIR/t1.root
253fi
254
255JDraw \
256 -f $WORKDIR/t1.root:\.\* \
257 -o $WORKDIR/t1.$FORMAT $BATCH
258
259let "YMAX = 0.0"
260
261for id in {1..31}; do
262
263 let "Y = $(JPrintResult -f $WORKDIR/f1\[${id}\].root:h0 -F GetMaximum)"
264
265 if (( $Y > $YMAX )); then
266 let "YMAX = $Y"
267 fi
268done
269
270let "YMAX = 10**($(printf "%1.0f" $((log10($YMAX) + 0.5))))"
271let "YMIN = $YMAX * 1.0e-4"
272
273for id in {1..31}; do
274
275 printf "Generating graphics PMT %2d\r" $id
276
277 JPlot1D \
278 -f $WORKDIR/"domino.root:h\[${id}\]" \
279 -f $WORKDIR/f1\[${id}\].root:h0 \
280 -f $WORKDIR/g1\[${id}\].root:h0 \
281 -y "$YMIN $YMAX" -Y \
282 -> "#Deltat [ns]" \
283 -\^ "dP/dt npe/ns" \
284 -T "$TITLE[$id]" \
285 -O "HIST][" \
286 -o $WORKDIR/pdf_${id}.$FORMAT $BATCH
287done
288
289printf "\n"
290
291montage \
292 -tile 8x4 \
293 -geometry +0+0 \
294 $WORKDIR/pdf_{1..31}.$FORMAT \
295 $WORKDIR/t1.$FORMAT \
296 pdf.$TYPE.$FORMAT >& /dev/null
297
298
299for id in {1..31}; do
300
301 printf "Generating graphics PMT %2d\r" $id
302
303 JSum1D \
304 -f $WORKDIR/"domino.root:h\[${id}\]" \
305 -o $WORKDIR/"Domino.root" \
306 -B
307
308 JSum1D \
309 -f $WORKDIR/f1\[${id}\].root:h0 \
310 -o $WORKDIR/F1.root \
311 -B
312
313 JSum1D \
314 -f $WORKDIR/g1\[${id}\].root:h0 \
315 -o $WORKDIR/G1.root \
316 -B
317
318 let "Y = $(JPrintResult -f $WORKDIR/F1.root:h0 -F GetMaximum)"
319 let "YMAX = 10**($(printf "%1.1f" $((log10($Y)))) + 0.2)"
320 let "YMIN = $YMAX * 1.0e-1"
321
322 JPlot1D \
323 -f $WORKDIR/"Domino.root:h\[${id}\]" \
324 -f $WORKDIR/F1.root:h0 \
325 -f $WORKDIR/G1.root:h0 \
326 -y "$YMIN $YMAX" -Y \
327 -> "#Deltat [ns]" \
328 -\^ "#int P" \
329 -T "$TITLE[$id]" \
330 -o $WORKDIR/PDF_${id}.$FORMAT $BATCH
331
332 rm -f $WORKDIR/Domino.root
333 rm -f $WORKDIR/F1.root
334 rm -f $WORKDIR/G1.root
335done
336
337printf "\n"
338
339montage \
340 -tile 8x4 \
341 -geometry +0+0 \
342 $WORKDIR/PDF_{1..31}.$FORMAT \
343 $WORKDIR/t1.$FORMAT \
344 PDF.$TYPE.$FORMAT >& /dev/null
Generated by doxygen 1.10.0