← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/sumres.py

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/sumres.py
5
5
#  Replaces the sumres.f and sumres2.f files
6
6
#  MZ, 2011-10-22
7
7
 
 
8
from __future__ import division
8
9
import math
9
10
import sys
10
11
import random
30
31
processes=[]
31
32
tot=0
32
33
err=0
 
34
totABS=0
 
35
errABS=0
33
36
 
34
37
# open the file containing the list of directories
35
38
file=open("dirs.txt")
37
40
file.close()
38
41
dirs.remove('')
39
42
 
40
 
 
41
 
for line in lines:
42
 
    list = line.split()
 
43
# The syntax of lines should be first the ABS cross section for the
 
44
# channel and the line after that the cross section for the same
 
45
# channel.
 
46
for line in range(0,len(lines),2):
 
47
    list = lines[line].split()
43
48
    if list:
44
49
        proc={}
45
50
        proc['folder'] = list[0].split('/')[0]
46
51
        proc['subproc'] = proc['folder'][0:proc['folder'].rfind('_')]
47
52
        proc['channel'] = list[0].split('/')[1]
48
53
        dirs.remove(os.path.join(proc['folder'], proc['channel']))
49
 
        if list[3] != '[ABS]:':
 
54
        proc['resultABS'] = float(list[4])
 
55
        proc['errorABS'] = float(list[6])
 
56
        proc['err_percABS'] = proc['errorABS']/proc['resultABS']*100.
 
57
        processes.append(proc)
 
58
        totABS+= proc['resultABS']
 
59
        errABS+= math.pow(proc['errorABS'],2)
 
60
        list = lines[line+1].split()
 
61
        if list:
50
62
            proc['result'] = float(list[3])
51
63
            proc['error'] = float(list[5])
52
 
        else:
53
 
            proc['result'] = float(list[4])
54
 
            proc['error'] = float(list[6])
55
 
        proc['err_perc'] = proc['error']/proc['result']*100.
56
 
        processes.append(proc)
57
 
        tot+= proc['result']
58
 
        err+= math.pow(proc['error'],2)
59
 
 
 
64
            proc['err_perc'] = proc['error']/proc['result']*100.
 
65
            tot+= proc['result']
 
66
            err+= math.pow(proc['error'],2)
60
67
if dirs:
61
68
    print "%d jobs did not terminate correctly " % len(dirs)
62
69
    print '\n'.join(dirs)
63
70
 
64
 
processes.sort(key = lambda proc: -proc['error'])
 
71
processes.sort(key = lambda proc: -proc['errorABS'])
65
72
 
66
73
correct = len(processes) == nexpected
67
74
print "Found %d correctly terminated jobs " %len(processes) 
77
84
for proc in processes:
78
85
    content+='%(folder)20s  %(channel)15s   %(result)10.8e    %(error)6.4e       %(err_perc)6.4f%%  \n' %  proc
79
86
 
 
87
content+='\n\nABS cross-section per integration channel:\n'
 
88
for proc in processes:
 
89
    content+='%(folder)20s  %(channel)15s   %(resultABS)10.8e    %(errorABS)6.4e       %(err_percABS)6.4f%%  \n' %  proc
 
90
 
80
91
content+='\n\nCross-section per subprocess:\n'
81
92
#for subpr in sorted(set(subprocs)):
82
93
subprocesses=[]
115
126
    content+=  '%(subproc)20s    %(xsect)10.8e   %(err)6.4e\n' % subpr
116
127
 
117
128
 
118
 
content+='\nTotal: \n                      %10.8e +- %6.4e  (%6.4e%%)\n' %\
119
 
        (tot, math.sqrt(err), math.sqrt(err)/tot *100.)
 
129
content+='\nTotal ABS and \nTotal: \n                      %10.8e +- %6.4e  (%6.4e%%)\n                      %10.8e +- %6.4e  (%6.4e%%)\n' %\
 
