Manual Setting of Matrix Element

Asked by James Cockburn

Hi,

I want to do something a bit odd within MadGraph - basically, what I want to do is have it such that the colour/helicity summed/averaged matrix element returns the value 1, no matter what momentum is fed to the function. I was told this could be done by going into dsample.f and commenting out the call to the matrix element (which I believe is this dsig function) but then MadGraph tells me that I have a zero cross section. This function is stored in a variable called fx, so instead I try manually setting fx=1 and this time I get an exception (IEEE_INVALID_FLAG). I then tried going into the auto_dsig.f file and manually have that return 1 (or actually, just the pdf convolution) but the results that come out of this don't seem right. I also generated the standalone code and it seems like, given some phase space point, I get different results back from the standalone and the full code. In the latter, I am simply printing out the momentum passed into the function defined in matrix1.f and the value it spits out and comparing this to what I get in the standalone. The different between them is not a constant factor, so there is no coupling factor differences (I am using a fixed scale) so I have no idea what on Earth is going on...any help appreciated.

Cheers,
James

p.s. For reference, I am running a g g > g g process.

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Solved by:
Olivier Mattelaer
Solved:
Last query:
Last reply:
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#1

Dear James,

Indeed it looks like the dsample.f is not the correct file to edit.
the correct file to edit is matrix1.f

you need to add the line:
ANS=DSIGN(1d0*4*81, ANS)
before the line
IF (MULTI_CHANNEL) THEN

(if you have matrix2.f /… you have to do that for each file)

> . I also generated the standalone code and it seems like, given some phase space point, I get different results back from the standalone and the full code.

Yes this is normal due to the multi-channel method of integration (see hep-ph/0208156)

Cheers,

Olivier

> On Dec 2, 2015, at 17:22, James Cockburn <email address hidden> wrote:
>
> New question #276439 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/276439
>
> Hi,
>
> I want to do something a bit odd within MadGraph - basically, what I want to do is have it such that the colour/helicity summed/averaged matrix element returns the value 1, no matter what momentum is fed to the function. I was told this could be done by going into dsample.f and commenting out the call to the matrix element (which I believe is this dsig function) but then MadGraph tells me that I have a zero cross section. This function is stored in a variable called fx, so instead I try manually setting fx=1 and this time I get an exception (IEEE_INVALID_FLAG). I then tried going into the auto_dsig.f file and manually have that return 1 (or actually, just the pdf convolution) but the results that come out of this don't seem right. I also generated the standalone code and it seems like, given some phase space point, I get different results back from the standalone and the full code. In the latter, I am simply printing out the momentum passed into the function defined in matrix1.f and the value it spits out and comparing this to what I get in the standalone. The different between them is not a constant factor, so there is no coupling factor differences (I am using a fixed scale) so I have no idea what on Earth is going on...any help appreciated.
>
> Cheers,
> James
>
> p.s. For reference, I am running a g g > g g process.
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
James Cockburn (j-d-cockburn) said :
#2

Hi Olivier,

Thanks...I'll add that line. If I want to check that this is absolutely, definitely, 100% giving me 1 for the ME in MadGraph, what is the best way to do that? I would assume that there is a variable somewhere in the code which stores the value of the ME at each point, so I'd like to print that out somewhere.

Also, this factor seems to just be a random set of numbers, what is the reason behind them? Is this to do with the MULTI_CHANNEL integration? (Thanks for the reference on that, btw). If i were to generate, say, g g > g g g or u d > u g d or anything else, would this factor stay the same, or is this a process-dependent factor you gave me?

Cheers,
James

Revision history for this message
James Cockburn (j-d-cockburn) said :
#3

Also, there are two places in matrix1.f where we have this 'IF MULTI_CHANNEL THEN' line, so where am I having to put this extra line precisely?

- James

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#4

Hi James,

> I would assume that there is a variable somewhere in the
> code which stores the value of the ME at each point, so I'd like to
> print that out somewhere.

I’m not aware of such variable.

