← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_129738/PROC_129738/bin/internal/merge.pl

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_129738/PROC_129738/bin/internal/merge.pl
 
1
#!/usr/bin/perl -w
 
2
 
 
3
################################################################################
 
4
# merge.pl
 
5
#  Richard Corke (richard.corke@thep.lu.se)
 
6
#
 
7
# Changed structure to not read entire files into memory
 
8
# 
 
9
# Based on merge.pl v1.2 by Michel Herquet (UCL-CP3)
 
10
# DESCRIPTION : script to merge to LHE events files
 
11
# USAGE :
 
12
# ./merge.pl eventfile1.lhe.gz eventfile2.lhe.gz ... result.lhe.gz banner.txt
 
13
# Note that result can be eventfile1, eventfile2, ...
 
14
# OUTPUT : merged file, banner
 
15
################################################################################
 
16
 
 
17
use Compress::Zlib;
 
18
use List::Util qw[min max];
 
19
 
 
20
 
 
21
###########################################################################
 
22
# Initialisation
 
23
###########################################################################
 
24
# Set tag names for later
 
25
my $begin_header = '<header>';
 
26
my $end_header   = '</header>';
 
27
my $begin_event  = '<event>';
 
28
my $end_event    = '</event>';
 
29
my $begin_init   = '<init>';
 
30
my $end_init     = '</init>';
 
31
 
 
32
# Parse the command line arguments
 
33
if ( $#ARGV < 2 ) {
 
34
  die "This script must be called with at least three filenames as arguments!\n";
 
35
}
 
36
my $bannerfile = pop(@ARGV);
 
37
my $outfile = pop(@ARGV);
 
38
 
 
39
 
 
40
###########################################################################
 
41
# First pass - extract number of events and cross-sections
 
42
###########################################################################
 
