← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_427003/PROC_427003/bin/internal/Gridpack/replace.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_427003/PROC_427003/bin/internal/Gridpack/replace.pl
 
1
#!/usr/bin/perl -w
 
2
 
 
3
################################################################################
 
4
# replace.pl
 
5
#  Johan Alwall (alwall@slac.stanford.edu)
 
6
#
 
7
# DESCRIPTION : script to replace particle codes in MG/ME event files
 
8
# USAGE :
 
9
# replace.pl infile.lhe outfile.lhe
 
10
# Asks for list of particles to replace. To replace electrons and neutrinos
 
11
# with electrons, muons and taus, give
 
12
# 11: 11 13 15
 
13
# -12: -12 -14 -16
 
14
# done
 
15
# To replace electrons and positrons in W-W+ pair production please run
 
16
# script twice, ones for the e-/ve~, once for e+/ve.
 
17
# Automatically increases cross section with correct factor.
 
18
################################################################################
 
19
 
 
20
use Compress::Zlib;
 
21
use List::Util qw[min max];
 
22
 
 
23
print "Running script replace.pl to replace particle codes in event file.\n";
 
24
 
 
25
###########################################################################
 
26
# Initialisation
 
27
###########################################################################
 
28
# Set tag names for later
 
29
my $begin_header = '<header>';
 
30
my $end_header   = '</header>';
 
31
my $begin_event  = '<event>';
 
32
my $end_event    = '</event>';
 
33
my $begin_init   = '<init>';
 
34
my $end_init     = '</init>';
 
35
 
 
36
# Parse the command line arguments
 
37
if ( $#ARGV != 1 ) {
 
38
  die "Usage: replace.pl infile.lhe outfile.lhe\n";
 
39
}
 
40
my $outfile = pop(@ARGV);
 
41
my $infile = pop(@ARGV);
 
42
 
 
43
open INFILE, "<$infile" or die ("Error: Couldn't open file $infile\n");
 
44
 
 
45
###########################################################################
 
46
# Read which particles to replace
 
47
###########################################################################
 
48
 
 
49
print "Enter particles to be replaced in syntax PID : PID1 PID2 PID3 ...\n";
 
50
print "(note that they will be replaced together, so must be same number of\n";
 
51
print " particles in each line)\n";
 
52
print "Enter \"done\" or <return> when ready\n";
 
53
$line=<STDIN>;
 
54
my %particles;
 
55
my $nparts=0;
 
56
my $i=0;
 
