EFT: generate only interfence term?

Dear Sebastien,

This has been introduce in the latest version of the code (2.2.0).
We have introduce a new syntax
in addition of the usual order restriction
QED(<)=2 QCD(<)=2 NP(<)=2
which apply at the diagram level.
you have a new syntax
QED^2
QCD^2
NP^2
which are constraint which apply on the diagram square level.
For those restriction you have a higher variety of operator allowed
NP^2<=2
NP^2=2 (which like in the previous syntax is equivalent to <= for historical reason)
NP^2==2 (which is the real equality)
NP^2 >2

Since you want to have only
> 2 Re(SM* operator).

This is typically
NP^2==2
that you want.

> As this interference term can be negative, some events would have an associated weight = -1

Yes that will be the case.

Cheers,

Olivier

On Oct 6, 2014, at 12:31 PM, Sébastien Wertz <email address hidden> wrote:

> New question #255413 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Hello!
>
> I'd like to study the TopEffTh effective theory model implemented in MG, up to order 1/Lambda^2.
>
> To do so, it would be very convenient if I could generate events, not according to |M|^2 = |SM|^2 + 2 Re(SM* operator) + |operator|^2, but only according to the interference term 2 Re(SM* operator).
>
> Doing this for one operator at a time, my signal would be completely linear in the coupling of this operator, which would make like much easier for me.
>
> As this interference term can be negative, some events would have an associated weight = -1. Alternatively, I could use a weighted distribution of events (skipping the unweighting step in the generation), where some of the weights would be negative. To achieve this, is there a way to fiddle with how MG writes out the (averaged, summed) matrix element?
>
> If not, another solution could be to compute, for each PS-point, the value of the full |M|^2, and to subtract the value of |SM|^2 in order to approximate the interference term (as computing both matrix elements using standard MG output seems more straightforward). Of course, if possible, it would be even better to also subtract the |operator|^2 term (I don't know if it possible to write an amplitude which has only operator-diagrams, without the SM ones).
>
> Looking forward to hearing from you,
>
> Sébastien
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

That is great, thank you for the quick answer!

Sébastien

Actually I still have some questions...

I'm trying to generate plain and simple p p > t t~, with, as you said, NP^2==2, but using fully-leptonic top decays (because I need the spin correlations and don't want to generate a hundred million events to keep only 20 000).

I then get the message:

Error detected in sub-command generate p p > t t~ NP^2==2 , ( t > b W+ , W+ > vl l+ ) , ( t~ > b~ W- , W- > vl~ l- )
write debug file MG5_debug
If you need help with this issue please contact us on https://answers.launchpad.net/madgraph5
MadGraph5Error : Decay processes cannot specify squared orders constraints.

Same thing if I try: p p > t t~ NP^2==2 , ( t > b W+ NP=0 , W+ > vl l+ ) , ( t~ > b~ W- NP=0 , W- > vl~ l- )

Weird thing is, p p > t t~ , ( t > b W+ NP=0 , W+ > vl l+ ) , ( t~ > b~ W- NP=0 , W- > vl~ l- ) NP^2==2
doesn't output any error message but produdes strange results, with diagrams of order NP=6, no SM diagrams (needed for the interference), and so on...

So, there is at the moment no workaround to specify the order at the matrix-element-squared-level if there are decays involved?

Cheers,

Sébastien

Hi,

No indeed this is not implemented (at least for the moment). If you do not want anything special in the decay then one workaround is probably to use MadSpin (arXiv:1212.3460).
I never think about such combination before, such you will need to make some dirty modification:

One way to make it work is
1) to first run
p p > t t~ NP^2==2
and generate your sample.
2) edit the lhe file and change the proc_card information to claim to madspin that is a
p p > t t~ sample (i.e. pass to p p > t t~ NP=2)
3) run
./bin/madevent
and type
decay_events
and use MadSpin as normal.

Cheers,

Olivier

On Oct 6, 2014, at 4:06 PM, Sébastien Wertz <email address hidden> wrote:

> Question #255413 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Status: Solved => Open
>
> Sébastien Wertz is still having a problem:
> Actually I still have some questions...
>
> I'm trying to generate plain and simple p p > t t~, with, as you said,
> NP^2==2, but using fully-leptonic top decays (because I need the spin
> correlations and don't want to generate a hundred million events to keep
> only 20 000).
>
> I then get the message:
>
> Error detected in sub-command generate p p > t t~ NP^2==2 , ( t > b W+ , W+ > vl l+ ) , ( t~ > b~ W- , W- > vl~ l- )
> write debug file MG5_debug
> If you need help with this issue please contact us on https://answers.launchpad.net/madgraph5
> MadGraph5Error : Decay processes cannot specify squared orders constraints.
>
> Same thing if I try: p p > t t~ NP^2==2 , ( t > b W+ NP=0 , W+ > vl l+ )
> , ( t~ > b~ W- NP=0 , W- > vl~ l- )
>
> Weird thing is, p p > t t~ , ( t > b W+ NP=0 , W+ > vl l+ ) , ( t~ > b~ W- NP=0 , W- > vl~ l- ) NP^2==2
> doesn't output any error message but produdes strange results, with diagrams of order NP=6, no SM diagrams (needed for the interference), and so on...
>
> So, there is at the moment no workaround to specify the order at the
> matrix-element-squared-level if there are decays involved?
>
> Cheers,
>
> Sébastien
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hello,

Yes, MadSpin would indeed solve my problem! (it this the kind of situation one could hope to see taken care of in the next update? :) )

No, I do not need anything special in the decays (just the regular SM top decays, ... for now!), except that I want to restrict to the fully-leptonic case.

Would there be a way to implement the solution you've described before launching a multi_run? I'm afraid I'll have to run MadSpin and Pythia on each sub-directory of events separately (Writing a script would do, but if I can avoid it...), since the multi_run would then only take care of the partonic level...

Thanks again for your great help,

Sébastien

Hi,

> Yes, MadSpin would indeed solve my problem! (it this the kind of
> situation one could hope to see taken care of in the next update? :) )

I’m not sure yet to be honest, such method can introduce small biais and if automate the method describe here, the user will not be
aware of such problem. So I’m not yet decide of what is the correct approach here.

> No, I do not need anything special in the decays (just the regular SM
> top decays, ... for now!), except that I want to restrict to the fully-
> leptonic case.

Yes no problem on that.

> Would there be a way to implement the solution you've described before
> launching a multi_run? I'm afraid I'll have to run MadSpin and Pythia on
> each sub-directory of events separately (Writing a script would do, but
> if I can avoid it...), since the multi_run would then only take care of
> the partonic level…

or a simple way to go around this is to modify the line 2529 of MadSpin/decay.py
        # 2. compute the production matrix element -----------------------------
        processes = [line[9:].strip() for line in self.banner.proc_card
                     if line.startswith('generate')]
        processes += [' '.join(line.split()[2:]) for line in self.banner.proc_card
                      if re.search(‘^\s*add\s+process', line)]
and adding after those the following line:
 processes = [l.replace(’NP^2==2’, ’NP2=2’) for l in processes]

Note that I have not tested the above proposal at all…

Cheers,

Olivier

On Oct 6, 2014, at 5:46 PM, Sébastien Wertz <email address hidden> wrote:

> Question #255413 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Status: Answered => Open
>
> Sébastien Wertz is still having a problem:
> Hello,
>
> Yes, MadSpin would indeed solve my problem! (it this the kind of
> situation one could hope to see taken care of in the next update? :) )
>
> No, I do not need anything special in the decays (just the regular SM
> top decays, ... for now!), except that I want to restrict to the fully-
> leptonic case.
>
> Would there be a way to implement the solution you've described before
> launching a multi_run? I'm afraid I'll have to run MadSpin and Pythia on
> each sub-directory of events separately (Writing a script would do, but
> if I can avoid it...), since the multi_run would then only take care of
> the partonic level...
>
> Thanks again for your great help,
>
> Sébastien
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

Following our e-mail conversation, I tried the following:

- Output the full matrix element, with decays, using the constrain NP=2
- For each subprocess, identify (using the matrix1.ps file) the diagrams which do not contain effective operators
- Edit each matrix1.f file in the following way:

Under "C LOCAL VARIABLES ":
* Declare INTEGERS X,Y
* Replace COMPLEX*16 JAMP(NCOLOR,NAMPSO,2) by COMPLEX*16 JAMP(NCOLOR,NAMPSO,2,2)

Under "C JAMPs contributing to orders ALL_ORDERS=1":
* Modify the definitions of the JAMP(.,.) so that JAMP(.,.,1) contains only the SM diagrams, and JAMP(.,.,2) the other ones. In the end, JAMP(.,.,1) + JAMP(.,.,2) = the old JAMP(.,.)

Then, in the loop computing MATRIX1:
* Also loop on X,Y = 0 ,1
* Replace ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M) by ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M,X)
* Replace MATRIX1=MATRIX1+ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
by:
IF ((X+Y).EQ.3) THEN
    MATRIX1=MATRIX1+ZTEMP*DCONJG(JAMP(I,N,Y))/DENOM(I)