> Also, this factor seems to just be a random set of numbers, what is the
> reason behind them?

The idea is to split the integral in a sum of simpler integral which are easier (i.e. faster) to integrate.
What you call random number are just the weight corresponding to the splitting over the phase-space between the various integral.

> Is this to do with the MULTI_CHANNEL integration?

yes

> (Thanks for the reference on that, btw). If i were to generate, say, g g
>> g g g or u d > u g d or anything else, would this factor stay the
> same, or is this a process-dependent factor you gave me?

If you look at equation 2.3 to 2.5 of the reference, you will see the factor that is used.
and yes this factor is process dependent.

> Also, there are two places in matrix1.f where we have this 'IF
> MULTI_CHANNEL THEN' line, so where am I having to put this extra line
> precisely?

before this one:

      IF (MULTI_CHANNEL) THEN
        XTOT=0D0
        DO I=1,NDIAGS
          XTOT=XTOT+AMP2(I)
        ENDDO
        IF (XTOT.NE.0D0) THEN
          ANS=ANS*AMP2(SUBDIAG(1))/XTOT
        ELSE
          ANS=0D0
        ENDIF
      ENDIF

You can see here the factor that are described in Equation 2.5.
You can also see that ANS is directly modified by the multi-channel weight and that we do not keep in memory the “real” value of the matrix element.

Cheers,

Olivier

> On Dec 3, 2015, at 10:52, James Cockburn <email address hidden> wrote:
>
> Question #276439 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/276439
>
> Status: Answered => Open
>
> James Cockburn is still having a problem:
> Hi Olivier,
>
> Thanks...I'll add that line. If I want to check that this is absolutely,
> definitely, 100% giving me 1 for the ME in MadGraph, what is the best
> way to do that? I would assume that there is a variable somewhere in the
> code which stores the value of the ME at each point, so I'd like to
> print that out somewhere.
>
> Also, this factor seems to just be a random set of numbers, what is the
> reason behind them? Is this to do with the MULTI_CHANNEL integration?
> (Thanks for the reference on that, btw). If i were to generate, say, g g
>> g g g or u d > u g d or anything else, would this factor stay the
> same, or is this a process-dependent factor you gave me?
>
> Cheers,
> James
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
James Cockburn (j-d-cockburn) said :
#5

Hi Olivier,

>> I would assume that there is a variable somewhere in the
>> code which stores the value of the ME at each point, so I'd like to
>> print that out somewhere.

>I’m not aware of such variable.

Alright...if I understand this multi-channel thing correctly, then my total amplitude should be equal to the sum of these fi. However, what I don't understand is where MG is calculating these fi... presumably it's doing it at ANS=ANS*AMP2(SUBDIAG(1))/XTOT, but then what is this SUBDIAG thing? And why are we only looking at SUBDIAG(1) and not all possible SUBDIAG(J) or something like this?

I guess this question boils down to one thing - this function called SMATRIX1 in matrix1.f, what exactly is it outputting? The comments underneath it say 'Returns amplitude squared summed/avg over colors and helicities' but clearly this is not true because then I could have just set ANS=1 and have been done with it.

> If you look at equation 2.3 to 2.5 of the reference, you will see the factor that is used.
and yes this factor is process dependent.

I see that SOME factor is used, yes...but I don't see why, for gg > gg for example, I can 'undo' it by multiplying by 4*81. Or am I misunderstanding? And how can I work out what factor to use for a general process with a well-defined (i.e.. not proton-proton) incoming state?

Cheers,
James

Revision history for this message
James Cockburn (j-d-cockburn) said :
#6

Has this 'suppress amplitude' bool got something to do with this? I found it in driver.f. I tried forcing it to return false, and I can generate things but also get a few 'STOP 1' messages printed to my screen

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#7

Hi,

> Alright...if I understand this multi-channel thing correctly, then my
> total amplitude should be equal to the sum of these fi.

Correct

> However, what I
> don't understand is where MG is calculating these fi... presumably it's
> doing it at ANS=ANS*AMP2(SUBDIAG(1))/XTOT,

