Iguana 0.0.0
Implementation Guardian of Analysis Algorithms
Loading...
Searching...
No Matches
iguana_ex_fortran_01_action_functions.f
Go to the documentation of this file.
1
11
12 use iso_c_binding
13 implicit none
14
15c ------------------------------------------------------------------
16c data declarations
17c ------------------------------------------------------------------
18c NOTE: using `iso_c_binding` types for input and output of C-bound
19c functions and subroutines, i.e., for HIPO and Iguana usage;
20c using standard F77 types may still work, but might be
21c compiler dependent in some cases
22
23c program parameters
24 integer*4 argc
25 character*1024 in_file ! HIPO file
26 integer num_events ! number of events to read
27 character*1024 config_file ! YAML configuration file
28 character*16 num_events_arg
29 logical config_file_set
30 character(kind=c_char, len=1024)
31 & in_file_cstr, config_file_cstr, etc_dir
32
33c HIPO and bank variables
34 integer counter ! event counter
35 integer(c_int) reader_status ! hipo event loop vars
36 integer(c_int) nrows, nrows_c ! number of rows
37 integer(c_int) nr ! unused
38 integer n_max ! max number of rows we can read
39 parameter(n_max=50)
40
41c REC::Particle columns
42 integer(c_int) pid(n_max)
43 real(c_float) px(n_max), py(n_max), pz(n_max)
44 real(c_float) vz(n_max)
45 integer(c_int) stat(n_max)
46 integer(c_int) sector(n_max)
47 real(c_float) torus(n_max)
48 integer(c_int) runnum(n_max)
49
50c iguana algorithm outputs
51 logical(c_bool) accept(n_max) ! filter
52 real(c_double) qx, qy, qz, qe ! q vector
53 real(c_double) q2, x, y, w, nu ! inclusive kinematics
54 real(c_double) beampz, targetm ! beam and target
55 integer(c_int) key_vz_filter ! key for Z-vertex filter
56 integer(c_int) key_inc_kin ! key for inclusive kinematics
57
58c iguana algorithm indices
59 integer(c_int) algo_eb_filter, algo_vz_filter,
60 & algo_inc_kin, algo_mom_cor
61
62c misc.
63 integer i
64 real p, p_ele
65 integer i_ele
66 logical found_ele
67
68c ------------------------------------------------------------------
69c open the input file
70c ------------------------------------------------------------------
71
72c parse arguments
73 num_events = 10
74 config_file_set = .false.
75 argc = iargc()
76 if(argc.lt.1) then
77 etc_dir = ''
78 call iguana_getconfiginstallationprefix(etc_dir)
79 print *, 'ERROR: please at least specify a HIPO_FILE'
80 print *, ''
81 print *, 'ARGS: ', 'HIPO_FILE', ' ', 'NUM_EVENTS'
82 print *, ''
83 print *, ' HIPO_FILE: ', 'the input HIPO file'
84 print *, ''
85 print *, ' NUM_EVENTS (optional): ',
86 & 'the number of events (0 for all)'
87 print *, ' default: ', num_events
88 print *, ''
89 print *, ' CONFIG_FILE (optional): ',
90 & 'algorithm configuration file'
91 print *, ' default: ', 'use the internal defaults'
92 print *, ' example config file: ',
93 & trim(etc_dir)//'/examples/my_z_vertex_cuts.yaml'
94 stop
95 else
96 call getarg(1, in_file)
97 in_file_cstr = trim(in_file)//c_null_char
98 end if
99 if(argc.ge.2) then
100 call getarg(2, num_events_arg)
101 read(num_events_arg,*) num_events
102 end if
103 if(argc.ge.3) then
104 call getarg(3, config_file)
105 config_file_cstr = trim(config_file)//c_null_char
106 config_file_set = .true.
107 endif
108
109c open the input HIPO file
110 call hipo_file_open(in_file_cstr)
111 reader_status = 0
112 counter = 0
113
114c ------------------------------------------------------------------
115c create iguana algorithms
116c ------------------------------------------------------------------
117
118c before anything for Iguana, call `iguana_create()`
119 call iguana_create()
120! call iguana_bindings_set_verbose() ! enable additional log print
121
122c then create the algorithm instances
123c - the 1st argument is an integer, the algorithm index, which
124c references the created instance of the algorithm; it will be
125c set after calling this subroutine and you will need it to call
126c other iguana subroutines (namely, "action functions")
127c - the 2nd argument is the algorithm name
128c IMPORTANT: any time you pass a string to iguana, be sure that
129c it has the correct termination by appending:
130c `//c_null_char`
131 call iguana_algo_create(
132 & algo_eb_filter,
133 & 'clas12::EventBuilderFilter'//c_null_char)
134 call iguana_algo_create(
135 & algo_vz_filter,
136 & 'clas12::ZVertexFilter'//c_null_char)
137 call iguana_algo_create(
138 & algo_inc_kin,
139 & 'physics::InclusiveKinematics'//c_null_char)
140 call iguana_algo_create(
141 & algo_mom_cor,
142 & 'clas12::MomentumCorrection'//c_null_char)
143
144c ------------------------------------------------------------------
145c configure and start iguana algorithms
146c ------------------------------------------------------------------
147
148c set log levels (don't forget `//c_null_char`)
149 call iguana_algo_set_log_level(
150 & algo_eb_filter, 'debug'//c_null_char)
151 call iguana_algo_set_log_level(
152 & algo_vz_filter, 'debug'//c_null_char)
153 call iguana_algo_set_log_level(
154 & algo_inc_kin, 'debug'//c_null_char)
155 call iguana_algo_set_log_level(
156 & algo_mom_cor, 'debug'//c_null_char)
157
158c configure algorithms with a configuration file
159 if(config_file_set) then
160 call iguana_set_config_file(config_file_cstr)
161 endif
162
163c start all algorithms, which "locks" their configuration
164 call iguana_start()
165
166c ------------------------------------------------------------------
167c event loop
168c ------------------------------------------------------------------
169
170 10 if(reader_status.eq.0 .and.
171 & (num_events.eq.0 .or. counter.lt.num_events)) then
172
173 print *, ''
174 print *, '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> event ', counter
175
176c read banks
177 call hipo_file_next(reader_status)
178 call hipo_read_bank('REC::Particle', nrows)
179 call hipo_read_bank('RUN::config', nrows_c)
180 call hipo_read_int('REC::Particle', 'pid', nr, pid, n_max)
181 call hipo_read_float('REC::Particle', 'px', nr, px, n_max)
182 call hipo_read_float('REC::Particle', 'py', nr, py, n_max)
183 call hipo_read_float('REC::Particle', 'pz', nr, pz, n_max)
184 call hipo_read_float('REC::Particle', 'vz', nr, vz, n_max)
185 call hipo_read_int('REC::Particle', 'status', nr, stat, n_max)
186 call hipo_read_float('RUN::config', 'torus', nr, torus, n_max)
187 call hipo_read_int('RUN::config', 'run', nr, runnum, n_max)
188
189c before using the Z-vertext filter, we must "prepare" the
190c algorithm's configuration for this event; the resulting
191c 'key_vz_filter' must be passed to the action function;
192 call iguana_clas12_zvertexfilter_prepareevent(
193 & algo_vz_filter, runnum(1), key_vz_filter)
194c similarly for the inclusive kinematics algorithm; use '-1'
195c for the beam energy, so RCDB is used to get the energy
196 call iguana_physics_inclusivekinematics_prepareevent(
197 & algo_inc_kin, runnum(1), -1, key_inc_kin)
198
199c call iguana filters
200c - the `logical` variable `accept` must be initialized to
201c `.true.`, since we will use it to "chain" the filters
202c - the event builder filter trivial: by default it accepts only
203c `REC::Particle::pid == 11 or -211` (simple example algorithm)
204c - the AND with the z-vertex filter is the final filter, `accept`
205 print *, '===> filter particles with iguana:'
206 do i=1, nrows
207 accept(i) = .true.
208 call iguana_clas12_eventbuilderfilter_filter(
209 & algo_eb_filter, pid(i), accept(i))
210 call iguana_clas12_zvertexfilter_filter(
211 & algo_vz_filter, vz(i), pid(i), stat(i),
212 & key_vz_filter, accept(i))
213 print *, ' i = ', i, ' pid = ', pid(i), ' vz = ', vz(i),
214 & ' status = ', stat(i), ' => accept = ', accept(i)
215 enddo
216
217c get sector number
218c FIXME: we have the algorithm `iguana::clas12::SectorFinder`,
219c but it is not yet compatible with Fortran bindings; until
220c then, assume the sector number is 1
221 do i=1, nrows
222 sector(i) = 1 ! FIXME
223 enddo
224
225c momentum corrections
226 if(nrows_c.lt.1) then
227 print *, '===> no RUN::config bank; cannot ',
228 & 'apply momentum corrections'
229 else
230 print *, '===> momentum corrections:'
231 do i=1, nrows
232 if(accept(i)) then
233 print *, ' i = ', i
234 print *, ' before: p = (', px(i), py(i), pz(i), ')'
235 call iguana_clas12_momentumcorrection_transform(
236 & algo_mom_cor,
237 & px(i), py(i), pz(i),
238 & sector(i), pid(i), torus(1),
239 & px(i), py(i), pz(i))
240 print *, ' after: p = (', px(i), py(i), pz(i), ')'
241 endif
242 enddo
243 endif
244
245c simple electron finder: trigger and highest |p|
246 p_ele = 0
247 found_ele = .false.
248 print *, '===> finding electron...'
249 do i=1, nrows
250 if(accept(i)) then
251 print *, ' i = ', i, ' status = ', stat(i)
252 if(pid(i).eq.11 .and. stat(i).lt.0) then
253 p = sqrt(px(i)**2 + py(i)**2 + pz(i)**2)
254 if(p.gt.p_ele) then
255 i_ele = i
256 p_ele = p
257 found_ele = .true.
258 endif
259 endif
260 endif
261 enddo
262 if(found_ele) then
263 print *, '===> found DIS electron:'
264 print *, ' i = ', i_ele
265 print *, ' p = ', p_ele
266 else
267 print *, '===> no DIS electron'
268 endif
269
270c compute DIS kinematics with iguana, if electron is found
271 if(found_ele) then
272 call iguana_physics_inclusivekinematics_computefromlepton(
273 & algo_inc_kin,
274 & px(i_ele), py(i_ele), pz(i_ele),
275 & key_inc_kin,
276 & qx, qy, qz, qe,
277 & q2, x, y, w, nu,
278 & beampz, targetm)
279 print *, '===> inclusive kinematics:'
280 print *, ' q = (', qx, qy, qz, qe, ')'
281 print *, ' Q2 = ', q2
282 print *, ' x = ', x
283 print *, ' y = ', y
284 print *, ' W = ', w
285 print *, ' nu = ', nu
286 endif
287
288 counter = counter + 1
289 goto 10
290 endif
291
292c ------------------------------------------------------------------
293c cleanup
294c ------------------------------------------------------------------
295
296c don't forget to call `iguana_stop()` to stop the algorithms
297c and free the allocated memory
298 call iguana_stop()
299
300 end program
program iguana_ex_fortran_01_action_functions
Fortran example program.