ENDIF

That way, I sum only the product terms containing one SM diagram and one operator-diagram, giving me a matrix-element-squared of order 1/Lambda².

Here comes the part I'm not sure about: There is a second loop computing JAMP2 from the different diagrams. I don't understand how this quantity is used. As the loop uses JAMP(I,N) I replaced it by JAMP(I,N,1)+JAMP(I,N,2) and that way it seems to work.

The events and the cross-section I get seem to make sense, although I'm a bit surprised that all my events have positive weight. I'll check that by using a negative operator coefficient...

Can you please confirm that this is the right way to go?

* There's a typo in my comment: I don't loop over X,Y = 0 ,1 but over X,Y = 1,2 or course.... (why can't I edit previous comments?)

Seems good!

JAMP2 defines the probability to assign an event to a given color flow (which is important for the shower)

In Valentin implementation they are computed like:
      DO I = 1, NCOLOR
        DO M = 1, NAMPSO
          DO N = 1, NAMPSO
            IF (CHOSEN_SO_CONFIGS(SQSOINDEX1(M,N))) THEN
              JAMP2(I)=JAMP2(I)+ABS(JAMP(I,M)*DCONJG(JAMP(I,N)))
            ENDIF
          ENDDO
        ENDDO
      ENDDO

Where M,N plays the role of your X,Y.

Cheers,

Olivier

On Oct 8, 2014, at 9:11 AM, Sébastien Wertz <email address hidden> wrote:

> Question #255413 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Sébastien Wertz posted a new comment:
> Hi Olivier,
>
> Following our e-mail conversation, I tried the following:
>
> - Output the full matrix element, with decays, using the constrain NP=2
> - For each subprocess, identify (using the matrix1.ps file) the diagrams which do not contain effective operators
> - Edit each matrix1.f file in the following way:
>
> Under "C LOCAL VARIABLES ":
> * Declare INTEGERS X,Y
> * Replace COMPLEX*16 JAMP(NCOLOR,NAMPSO,2) by COMPLEX*16 JAMP(NCOLOR,NAMPSO,2,2)
>
> Under "C JAMPs contributing to orders ALL_ORDERS=1":
> * Modify the definitions of the JAMP(.,.) so that JAMP(.,.,1) contains only the SM diagrams, and JAMP(.,.,2) the other ones. In the end, JAMP(.,.,1) + JAMP(.,.,2) = the old JAMP(.,.)
>
> Then, in the loop computing MATRIX1:
> * Also loop on X,Y = 0 ,1
> * Replace ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M) by ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M,X)
> * Replace MATRIX1=MATRIX1+ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
> by:
> IF ((X+Y).EQ.3) THEN
> MATRIX1=MATRIX1+ZTEMP*DCONJG(JAMP(I,N,Y))/DENOM(I)
> ENDIF
>
> That way, I sum only the product terms containing one SM diagram and one
> operator-diagram, giving me a matrix-element-squared of order 1/Lambda².
>
> Here comes the part I'm not sure about: There is a second loop computing
> JAMP2 from the different diagrams. I don't understand how this quantity
> is used. As the loop uses JAMP(I,N) I replaced it by
> JAMP(I,N,1)+JAMP(I,N,2) and that way it seems to work.
>
> The events and the cross-section I get seem to make sense, although I'm
> a bit surprised that all my events have positive weight. I'll check that
> by using a negative operator coefficient...
>
> Can you please confirm that this is the right way to go?
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi,

I'm sorry but I'm not sure how I should modify this section for my particular case.

Given an amplitude A = A1 + A2 + A3, where the Ai have different color structures (ie JAMP(i,.,.) here), JAMP2(i) is defined as |Ai|^2, is that it?

So, I should write:
DO X = 1, 2
  DO Y = 1, 2
    DO I = 1, NCOLOR
      DO M = 1, NAMPSO
         DO N = 1, NAMPSO
          IF (CHOSEN_SO_CONFIGS(SQSOINDEX1(M,N))) THEN
            IF ((X+Y).EQ.3) THEN
              JAMP2(I)=JAMP2(I)+ABS(JAMP(I,M,X)*DCONJG(JAMP(I,N,Y)))
           ENDIF
         ENDIF
       ENDDO
     ENDDO
   ENDDO
  ENDDO
ENDDO