Correct

> but then what is this SUBDIAG
> thing? And why are we only looking at SUBDIAG(1) and not all possible
> SUBDIAG(J) or something like this?

SUBDIAG(1) correspond to the variable i of the sum.

> and not all possible
> SUBDIAG(J) or something like this?

This is (I guess) an historical reason.
The current way to run the code is to run N times the executable and for each of those initialize SUBDIAG(1) to a value between 1 and N.
I guess that in the past, only one executable was used and then the code was selecting in a (smart) random way the value between 1 and N.
and therefore at that time it was useful to have SUBDIAG as a list.

> I guess this question boils down to one thing - this function called
> SMATRIX1 in matrix1.f, what exactly is it outputting? The comments
> underneath it say 'Returns amplitude squared summed/avg over colors and
> helicities' but clearly this is not true because then I could have just
> set ANS=1 and have been done with it.

The real output is the function to integrate.
The function to integrate depends of the run_option
It can contain the “suppress amplitude factor” (I prefer to call this factor the weight associate to the multi-channel)
It can also do (or not) the average over helicity/… depending of various options mode.

> but I don't see why, for gg > gg
> for example, I can ‘undo' it by multiplying by 4*81.

Oh that factor is just the factor for the helicity and color factor (should be 4*64 actually, sorry for the confusion).
You can include it or not depending if you want the matrix element or the average matrix element to be normalized to one.

> Has this 'suppress amplitude' bool got something to do with this? I
> found it in driver.f. I tried forcing it to return false, and I can
> generate things but also get a few 'STOP 1' messages printed to my
> screen

I will not do that since then the computation is clearly wrong by default.
The only reason of this option is to allow to compare the total cross-section for a given Feynman diagram when debugging the code and running each integral by hand.
So we basically never run with this mode on False and can not really comment on it.

Cheers,

Olivier

> On Dec 3, 2015, at 11:42, James Cockburn <email address hidden> wrote:
>
> Question #276439 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/276439
>
> Status: Answered => Open
>
> James Cockburn is still having a problem:
> Hi Olivier,
>
>>> I would assume that there is a variable somewhere in the
>>> code which stores the value of the ME at each point, so I'd like to
>>> print that out somewhere.
>
>> I’m not aware of such variable.
>
> Alright...if I understand this multi-channel thing correctly, then my
> total amplitude should be equal to the sum of these fi. However, what I
> don't understand is where MG is calculating these fi... presumably it's
> doing it at ANS=ANS*AMP2(SUBDIAG(1))/XTOT, but then what is this SUBDIAG
> thing? And why are we only looking at SUBDIAG(1) and not all possible
> SUBDIAG(J) or something like this?
>
> I guess this question boils down to one thing - this function called
> SMATRIX1 in matrix1.f, what exactly is it outputting? The comments
> underneath it say 'Returns amplitude squared summed/avg over colors and
> helicities' but clearly this is not true because then I could have just
> set ANS=1 and have been done with it.
>
>> If you look at equation 2.3 to 2.5 of the reference, you will see the factor that is used.
> and yes this factor is process dependent.
>
> I see that SOME factor is used, yes...but I don't see why, for gg > gg
> for example, I can 'undo' it by multiplying by 4*81. Or am I
> misunderstanding? And how can I work out what factor to use for a
> general process with a well-defined (i.e.. not proton-proton) incoming
> state?
>
> Cheers,
> James
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
James Cockburn (j-d-cockburn) said :
#8

Okay so now I feel like I've gone around in a circle a bit now -

If that factor is just the helicity average and colour factor, and I want is such that the matrix element summed/averaged over helicity/colour =1, then I should just manually set ANS=1 before the 'if multi_channel" bit, right?

But then this multi_channel code multiplies by a factor, which means that we instead output one of these fi things which is what we integrate. If we do this for all fis, then we get the integral of the actual function we wanted, i.e. the amplitude.

