Merge lp:~maddevelopers/mg5amcnlo/1.4.5 into lp:~madteam/mg5amcnlo/trunk

[Olivier Mattelaer]

mainly cluster support.
This might be useful to limit the problem on the ucl site.
(mainly allow up to 5 minutes of unresponsive condor)
and wait up to 5 minutes if the files results.dat didn't exits.

Cheers,

Olivier

[Johan Alwall (johan-alwall)]

Hello Olivier,

Very good. Comments/questions below:

> This might be useful to limit the problem on the ucl site.
> (mainly allow up to 5 minutes of unresponsive condor)
> and wait up to 5 minutes if the files results.dat didn't exits.
>
> === modified file 'UpdateNotes.txt'
> +1.4.5 (11/04/12)  OM: correct pointless bug #975647 (SLAH convention problem)
> +                      Thanks to Sho Iwamoto

Replace "pointless" with "minor".

> === modified file 'madgraph/various/cluster.py'
> --- madgraph/various/cluster.py 2012-03-20 19:55:40 +0000
> +++ madgraph/various/cluster.py 2012-04-11 23:06:24 +0000
> @@ -31,7 +31,7 @@
>  class NotImplemented(MadGraph5Error):
>     pass
>
> -def multiple_try(nb_try=5, sleep=1):
> +def multiple_try(nb_try=5, sleep=20):

This only affects waiting time if first attempt fails, right?

>     def deco_retry(f):
>         def deco_f_retry(*args, **opt):
> @@ -41,7 +41,7 @@
>                 except KeyboardInterrupt:
>                     raise
>                 except:
> -                    time.sleep(sleep)
> +                    time.sleep(sleep * nb_try)

Is this sleep * number of attempts already made, so that the time increases? That would make sense.

Cheers,
Johan

[lp:~maddevelopers/mg5amcnlo/1.4.5]

212. By Olivier Mattelaer on 2012-04-12

make the waiting time increasing (thanks Johan to have pointed my mistake on that)
and modify the updateNote

213. By Olivier Mattelaer on 2012-04-12

Add waiting for results.dat if not present (up to one minut).
Add the possibility to clean to cluster in case of error. (done by default)
-refactoring the '-f' option to pass it to a global variable in madevent cmd interface

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

>> === modified file 'madgraph/various/cluster.py'
>> --- madgraph/various/cluster.py 2012-03-20 19:55:40 +0000
>> +++ madgraph/various/cluster.py 2012-04-11 23:06:24 +0000
>> @@ -31,7 +31,7 @@
>> class NotImplemented(MadGraph5Error):
>> pass
>>
>> -def multiple_try(nb_try=5, sleep=1):
>> +def multiple_try(nb_try=5, sleep=20):
>
> This only affects waiting time if first attempt fails, right?

Yes.

>> def deco_retry(f):
>> def deco_f_retry(*args, **opt):
>> @@ -41,7 +41,7 @@
>> except KeyboardInterrupt:
>> raise
>> except:
>> - time.sleep(sleep)
>> + time.sleep(sleep * nb_try)
>
> Is this sleep * number of attempts already made, so that the time
> increases? That would make sense.
>

That was the idea, but looks like that I put the wrong variable.
Thanks to have point this.

So now the waiting time will be
20
40
60
80
100

So a total of 5 minuts.

I also note that one of my modification was not done.
It was adding a multiple try to the results reading since it looks
like that some time, they are filesystem problem
and that the file is not yet on disk.

I have also added that we propose the clean the queue is something
goes wrong.
(This is done by default on the cluster)
The reason not always doing it is that this might remove the reason of
the problem.

Cheers,

Olivier

> Cheers,
> Johan
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/1.4.5/+merge/
> 101662
> You proposed lp:~maddevelopers/madgraph5/1.4.5 for merging.

[lp:~maddevelopers/mg5amcnlo/1.4.5]

214. By Johan Alwall on 2012-04-13

Proper reading of cross section and number of events from Pythia log

215. By Johan Alwall on 2012-04-13

Proper reading of cross section and number of events from Pythia log

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

I've fixed the error that you pointed to me.

I've also changed (again) the way to read the xsec for pythia. Just trying to factorize the
step for more clarity. I take an "old" code that I was using for MW. This reads a file line by line starting at the end of the file. Feel free to modify this class if you want.
(if you don't like the regexp line by line stuff for example).

I've also modified the function get_pythia_error
1) I've removed the argument cross which is pointless now (since you pass nacc)
2) I've changed all the call associate and fixed two of them where the argument were not
   corrected (n_acc = n_gen for two of those)

Otherwise, you change the formula for the computation of the pythia error. Why?
This one is obviously wrong.

Cheers,

Olivier

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/1.4.5]

216. By Olivier Mattelaer on 2012-04-15

fix few small bug
+ review Johan modifications

217. By Johan Alwall on 2012-04-16

Increased number of tries for non-zero helicities to 100 in genps.inc

[Johan Alwall (johan-alwall)]

Hello Olivier,