2 questions regarding this (we're almost done!):
- what are the AMP2(.) defined just before this loop? They're not used to compute JAMP2
- what are CHOSEN_SO_CONFIGS() and SQSOINDEX1() ? If I understand right they're used to restrict the order of the |M|^2 (ie apply the NP^2==2 rule)? So I could just remove everyting related to that, that is, the loops over M and N in the computation of JAMP and JAMP2?

Thanks again!

Sébastien

Hi,

No your “X” is the “N” one in the code that I quote
and your “Y” is the “M” one in the code that I quote

So you should write

> DO X = 1, 2
> DO Y = 1, 2
> DO I = 1, NCOLOR
> IF ((X+Y).EQ.3) THEN
> JAMP2(I)=JAMP2(I)+ABS(JAMP(I,X)*DCONJG(JAMP(I,Y)))
> ENDIF
> ENDDO
> ENDDO
> ENDDO

> - what are the AMP2(.) defined just before this loop? They’re not used to compute JAMP2

This are number for the phase-space integration. You can change to any value (as long as the sum is different of zero) and it will only change the efficiency of the computation (i.e. the numerical error on cross-section/time) but not the result.
They are not really good value to put here for the interference computation but the current seems to works well.

> - what are CHOSEN_SO_CONFIGS() and SQSOINDEX1() ? If I understand right they're used to restrict the order of the |M|^2 (ie apply the NP^2==2 rule)? So I could just remove everyting related to that, that is, the loops over M and N in the computation of JAMP and JAMP2?

I do not know those details, but they play the role of your " IF ((X+Y).EQ.3) THEN”

Cheers,

Olivier

On Oct 8, 2014, at 10:31 AM, Sébastien Wertz <email address hidden> wrote:

> Question #255413 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Status: Answered => Open
>
> Sébastien Wertz is still having a problem:
> Hi,
>
> I'm sorry but I'm not sure how I should modify this section for my
> particular case.
>
> Given an amplitude A = A1 + A2 + A3, where the Ai have different color
> structures (ie JAMP(i,.,.) here), JAMP2(i) is defined as |Ai|^2, is that
> it?
>
> So, I should write:
> DO X = 1, 2
> DO Y = 1, 2
> DO I = 1, NCOLOR
> DO M = 1, NAMPSO
> DO N = 1, NAMPSO
> IF (CHOSEN_SO_CONFIGS(SQSOINDEX1(M,N))) THEN
> IF ((X+Y).EQ.3) THEN
> JAMP2(I)=JAMP2(I)+ABS(JAMP(I,M,X)*DCONJG(JAMP(I,N,Y)))
> ENDIF
> ENDIF
> ENDDO
> ENDDO
> ENDDO
> ENDDO
> ENDDO
>
> 2 questions regarding this (we're almost done!):
> - what are the AMP2(.) defined just before this loop? They're not used to compute JAMP2
> - what are CHOSEN_SO_CONFIGS() and SQSOINDEX1() ? If I understand right they're used to restrict the order of the |M|^2 (ie apply the NP^2==2 rule)? So I could just remove everyting related to that, that is, the loops over M and N in the computation of JAMP and JAMP2?
>
> Thanks again!
>
> Sébastien
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

I did some cross-checks and it seems to work: generating a sample of events with just one operator set to 1 gives me cross-section X, and setting the operator to -1 gives a cross-section -X and negative event weights!

Of course these are only preliminary studies and I will really need to find a way to automatise this procedure (Modifying the ME by hand is fine when there are no more than 10 diagrams per sub-process, and 5 or 6 subprocesses... But for instance in the case of p p > t t~ j there are already about 900 diagrams at LO. And we're not event talking about NLO...). I'll see what is possible to do with the Louvain MG experts.

Thank you for patience,

Sébastien

Hi,

> I’ll see what is possible to do with the Louvain MG experts.

The real expert on this is Valentin Hirschi (currently at SLAC)

Cheers,

Olivier

On Oct 9, 2014, at 2:46 PM, Sébastien Wertz <email address hidden> wrote:

> Question #255413 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/255413
>
> Status: Answered => Solved
>
> Sébastien Wertz confirmed that the question is solved:
> Hi Olivier,
>
> I did some cross-checks and it seems to work: generating a sample of
> events with just one operator set to 1 gives me cross-section X, and
> setting the operator to -1 gives a cross-section -X and negative event
> weights!
>
> Of course these are only preliminary studies and I will really need to
> find a way to automatise this procedure (Modifying the ME by hand is
> fine when there are no more than 10 diagrams per sub-process, and 5 or 6
> subprocesses... But for instance in the case of p p > t t~ j there are
> already about 900 diagrams at LO. And we're not event talking about
> NLO...). I'll see what is possible to do with the Louvain MG experts.
>
> Thank you for patience,
>
> Sébastien
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.