So this means that actually I should set the output of the smatrix function to be just 1/N, where N is the number of fi functions we want, right? Otherwise, we are weighting by the amplitude of a given Feynman diagram which we haven't been able to manually change meaning that the sum of the fi will not be 1.

You know, it would be really really really useful if there was a way of running MadGraph without this multi-channel thing. Is this possible in, say, a previous version?

Cheers,
James

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#9

Hi,

> If that factor is just the helicity average and colour factor, and I
> want is such that the matrix element summed/averaged over
> helicity/colour =1, then I should just manually set ANS=1 before the 'if
> multi_channel” bit, right?

No because we divide by those number after the multi-channel via the line:
      ANS=ANS/DBLE(IDEN)
This line also contains the symmetry factor.

> But then this multi_channel code multiplies by a factor, which means
> that we instead output one of these fi things which is what we
> integrate. If we do this for all fis, then we get the integral of the
> actual function we wanted, i.e. the amplitude.

correct

> So this means that actually I should set the output of the smatrix
> function to be just 1/N, where N is the number of fi functions we want,
> right? Otherwise, we are weighting by the amplitude of a given Feynman
> diagram which we haven't been able to manually change meaning that the
> sum of the fi will not be 1.

I do not agree on that statement.
The only reason to return 1/N is if you do not include the weight of the multi-channel.
Since then you perform N times the same integral (using N different parametrization of the phase-space).
but if you include the weight of the multi-channel you certainly do not need the factor 1/N since those weight are automatically normalize to one by construction.

> You know, it would be really really really useful if there was a way of
> running MadGraph without this multi-channel thing. Is this possible in,
> say, a previous version?

As I said we can do that in debug mode. Here is how to do it:
Go to one of the Subprocesses/PXXXX/GXXX directory (Please do not use a GXXX which as a dot in his name)
and change the input_app.txt and setting suppress amplitude to False.
then run ../madevent < input_app.txt
Then you run for one given parametrization of the phase-space without this multi-channel.
and will have the associate results.
Off course this is not testing the full phase-space Volume but only testing one of the (many) channel of integration so not really what you ask for.

I will never make this available in a full automated way, since using only one change of variable is a terrible idea and you will have a lot of trouble of convergence if you do so.

Cheers,

Olivier

> On Dec 3, 2015, at 17:03, James Cockburn <email address hidden> wrote:
>
> Question #276439 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/276439
>
> Status: Answered => Open
>
> James Cockburn is still having a problem:
> Okay so now I feel like I've gone around in a circle a bit now -
>
> If that factor is just the helicity average and colour factor, and I
> want is such that the matrix element summed/averaged over
> helicity/colour =1, then I should just manually set ANS=1 before the 'if
> multi_channel" bit, right?
>
> But then this multi_channel code multiplies by a factor, which means
> that we instead output one of these fi things which is what we
> integrate. If we do this for all fis, then we get the integral of the
> actual function we wanted, i.e. the amplitude.
>
> So this means that actually I should set the output of the smatrix
> function to be just 1/N, where N is the number of fi functions we want,
> right? Otherwise, we are weighting by the amplitude of a given Feynman
> diagram which we haven't been able to manually change meaning that the
> sum of the fi will not be 1.
>
> You know, it would be really really really useful if there was a way of
> running MadGraph without this multi-channel thing. Is this possible in,
> say, a previous version?
>
> Cheers,
> James
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
James Cockburn (j-d-cockburn) said :
#10

Hi Olivier,

>> So this means that actually I should set the output of the smatrix
>> function to be just 1/N, where N is the number of fi functions we want,
>> right? Otherwise, we are weighting by the amplitude of a given Feynman
>> diagram which we haven't been able to manually change meaning that the
>> sum of the fi will not be 1.

>I do not agree on that statement.
>The only reason to return 1/N is if you do not include the weight of the multi-channel.
>Since then you perform N times the same integral (using N different parametrization of the phase-space).
>but if you include the weight of the multi-channel you certainly do not need the factor 1/N since those weight are automatically >normalize to one by construction.

This was the key point that I was misunderstanding. I see now (eq 2.5 in the reference) that what you say is true. Apologies for that.

