Bias Re-weighting in DY+Jets

Asked by Jonas Roemer on 2017-10-26

Dear experts,

I have a question regarding the bias re-weighting. I want to produce a DY sample with a smooth mass tail up to masses of a few TeV for the excited lepton analysis in CMS. In order to avoid producing several mass-binned tail samples I want to use the bias re-weighting technique.

In the end I want to select events in the 2l2j final state and I need the NLO shape as the LO shape does not fit the observed spectrum well. For testing purposes I also tried the production at LO and I found issues even though different to the ones in NLO.

For LO I wrote a bias re-weighting module as explained in the wiki (https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias). I calculate the mass of the Z/gamma and use it in the formula to calculate the bias weight (first power-law then exponential for the high mass tails.) The code (besides variable) definitions is:

C parameters = {'bias_normalization': 1000.0,
C 'bias_exponent': 4.29171,
C 'bias_switch': 1.73642d3,
C 'bias_exponential': 1.41843d-3,
C 'bias_cutoff': 1d4}
C

C --------------------
C BEGIN IMPLEMENTATION
C --------------------

do i=1,nexternal
if (is_a_l(i)) then
sump(0)=sump(0)+p(0,i)
sump(1)=sump(1)+p(1,i)
sump(2)=sump(2)+p(2,i)
sump(3)=sump(3)+p(3,i)
endif
enddo

mass = sqrt(sump(0)**2 - sump(1)**2 - sump(2)**2 - sump(3)**2)

bias_weight=(mass/bias_normalization)**bias_exponent
bias_weight=(mass**bias_exponent/bias_normalization)
if (mass.gt.bias_switch) then
bias_weight=bias_weight*
& exp((mass-bias_switch)*bias_exponential)
endif

return

If I generate p p > e- e+ I get a really nice mass spectrum which is smooth up to masses of ~10TeV. I can also add jets (p p > e- e+ j) and it still works fine. However, combining both processes into one job does not succeed:

generate p p > e+ e- j / h @0
add process p p > e+ e- j / h @1

I only get events containing zero jet and no events containing one jet. What am I doing wrong here? If I do not use my bias module everything works fine (i.e. I get a fraction of 0.175 containing one jet).

For NLO I use the same code in the cuts.f with the change of is_a_l(i) -> abs(ipdg(i)).eq.11 .or. abs(ipdg(i)).eq.13 .or. abs(ipdg(i)).eq.15.

Here, the production without jets works fine (generate p p > e+ e- [QCD]) and I get again a nice smooth spectrum. If I add a jet (generate p p > e+ e- j [QCD]) the generation fails "during the collection of results". It gives me two log files to check and in both log files I get the message:

Error 3 in kinematics_driver: negtive sqrt
-Infinity 100

It seems that after receiving this error 100 times the generation stops.

Concluding, it seems that my bias re-weighting code works fine for Drell-Yan without jets but have some issues when adding jets. Can you point me to the mistakes I am making and how a possible solution could look like. If needed I am happy to provide more code/files/resulting plots/etc.

best and thank you,
Jonas

Question information

Language:
English Edit question
Status:
For:
Assignee:
Paolo Torrielli Edit question
Last query:
2017-11-23
2017-11-27
 Olivier Mattelaer (olivier-mattelaer) said on 2017-10-26: #1

Hi,

This reminds me an old bug which was occuring only in presence of matching/merging.

For NLO, are you sure that your bias is infrared save? This might be the problem.

Cheers,

Olivier

 Jonas Roemer (jonasroemer) said on 2017-10-26: #2

Hi,

Oh sorry, forgot to mention it. I am using madgraph in the version 2.6.0. I guess that is the latest version, right?

For NLO: I hoped that taking the invariant mass of the two leptons would be infrared safe. But I am sadly far from being a generator or theory expert so I am not entirely sure about it.

best,
Jonas

 Jonas Roemer (jonasroemer) said on 2017-11-06: #3

Hi,

I tested some things like adding a constant term in the sqrt of the mass calculation to see if that sqrt gets negative but it did not solve the problem add all.

Do you have an idea why it does not work?

best,
Jonas

 Olivier Mattelaer (olivier-mattelaer) said on 2017-11-07: #4