43
foreach $infile (@ARGV) {
 
44
  $gzin = gzopen($infile, "r") || die ("Couldn't open file $infile\n");
 
45
 
 
46
  # No. events and cross-section from current file
 
47
  $noevents = $xsec = -1;
 
48
 
 
49
  # Storage for current file's init block
 
50
  @gzinit = ();
 
51
 
 
52
  # Keep track if we are in the init block or not
 
53
  $initblock = 0;
 
54
 
 
55
  while (1) {
 
56
    $gzbytes = $gzin->gzreadline($gzline);
 
57
    if ($gzbytes == -1) { die("Error reading from file $infile\n"); }
 
58
    if ($gzbytes == 0)  { last; }
 
59
 
 
60
    # Extract <MGGenerationInfo> information
 
61
    if ($initblock == 0) {
 
62
      if ($gzline =~ s/#  Number of Events\s*:\s*(.*)\n/$1/) { $noevents = $gzline; }
 
63
      if ($gzline =~ s/#  Integrated weight \(pb\)\s*:\s*(.*)\n/$1/) { $xsec = $gzline; }
 
64
 
 
65
      # Check if we enter <init> block
 
66
      if ($gzline =~ m/$begin_init/) { $initblock++; next; }
 
67
 
 
68
    # Extract <init> information
 
69
    } else {
 
70
 
 
71
      # Check for end of <init> block
 
72
      if ($gzline =~ m/$end_init/) { last; }
 
73
 
 
74
      # Remove leading whitespace and split
 
75
      $gzline  =~ s/^\s+//;
 
76
      @gzparam = split(/\s+/, $gzline);
 
77
 
 
78
      # Check <init> block line has right no. of entries
 
79
      if ($initblock == 1) {
 
80
        if ($#gzparam != 9) { die "Not right number of params in init"; }
 
81
      } else {
 
82
        if ($#gzparam != 3) { die "Not right number of params in init"; }
 
83
      }
 
84
 
 
85
      push(@gzinit, [ @gzparam ]);
 
86
      $initblock++;
 
87
 
 
88
    }
 
89
  }
 
90
  $gzin->gzclose();
 
91
 
 
92
  # Check the file contained all the information we need
 
93
  if ($noevents == -1 || $xsec == -1) {
 
94
    die("Couldn't extract No. Events / Cross-section from $infile")
 
95
  }
 
96
 
 
97
  # Store information for later
 
98
  push(@infiles, [ $infile, $noevents, $xsec, [ @gzinit ] ]);
 
99
}
 
100
 
 
101
 
 
102
###########################################################################
 
103
# Check/compute cross-sections, totals and unit weight
 
104
###########################################################################
 
105
$oldxsec        = $infiles[0][2];
 
106
@oldinit        = @{$infiles[0][3]};
 
107
$totevents = 0;  $totxsec = 0.0;
 
108
foreach $infile (@infiles) {
 
109
  print "Input file: $infile->[0]\n";
 
110
  print " No. Events = $infile->[1], Cross-section = $infile->[2]\n";
 
111
 
 
112
  # Check that cross sections do not differ too much
 
113
  $newxsec = $infile->[2];
 
114
  if (abs(($newxsec - $oldxsec) / $newxsec) > 0.05 ) {
 
115
    print " WARNING Cross sections do not agree with a 5\% precision!\n";
 
116
  }
 
117
  $oldxsec = $newxsec;
 
118
 
 
119
  @currinit = @{$infile->[3]};
 
120
  # Same number of lines in <init> block?
 
121
  if ($#oldinit != $#currinit)             { die("Init blocks do not match"); }
 
122
 
 
123
  # Same number of entries on first line of <init> block?
 
124
  if ($#{$oldinit[0]} != $#{$currinit[0]}) { die("Init blocks do not match"); }
 
125
 
 
126
  # All entries the same on first line of <init> block?
 
127
  for ($i = 0; $i <= $#{$oldinit[0]}; $i++) {
 
128
    if ($oldinit[0][$i] != $currinit[0][$i])
 
129
      { die("Init blocks do not match"); }
 
130
  }
 
131
 
 
132
  # Create new init block (overwrite first file's init block data)
 
133
  for ($i = 1; $i <= $#oldinit; $i++) {
 
134
    if ($oldinit[$i][3] != $currinit[$i][3]) { die("Init blocks do not match"); }
 
135
 
 
136
    print " xsecup = $currinit[$i][0], xerrup = $currinit[$i][1]\n";
 
137
    print " xmaxup = $currinit[$i][2], lprup = $currinit[$i][3]\n";
 
138
 
 
139
    # XSECUP = sum(xsecup * no.events) / tot.events
 
140
    # XERRUP = sqrt( sum(sigma^2 * no.events^2) ) / tot.events
 
141
 
 
142
    # Here we temporarily store:
 
143
    #  sum(xsecup * no.events)
 
144
    #  sum(sigma^2 * no.events^2)
 
145
    if (\$oldinit == \$currinit) {
 
146
      $oldinit[$i][0] *= $infile->[1];
 
147
      $oldinit[$i][1] *= $oldinit[$i][1] * $infile->[1]**2;
 
148
 
 
149
    } else {
 
150
      $oldinit[$i][0] += ($currinit[$i][0] * $infile->[1]);
 
151
      $oldinit[$i][1] += $currinit[$i][1]**2 * $infile->[1]**2;
 
152
 
 
153
    }
 
154
 
 
155
    # XMAXUP = max(xmaxup)
 
156
    $oldinit[$i][2] = max($oldinit[$i][2], $currinit[$i][2]);
 
157
  }
 
158
 
 
159
  # Total number of events and total cross-section
 
160
  $totevents += $infile->[1];
 
161
  $totxsec += ($infile->[1] * $infile->[2]);
 
162
 
 
163
  print "\n";
 
164
}
 
165
print "\n";
 
166
 
 
167
# Finish calculation of XSECUP and XERRUP
 
168
for ($i = 1; $i <= $#oldinit; $i++) {
 
169
  $oldinit[$i][0] /= $totevents;
 
170
  $oldinit[$i][0] = sprintf('%0.5E', $oldinit[$i][0]);
 
171
 
 
172
  $oldinit[$i][1] = sqrt($oldinit[$i][1]);
 
173
  $oldinit[$i][1] /= $totevents;
 
174
  $oldinit[$i][1] = sprintf('%0.5E', $oldinit[$i][1]);
 
175
}
 
176
 
 
177
# Finish calculation of total xsec and new unit weight
 
178
$totxsec /= $totevents;
 
179
$dispxsec = sprintf('%0.5E', $totxsec);
 
180
$uwgt = sprintf('%0.5E', $totxsec / $totevents);
 
181
 
 
182
# Display new information
 
183
print "Banner file: $bannerfile\n";
 
184
print "Output file: $outfile\n";
 
185
print " No. Events = $totevents, Cross-section = $dispxsec\n";
 
186
for ($i = 1; $i <= $#oldinit; $i++) {
 
187
  print " xsecup = $oldinit[$i][0], xerrup = $oldinit[$i][1]\n";
 
188
  print " xmaxup = $oldinit[$i][2], lprup = $oldinit[$i][3]\n";
 
189
}
 
190
print "\n";
 
191
print " Unit weight = $uwgt\n";
 
192
 
 
193
 
 
194
###########################################################################
 
195
# Second pass - output file with new information
 
196
###########################################################################
 
197
# The first file is written out with changes to the header, init block
 
198
# and unit weight of each event.
 
199
# All events from other files are then written out with their unit weight
 
200
# changed.
 
201
###########################################################################
 
202
 
 
203
$gzout    = gzopen($outfile, "w") || die ("Couldn't open file $outfile\n");
 
204
open(BANNER, ">$bannerfile")      || die ("Couldn't open file $outfile\n");
 
205
 
 
206
$stage = 0;
 
207
$eventcount = 0;
 
208
foreach $infile (@infiles) {
 
209
  $gzin = gzopen($infile->[0], "r") || die ("Couldn't open file $infile\n");
 
210
 
 
211
  while (1) {
 
212
    $gzbytes = $gzin->gzreadline($gzline);
 
213
    if ($gzbytes == -1) { die("Error reading from file $infile\n"); }
 
214
    if ($gzbytes == 0)  { last; }
 
215
 
 
216
    # Pre-header
 
217
    if ($stage == 0) {
 
218
      if ($gzline =~ m/$begin_header/) { $stage++; }
 
219
 
 
220
    # Header (and output to banner)
 
221
    } elsif ($stage == 1) {
 
222
      $gzline =~ s/#  Integrated weight \(pb\)\s*:(.*)\n/#  Integrated weight (pb)  :  $dispxsec\n/;
 
223
      $gzline =~ s/#  Number of Events\s*:(.*)\n/#  Number of Events        :  $totevents\n/;
 
224
      $gzline =~ s/#  Unit wgt\s*:(.*)\n/#  Unit wgt                :  $uwgt\n/;
 
225
 
 
226
      if ($gzline =~ m/$end_header/)   { $stage++; }
 
227
      else { print BANNER $gzline; }
 
228
 
 
229
    # Init block
 
230
    } elsif ($stage == 2) {
 
231
      if ($gzline =~ m/$end_init/) {
 
232
        $gzline = "<init>\n";
 
233
 
 
234
        for ($i = 0; $i <= $#oldinit; $i++) {
 
235
          $gzline .= "  ";
 
236
          for ($j = 0; $j <= $#{$oldinit[$i]}; $j++) {
 
237
            $gzline .= "$oldinit[$i][$j] ";
 
238
          }
 
239
          $gzline .= "\n";
 
240
        }
 
241
 
 
242
        $gzline .= "</init>\n";
 
243
        
 
244
        $stage++;
 
245
      } else { next; }
 
246
 
 
247
    # Pre-event
 
248
    } elsif ($stage == 3) {
 
249
      if ($gzline =~ m/$begin_event/)  { $stage++; } else { next; }
 
250
 
 
251
    # Event information
 
252
    } elsif ($stage == 4) {
 
253
      $gzline  =~ s/^\s+//;
 
254
      @gzparam = split(/\s+/, $gzline);
 
255
      if ($#gzparam != 5) { die "Not right number of param in first line of event"; }
 
256
      $gzline = " $gzparam[0] $gzparam[1] $uwgt $gzparam[3] $gzparam[4] $gzparam[5]\n";
 
257
 
 
258
      $stage++;
 
259
 
 
260
    # Event particles
 
261
    } elsif ($stage == 5) {
 
262
      if ($gzline =~ m/$end_event/)    { $stage = 3; $eventcount++ }
 
263
 
 
264
    }
 
265
 
 
266
    # Write out the line
 
267
    $gzout->gzwrite($gzline);
 
268
  }
 
269
}
 
270
 
 
271
# Write out closing tag and close
 
272
$gzout->gzwrite("</LesHouchesEvents>\n");
 
273
$gzout->gzclose();
 
274
 
 
275
print "Wrote $eventcount events\n";
 
276
exit(0);
 
277