Thank you so much for all of your help so far, but I have one final question if I may. I see that towards the end of section 2 of the reference you gave me, there is a sentence '...when some diagrams do not display any peaked behaviour and are dropped from equation 2.5'. Does this mean that MG actually approximates the matrix element to some extent? The way it reads makes it sound like you generate all of these 'channels' from the Feynman diagrams (which I agree can be fewer than the total number of diagrams), check for peaked behaviour in the channels and disregard the channels that don't have it. I think this might be an okay thing to do if the Atot in equation 2.5 remains the actual, physical Atot formed by consideration of all the diagrams, but the notation of 2.5 suggests that this is not the case and it is actually formed from the sum of all the channels. So if you neglect one channel, you don't integrate the full Atot, only a portion of it. Do I have this right?

Cheers,
James

Revision history for this message
Best Olivier Mattelaer (olivier-mattelaer) said :
#11

> Thank you so much for all of your help so far, but I have one final
> question if I may. I see that towards the end of section 2 of the
> reference you gave me, there is a sentence '...when some diagrams do not
> display any peaked behaviour and are dropped from equation 2.5'. Does
> this mean that MG actually approximates the matrix element to some
> extent? The way it reads makes it sound like you generate all of these
> 'channels' from the Feynman diagrams (which I agree can be fewer than
> the total number of diagrams), check for peaked behaviour in the
> channels and disregard the channels that don't have it. I think this
> might be an okay thing to do if the Atot in equation 2.5 remains the
> actual, physical Atot formed by consideration of all the diagrams,

You are correct the Atot in equation 2.5 is not (can not be) touched and always remains the full (physical) amplitude square.

> but
> the notation of 2.5 suggests that this is not the case and it is
> actually formed from the sum of all the channels. So if you neglect one
> channel, you don't integrate the full Atot, only a portion of it. Do I
> have this right?

Looks like the notation is confusing then. Modify Atot will be terrible since it will induce bias since equation 2.3 will not be valid anymore.

Cheers,

Olivier

> On Dec 4, 2015, at 11:27, James Cockburn <email address hidden> wrote:
>
> Question #276439 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/276439
>
> Status: Answered => Open
>
> James Cockburn is still having a problem:
> Hi Olivier,
>
>>> So this means that actually I should set the output of the smatrix
>>> function to be just 1/N, where N is the number of fi functions we want,
>>> right? Otherwise, we are weighting by the amplitude of a given Feynman
>>> diagram which we haven't been able to manually change meaning that the
>>> sum of the fi will not be 1.
>
>> I do not agree on that statement.
>> The only reason to return 1/N is if you do not include the weight of the multi-channel.
>> Since then you perform N times the same integral (using N different parametrization of the phase-space).
>> but if you include the weight of the multi-channel you certainly do not need the factor 1/N since those weight are automatically >normalize to one by construction.
>
> This was the key point that I was misunderstanding. I see now (eq 2.5 in
> the reference) that what you say is true. Apologies for that.
>
> Thank you so much for all of your help so far, but I have one final
> question if I may. I see that towards the end of section 2 of the
> reference you gave me, there is a sentence '...when some diagrams do not
> display any peaked behaviour and are dropped from equation 2.5'. Does
> this mean that MG actually approximates the matrix element to some
> extent? The way it reads makes it sound like you generate all of these
> 'channels' from the Feynman diagrams (which I agree can be fewer than
> the total number of diagrams), check for peaked behaviour in the
> channels and disregard the channels that don't have it. I think this
> might be an okay thing to do if the Atot in equation 2.5 remains the
> actual, physical Atot formed by consideration of all the diagrams, but
> the notation of 2.5 suggests that this is not the case and it is
> actually formed from the sum of all the channels. So if you neglect one
> channel, you don't integrate the full Atot, only a portion of it. Do I
> have this right?
>
> Cheers,
> James
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
James Cockburn (j-d-cockburn) said :
#12

Thanks Olivier Mattelaer, that solved my question.