when I print the bias_weight value from the function that you wrote, I have the following value:

4.4691152697180800E+306
4.7975613598866533E+306
5.1501265095405160E+306
5.5285804658568964E+306
5.9348227290446513E+306
Note: The following floating-point exceptions are signalling: IEEE_OVERFLOW_FLAG
STOP 1
6.3708920579627242E+306
6.8389766715568199E+306
7.3414251970146957E+306
7.8807584192600043E+306
8.4596818903892227E+306
9.0810994619334662E+306
9.7481278074355601E+306
1.0464112007771812E+307
1.1232642276896887E+307
1.2057571911424455E+307
1.2943036553496319E+307

One weird stuff is that such number are always increasing and so large.
Did you try to reset sump to 0 at the beginning of each call?

Cheers,

Olivier

> On Nov 6, 2017, at 12:44, Jonas Roemer <email address hidden> wrote:
>
> Question #660025 on MadGraph5_aMC@NLO changed:
>
> Jonas Roemer posted a new comment:
> Hi,
>
> I tested some things like adding a constant term in the sqrt of the mass
> calculation to see if that sqrt gets negative but it did not solve the
>
> Do you have an idea why it does not work?
>
> best,
> Jonas
>
> --

 Jonas Roemer (jonasroemer) said on 2017-11-08: #5

Hi Olivier,

Are these numbers from LO or NLO calculations. In fact, for LO, I did not set sump to 0 at the start, for NLO I fixed that in the meantime and set them to 0. I also added some debug output and found out a few things.

For LO: If I generate 0 jets and 1 jets (generate p p > e+ e- j / h @0, add process p p > e+ e- j / h @1) I find 2 leptons in the loop over all particles (do i=1,nexternal) for the process 0 (so the pp>ee) but 0 leptons for the process 1 (pp>eej). Also nexternals is 4 for the process 1 while it should be 5 (if i generate only pp>eej it is 5). Checking the SubProcesses/P1_qq_llg/nexternal.inc it also says 5 but for some reason this value is not imported. Looking into the bias code i find the line: include '../../nexternal.inc' which seems to include the nexternal.inc in the Source folder but not the SubProcesses/P* folder. May the wrong import be the reason and how do I provide the correct import to the nexternal.inc in the SubProcesses folder? I just copied the code from the ptj_bias code and did not think about it.

For NLO: I print the weight and the mass every time the bias_weight_function is called and look when the kinematics driver has negative a sqrt. The event generation is still running and I have only found a few instances of negative sqrts but so far nothing unusual seems to happen. The masses are in the range of 500-800 GeV and the weights are in the order of magnitude of 1e-2. I also find a lot of events in this mass/weight range without error.

I also calculated the expected range of weights. For 50GeV I would get 5.4e-07 and for 10TeV 5.0e8. At least after the addition of sump=0 this is also what I see in the debug output.

Thank You,
Jonas

 Olivier Mattelaer (olivier-mattelaer) said on 2017-11-08: #6

Hi,

I did not test this for LO, since my understanding is that you had it working in that case.
But indeed have different multiplicity would be problematic.
I thought that the nexternal.inc in Source was suppose to be set to the maximal value to avoid such kind of issue.
I need to check this, but this is likely to fix the problem.

So do I understand correctly that your NLO bias is working now?
If not do you know which sqrt is not working? is it the one that you define in your bias function or another one?

Cheers,

Olivier

 Jonas Roemer (jonasroemer) said on 2017-11-09: #7

Hi,

It sadly failed again (every run takes several hours). I do not know which sqrt fails but it is definitely not the one in the bias weight function. So I guess it is some internal sqrt. I also see no clear correlation between the bias_weight/mass and the failure of a sqrt. It also happens rarely - most events are generated without error. If it helps I can send you the log file for one job that failed. I append the exact bias code I am using now:

subroutine bias_weight_function(p,ipdg,bias_wgt)

implicit none
include 'nexternal.inc'
double precision bias_wgt,p(0:3,nexternal)
integer ipdg(nexternal),i
double precision sump(0:3)
double precision mass
double precision bias_normalization
double precision bias_exponent
double precision bias_switch
double precision bias_exponential
double precision bias_cutoff
logical do_mass