130
        (totABS, math.sqrt(errABS), math.sqrt(errABS)/totABS *100.,tot, math.sqrt(err), math.sqrt(err)/tot *100.)
120
131
 
121
132
if not correct:
122
133
    sys.exit('ERROR: not all jobs terminated correctly\n')
143
154
        proc['lhefile'] = os.path.join(proc['folder'], proc['channel'], 'events.lhe')
144
155
        proc['nevents'] = 0
145
156
    while totevts :
146
 
        target = random.random() * tot
 
157
        target = random.random() * totABS
147
158
        crosssum = 0.
148
159
        i = 0
149
160
        while i<len(processes) and crosssum < target:
150
161
            proc = processes[i]
151
 
            crosssum += proc['result']
 
162
            crosssum += proc['resultABS']
152
163
            i += 1            
153
164
        totevts -= 1
154
165
        i -= 1
162
173
 
163
174
    content_evts = ''
164
175
    for proc in processes:
165
 
        content_evts+= ' '+proc['lhefile']+'              %(nevents)10d            %(result)10.8e        1.0 \n' %  proc
 
176
        content_evts+= ' '+proc['lhefile']+'              %(nevents)10d            %(resultABS)10.8e        1.0 \n' %  proc
166
177
        nevts_file = open(os.path.join(proc['folder'], proc['channel'], 'nevts'),'w')
167
178
        nevts_file.write('%10d\n' % proc['nevents'])
168
179
        nevts_file.close()
179
190
        for line in fileinputs:
180
191
            i += 1
181
192
            if i == 2:
182
 
                accuracy=min(math.sqrt(tot/(req_acc2_inv*proc['result'])),0.2)
 
193
                accuracy=min(math.sqrt(totABS/(req_acc2_inv*proc['resultABS'])),0.2)
183
194
                fileinputschannel.write('%10.8e\n' % accuracy)
184
195
            elif i == 8:
185
196
                fileinputschannel.write('1        ! MINT mode\n')
191
202
    evts_file = open('nevents_unweighted', 'w')
192
203
    evts_file.write(content_evts)
193
204
    evts_file.close()
 
205
 
 
206
# if nevents = -1 and req_acc >= 0, we need to determine the required
 
207
# accuracy in each of the channels: this is for fixed order running!
 
208
elif req_acc>=0 and nevents==-1:
 
209
    req_accABS=req_acc*abs(tot)/totABS
 
210
    content_evts = ''
 
211
    for proc in processes:
 
212
        if proc['channel'][0:3] == 'all':
 
213
            fileinputs = open("madin.all")
 
214
        elif proc['channel'][0:4] == 'novB':
 
215
            fileinputs = open("madin.novB")
 
216
        elif proc['channel'][0:4] == 'born':
 
217
            fileinputs = open("madin.born")
 
218
        elif proc['channel'][0:4] == 'grid':
 
219
            fileinputs = open("madin.grid")
 
220
        elif proc['channel'][0:4] == 'viSB':
 
221
            fileinputs = open("madin.viSB")
 
222
        elif proc['channel'][0:4] == 'virt':
 
223
            fileinputs = open("madin.virt")
 
224
        elif proc['channel'][0:4] == 'novi':
 
225
            fileinputs = open("madin.novi")
 
226
        else:
 
227
            sys.exit("ERROR, DONT KNOW WHICH INPUTS TO USE")
 
228
        fileinputschannel = open(os.path.join(proc['folder'], proc['channel'], 'madinM1'),'w')
 
229
        i=0
 
230
        for line in fileinputs:
 
231
            i += 1
 
232
            if i == 2:
 
233
                accuracy=req_accABS*math.sqrt(totABS*proc['resultABS'])
 
234
                fileinputschannel.write('%10.8e\n' % accuracy)
 
235
            elif i == 9:
 
236
                fileinputschannel.write('-1        ! restart from existing grids\n')
 
237
            else:
 
238
                fileinputschannel.write(line)
 
239
        fileinputschannel.close()
 
240
        fileinputs.close()