57
 
 
58
while ($line !~ m/^\s*done\s*$/i && $line !~ m/^\s*$/){
 
59
    # Remove leading whitespace
 
60
    $line  =~ s/^\s+//;
 
61
    if($line =~ m/^(-?\d+)\s*\:\s*((-?\d+\s+)+)/){
 
62
        @tmp = split (/\s+/,$2);
 
63
        if($i == 0) { $nparts = $#tmp+1; }
 
64
        elsif($#tmp+1 != $nparts){
 
65
            die "Error: Not same number of particles in each replace: ",$#tmp+1," $nparts\n";
 
66
        }
 
67
        $i = $1;
 
68
        $particles{$i}=[ @tmp ];
 
69
        print "Replacing $i with @{$particles{$i}}\n";
 
70
    }
 
71
    elsif($line !~ m/^#/){
 
72
        print "Please use syntax PID1 : PID2 PID3 ... or \"done\" or <return> when ready\n";
 
73
    }
 
74
    $line=<STDIN>;
 
75
}
 
76
@keys = keys %particles;
 
77
if($#keys < 0) {die "Error: No particles to replace\n";}
 
78
print "Multiplying cross section and weight by $nparts\n";
 
79
 
 
80
###########################################################################
 
81
# Go through file and perform replacements
 
82
###########################################################################
 
83
 
 
84
open OUTFILE, ">$outfile" or die ("Error: Couldn't open file $outfile for writing\n");
 
85
 
 
86
# No. events and cross-section numbers file
 
87
$nevents = $xsec = $trunc = $unitwgt = -1;
 
88
 
 
89
# Keep track in which block we are
 
90
$initblock = 0;
 
91
$headerblock = 0;
 
92
$eventblock = 0;
 
93
$eventcount = 0;
 
94
$rdnseed = 0;
 
95
 
 
96
while (my $line = <INFILE>) {
 
97
 
 
98
    # Extract <MGGenerationInfo> information
 
99
    if ($line =~ m/$end_header/) { 
 
100
        if($nevents == -1 || $xsec == -1 || $trunc == -1 || $unitwgt == -1) {
 
101
            die "Error: Couldn't find number of events and cross section in $infile.\n";
 
102
        }
 
103
        $headerblock = 0; 
 
104
        print OUTFILE "<ReplaceParticleInfo>\n";
 
105
        printf OUTFILE "#  Number of Events        : %11i\n",$nevents;
 
106
        printf OUTFILE "#  Integrated weight (pb)  : %11.4E\n",$xsec;
 
107
        printf OUTFILE "#  Truncated wgt (pb)      : %11.4E\n",$trunc;
 
108
        printf OUTFILE "#  Unit wgt                : %11.4E\n",$unitwgt;
 
109
        if($rdnseed > 0){
 
110
            print OUTFILE " ", $rdnseed+1, "  = gseed ! Random seed for next iteration of replace.pl\n";}
 
111
        print OUTFILE "</ReplaceParticleInfo>\n";
 
112
        if($rdnseed > 0) {print "Initialize random seed with $rdnseed\n";srand($rdnseed);}
 
113
        else {print "Warning: Random seed 0, use default random seed (unreproducible)\n";}
 
114
    }
 
115
    if ($line =~ m/$end_init/) { $initblock=0; }
 
116
    if ($line =~ m/$end_event/) { 
 
117
        $eventcount++;
 
118
        $eventblock=0; 
 
119
    }
 
120
    
 
121
    if ($headerblock == 1) {
 
122
        # In header
 
123
        if ($line =~ m/^#\s+Number of Events\s*:\s*(.*)\n/) { $nevents = $1; }
 
124
        if ($line =~ m/^#\s+Integrated weight \(pb\)\s*:\s*(.*)\n/) { $xsec = $1*$nparts; }
 
125
        if ($line =~ m/^#\s+Truncated wgt \(pb\)\s*:\s*(.*)\n/) { $trunc = $1*$nparts; }
 
126
        if ($line =~ m/^#\s+Unit wgt\s*:\s*(.*)\n/) { $unitwgt = $1*$nparts; }
 
127
        if ($line =~ m/^\s*(\d+)\s*=\s*gseed/) { $rdnseed = $1; }
 
128
    } elsif ($initblock > 0) {
 
129
        # In <init> block
 
130
 
 
131
        if($initblock > 1) {
 
132
        # Remove leading whitespace and split
 
133
            $line  =~ s/^\s+//;
 
134
            @param = split(/\s+/, $line);
 
135
            if ($#param != 3) { die "Error: Wrong number of params in init ($#param)"; }
 
136
            $param[0]*=$nparts;
 
137
            $param[1]*=$nparts;
 
138
            $param[2]*=$nparts;
 
139
            $line = sprintf " %18.11E %18.11E %18.11E %3i\n",$param[0],$param[1],$param[2],$param[3];
 
140
        }
 
141
        $initblock++;
 
142
 
 
143
    } elsif ($eventblock > 0) {
 
144
        # In <event> block
 
145
        # Remove leading whitespace and split
 
146
        $line  =~ m/^\s*(.*)/;
 
147
        if($line=~/^\#/){
 
148
        #    print OUTFILE $line;
 
149
        }
 
150
        else{
 
151
            @param = split(/\s+/, $1);
 
152
            if($eventblock == 1){
 
153
                if ($#param != 5) { die "Error: Wrong number of params in event $eventcount \($#param / 5\)"; }
 
154
                $param[2]*=$nparts;
 
155
                $line = sprintf "%2i %3i %13.6E %13.6E %13.6E %13.6E\n",$param[0],$param[1],$param[2],$param[3],$param[4],$param[5];
 
156
                
 
157
                # Randomly choose a particle from the lists
 
158
                $rnumber = int(rand($nparts));
 
159
            }
 
160
            else {
 
161
#           @param = split(/\s+/, $1);
 
162
                if ($#param != 12) { die "Error: Wrong number of params in event $eventcount \($#param / 12\)"; }
 
163
                if($particles{$param[0]}){
 
164
                    $line = sprintf "%9i %4i %4i %4i %4i %4i %18.10E %18.10E %18.10E %18.10E %18.10E %1.0f. %2.0f.\n",
 
165
                    $particles{$param[0]}[$rnumber],$param[1],$param[2],$param[3],$param[4],$param[5],$param[6],$param[7],$param[8],$param[9],$param[10],$param[11],$param[12];
 
166
                }
 
167
            }
 
168
        }
 
169
        $eventblock++;
 
170
    }
 
171
            
 
172
    if ($line =~ m/$begin_header/) { $headerblock = 1; }
 
173
    if ($line =~ m/$begin_init/) { $initblock=1; }
 
174
    if ($line =~ m/$begin_event/) { $eventblock=1; }
 
175
 
 
176
    print OUTFILE "$line";
 
177
}
 
178
 
 
179
close OUTFILE;
 
180
close INFILE;
 
181
 
 
182
print "Wrote $eventcount events\n";
 
183
if( $eventcount < $nevents ) { print "Warning: $infile ended early\n"; }
 
184
 
 
185
exit(0);
 
186
 
 
187