bias_normalization=3.64711d13
bias_exponent=4.29171d0
bias_switch=1.73642d3
bias_exponential=1.41843d-3
bias_cutoff=1d4

sump(0)=0d0
sump(1)=0d0
sump(2)=0d0
sump(3)=0d0
do_mass=.false.
do i=1,nexternal
if (abs(ipdg(i)).eq.11 .or. abs(ipdg(i)).eq.13 .or.
& abs(ipdg(i)).eq.15) then
do_mass=.true.
sump(0)=sump(0)+p(0,i)
sump(1)=sump(1)+p(1,i)
sump(2)=sump(2)+p(2,i)
sump(3)=sump(3)+p(3,i)
endif
enddo

if (.not.do_mass) then
bias_wgt=1d0
write(*,*) bias_wgt
return
endif

mass = sump(0)**2-sump(1)**2-sump(2)**2-sump(3)**2
if (mass.gt.0d0) then
mass=sqrt(mass)
else
mass=0d0
endif

bias_wgt=(mass**bias_exponent/bias_normalization)
if (mass.gt.bias_switch) then
bias_wgt=bias_wgt*exp((mass-bias_switch)*bias_exponential)
endif

write(*,*) bias_wgt, mass

return
end

Neither is do_mass false nor is mass=0d0 so both cases never seem to happen.

The commands to execute are:

generate p p > e+ e- j [QCD]
output DY_EE_1jhet_biased_10k_events -nojpeg
launch
done
set event_norm bias
set mll_sf 50 (did not work without this as well)
set nevents 10000
done

Thank You,
Jonas

 Olivier Mattelaer (olivier-mattelaer) said on 2017-11-23: #8

Looks like that the sqrt which goes negative is the one in mountecarlocounter.f (line 1469)
I will assign Paolo to this thread since he is the expert in counterterm.

Cheers,

Olivier

 Jonas Roemer (jonasroemer) said on 2017-11-23: #9

Hi,

I also looked in the last days into this issue again and I think I found a "fix" (run just completed yesterday evening). I now use pythia8 for showering (it was herwig before I guess since I did not change anything) and it went through without crashing. However, I was not able to check the distributions since madanalysis refuses to analyze the resulting lhe file.

I see the issue in the same file in the if(shower_mc.eq.'HERWIG6'.or.shower_mc.eq.'HERWIGPP') block, so switching to Pythia "fixed" the issue as the code is not called anymore. But I doubt that this is the desired solution for the problem.

Thank You,
Jonas

 Olivier Mattelaer (olivier-mattelaer) said on 2017-11-23: #10

Hi,

Depend actually of what the problem in that block is.
Paolo should be able to comment more (I hope)

Olivier

> On Nov 23, 2017, at 12:53, Jonas Roemer <email address hidden> wrote:
>
> Question #660025 on MadGraph5_aMC@NLO changed:
>
> Hi,
>
> I also looked in the last days into this issue again and I think I found
> a "fix" (run just completed yesterday evening). I now use pythia8 for
> showering (it was herwig before I guess since I did not change anything)
> and it went through without crashing. However, I was not able to check
> the distributions since madanalysis refuses to analyze the resulting lhe
> file.
>
> I see the issue in the same file in the
> if(shower_mc.eq.'HERWIG6'.or.shower_mc.eq.'HERWIGPP') block, so
> switching to Pythia "fixed" the issue as the code is not called anymore.
> But I doubt that this is the desired solution for the problem.
>
> Thank You,
> Jonas
>
> --

 Paolo Torrielli (paolotorriell) said on 2017-11-27: #11

Dear all,
as Olivier correctly pointed out, the error is caused by line 1469 in montecarlocounter.f, which represents the definition of the shower variable for Herwig. It should always be positive, so one should understand what is the quantity that drives qMCarg negative. For example, could you print the various values of 'w2' and of 'zeta2'? Once you know who is negative (or if it is zeta2 > 1), could you for example print the various quantities that enter function get_zeta, in order to understand which is the one responsible for a negative square root?
Thanks.
Cheers.
Paolo