~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to madgraph/various/combine_runs.py

  • Committer: olivier Mattelaer
  • Date: 2015-03-05 00:14:16 UTC
  • mfrom: (258.1.9 2.3)
  • mto: (258.8.1 2.3)
  • mto: This revision was merged to the branch mainline in revision 259.
  • Revision ID: olivier.mattelaer@uclouvain.be-20150305001416-y9mzeykfzwnl9t0j
partial merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
################################################################################
2
 
#
3
 
# Copyright (c) 2012 The MadGraph5_aMC@NLO Development team and Contributors
4
 
#
5
 
# This file is a part of the MadGraph5_aMC@NLO project, an application which 
6
 
# automatically generates Feynman diagrams and matrix elements for arbitrary
7
 
# high-energy processes in the Standard Model and beyond.
8
 
#
9
 
# It is subject to the MadGraph5_aMC@NLO license which should accompany this 
10
 
# distribution.
11
 
#
12
 
# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
13
 
#
14
 
################################################################################
15
 
"""Program to combine results from channels that have been
16
 
     split into multiple jobs. Multi-job channels are identified
17
 
     by local file mjobs.dat in the channel directory.
18
 
"""
19
 
from __future__ import division
20
 
import math
21
 
import os
22
 
import re
23
 
import logging
24
 
 
25
 
try:
26
 
    import madgraph.various.sum_html as sum_html
27
 
    from madgraph import InvalidCmd, MadGraph5Error, MG5DIR
28
 
except ImportError:
29
 
    import internal.sum_html as sum_html
30
 
    from internal import InvalidCmd, MadGraph5Error
31
 
    
32
 
logger = logging.getLogger('madevent.combine_run') # -> stdout
33
 
 
34
 
#usefull shortcut
35
 
pjoin = os.path.join
36
 
 
37
 
   
38
 
def get_inc_file(path):
39
 
    """read the information of fortran inc files and returns
40
 
       the definition in a dictionary format.
41
 
       This catch PARAMETER (NAME = VALUE)"""
42
 
       
43
 
    pat = re.compile(r'''PARAMETER\s*\((?P<name>[_\w]*)\s*=\s*(?P<value>[\+\-\ded]*)\)''',
44
 
                     re.I)
45
 
        
46
 
    out = {}   
47
 
    for name, value in pat.findall(open(path).read()):
48
 
        orig_value = str(value)
49
 
        try:
50
 
            out[name.lower()] = float(value.replace('d','e'))
51
 
        except ValueError:
52
 
            out[name] = orig_value
53
 
    return out
54
 
 
55
 
class CombineRuns(object):
56
 
    
57
 
    def __init__(self, me_dir, subproc=None):
58
 
        
59
 
        self.me_dir = me_dir
60
 
        
61
 
        if not subproc:
62
 
            subproc = [l.strip() for l in open(pjoin(self.me_dir,'SubProcesses', 
63
 
                                                                 'subproc.mg'))]
64
 
        self.subproc = subproc
65
 
        maxpart = get_inc_file(pjoin(me_dir, 'Source', 'maxparticles.inc'))
66
 
        self.maxparticles = maxpart['max_particles']
67
 
    
68
 
    
69
 
        for procname in self.subproc:
70
 
            path = pjoin(self.me_dir,'SubProcesses', procname)
71
 
            channels = self.get_channels(path)
72
 
            for channel in channels:
73
 
                self.sum_multichannel(channel)
74
 
    
75
 
    def sum_multichannel(self, channel):
76
 
        """Looks in channel to see if there are multiple runs that
77
 
        need to be combined. If so combines them into single run"""
78
 
       
79
 
        alphabet = "abcdefghijklmnopqrstuvwxyz"
80
 
 
81
 
        if os.path.exists(pjoin(channel, 'multijob.dat')):
82
 
            njobs = int(open(pjoin(channel, 'multijob.dat')).read())
83
 
        else:
84
 
            return
85
 
        results = sum_html.Combine_results(channel)
86
 
        if njobs:
87
 
            logger.debug('find %s multijob in %s' % (njobs, channel))
88
 
        else:
89
 
            return
90
 
        for i in range(njobs):
91
 
            if channel.endswith(os.path.pathsep):
92
 
                path = channel[:-1] + alphabet[i] 
93
 
            else:
94
 
                path = channel + alphabet[i] 
95
 
            results.add_results(name=alphabet[i], 
96
 
                                filepath=pjoin(path, 'results.dat'))
97
 
        
98
 
        results.compute_average()
99
 
        if results.xsec:
100
 
            results.write_results_dat(pjoin(channel, 'results.dat'))
101
 
        else:
102
 
            return
103
 
        ### Adding information in the log file
104
 
        fsock = open(pjoin(channel, 'log.txt'), 'a')
105
 
        fsock.write('--------------------- Multi run with %s jobs. ---------------------\n'
106
 
                    % njobs)
107
 
        for r in results:
108
 
            fsock.write('job %s : %s %s %s\n' % (r.name, r.xsec, r.axsec, r.nunwgt))  
109
 
            
110
 
        #Now read in all of the events and write them
111
 
        #back out with the appropriate scaled weight
112
 
        fsock = open(pjoin(channel, 'events.lhe'), 'w')
113
 
        wgt = results.axsec / results.nunwgt
114
 
        tot_nevents=0
115
 
        for result in results:  
116
 
            i = result.name
117
 
            if channel.endswith(os.path.pathsep):
118
 
                path = channel[:-1] + i 
119
 
            else:
120
 
                path = channel + i 
121
 
            nw = self.copy_events(fsock, pjoin(path,'events.lhe'), wgt)
122
 
            #tot_events += nw
123
 
        #logger.debug("Combined %s events to %s " % (tot_events, channel))
124
 
 
125
 
 
126
 
    def copy_events(self, fsock, input, new_wgt):
127
 
        """ Copy events from separate runs into one file w/ appropriate wgts"""
128
 
        
129
 
        def get_fortran_str(nb):
130
 
            data = '%E' % nb
131
 
            nb, power = data.split('E')
132
 
            nb = abs(float(nb)) /10
133
 
            power = int(power) + 1
134
 
            return '%.7fE%+03i' %(nb,power)
135
 
        new_wgt = get_fortran_str(new_wgt)
136
 
        old_line = ""
137
 
        for line in open(input):
138
 
            if old_line.startswith("<event>"):
139
 
                data = line.split()
140
 
                if not len(data) == 6:
141
 
                    raise MadGraph5Error, "Line after <event> should have 6 entries"
142
 
                if float(data[2]) > 0:
143
 
                    sign = ''
144
 
                else:
145
 
                    sign = '-'
146
 
                    
147
 
                line= ' %s  %s%s  %s\n' % ('   '.join(data[:2]), sign,
148
 
                                           new_wgt, '  '.join(data[3:]))
149
 
            fsock.write(line)
150
 
            old_line = line
151
 
           
152
 
    def get_channels(self, proc_path):
153
 
        """Opens file symfact.dat to determine all channels"""
154
 
        sympath = os.path.join(proc_path, 'symfact.dat')
155
 
        
156
 
        #ncode is number of digits needed for the bw coding
157
 
        
158
 
        ncode = int(math.log10(3)*(self.maxparticles-3))+1
159
 
        channels = []
160
 
        for line in open(sympath):
161
 
            try:
162
 
                xi, j = line.split()
163
 
            except Exception:
164
 
                break
165
 
            xi, j  = float(xi), int(j)
166
 
            
167
 
            if j > 0:
168
 
                k = int(xi) 
169
 
                npos = int(math.log10(k))+1
170
 
                #Write with correct number of digits
171
 
                if xi == k:
172
 
                    dirname = 'G%i' % k
173
 
                else:
174
 
                    dirname = 'G%.{0}f'.format(ncode) % xi
175
 
                channels.append(os.path.join(proc_path,dirname))
176
 
        return channels
177
 
    
178
 
        
179
 
        
180