« back to all changes in this revision
Viewing changes to Template/Source/combine_runs.f
- browse files at revision 147.7.143
- compare with another revision
- download diff
- download tarball
- view history from revision 147.7.143
- Committer: olivier Mattelaer
- Date: 2012-09-27 21:47:26 UTC
- mfrom: (147.13.15 transfer_file_on_exit)
- mto: (147.16.5 new_history)
- mto: This revision was merged to the branch mainline in revision 223.
- Revision ID: olivier.mattelaer@uclouvain.be-20120927214726-8jn3hdm0jiulk14r
merge lp:~maddevelopers/madgraph5/transfer_on_exit
- files added:
- madgraph/various/combine_runs.py
- files removed:
- Template/Source/combine_runs.f
- files modified:
- Template/Source/PDF/opendata.f
added
removed
Template/Source/combine_runs.f
1
program test
2
c*****************************************************************
3
c Program to combine results from channels that have been
4
c split into multiple jobs. Multi-job channels are identified
5
c by local file mjobs.dat in the channel directory.
6
c****************************************************************
7
implicit none
8
c
9
c Constants
10
c
11
include 'run_config.inc'
12
integer maxsubprocesses
13
parameter (maxsubprocesses=999)
14
c
15
c Local
16
c
17
character*300 procname(maxsubprocesses)
18
character*300 channame(maxsubprocesses)
19
character*300 pathname
20
integer iproc,nproc
21
integer ichan,nchan
22
integer jc
23
24
c-----
25
c Begin Code
26
c-----
27
28
c
29
c Get list of subprocesses
30
c
31
call get_subprocess(procname,nproc)
32
33
do iproc = 1, nproc !Loop over each subprocess
34
jc = index(procname(iproc),' ')
35
pathname = procname(iproc)(1:jc)
36
c
37
c get list of integration channels for this subprocess
38
c
39
call get_channels(pathname,channame,nchan)
40
do ichan = 1, nchan
41
call sum_multichannel(channame(ichan))
42
enddo
43
enddo
44
45
end
46
47
subroutine get_subprocess(subname,ns)
48
c*****************************************************************
49
c Opens file subproc.mg to determine all subprocesses
50
c****************************************************************
51
implicit none
52
c
53
c Constants
54
c
55
character*(*) plist
56
parameter (plist='subproc.mg')
57
integer maxsubprocesses
58
parameter (maxsubprocesses=999)
59
c
60
c Arguments
61
c
62
character*(*) subname(maxsubprocesses)
63
integer ns
64
c-----
65
c Begin Code
66
c-----
67
ns = 1
68
open(unit=15, file=plist,status='old',err=99)
69
do while (.true.)
70
read(15,*,err=999,end=999) subname(ns)
71
ns=ns+1
72
enddo
73
99 subname(ns) = './'
74
write(*,*) "Did not find ", plist
75
return
76
999 ns = ns-1
77
write(*,*) "Found ", ns," subprocesses"
78
close(15)
79
end
80
81
subroutine get_channels(pathname,channame,nchan)
82
c*****************************************************************
83
c Opens file symfact.dat to determine all channels
84
c****************************************************************
85
implicit none
86
c
87
c Constants
88
c
89
character*(*) symfile
90
parameter (symfile='symfact.dat')
91
include 'maxparticles.inc'
92
integer maxsubprocesses
93
parameter (maxsubprocesses=999)
94
c
95
c Arguments
96
c
97
character*(*) channame(maxsubprocesses)
98
character*(*) pathname
99
integer nchan
100
c
101
c Local
102
c
103
character*300 fname
104
character*50 dirname
105
integer jc,ip
106
double precision xi
107
integer j,k
108
integer ncode,npos
109
character*20 formstr
110
c-----
111
c Begin Code
112
c-----
113
jc = index(pathname," ")
114
c ncode is number of digits needed for the bw coding
115
ncode=int(dlog10(3d0)*(max_particles-3))+1
116
fname = pathname(1:jc-1) // "/" // symfile
117
nchan = 0
118
open(unit=35, file=fname,status='old',err=99)
119
do while (.true.)
120
read(35,*,err=99,end=99) xi,j
121
if (j .gt. 0) then
122
k = int(xi*(1+10**(-ncode)))
123
npos=int(dlog10(dble(k)))+1
124
if ( (xi-k) .eq. 0) then
125
c Write with correct number of digits
126
write(formstr,'(a,i1,a)') '(a,i',npos,',a)'
127
write(dirname, formstr) 'G',k,'/'
128
else !Handle B.W.
129
c Write with correct number of digits
130
write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
131
$ '.',ncode,',a)'
132
write(dirname,formstr) 'G',xi,'/'
133
endif
134
ip = index(dirname,'/')
135
nchan=nchan+1
136
channame(nchan) = fname(1:jc-1)// "/" //dirname(1:ip)
137
endif
138
98 enddo
139
99 close(35)
140
write(*,*) "Found ", nchan, " channels."
141
end
142
143
subroutine sum_multichannel(channame)
144
c*****************************************************************
145
c Looks in channel to see if there are multiple runs that
146
c need to be combined. If so combines them into single run
147
c****************************************************************
148
implicit none
149
c
150
c Constants
151
c
152
character*(*) mfile
153
parameter (mfile='multijob.dat')
154
integer maxjobs
155
parameter (maxjobs=26)
156
c
157
c Arguments
158
c
159
character*(*) channame
160
c
161
c Local
162
c
163
character*300 fname
164
character*300 pathname
165
character*26 alphabet
166
integer jc,i
167
integer njobs
168
double precision xsec(maxjobs),xerr(maxjobs),rxsec(maxjobs)
169
double precision xsec_it(9,maxjobs),xerr_it(9,maxjobs)
170
double precision rxsec_it(9,maxjobs)
171
double precision eff(9,maxjobs),wmax(9,maxjobs)
172
integer nevents(maxjobs),ntry(maxjobs)
173
double precision xsec_tot, rxsec_tot, xerr_tot
174
integer tot_nevents, tot_ntry
175
integer lunw,nw
176
double precision wgt
177
c-----
178
c Begin Code
179
c-----
180
alphabet="abcdefghijklmnopqrstuvwxyz"
181
jc = index(channame," ")
182
fname = channame(1:jc-1) // mfile
183
njobs = 0
184
open(unit=35, file=fname,status='old',err=99)
185
read(35,*,err=99,end=99) njobs
186
99 close(35)
187
if(njobs .gt. 0) then
188
write(*,*) "Found ", njobs, " jobs in ",fname
189
c
190
c Read in results from each of the jobs
191
c
192
do i = 1, njobs
193
pathname = channame(1:jc-2) // alphabet(i:i)
194
call get_results(pathname,xsec(i),rxsec(i),xerr(i),ntry(i),nevents(i),
195
$ xsec_it(1,i),xerr_it(1,i),rxsec_it(1,i),eff(1,i),wmax(1,i))
196
enddo
197
c
198
c Process results from each job to get final statistics
199
c
200
xsec_tot = 0d0
201
rxsec_tot = 0d0
202
xerr_tot = 0d0
203
tot_ntry = 0
204
tot_nevents = 0
205
do i = 1, njobs
206
xsec_tot = xsec_tot+xsec(i)
207
rxsec_tot = rxsec_tot+rxsec(i)
208
xerr_tot = xerr_tot+xerr(i)**2
209
tot_ntry = tot_ntry+ntry(i)
210
tot_nevents = tot_nevents+nevents(i)
211
enddo
212
xsec_tot = xsec_tot / njobs
213
rxsec_tot = rxsec_tot / njobs
214
xerr_tot = sqrt(xerr_tot)/njobs
215
216
pathname = channame(1:jc-2)
217
i=1
218
call put_results(pathname,xsec_tot,rxsec_tot,xerr_tot,tot_ntry,
219
$ tot_nevents,xsec_it(1,i),rxsec_it(1,i),xerr_it(1,i),eff(1,i),
220
$ wmax(1,i))
221
call write_logfile(pathname,xsec,rxsec,xerr,ntry,nevents,
222
$ xsec_it,rxsec_it,xerr_it,eff,wmax,njobs)
223
c
224
c Now read in all of the events and write them
225
c back out with the appropriate scaled weight
226
c
227
lunw=15
228
jc = index(pathname," ")
229
open(unit=lunw,file=pathname(1:jc-1)//"/"//"events.lhe",
230
$ status="unknown",err=999)
231
write(*,*) "Placing combined events in ",
232
$ pathname(1:jc-1)//"/"//"events.lhe"
233
wgt = xsec_tot / tot_nevents
234
tot_nevents=0
235
do i = 1, njobs
236
jc = index(channame," ")
237
pathname = channame(1:jc-2) // alphabet(i:i)
238
call copy_events(pathname,lunw,wgt,nw)
239
tot_nevents=tot_nevents+nw
240
write(*,*) "Added ",nw," events from run ", alphabet(i:i)
241
if (nw .ne. nevents(i)) then
242
write(*,*) "Error writing events ",i,nw,nevents(i),
243
$ pathname
244
else
245
endif
246
enddo
247
write(*,*) "Combined ",tot_nevents," to ",channame(1:jc-2)
248
close(lunw)
249
endif
250
return
251
999 close(lunw)
252
write(*,*) "Error, unable to open events.lhe file for output"
253
$ ,pathname
254
end
255
256
subroutine copy_events(pathname,lunw,wgt,nw)
257
c*********************************************************************
258
c Copy events from separate runs into one file w/ appropriate wgts
259
c*********************************************************************
260
implicit none
261
c
262
c Constants
263
c
264
character*(*) eventfname
265
parameter (eventfname="events.lhe")
266
include 'maxparticles.inc'
267
integer maxexternal !Max number external momenta
268
parameter (maxexternal=2*max_particles-3)
269
c
270
c Arguments
271
c
272
character*(*) pathname
273
double precision wgt
274
integer lunw,nw
275
c
276
c Local
277
c
278
character*300 fname
279
integer jc
280
double precision P(0:4,maxexternal),xwgt
281
integer n,ic(7,maxexternal),ievent
282
double precision scale,aqcd,aqed
283
character*300 buff
284
logical done
285
c-----
286
c Begin Code
287
c-----
288
jc = index(pathname," ")
289
fname = pathname(1:jc-1) // "/" // eventfname
290
nw = 0
291
write(*,*) "Copy from ",fname(1:jc+10)
292
open(unit=35, file=fname,status='old',err=999)
293
done = .false.
294
do while (.not. done .and. nw < 999999)
295
done = .true.
296
call read_event(35,P,xwgt,n,ic,ievent,scale,
297
$ aqcd,aqed,buff,done)
298
if (.not. done) then
299
wgt = dsign(wgt,xwgt)
300
call write_event(lunw,P,wgt,n,ic,ievent,scale,
301
$ aqcd,aqed,buff)
302
nw=nw+1
303
endif
304
enddo
305
99 close(35)
306
return
307
999 close(35)
308
write(*,*) "Error unable to open ",fname
309
end
310
311
subroutine get_results(pathname,xsec,rxsec,xerr,ntry,nevents,
312
$ xsec_it,xerr_it,rxsec_it,eff,wmax)
313
c*********************************************************************
314
c Read run results from file results.dat
315
c*********************************************************************
316
implicit none
317
c
318
c Constants
319
c
320
character*(*) resultfname
321
parameter (resultfname="results.dat")
322
c
323
c Arguments
324
c
325
character*(*) pathname
326
double precision xsec,rxsec,xerr
327
double precision xsec_it(9),rxsec_it(9),xerr_it(9)
328
double precision eff(9),wmax(9)
329
integer nevents,ntry
330
c
331
c Local
332
c
333
character*300 fname
334
integer jc,i,j
335
double precision x1,x2,x3,x4
336
c-----
337
c Begin Code
338
c-----
339
jc = index(pathname," ")
340
fname = pathname(1:jc-1) // "/" // resultfname
341
nevents = 0
342
i=1
343
open(unit=35, file=fname,status='old',err=99)
344
read(35,*,err=99,end=99) xsec,xerr, x1, ntry,x3,x4,nevents,x1,x1,rxsec
345
do while (.true.)
346
read(35,*,end=99,err=99) x1,xsec_it(i),xerr_it(i),
347
$ eff(i),wmax(i),rxsec_it(i)
348
i=i+1
349
enddo
350
99 close(35)
351
i=i-1
352
xsec_it(i+1) = -1d0
353
if (nevents .gt. 0) then
354
write(*,*) "Found ",nevents, " events in ", pathname(1:jc-1)
355
i=1
356
do while (xsec_it(i) .gt. 0d0)
357
write(*,*) i,rxsec_it(i),xsec_it(i),xerr_it(i),eff(i),wmax(i)
358
i=i+1
359
enddo
360
endif
361
end
362
363
subroutine write_logfile(pathname,xsec,rxsec,xerr,ntry,nevents,
364
$ xsec_it,rxsec_it,xerr_it,eff,wmax,njobs)
365
c*********************************************************************
366
c Read run results from file results.dat
367
c*********************************************************************
368
implicit none
369
c
370
c Constants
371
c
372
integer maxjobs
373
parameter (maxjobs=26)
374
character*(*) logfname
375
parameter (logfname="log.txt")
376
c
377
c Arguments
378
c
379
character*(*) pathname
380
double precision xsec(maxjobs),rxsec(maxjobs),xerr(maxjobs)
381
double precision xsec_it(9,maxjobs),xerr_it(9,maxjobs)
382
double precision rxsec_it(9,maxjobs)
383
double precision eff(9,maxjobs),wmax(9,maxjobs)
384
integer nevents(maxjobs),ntry(maxjobs)
385
double precision xsec_tot, rxsec_tot, xerr_tot
386
integer tot_nevents, tot_ntry
387
integer njobs
388
c
389
c Local
390
c
391
character*300 fname
392
character*26 alphabet
393
integer jc,i,j
394
c-----
395
c Begin Code
396
c-----
397
alphabet="abcdefghijklmnopqrstuvwxyz"
398
jc = index(pathname," ")
399
fname = pathname(1:jc-1) // "/" // logfname
400
open(unit=35,file=fname,status="old",access="append",err=999)
401
402
write(35,*) '--------------------- Multi run with ',njobs,
403
$ ' jobs. ','---------------------'
404
do i=1,njobs
405
write(35,*) "job ",alphabet(i:i),":",rxsec(i),xsec(i),nevents(i)
406
enddo
407
close(35)
408
return
409
999 write(*,*) "Error opening ",fname(1:jc+10), "for append"
410
close(35)
411
end
412
413
subroutine put_results(pathname,xsec,rxsec,xerr,ntry,nevents,
414
$ xsec_it,rxsec_it,xerr_it,eff,wmax)
415
c*********************************************************************
416
c Read run results from file results.dat
417
c*********************************************************************
418
implicit none
419
c
420
c Constants
421
c
422
character*(*) resultfname
423
parameter (resultfname="results.dat")
424
c
425
c Arguments
426
c
427
character*(*) pathname
428
double precision xsec,rxsec,xerr
429
double precision xsec_it(9),rxsec_it(9),xerr_it(9)
430
double precision eff(9),wmax(9)
431
integer nevents,ntry
432
c
433
c Local
434
c
435
character*300 fname
436
integer jc,i,j
437
double precision x1,x2
438
integer i3,i4
439
c-----
440
c Begin Code
441
c-----
442
jc = index(pathname," ")
443
fname = pathname(1:jc-1) // "/" // resultfname
444
x1 = 0d0
445
x2 = 0d0
446
i3 = 0
447
i4 = 0
448
open(unit=35, file=fname,status='unknown',err=99)
449
write(35,'(3e12.5,2i9,i5,i9,e10.3,e11.3,e13.5)') xsec,xerr, x1, ntry,i3,i4
450
$ ,nevents,(1.00*nevents)/xsec,0d0,rxsec
451
i=1
452
do while (xsec_it(i) .gt. 0 .and. i .lt. 10)
453
write(35,'(i4,5e15.5)') i,xsec_it(i),xerr_it(i),
454
$ eff(i),wmax(i),rxsec_it(i)
455
i=i+1
456
enddo
457
99 close(35)
458
end
Loggerhead is a web-based interface for Breezy
Version: 2.0.1