> I've also changed (again) the way to read the xsec for pythia. Just trying to
> factorize the
> step for more clarity. I take an "old" code that I was using for MW. This
> reads a file line by line starting at the end of the file. Feel free to modify
> this class if you want.
> (if you don't like the regexp line by line stuff for example).

Excellent.

> Otherwise, you change the formula for the computation of the pythia error.
> Why?
> This one is obviously wrong.

The other one was wrong when you ask for fewer events, as well as when you have low matching efficiency. For small number of events, the uncertainty should be Poissonian, and not related to the total xsec uncertainty. On the other hand, if we add the two uncertainties for 100% acceptance of events, we clearly get twice the actual uncertainty. So this was the best approximation I could think of, which gives a reasonable result in both limits. Please let me know if you have some other suggestion.

Thanks,
Johan

>
>
> Cheers,
>
> Olivier

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

> The other one was wrong when you ask for fewer events, as well as
> when you have low matching efficiency. For small number of events,
> the uncertainty should be Poissonian, and not related to the total
> xsec uncertainty.

In general you have the following formula:
\sigma_m = \sigma * N_acc / N_gen
this means that the associated error should be (supposing full
correlation)
error_m = error * N_acc / N_gen + \sigma * error(N_acc) / N_gen

If we take poissonian error for N_acc we get the following formula:

error_m = error * N_acc / N_gen + \sigma * sqrt(N_acc) / N_gen

I don't see any reason to drop the first term if N_acc is small.

> On the other hand, if we add the two uncertainties for 100%
> acceptance of events, we clearly get twice the actual uncertainty.

That's assuming that error = \sigma / sqrt(N) which is the error that
we would like to have, but in principle we don't have any relation
between the two.
In principle, the above formula still stands. So we should have this
factor two (or sqrt(2) if we use uncorellated error) if the relation
that you pre-suppose between the error and sigma is valid.
The only solution to not have it, is to not use Poissonian
distribution...
So I think that indeed this factor 2 (or sqrt(2)) is in fact real.

Concerning the point of choosing between correlated error and
uncorrelated error. I'm starting to change my mind in favor of the
second one.

Cheers,

Olivier

On Apr 16, 2012, at 4:00 AM, Johan Alwall wrote:

> Hello Olivier,
>
>> I've also changed (again) the way to read the xsec for pythia. Just
>> trying to
>> factorize the
>> step for more clarity. I take an "old" code that I was using for
>> MW. This
>> reads a file line by line starting at the end of the file. Feel
>> free to modify
>> this class if you want.
>> (if you don't like the regexp line by line stuff for example).
>
> Excellent.
>
>> Otherwise, you change the formula for the computation of the pythia
>> error.
>> Why?
>> This one is obviously wrong.
>
> The other one was wrong when you ask for fewer events, as well as
> when you have low matching efficiency. For small number of events,
> the uncertainty should be Poissonian, and not related to the total
> xsec uncertainty. On the other hand, if we add the two uncertainties
> for 100% acceptance of events, we clearly get twice the actual
> uncertainty. So this was the best approximation I could think of,
> which gives a reasonable result in both limits. Please let me know
> if you have some other suggestion.
>
> Thanks,
> Johan
>
>>
>>
>> Cheers,
>>
>> Olivier
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/1.4.5/+merge/
> 101662
> You proposed lp:~maddevelopers/madgraph5/1.4.5 for merging.

[Johan Alwall (johan-alwall)]

Hi Olivier,

> In general you have the following formula:
> \sigma_m = \sigma * N_acc / N_gen
> this means that the associated error should be (supposing full
> correlation)
> error_m = error * N_acc / N_gen + \sigma * error(N_acc) / N_gen
>
> If we take poissonian error for N_acc we get the following formula:
>
> error_m = error * N_acc / N_gen + \sigma * sqrt(N_acc) / N_gen

This all sounds reasonable, except that it can't be correct. If N_acc is small and N_gen large, say e.g. N_acc = 100 and N_gen = 1M, you get
error_m = (small) + \sigma * 10 / 1M,
so an uncertainty in the cross section of 1e-5. However, if I pick events from 1M events with a probability p=1e-4 (which is what we are doing here), the Poisson uncertainty is 1/sqrt(\mu)=1/sqrt(100) = 10%.

On the other hand, if we have 10,000 events and all pass Pythia, the uncertainty on the cross section can never be different from the uncertainty from MadGraph.

In fact, Rikkert went through the details of this scenario (i.e., the uncertainty on the number of events picked from a finite pool) in a wiki page: http://cp3wks05.fynu.ucl.ac.be/twiki/bin/view/Main/NumberOfEvents

Please refer to the corresponding page in the new wiki. Let me know if you can't find it.

All the best,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

This is the link that you look for.
https://server06.fynu.ucl.ac.be/projects/madgraph/wiki/NumberOfEvents
(all the pages are transfered on this wiki now).
But I don't think that this is the same problem.

> This all sounds reasonable, except that it can't be correct. If
> N_acc is small and N_gen large, say e.g. N_acc = 100 and N_gen = 1M,
> you get
> error_m = (small) + \sigma * 10 / 1M,
> so an uncertainty in the cross section of 1e-5. However, if I pick
> events from 1M events with a probability p=1e-4 (which is what we
> are doing here), the Poisson uncertainty is 1/sqrt(\mu)=1/sqrt(100)
> = 10%.

if you compute the relative error, you have the following formula:

error_m / sigma_m = error/sigma + 1 / sqrt(Nacc)

and indeed you have a relative error of 10%. So this is correct. (I
think you were confuse between sigma and sigma_m)

Cheers,

Olivier

On 16-avr.-12, at 21:34, Johan Alwall wrote:

> Hi Olivier,
>
>> In general you have the following formula:
>> \sigma_m = \sigma * N_acc / N_gen
>> this means that the associated error should be (supposing full
>> correlation)
>> error_m = error * N_acc / N_gen + \sigma * error(N_acc) / N_gen
>>
>> If we take poissonian error for N_acc we get the following formula:
>>
>> error_m = error * N_acc / N_gen + \sigma * sqrt(N_acc) / N_gen
>
> This all sounds reasonable, except that it can't be correct. If
> N_acc is small and N_gen large, say e.g. N_acc = 100 and N_gen = 1M,
> you get
> error_m = (small) + \sigma * 10 / 1M,
> so an uncertainty in the cross section of 1e-5. However, if I pick
> events from 1M events with a probability p=1e-4 (which is what we
> are doing here), the Poisson uncertainty is 1/sqrt(\mu)=1/sqrt(100)
> = 10%.
>
> On the other hand, if we have 10,000 events and all pass Pythia, the
> uncertainty on the cross section can never be different from the
> uncertainty from MadGraph.
>
> In fact, Rikkert went through the details of this scenario (i.e.,
> the uncertainty on the number of events picked from a finite pool)
> in a wiki page: http://cp3wks05.fynu.ucl.ac.be/twiki/bin/view/Main/NumberOfEvents
>
> Please refer to the corresponding page in the new wiki. Let me know
> if you can't find it.
>
> All the best,
> Johan
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/1.4.5/+merge/
> 101662
> You proposed lp:~maddevelopers/madgraph5/1.4.5 for merging.

[Johan Alwall (johan-alwall)]

Hello Olivier,

> if you compute the relative error, you have the following formula:
>
> error_m / sigma_m = error/sigma + 1 / sqrt(Nacc)
>
> and indeed you have a relative error of 10%. So this is correct. (I
> think you were confuse between sigma and sigma_m)

Indeed, I was using sigma_m in my formula, since that's what's read from the Pythia log file. So that part of our formulas are identical. Note that in general sigma_m is NOT sigma*n/N_gen, since n can be chosen using niter=n in the pythia_card. In fact, sigma should only be used to calculate the relative error of the full sample.

> This is the link that you look for.
> https://server06.fynu.ucl.ac.be/projects/madgraph/wiki/NumberOfEvents
> (all the pages are transfered on this wiki now).
> But I don't think that this is the same problem.

It should be quite similar. The relative uncertainty must go from the MadGraph error (when all events are accepted) to 1/sqrt(N) (when few events are accepted). This is exactly what you get for the statistics described there. However, my formula is clearly wrong, since you get too large error when you ask for a subset of the events using the niter flag in the pythia_card. So in fact, both tried events and accepted events are needed to calculate the error (the total number of events should be the tried events rather than Ngen, which is not relevant in the niter < Ngen case).

Johan

[Olivier Mattelaer (olivier-mattelaer)]

On 17-avr.-12, at 00:33, Johan Alwall wrote:

> Hello Olivier,
>
>> if you compute the relative error, you have the following formula:
>>
>> error_m / sigma_m = error/sigma + 1 / sqrt(Nacc)
>>
>> and indeed you have a relative error of 10%. So this is correct. (I
>> think you were confuse between sigma and sigma_m)
>
> Indeed, I was using sigma_m in my formula, since that's what's read
> from the Pythia log file. So that part of our formulas are
> identical. Note that in general sigma_m is NOT sigma*n/N_gen, since
> n can be chosen using niter=n in the pythia_card. In fact, sigma
> should only be used to calculate the relative error of the full
> sample.

Ok this is problematic (if you fix niter). And then the error depends
on how pythia calculates the cross-section.
Do you how how this is done?
is it sigma_m = n / n_gen
where this time n_gen is not fixed anymore (i.e. that's the number of
events required to have n accepted events)?

>> This is the link that you look for.
>> https://server06.fynu.ucl.ac.be/projects/madgraph/wiki/NumberOfEvents
>> (all the pages are transfered on this wiki now).
>> But I don't think that this is the same problem.
>
> It should be quite similar. The relative uncertainty must go from
> the MadGraph error (when all events are accepted) to 1/sqrt(N) (when
> few events are accepted). This is exactly what you get for the
> statistics described there. However, my formula is clearly wrong,
> since you get too large error when you ask for a subset of the
> events using the niter flag in the pythia_card. So in fact, both
> tried events and accepted events are needed to calculate the error
> (the total number of events should be the tried events rather than
> Ngen, which is not relevant in the niter < Ngen case).

In fact the real point, is that we don't have a poisson distribution
but a binomial distributions (which is the correct distribution for
sucess/failure experiment):
(from wikipedia: http://en.wikipedia.org/wiki/Binomial_distribution):
 >In probability theory and statistics, the binomial distribution is
the discrete probability distribution of the number of successes in a
sequence of n independent yes/no experiments, each of which yields
 >success with probability p.
and the variance square associate is n * p* (1-p)

and therefore error associate to it is (supposing niter=-1)
error(N_acc) = sqrt(N_gen*(N_acc/N_gen)*(1-N_acc/N_gen))

The binomial can be approximate to a Poisson if p <<1 but not in the
case where p is close to one.

Now the case where niter is not -1 requires a Negative binomial
distribution:
(from wikipedia : http://en.wikipedia.org/wiki/Negative_binomial_distribution)
 >In probability theory and statistics, the negative binomial
distribution is a discrete probability distribution of the number of
successes in a sequence of Bernoulli trials before a specified (non-
random) number >of failures (denoted r) occur.

I would say that the error on N_gen is the same of the number of
"suceesses" (which is here not passing the matching)
in that case the error on N_gen should be
error(N_gen) = sqrt( p * r / (1-p)**2)

i.e.

error(N_gen) = sqrt( (1- N_acc/...

[Johan Alwall (johan-alwall)]

Hello Olivier,

> Ok this is problematic (if you fix niter). And then the error depends
> on how pythia calculates the cross-section.
> Do you how how this is done?
> is it sigma_m = n / n_gen
> where this time n_gen is not fixed anymore (i.e. that's the number of
> events required to have n accepted events)?

That sounds right.

> In fact the real point, is that we don't have a poisson distribution
> but a binomial distributions (which is the correct distribution for
> sucess/failure experiment):
> (from wikipedia: http://en.wikipedia.org/wiki/Binomial_distribution):
> >In probability theory and statistics, the binomial distribution is
> the discrete probability distribution of the number of successes in a
> sequence of n independent yes/no experiments, each of which yields
> >success with probability p.
> and the variance square associate is n * p* (1-p)
>
> and therefore error associate to it is (supposing niter=-1)
> error(N_acc) = sqrt(N_gen*(N_acc/N_gen)*(1-N_acc/N_gen))

Exactly. This is exactly what Rikkert wrote about.

> So in the first case (niter=-1) the error is
>
> error * N_acc/N_gen + sigma * sqrt(N_gen*(N_acc/N_gen)*(1-N_acc/
> N_gen)) / N_gen

Exactly. To write it a bit shorter and easier to read:
error * N_acc/N_gen + sigma_m * sqrt((1-N_acc/N_gen)/N_acc)

> while in the second (niter=nacc) the error is
>
> error * N_acc/N_gen + sigma * N_acc / N_gen**2 * (N_gen*sqrt((1-N_acc/
> N_gen)/N_acc)
> i.e.
> error * N_acc/N_gen + sigma * N_acc / N_gen *sqrt((1-N_acc/N_gen)/N_acc)

I don't think this is completely correct however. Note that N_try (which you call "N_gen" above, if I understand you correctly) is here restricted by the actual number N_gen of events in the event file. This is not reflected in the formula.

I would suggest to simply use the first formula also for the second case. I don't think anybody expects us to have the perfect treatment in this case, and that should give a close enough approximation (it has the same limiting behavior when N_acc/N_try = 1 and N_acc/N_try -> 0).

> (or off course their uncorellated equivalent.)

That's right, the errors should be uncorrelated.

Cheers,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

>> while in the second (niter=nacc) the error is
>>
>> error * N_acc/N_gen + sigma * N_acc / N_gen**2 * (N_gen*sqrt((1-
>> N_acc/
>> N_gen)/N_acc)
>> i.e.
>> error * N_acc/N_gen + sigma * N_acc / N_gen *sqrt((1-N_acc/N_gen)/
>> N_acc)
>
> I don't think this is completely correct however. Note that N_try
> (which you call "N_gen" above, if I understand you correctly) is
> here restricted by the actual number N_gen of events in the event
> file. This is not reflected in the formula.

This formula is correct if niter=n_acc (i.e. if you reached the
requested number of events) Otherwise the previous formula is the
valid one.

> I would suggest to simply use the first formula also for the second
> case. I don't think anybody expects us to have the perfect treatment
> in this case, and that should give a close enough approximation (it
> has the same limiting behavior when N_acc/N_try = 1 and N_acc/N_try -
> > 0).

> I don't think anybody expects us to have the perfect treatment in
> this case

That's a wrong argument, off course that if we quote a number they
will expect that we give the correct number (always better to say
nothing that saying smtg wrong)
In fact, including it is not that difficult, since the information is
on the same line that we parse. The only point is that we need to
store the value of the matched error, which is fine.
I can do it if you want.

Cheers,

Olivier

On Apr 17, 2012, at 2:19 AM, Johan Alwall wrote:

> Hello Olivier,
>
>> Ok this is problematic (if you fix niter). And then the error depends
>> on how pythia calculates the cross-section.
>> Do you how how this is done?
>> is it sigma_m = n / n_gen
>> where this time n_gen is not fixed anymore (i.e. that's the number of
>> events required to have n accepted events)?
>
> That sounds right.
>
>> In fact the real point, is that we don't have a poisson distribution
>> but a binomial distributions (which is the correct distribution for
>> sucess/failure experiment):
>> (from wikipedia: http://en.wikipedia.org/wiki/
>> Binomial_distribution):
>>> In probability theory and statistics, the binomial distribution is
>> the discrete probability distribution of the number of successes in a
>> sequence of n independent yes/no experiments, each of which yields
>>> success with probability p.
>> and the variance square associate is n * p* (1-p)
>>
>> and therefore error associate to it is (supposing niter=-1)
>> error(N_acc) = sqrt(N_gen*(N_acc/N_gen)*(1-N_acc/N_gen))
>
> Exactly. This is exactly what Rikkert wrote about.
>
>> So in the first case (niter=-1) the error is
>>
>> error * N_acc/N_gen + sigma * sqrt(N_gen*(N_acc/N_gen)*(1-N_acc/
>> N_gen)) / N_gen
>
> Exactly. To write it a bit shorter and easier to read:
> error * N_acc/N_gen + sigma_m * sqrt((1-N_acc/N_gen)/N_acc)
>
>> while in the second (niter=nacc) the error is
>>
>> error * N_acc/N_gen + sigma * N_acc / N_gen**2 * (N_gen*sqrt((1-
>> N_acc/
>> N_gen)/N_acc)
>> i.e.
>> error * N_acc/N_gen + sigma * N_acc / N_gen *sqrt((1-N_acc/N_gen)/
>> N_acc)
>
> I don't think this is completely correct however. Note that N_try
> (which you call "N_gen" above, if I understand you corr...

[lp:~maddevelopers/mg5amcnlo/1.4.5]

218. By mattelaer-olivier on 2012-04-17

add a warning if the user specify a order larger than the one allowed in couplings_order

[Johan Alwall (johan-alwall)]

Hello Olivier,

> >> while in the second (niter=nacc) the error is
> >>
> >> error * N_acc/N_gen + sigma * N_acc / N_gen**2 * (N_gen*sqrt((1-
> >> N_acc/
> >> N_gen)/N_acc)
> >> i.e.
> >> error * N_acc/N_gen + sigma * N_acc / N_gen *sqrt((1-N_acc/N_gen)/
> >> N_acc)
> >
> > I don't think this is completely correct however. Note that N_try
> > (which you call "N_gen" above, if I understand you correctly) is
> > here restricted by the actual number N_gen of events in the event
> > file. This is not reflected in the formula.
>
> This formula is correct if niter=n_acc (i.e. if you reached the
> requested number of events) Otherwise the previous formula is the
> valid one.
>
> > I would suggest to simply use the first formula also for the second
> > case. I don't think anybody expects us to have the perfect treatment
> > in this case, and that should give a close enough approximation (it
> > has the same limiting behavior when N_acc/N_try = 1 and N_acc/N_try -
> > > 0).
>
> > I don't think anybody expects us to have the perfect treatment in
> > this case
>
> That's a wrong argument, off course that if we quote a number they
> will expect that we give the correct number (always better to say
> nothing that saying smtg wrong)
> In fact, including it is not that difficult, since the information is
> on the same line that we parse. The only point is that we need to
> store the value of the matched error, which is fine.
> I can do it if you want.

Ok, sure. So, basically: If N_try = N_gen, we use the first formula. If N_try < N_gen, we use the second.

Excellent.

Johan

review: Approve

[lp:~maddevelopers/mg5amcnlo/1.4.5]

219. By omatt on 2012-04-18

raise an error if condor_q is unresponsive

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

> Ok, sure. So, basically: If N_try = N_gen, we use the first formula.
> If N_try < N_gen, we use the second.

Yes but the funny part is that the formula is the same (up to Ntry ->
Ngen):
# error for a Negative binomial distribution:
error_m = sqrt((error * Nacc/Ntry)**2 + sigma_m**2 *(1-Nacc/Ntry)/Nacc)

# error for a binomial distributions:
error_m = sqrt((error * Nacc/Ngen)**2 + sigma_m**2 *(1-Nacc/Ngen)/ Nacc)

even if they are compute via different formula...
So we don't have to split the formula at the end.

This makes me slightly uncomfortable but I think this makes some sense
if you take the case where in
order to have N_acc events you reach the points where Ntry=Ngen.
Having the same formula allows to have a perfectly smooth transition
between the two cases.

Do you think it makes sense?

Cheers,

Olivier

On Apr 18, 2012, at 12:03 AM, Johan Alwall wrote:

> Review: Approve
>
> Hello Olivier,
>
>>>> while in the second (niter=nacc) the error is
>>>>
>>>> error * N_acc/N_gen + sigma * N_acc / N_gen**2 * (N_gen*sqrt((1-
>>>> N_acc/
>>>> N_gen)/N_acc)
>>>> i.e.
>>>> error * N_acc/N_gen + sigma * N_acc / N_gen *sqrt((1-N_acc/N_gen)/
>>>> N_acc)
>>>
>>> I don't think this is completely correct however. Note that N_try
>>> (which you call "N_gen" above, if I understand you correctly) is
>>> here restricted by the actual number N_gen of events in the event
>>> file. This is not reflected in the formula.
>>
>> This formula is correct if niter=n_acc (i.e. if you reached the
>> requested number of events) Otherwise the previous formula is the
>> valid one.
>>
>>> I would suggest to simply use the first formula also for the second
>>> case. I don't think anybody expects us to have the perfect treatment
>>> in this case, and that should give a close enough approximation (it
>>> has the same limiting behavior when N_acc/N_try = 1 and N_acc/
>>> N_try -
>>>> 0).
>>
>>> I don't think anybody expects us to have the perfect treatment in
>>> this case
>>
>> That's a wrong argument, off course that if we quote a number they
>> will expect that we give the correct number (always better to say
>> nothing that saying smtg wrong)
>> In fact, including it is not that difficult, since the information is
>> on the same line that we parse. The only point is that we need to
>> store the value of the matched error, which is fine.
>> I can do it if you want.
>
> Ok, sure. So, basically: If N_try = N_gen, we use the first formula.
> If N_try < N_gen, we use the second.
>
> Excellent.
>
> Johan
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/1.4.5/+merge/
> 101662
> You proposed lp:~maddevelopers/madgraph5/1.4.5 for merging.

[Johan Alwall (johan-alwall)]

Very good. Since the limits should be the same, in fact I don't see how the formulas could differ. So this is just a confirmation that indeed the statistical frameworks are correct.

All the best,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Ok Excellent,

I've just fixed a series of small bug.

And add the fact that multi-run ALWAYS change the seed between the run (so even if the run_card is for a fixed value)

I will release 1.4.5 tomorrow.

Cheers and thanks,

Olivier

review: Approve

[lp:~maddevelopers/mg5amcnlo/1.4.5]

220. By mattelaer-olivier on 2012-04-18

change date

221. By mattelaer-olivier on 2012-04-19

put the correct formula for the pythia error

222. By mattelaer-olivier on 2012-04-19

call the function quit when MadEvent is call like ./bin/madevent cmd.cmd

223. By mattelaer-olivier on 2012-04-19

fix an html output problem (when launching generate_events RUN twice)

224. By mattelaer-olivier on 2012-04-19

make stdout/stderr different for pthia/pgs/Delphes
update the seed between the multi-run processes even if
the run_card has a fixed value.

[Johan Alwall (johan-alwall)]

Hello Olivier,

I think all your modifs are ok, just remove any debug print statements.
Oh, also please change the following statement in the UpdateNotes (which can be confusing and worrying for a user):
                  OM+JA: Correct the way the events/error are calculated for matched
                      cross-sections.
to:
                  OM+JA: Correct the display of number of events and error for Pythia
                      in the html files.

Cheers,
Johan

[lp:~maddevelopers/mg5amcnlo/1.4.5]

225. By Olivier Mattelaer on 2012-04-20

remove pointless print
update UpdateNotes
Fix a compilation problems for MadAnalysis

=== modified file 'Template/Source/genps.inc'
--- Template/Source/genps.inc	2011-05-29 07:23:58 +0000
+++ Template/Source/genps.inc	2012-04-20 05:15:24 +0000
@@ -32,4 +32,4 @@
       REAL*8 LIMHEL
       PARAMETER(LIMHEL=1e-6)
       INTEGER MAXTRIES
-      PARAMETER(MAXTRIES=10)
+      PARAMETER(MAXTRIES=100)

=== modified file 'Template/bin/internal/run_pythia'
--- Template/bin/internal/run_pythia	2011-10-18 03:42:31 +0000
+++ Template/bin/internal/run_pythia	2012-04-20 05:15:24 +0000
@@ -19,4 +19,5 @@
 echo "      LHAPATH=$pydir/PDFsets" >> ../Cards/pythia_card.dat
 export PDG_MASS_TBL=$pydir/mass_width_2004.mc
 
-$pydir/pythia && touch pythia.done
\ No newline at end of file
+$pydir/pythia && touch pythia.done
+

=== modified file 'Template/bin/madevent'
--- Template/bin/madevent	2012-02-08 05:31:40 +0000
+++ Template/bin/madevent	2012-04-20 05:15:24 +0000
@@ -126,11 +126,13 @@
             cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')
             cmd_line.use_rawinput = False
             cmd_line.run_cmd('import command ' + input_file)
+            cmd_line.run_cmd('quit')
             sys.exit()
         else:
             cmd_line = cmd_interface.MadEventCmdShell()
             cmd_line.use_rawinput = False
             cmd_line.run_cmd('import command ' + input_file)
+            cmd_line.run_cmd('quit')
             sys.exit()
     else:
         # Interactive mode

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2012-03-29 09:08:22 +0000
+++ UpdateNotes.txt	2012-04-20 05:15:24 +0000
@@ -1,7 +1,22 @@
 Update notes for MadGraph 5 (in reverse time order)
 
+1.4.5 (11/04/12)  OM: Change the seed automatically in multi_run. (Even if the seed
+                      was set to a non automatic value in the card.)
+                  OM: correct a minor bug #975647 (SLAH convention problem) 
+                      Thanks to Sho Iwamoto
+                  OM: Improve cluster support (more secure and complete version)
+                  JA: Increased the number of events tested for non-zero helicity
+                      configurations (needed for goldstino processes).
+                  OM: Add a command to remove the file RunWeb which were not always
+                      deleted correctly
+                  OM+JA: Correct the display of number of events and error for Pythia 
+                     in the html files.
+                  OM: Changed the way the stdout/stderr are treated on the cluster
+                      since some cluster cann't support to have the same output file
+                      for both. (thanks abhishek)
+
 1.4.4 (29/03/12)  OM: Added a command: "output aloha" which allows to creates a 
-       		      subset (or all) of the aloha routines linked to the
+                      subset (or all) of the aloha routines linked to the
                       current model
                   OM: allow to choose the duration of the timer for the questions.
                       (via ./input/mg5_configuration.txt)

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2012-03-29 09:08:22 +0000
+++ madgraph/VERSION	2012-04-20 05:15:24 +0000
@@ -1,3 +1,3 @@
-version = 1.4.4
-date = 2012-03-29
+version = 1.4.5.beta1
+date = 2012-04-18
 

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2012-03-28 03:42:34 +0000
+++ madgraph/core/base_objects.py	2012-04-20 05:15:24 +0000
@@ -1989,7 +1989,11 @@
         tmp = [(k,v) for (k,v) in expansion_orders.items() if 0 < v < 99]
         for (k,v) in tmp:  
             if k in orders:
-                orders[k] = min(orders[k], v)
+                if v < orders[k]:
+                    logger.warning('''The coupling order (%s=%s) specified is larger than the one allowed 
+             by the model builder. The maximal value allowed is %s. 
+             We set the %s order to this value''' % (k,orders[k],v,k))
+                    orders[k] = v
             else:
                 orders[k] = v
 

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2012-03-28 05:07:30 +0000
+++ madgraph/interface/madevent_interface.py	2012-04-20 05:15:24 +0000
@@ -108,10 +108,12 @@
     InvalidCmd = InvalidCmd
     ConfigurationError = MadGraph5Error
 
-    
     def __init__(self, *arg, **opt):
         """Init history and line continuation"""
         
+        # Tag allowing/forbiding question
+        self.force = False
+        
         # If possible, build an info line with current version number 
         # and date, from the VERSION text file
         info = misc.get_pkg_info()
@@ -196,7 +198,10 @@
     
     def postcmd(self, stop, line):
         """ Update the status of  the run for finishing interactive command """
-                        
+        
+        # relaxing the tag forbidding question
+        self.force = False
+        
         if not self.use_rawinput:
             return stop
         
@@ -247,7 +252,7 @@
         """If a ME run is currently running add a link in the html output"""
 
         self.add_error_log_in_html()
-        cmd.Cmd.nice_user_error(self, error, line)
+        cmd.Cmd.nice_user_error(self, error, line)            
         
     def nice_config_error(self, error, line):
         """If a ME run is currently running add a link in the html output"""
@@ -585,11 +590,7 @@
     def check_generate_events(self, args):
         """check that the argument for generate_events are valid"""
         
-        force = False
         run = None
-        if '-f' in args:
-            force = True
-            args.remove('-f')
         if args and args[-1].startswith('--laststep='):
             run = args[-1].split('=')[-1]
             if run not in ['auto','parton', 'pythia', 'pgs', 'delphes']:
@@ -607,19 +608,16 @@
             self.help_generate_events()
             raise self.InvalidCmd('Too many argument for generate_events command: %s' % cmd)
                     
-        return force, run
+        return run
 
     def check_calculate_decay_widths(self, args):
         """check that the argument for calculate_decay_widths are valid"""
         
         if self.ninitial != 1:
             raise self.InvalidCmd('Can only calculate decay widths for decay processes A > B C ...')
-        force = False
+
         accuracy = 0.01
         run = None
-        if '-f' in args:
-            force = True
-            args.remove('-f')
         if args and args[-1].startswith('--accuracy='):
             try:
                 accuracy = float(args[-1].split('=')[-1])
@@ -630,18 +628,14 @@
             self.help_calculate_decay_widths()
             raise self.InvalidCmd('Too many argument for calculate_decay_widths command: %s' % cmd)
                     
-        return force, accuracy
+        return accuracy
 
 
 
     def check_multi_run(self, args):
         """check that the argument for survey are valid"""
 
-        force = False
         run = None
-        if '-f' in args:
-            force = True
-            args.remove('-f')
         
         if not len(args):
             self.help_multi_run()
@@ -669,7 +663,7 @@
         self.check_survey(args, cmd='multi_run')
         args.insert(0, int(nb_run))
         
-        return force, run
+        return run
 
     def check_refine(self, args):
         """check that the argument for survey are valid"""
@@ -747,6 +741,20 @@
         elif laststep:
             raise self.InvalidCmd('only one laststep argument is allowed')
      
+        # If not pythia-pgs path
+        if not self.options['pythia-pgs_path']:
+            logger.info('Retry to read configuration file to find pythia-pgs path')
+            self.set_configuration()
+            
+        if not self.options['pythia-pgs_path'] or not \
+            os.path.exists(pjoin(self.options['pythia-pgs_path'],'src')):
+            error_msg = 'No pythia-pgs path correctly set.'
+            error_msg += 'Please use the set command to define the path and retry.'
+            error_msg += 'You can also define it in the configuration file.'
+            raise self.InvalidCmd(error_msg)
+     
+     
+     
         tag = [a for a in args if a.startswith('--tag=')]
         if tag: 
             args.remove(tag[0])
@@ -775,19 +783,6 @@
         input_file = pjoin(self.me_dir,'Events',self.run_name, 'unweighted_events.lhe')
         output_file = pjoin(self.me_dir, 'Events', 'unweighted_events.lhe')
         os.system('gunzip -c %s > %s' % (input_file, output_file))
-            
-
-        # If not pythia-pgs path
-        if not self.options['pythia-pgs_path']:
-            logger.info('Retry to read configuration file to find pythia-pgs path')
-            self.set_configuration()
-            
-        if not self.options['pythia-pgs_path'] or not \
-            os.path.exists(pjoin(self.options['pythia-pgs_path'],'src')):
-            error_msg = 'No pythia-pgs path correctly set.'
-            error_msg += 'Please use the set command to define the path and retry.'
-            error_msg += 'You can also define it in the configuration file.'
-            raise self.InvalidCmd(error_msg)
         
         args.append(mode)
     
@@ -801,14 +796,6 @@
             tag = tag[0]
             tmp_args.remove('--tag=%s' % tag)
 
-        try:
-            tmp_args.remove('-f')
-        except:
-            pass
-        else:
-            if args[0] == '-f':
-                args.pop(0)
-                args.append('-f')
 
         if len(tmp_args) == 0:
             self.help_remove()
@@ -817,7 +804,7 @@
             return tmp_args[0], tag, ['all']
         else:
             for arg in tmp_args[1:]:
-                if arg not in self._clean_mode and arg != '-f':
+                if arg not in self._clean_mode:
                     self.help_clean()
                     raise self.InvalidCmd('%s is not a valid options for clean command'\
                                               % arg)
@@ -855,10 +842,6 @@
             args.append('all')
             return
 
-        force = False
-        if '-f' in args:
-            args.remove('-f')
-            force = True
         
         if args[0] not in self._plot_mode:
             self.set_run_name(args[0], level='plot')
@@ -874,9 +857,6 @@
                  self.help_plot()
                  raise self.InvalidCmd('unknown options %s' % arg)        
     
-        if force:
-            args.append('-f')
-    
     
     def check_pgs(self, arg):
         """Check the argument for pythia command
@@ -904,7 +884,7 @@
         
         if len(arg) == 0 and not self.run_name:
             if self.results.lastrun:
-                args.insert(0, self.results.lastrun)
+                arg.insert(0, self.results.lastrun)
             else:
                 raise self.InvalidCmd('No run name currently define. Please add this information.')             
         
@@ -1402,7 +1382,7 @@
             elif arg == '-m':
                 self.cluster_mode = 2
             elif arg == '-f':
-                continue
+                self.force = True
             elif not arg.startswith('--'):
                 if error:
                     raise self.InvalidCmd('%s argument cannot start with - symbol' % arg)
@@ -1419,8 +1399,7 @@
                 self.web = True
                 self.cluster_mode = 1
                 self.results.def_web_mode(True)
-                if not '-f' in args:
-                    args.append('-f')
+                self.force = True
             else:
                 continue
             args.remove(arg)
@@ -1544,11 +1523,7 @@
         
         args = self.split_arg(line)
         #check the validity of the arguments
-        self.check_banner_run(args)
-        if '-f' in args:
-            force = True
-        else:
-            force = False     
+        self.check_banner_run(args)    
              
         banner_mod.split_banner(args[0], self.me_dir, proc_card=False)
         
@@ -1562,13 +1537,13 @@
         
         
         # Check if we want to modify the run
-        if not force:
+        if not self.force:
             ans = self.ask('Do you want to modify the Cards?', 'n', ['y','n'])
             if ans == 'n':
-                force = True
+                self.force = True
         
         # Call Generate events
-        self.exec_cmd('generate_events %s %s' % (self.run_name, force and '-f' or ''))
+        self.exec_cmd('generate_events %s %s' % (self.run_name, self.force and '-f' or ''))
  
  
  
@@ -1704,8 +1679,8 @@
         
         args = self.split_arg(line)
         # Check argument's validity
-        force, mode = self.check_generate_events(args)
-        self.ask_run_configuration(mode, force)
+        mode = self.check_generate_events(args)
+        self.ask_run_configuration(mode)
         if not args:
             # No run name assigned -> assigned one automaticaly 
             self.set_run_name(self.find_available_run_name(self.me_dir), None, 'parton')
@@ -1748,14 +1723,13 @@
             self.exec_cmd('refine %s' % nb_event, postcmd=False)
             self.exec_cmd('refine %s' % nb_event, postcmd=False)
             self.exec_cmd('combine_events', postcmd=False)
+            self.print_results_in_shell(self.results.current)
             self.create_plot('parton')
             self.exec_cmd('store_events', postcmd=False)
             self.exec_cmd('pythia --no_default', postcmd=False, printcmd=False)
             # pythia launches pgs/delphes if needed
             self.store_result()
         
-        self.print_results_in_shell(self.results.current)
-            
             
     def print_results_in_shell(self, data):
         """Have a nice results prints in the shell,
@@ -1767,18 +1741,12 @@
         else:
             logger.info("     Cross-section :   %.4g +- %.4g pb" % (data['cross'], data['error']))
         logger.info("     Nb of events :  %s" % data['nb_event'] )
-        if data['cross_pythia']:
-            error = data.get_pythia_error(data['cross'], data['error'], 
-                                        data['cross_pythia'], data['nb_event'])
-            nb_event = 0
-            if data['cross']:
-                nb_event = int(0.5+(data['nb_event'] * data['cross_pythia'] /data['cross']))
-            
+        if data['cross_pythia'] and data['nb_event_pythia']:
             if self.ninitial == 1:
-                logger.info("     Matched Width :   %.4g +- %.4g GeV" % (data['cross_pythia'], error))
+                logger.info("     Matched Width :   %.4g +- %.4g GeV" % (data['cross_pythia'], data['error_pythia']))
             else:
-                logger.info("     Matched Cross-section :   %.4g +- %.4g pb" % (data['cross'], error))            
-            logger.info("     Nb of events after Matching :  %s" % nb_event)
+                logger.info("     Matched Cross-section :   %.4g +- %.4g pb" % (data['cross_pythia'], data['error_pythia']))            
+            logger.info("     Nb of events after Matching :  %s" % data['nb_event_pythia'])
         logger.info(" " )
         
             
@@ -1790,8 +1758,8 @@
         
         args = self.split_arg(line)
         # Check argument's validity
-        force, accuracy = self.check_calculate_decay_widths(args)
-        self.ask_run_configuration('parton', force)
+        accuracy = self.check_calculate_decay_widths(args)
+        self.ask_run_configuration('parton')
         if not args:
             # No run name assigned -> assigned one automaticaly 
             self.set_run_name(self.find_available_run_name(self.me_dir))
@@ -1926,8 +1894,8 @@
         
         args = self.split_arg(line)
         # Check argument's validity
-        force, mode = self.check_multi_run(args)
-        self.ask_run_configuration(mode, force)
+        mode = self.check_multi_run(args)
+        self.ask_run_configuration(mode)
         main_name = self.run_name
         nb_run = args.pop(0)
         crossoversig = 0
@@ -1943,8 +1911,12 @@
             crossoversig+=cross/error**2
             inv_sq_err+=1.0/error**2
             self.results[main_name][-1]['cross'] = crossoversig/inv_sq_err
-            self.results[main_name][-1]['error'] = math.sqrt(1.0/inv_sq_err) 
-
+            self.results[main_name][-1]['error'] = math.sqrt(1.0/inv_sq_err)
+            if self.run_card['iseed'] != 0:
+                seed = self.random + 1
+                text = open(pjoin(self.me_dir, 'Cards','run_card.dat')).read()
+                (t,n) = re.subn(r'\s\d+\s*= iseed','   %s       = iseed' % seed,text)
+                open(pjoin(self.me_dir, 'Cards','run_card.dat'),'w').write(t)
         
         self.run_name = main_name
         self.results.def_current(main_name)
@@ -2281,24 +2253,18 @@
         
         # Check argument's validity
         args = self.split_arg(line)
-        if '-f' in args:
-            force = True
-            args.remove('-f')
-        else:
-            force = False
         if '--no_default' in args:
             if not os.path.exists(pjoin(self.me_dir, 'Cards', 'pythia_card.dat')):
                 return
-            force = True
             no_default = True
             args.remove('--no_default')
         else:
             no_default = False
                                     
         self.check_pythia(args)        
-        # the args are modify and the last arg is always the mode
-        
-        self.ask_pythia_run_configuration(args[-1], force)
+        # the args are modify and the last arg is always the mode 
+        if not no_default:
+            self.ask_pythia_run_configuration(args[-1])
 
         # Update the banner with the pythia card
         if not self.banner:
@@ -2307,15 +2273,6 @@
         # initialize / remove lhapdf mode        
         self.configure_directory()
 
-        #if not force:
-        #    if os.path.exists(pjoin(self.me_dir, 'Events', self.run_name, '%s_pythia.log' % tag)):
-        #        question = 'Previous run of pythia detected. Do you want to remove it?'
-        #        ans = self.ask(question, 'y', choices=['y','n'], timeout = 20)
-        #        if ans == 'n':
-        #            return
-         
-        #self.exec_cmd('remove %s pythia -f' % self.run_name)
-
         pythia_src = pjoin(self.options['pythia-pgs_path'],'src')
         
         self.update_status('Running Pythia', 'pythia')
@@ -2352,12 +2309,36 @@
         self.to_store.append('pythia')
         
         # Find the matched cross-section
-        if int(self.run_card['ickkw']):    
-            pythia_log = open(pjoin(self.me_dir,'Events', self.run_name, '%s_pythia.log' % tag))
-            cs_info = misc.get_last_line(pythia_log)
-            # line should be of type: Cross section (pb):    1840.20000000006 
-            cs_info = cs_info.split(':')[1]
-            self.results.add_detail('cross_pythia', cs_info)
+        if int(self.run_card['ickkw']):
+            # read the line from the bottom of the file
+            pythia_log = misc.BackRead(pjoin(self.me_dir,'Events', self.run_name, 
+                                                         '%s_pythia.log' % tag))
+            pythiare = re.compile("\s*I\s+0 All included subprocesses\s+I\s+(?P<generated>\d+)\s+(?P<tried>\d+)\s+I\s+(?P<xsec>[\d\.D\-+]+)\s+I")            
+            for line in pythia_log:
+                info = pythiare.search(line)
+                if not info:
+                    continue
+                try:
+                    # Pythia cross section in mb, we want pb
+                    sigma_m = float(info.group('xsec').replace('D','E')) *1e9
+                    Nacc = int(info.group('generated'))
+                    Ntry = int(info.group('tried'))
+                except ValueError:
+                    # xsec is not float - this should not happen
+                    self.results.add_detail('cross_pythia', 0)
+                    self.results.add_detail('nb_event_pythia', 0)
+                    self.results.add_detail('error_pythia', 0)
+                else:
+                    self.results.add_detail('cross_pythia', sigma_m)
+                    self.results.add_detail('nb_event_pythia', Nacc)
+                    #compute pythia error
+                    error = self.results[self.run_name].return_tag(self.run_tag)['error']                    
+                    error_m = math.sqrt((error * Nacc/Ntry)**2 + sigma_m**2 *(1-Nacc/Ntry)/Nacc)
+                    # works both for fixed number of generated events and fixed accepted events
+                    self.results.add_detail('error_pythia', error_m)
+                break                 
+
+            pythia_log.close()
         
         pydir = pjoin(self.options['pythia-pgs_path'], 'src')
         eradir = self.options['exrootanalysis_path']
@@ -2431,6 +2412,8 @@
         self.exec_cmd('pgs --no_default', postcmd=False, printcmd=False)
         if self.options['delphes_path']:
             self.exec_cmd('delphes --no_default', postcmd=False, printcmd=False)
+
+        self.print_results_in_shell(self.results.current)
     
     def get_available_tag(self):
         """create automatically a tag"""
@@ -2510,7 +2493,7 @@
                 to_suppress = [f for f in to_suppress if 'delphes' in f 
                                                       or 'pgs' in f 
                                                       or 'pythia' in f]
-        if '-f' not in args and len(to_suppress):
+        if not self.force and len(to_suppress):
             question = 'Do you want to suppress the following files?\n     %s' % \
                                '\n    '.join(to_suppress)
             ans = self.ask(question, 'y', choices=['y','n'])
@@ -2544,7 +2527,7 @@
                 to_suppress += glob.glob(pjoin(self.me_dir, 'SubProcesses', '*','%s*' % run))
                 to_suppress += glob.glob(pjoin(self.me_dir, 'SubProcesses', '*','*','%s*' % run))
 
-                if '-f' in args or len(to_suppress) == 0:
+                if self.force or len(to_suppress) == 0:
                     ans = 'y'
                 else:
                     question = 'Do you want to suppress the following files?\n     %s' % \
@@ -2685,11 +2668,6 @@
         
         args = self.split_arg(line)
         # Check argument's validity
-        if '-f' in args:
-            force = True
-            args.remove('-f')
-        else:
-            force = False
         if '--no_default' in args:
             no_default = True
             args.remove('--no_default')
@@ -2712,7 +2690,7 @@
                      pjoin(self.me_dir, 'Cards', 'pgs_card.dat'))
             logger.info('No pgs card found. Take the default one.')        
         
-        if not (no_default or force):
+        if not (no_default or self.force):
             self.ask_edit_cards(['pgs'], args)
             
         self.update_status('prepare PGS run', level=None)  
@@ -2807,11 +2785,6 @@
  
         args = self.split_arg(line)
         # Check argument's validity
-        if '-f' in args:
-            force = True
-            args.remove('-f')
-        else:
-            force = False
         if '--no_default' in args:
             no_default = True
             args.remove('--no_default')
@@ -2833,7 +2806,7 @@
         if not os.path.exists(pjoin(self.me_dir, 'Cards', 'delphes_trigger.dat')):    
             files.cp(pjoin(self.me_dir, 'Cards', 'delphes_trigger_default.dat'),
                      pjoin(self.me_dir, 'Cards', 'delphes_trigger.dat'))
-        if not (no_default or force):
+        if not (no_default or self.force):
             self.ask_edit_cards(['delphes', 'trigger'], args)
             
         self.update_status('Running Delphes', level=None)  
@@ -2984,8 +2957,19 @@
                     self.update_status((idle, run, finish, run_type), level=None)
             else:
                 update_status = lambda idle, run, finish: None
-            self.cluster.wait(self.me_dir, update_status)            
-
+            try:    
+                self.cluster.wait(self.me_dir, update_status)            
+            except Exception, error:
+                logger.info(error)
+                if not self.force:
+                    ans = self.ask('Cluster Error detected. Do you want to clean the queue?',
+                             default = 'y', answers=['y','n'])
+                else:
+                    ans = 'y'
+                if ans:
+                    self.cluster.remove()
+                raise
+                    
         if mode == 2:
             # Wait that all thread finish
             if not self.control_thread[2]:
@@ -3404,7 +3388,10 @@
             
             
         if not os.path.exists(event_path):
-            raise self.InvalidCmd, 'Events file %s does not exits' % event_path
+            if os.path.exists(event_path+'.gz'):
+                os.system('gzip -f %s.gz ' % event_path)
+            else:
+                raise self.InvalidCmd, 'Events file %s does not exits' % event_path
         
         self.update_status('Creating Plots for %s level' % mode, level = mode.lower())
                
@@ -3462,7 +3449,7 @@
 
 
     ############################################################################
-    def ask_run_configuration(self, mode=None, force=False):
+    def ask_run_configuration(self, mode=None):
         """Ask the question when launching generate_events/multi_run"""
         
         available_mode = ['0', '1']
@@ -3487,7 +3474,7 @@
             if '4' in available_mode:
                 question += """  4 / delphes :  MadEvent + Pythia + Delphes.\n"""
 
-            if not force:
+            if not self.force:
                 if not mode:
                     mode = self.ask(question, '0', options)
             elif not mode:
@@ -3524,7 +3511,7 @@
             self.add_card_to_run('trigger')
             cards.append('delphes_card.dat')
 
-        if force:
+        if self.force:
             return
 
         def get_question(mode):
@@ -3599,7 +3586,7 @@
                     
                     
     ############################################################################
-    def ask_pythia_run_configuration(self, mode=None, force=False):
+    def ask_pythia_run_configuration(self, mode=None):
         """Ask the question when launching pythia"""
         
         available_mode = ['0', '1', '2']
@@ -3614,7 +3601,7 @@
         if '3' in available_mode:
             question += """  3 / delphes  : Pythia + Delphes.\n"""
 
-        if not force:
+        if not self.force:
             if not mode:
                 mode = self.ask(question, '0', options)
         elif not mode:
@@ -3647,7 +3634,7 @@
             self.add_card_to_run('trigger')
             cards.append('delphes_card.dat')
         
-        if force:
+        if self.force:
             return mode
         
         # Ask the user if he wants to edit any of the files
@@ -3708,7 +3695,7 @@
     def ask_edit_cards(self, cards, fct_args):
         """Question for cards editions (used for pgs/delphes)"""
 
-        if '-f' in fct_args or '--no_default' in fct_args:
+        if self.force or '--no_default' in fct_args:
             return
         
         card_name = {'pgs': 'pgs_card.dat',

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2012-03-28 05:07:30 +0000
+++ madgraph/interface/madgraph_interface.py	2012-04-20 05:15:24 +0000
@@ -1896,7 +1896,6 @@
                 # check if a particle is asked more than once
                 if len(request_part) > len(set(request_part)):
                     for p in request_part:
-                        print p, request_part.count(p),present_part.count(p)
                         if request_part.count(p) > present_part.count(p):
                             continue
                         
@@ -2820,14 +2819,16 @@
                 path = os.path.join(MG5DIR, 'pythia-pgs', 'src', 'make_opts')
                 text = open(path).read()
                 text = text.replace('FC=g77','FC=gfortran')
-                open(path, 'w').writelines(text)            
+                open(path, 'w').writelines(text)    
             elif compiler == 'gfortran' and args[0] == 'MadAnalysis':
                 path = os.path.join(MG5DIR, 'MadAnalysis', 'makefile')
                 text = open(path).read()
-                text = text.replace('F77 = g77','F77 = gfortran')
-                open(path, 'w').writelines(text)            
-        if logger.level <= logging.INFO: 
-            misc.call(['make', 'clean'], )
+                text = text.replace('FC=g77','FC=gfortran')
+                open(path, 'w').writelines(text)
+                            
+        if logger.level <= logging.INFO:
+            devnull = open(os.devnull,'w') 
+            misc.call(['make', 'clean'], stdout=devnull, stderr=-2)
             status = misc.call(['make'], cwd = os.path.join(MG5DIR, name))
         else:
             misc.compile(['clean'], mode='', cwd = os.path.join(MG5DIR, name))

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2011-05-30 17:37:49 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2012-04-20 05:15:24 +0000
@@ -117,7 +117,7 @@
                GOODHEL(I,IMIRROR)=.TRUE.
                NGOOD(IMIRROR) = NGOOD(IMIRROR) +1
                IGOOD(NGOOD(IMIRROR),IMIRROR) = I
-               PRINT *,'Added good helicity ',I,TS(I)*NCOMB/ANS
+               PRINT *,'Added good helicity ',I,TS(I)*NCOMB/ANS,' in event ',NTRY(IMIRROR)
             ENDIF
          ENDDO
       ENDIF

=== modified file 'madgraph/various/cluster.py'
--- madgraph/various/cluster.py	2012-03-20 19:55:40 +0000
+++ madgraph/various/cluster.py	2012-04-20 05:15:24 +0000
@@ -31,7 +31,7 @@
 class NotImplemented(MadGraph5Error):
     pass
 
-def multiple_try(nb_try=5, sleep=1):
+def multiple_try(nb_try=5, sleep=20):
 
     def deco_retry(f):
         def deco_f_retry(*args, **opt):
@@ -41,7 +41,7 @@
                 except KeyboardInterrupt:
                     raise
                 except:
-                    time.sleep(sleep)
+                    time.sleep(sleep * (i+1))
             raise
         return deco_f_retry
     return deco_retry
@@ -118,6 +118,12 @@
     def launch_and_wait(self, prog, argument=[], cwd=None, stdout=None, 
                                                          stderr=None, log=None):
         """launch one job on the cluster and wait for it"""
+        
+        special_output = False # tag for concatanate the error with the output.
+        if stderr == -2 and stdout: 
+            #We are suppose to send the output to stdout
+            special_output = True
+            stderr = stdout + '.err'
         id = self.submit(prog, argument, cwd, stdout, stderr, log)
         while 1:        
             status = self.control_one_job(id)
@@ -125,7 +131,25 @@
                 time.sleep(20) #security to ensure that the file are really written on the disk
                 break
             time.sleep(30)
-            
+        
+        if special_output:
+            # combine the stdout and the stderr
+            #wait up to 50 s to see if those files exists
+            for i in range(5):
+                if os.path.exists(stdout):
+                    if not os.path.exists(stderr):
+                        time.sleep(5)
+                    if os.path.exists(stderr):
+                        err_text = open(stderr).read()
+                        if not err_text:
+                            return
+                        logger.warning(err_txt)                        
+                        text = open(stdout).read()
+                        open(stdout,'w').write(text + err_txt)
+                    else:
+                        return
+                time.sleep(10)
+                        
     def remove(self, *args):
         """ """
         logger.warning("""This cluster didn't support job removal, 
@@ -198,7 +222,12 @@
     def control_one_job(self, id):
         """ control the status of a single job with it's cluster id """
         cmd = 'condor_q '+str(id)+" -format \'%-2s \\n\' \'ifThenElse(JobStatus==0,\"U\",ifThenElse(JobStatus==1,\"I\",ifThenElse(JobStatus==2,\"R\",ifThenElse(JobStatus==3,\"X\",ifThenElse(JobStatus==4,\"C\",ifThenElse(JobStatus==5,\"H\",ifThenElse(JobStatus==6,\"E\",string(JobStatus))))))))\'"
-        status = subprocess.Popen([cmd], shell=True, stdout=subprocess.PIPE)
+        status = subprocess.Popen([cmd], shell=True, stdout=subprocess.PIPE, 
+                                                         stderr=subprocess.PIPE)
+        if status.returncode:
+            raise ClusterManagmentError, 'condor_q returns error: %s' % \
+                                                            status.stderr.read()
+
         return status.stdout.readline().strip()
     
     @check_interupt()
@@ -210,16 +239,22 @@
             return 0, 0, 0, 0
         
         cmd = "condor_q " + ' '.join(self.submitted_ids) + " -format \'%-2s \\n\' \'ifThenElse(JobStatus==0,\"U\",ifThenElse(JobStatus==1,\"I\",ifThenElse(JobStatus==2,\"R\",ifThenElse(JobStatus==3,\"X\",ifThenElse(JobStatus==4,\"C\",ifThenElse(JobStatus==5,\"H\",ifThenElse(JobStatus==6,\"E\",string(JobStatus))))))))\'"
-        status = subprocess.Popen([cmd], shell=True, stdout=subprocess.PIPE)
+        status = subprocess.Popen([cmd], shell=True, stdout=subprocess.PIPE, 
+                                                         stderr=subprocess.PIPE)
 
+        if status.returncode:
+            raise ClusterManagmentError, 'condor_q returns error: %s' % \
+                                                            status.stderr.read()
+            
+            
         idle, run, fail = 0, 0, 0
         for line in status.stdout:
             status = line.strip()
-            if status == 'I':
+            if status in ['I','U']:
                 idle += 1
             elif status == 'R':
                 run += 1
-            else:
+            elif status != 'C':
                 fail += 1
 
         return idle, run, self.submitted - (idle+run+fail), fail
@@ -240,6 +275,7 @@
     name = 'pbs'
     idle_tag = ['Q']
     running_tag = ['T','E','R']
+    complete_tag = ['C']
 
     @multiple_try()
     def submit(self, prog, argument=[], cwd=None, stdout=None, stderr=None, log=None):
@@ -324,6 +360,8 @@
                     idle += 1
                 elif status in self.running_tag:
                     run += 1
+                elif status in self.complete_tag:
+                    continue
                 else:
                     fail += 1
 

=== modified file 'madgraph/various/gen_crossxhtml.py'
--- madgraph/various/gen_crossxhtml.py	2012-03-24 03:06:44 +0000
+++ madgraph/various/gen_crossxhtml.py	2012-04-20 05:15:24 +0000
@@ -278,7 +278,8 @@
 
     def add_detail(self, name, value, run=None, tag=None):
         """ add information to current run (cross/error/event)"""
-        assert name in ['cross', 'error', 'nb_event', 'cross_pythia']
+        assert name in ['cross', 'error', 'nb_event', 'cross_pythia',
+                        'nb_event_pythia','error_pythia']
 
         if not run and not self.current:
             return
@@ -292,6 +293,8 @@
             run['cross_pythia'] = float(value)
         elif name == 'nb_event':
             run[name] = int(value)
+        elif name == 'nb_event_pythia':
+            run[name] = int(value)
         else:    
             run[name] = float(value)    
     
@@ -411,6 +414,13 @@
         except:
             self.web = False
 
+        # check if more than one parton output
+        parton = [r for r in self if r.parton]
+        # clean wrong previous run link
+        if len(parton)>1:
+            for p in parton[:-1]:
+                p.parton = []
+
         dico = self.info
         dico['run_span'] = sum([tag.get_nb_line() for tag in self], 1) -1
         dico['tag_data'] = '\n'.join([tag.get_html(self) for tag in self])
@@ -465,11 +475,10 @@
             output['cross'] = self[-2]['cross']
             output['error'] = self[-2]['error']
         elif (current.pgs or current.delphes) and not current['nb_event'] and len(self) > 1:
-            if self[-2]['cross_pythia']:
+            if self[-2]['cross_pythia'] and self[-2]['nb_event_pythia']:
                 output['cross'] = self[-2]['cross_pythia']
-                output['nb_event'] = int(0.5+(self[-2]['nb_event'] * current['cross'] /self[-2]['cross']))                           
-                output['error'] = current.get_pythia_error(self[-2]['cross'], 
-                       self[-2]['error'], current['cross'], current['nb_event'])
+                output['nb_event'] = self[-2]['nb_event_pythia']
+                output['error'] = self[-2]['error_pythia']
             else:
                 output['nb_event'] = self[-2]['nb_event']
                 output['cross'] = self[-2]['cross']
@@ -514,6 +523,7 @@
         self['nb_event'] = 0
         self['cross'] = 0
         self['cross_pythia'] = ''
+        self['nb_event_pythia'] = 0
         self['error'] = 0
         self.parton = [] 
         self.pythia = []
@@ -638,20 +648,7 @@
         
          return " <a  id='%(id)s' href='%(link1)s' onClick=\"check_link('%(link1)s','%(link2)s','%(id)s')\">%(name)s</a>" \
               % {'link1': link1, 'link2':link2, 'id': id, 'name':name}       
-    
-    def get_pythia_error(self, cross, error, pythia_cross, nb_event):
-        """compute the error associate to pythie"""
-        # pythia_cross = cross * n_acc / n_gen
-        # error_pythia = error * n_acc /n_gen + cross * sqrt(n_acc) / n_gen
         
-        if cross and nb_event: 
-            n_acc = int(0.5 + pythia_cross / cross * nb_event)
-            error_pythia = error * n_acc / nb_event 
-            error_pythia += cross * math.sqrt(n_acc) / nb_event
-        else:
-            error_pythia = 0
-        return error_pythia
-    
     def get_links(self, level):
         """ Get the links for a given level"""
         
@@ -776,9 +773,8 @@
         elif (self.pgs or self.delphes) and not self['nb_event']:
             if runresults[-2]['cross_pythia'] and runresults[-2]['cross']:
                 self['cross'] = runresults[-2]['cross_pythia']
-                self['nb_event'] = int(0.5+(runresults[-2]['nb_event'] * self['cross'] /runresults[-2]['cross']))                           
-                self['error'] = self.get_pythia_error(runresults[-2]['cross'], 
-                       runresults[-2]['error'], self['cross'], self['nb_event'])
+                self['error'] = runresults[-2]['error_pythia']
+                self['nb_event'] = runresults[-2]['nb_event_pythia']                           
             else:
                 self['nb_event'] = runresults[-2]['nb_event']
                 self['cross'] = runresults[-2]['cross']
@@ -807,13 +803,12 @@
                         local_dico['cross_span'] = nb_line -1
                     else:
                         local_dico['cross_span'] = nb_line
-                    if self['cross']:
-                        local_dico['nb_event'] = int(0.5+(self['nb_event'] * self['cross_pythia'] /self['cross']))
+                    if self['nb_event_pythia']:
+                        local_dico['nb_event'] = self['nb_event_pythia']
                     else:
                         local_dico['nb_event'] = 0
                     local_dico['cross'] = self['cross_pythia']
-                    local_dico['err'] = self.get_pythia_error(self['cross'],
-                           self['error'],self['cross_pythia'], self['nb_event']) 
+                    local_dico['err'] = self['error_pythia']
                 else:
                     local_dico['cross_span'] = nb_line
                     local_dico['cross'] = self['cross']
@@ -824,16 +819,15 @@
                 template = sub_part_template_parton
                 if self.parton:           
                     local_dico['cross_span'] = nb_line - 1
-                    if self['cross']:
-                        local_dico['nb_event'] = int(0.5+(self['nb_event'] * self['cross_pythia'] /self['cross']))
+                    if self['nb_event_pythia']:
+                        local_dico['nb_event'] = self['nb_event_pythia']
                     else:
                         local_dico['nb_event'] = 0
                 else:
                     local_dico['cross_span'] = nb_line
                     local_dico['nb_event'] = self['nb_event']
                 local_dico['cross'] = self['cross_pythia']
-                local_dico['err'] = self.get_pythia_error(self['cross'],
-                           self['error'],self['cross_pythia'], self['nb_event'])
+                local_dico['err'] = self['error_pythia']
             else:
                template = sub_part_template_pgs             
             

=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2012-03-20 07:16:18 +0000
+++ madgraph/various/misc.py	2012-04-20 05:15:24 +0000
@@ -261,6 +261,7 @@
         if len(lines) >= to_read or pos == 0:
             return lines[-to_read:offset and -offset or None]
         avg_line_length *= 1.3
+        avg_line_length = int(avg_line_length)
 
 ################################################################################
 # LAST LINE FUNCTION
@@ -270,6 +271,66 @@
     
     return tail(fsock, 1)[0]
     
+class BackRead(file):
+    """read a file returning the lines in reverse order for each call of readline()
+This actually just reads blocks (4096 bytes by default) of data from the end of
+the file and returns last line in an internal buffer."""
+
+
+    def readline(self):
+        """ readline in a backward way """
+        
+        while len(self.data) == 1 and ((self.blkcount * self.blksize) < self.size):
+          self.blkcount = self.blkcount + 1
+          line = self.data[0]
+          try:
+            self.seek(-self.blksize * self.blkcount, 2) # read from end of file
+            self.data = (self.read(self.blksize) + line).split('\n')
+          except IOError:  # can't seek before the beginning of the file
+            self.seek(0)
+            self.data = string.split(self.read(self.size - (self.blksize * (self.blkcount-1))) + line, '\n')
+    
+        if len(self.data) == 0:
+          return ""
+    
+        line = self.data.pop()
+        return line + '\n'
+
+    def __init__(self, filepos, blksize=4096):
+        """initialize the internal structures"""
+
+        # get the file size
+        self.size = os.stat(filepos)[6]
+        # how big of a block to read from the file...
+        self.blksize = blksize
+        # how many blocks we've read
+        self.blkcount = 1
+        file.__init__(self, filepos, 'rb')
+        # if the file is smaller than the blocksize, read a block,
+        # otherwise, read the whole thing...
+        if self.size > self.blksize:
+          self.seek(-self.blksize * self.blkcount, 2) # read from end of file
+        self.data = self.read(self.blksize).split('\n')
+        # strip the last item if it's empty...  a byproduct of the last line having
+        # a newline at the end of it
+        if not self.data[-1]:
+          self.data.pop()
+        
+    def next(self):
+        line = self.readline()
+        if line:
+            return line
+        else:
+            raise StopIteration
+        
+    def close(self):
+        """ close correctly file """
+        try:
+            self.close()
+        except:
+           pass
+
+
 
 
 #

=== modified file 'madgraph/various/sum_html.py'
--- madgraph/various/sum_html.py	2012-02-27 22:22:51 +0000
+++ madgraph/various/sum_html.py	2012-04-20 05:15:24 +0000
@@ -20,7 +20,10 @@
 logger = logging.getLogger('madevent.stdout') # -> stdout
 
 pjoin = os.path.join
-
+try:
+    import madgraph.various.cluster as cluster
+except:
+    import internal.cluster as cluster
 
 class OneResult(object):
     
@@ -41,6 +44,7 @@
         self.yerr_iter = []
         return
     
+    @cluster.multiple_try(nb_try=4,sleep=5)
     def read_results(self, filepath):
         """read results.dat and fullfill information"""
         

=== modified file 'models/mssm/parameters.py'
--- models/mssm/parameters.py	2011-03-30 05:36:25 +0000
+++ models/mssm/parameters.py	2012-04-20 05:15:24 +0000
@@ -215,7 +215,7 @@
                  value = 32337.4943,
                  texname = 'm_{H_u}^2',
                  lhablock = 'MSOFT',
-                 lhacode = [ 21 ])
+                 lhacode = [ 22 ])
 
 mHd2 = Parameter(name = 'mHd2',
                  nature = 'external',
@@ -223,7 +223,7 @@
                  value = -128800.134,
                  texname = 'm_{H_d}^2',
                  lhablock = 'MSOFT',
-                 lhacode = [ 22 ])
+                 lhacode = [ 21 ])
 
 RmQ211 = Parameter(name = 'RmQ211',
                    nature = 'external',

=== modified file 'tests/unit_tests/iolibs/test_export_v4.py'
--- tests/unit_tests/iolibs/test_export_v4.py	2012-03-28 04:08:18 +0000
+++ tests/unit_tests/iolibs/test_export_v4.py	2012-04-20 05:15:24 +0000
@@ -899,6 +899,7 @@
               NGOOD(IMIRROR) = NGOOD(IMIRROR) +1
               IGOOD(NGOOD(IMIRROR),IMIRROR) = I
               PRINT *,'Added good helicity ',I,TS(I)*NCOMB/ANS
+     $         ,' in event ',NTRY(IMIRROR)
             ENDIF
           ENDDO
         ENDIF