Use particle id in cuts.f

Hi,

First if you use a recent enough version of MG5aMC,
You actually have a couple of cut that you can define by pdg directly inside the run_card
(pt, rapidity and invariant mass)

> {} = pt_min_pdg ! pt cut for other particles (use pdg code). Applied on particle and anti-particle

> {} = pt_max_pdg ! pt cut for other particles (syntax e.g. {6: 100, 25: 50})
> {} = eta_min_pdg ! rap cut for other particles (use pdg code). Applied on particle and anti-particle
> {} = eta_max_pdg ! rap cut for other particles (syntax e.g. {6: 2.5, 23: 5})
> {} = mxx_min_pdg ! min invariant mass of a pair of particles X/X~ (e.g. {6:250})
> {'default': False} = mxx_only_part_antipart ! if True the invariant mass is applied only
> ! to pairs of particle/antiparticle and not to pairs of the same pdg codes.

> I notice that you use the a 3-inputs function idup(i, 1, iproc) by a loop of iproc. I guess "iproc" equals to 1 for the first process like P1_* in SubProcesses, and 2 for P2_*. And actually I don't know the function well.

No the iproc, corresponds to which Y in matrixY.f is inside inside a particular PX_*
Indeed we compute simultaneously multiple process who can differ by their PDG (but not by their mass/...).
When we merge two of such processes we have two cases:

1) they have the same matrix-element than another process
in that case we keep a single matrixY.f for both of them and cuts.f will ALWAYS see the pdg code of the first process. When writing the event, we will randomly choose which of the two (or more) pdg code we will write on disk (potentially pdf weighted choice if they differ in the initial state).
So in that case, it is not possible to do different cut between the two set of pdg code.

2) if they do not have the same matrix-element, then we create a new matrixY.f

So this is why idup has 3 entries:
1) the number of the particle
2) which flavour it is for this matrix-element
3) which matrix-element is currently consider

Cheers,

Olivier

PS: Note that you SHOULD not modify cuts.f, we have a dedicated file that you can modify called
dummy_functions.f which is designed for user/plugin-hook. please use this one.

> On 20 Feb 2019, at 08:53, Rui Zhang <email address hidden> wrote:
>
> New question #678667 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/678667
>
> Hello,
>
>
> I am trying to add my special cuts in madgraph. But I'm confused about how to use particle id in cuts.f(I hope I choose the correct file).
>
> I notice that you use the a 3-inputs function idup(i, 1, iproc) by a loop of iproc. I guess "iproc" equals to 1 for the first process like P1_* in SubProcesses, and 2 for P2_*. And actually I don't know the function well.
>
> I try to extract particles by adding an extra loop of iproc. And if the idup(i, 1, iproc) is what I want, I define the 4-momentum of the chosen particle as p(0:3,i)? Is it correct? Will I extract duplicate or wrong particles?
>
> I add this after your annotation "SPECIAL CUTS" and if it doesn't pass the cuts, I set "passcuts=.false.". But when I apply two complementary cuts, the two cross section don't change either. Of course I have done something wrong.
>
> I hope I asked the proper question. And any help will be appreciated.
>
>
> Thank you!
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thank you, Olivier!

I'm using 2.6.3.2. These cuts are related to all external particles so I think I can't use the basic cut in run_card.dat.

Actually, I tried to change P(0:3,3) to zero in cuts.f but the momentum in lhe file is not zero. So I just hesitate think if I change the wrong file.

Let me try to use dummy_functions.f... But where is it? I use "find" but shell tell me nothing.

It is in SubProcesses directory:

[SubProcesses]$ ls dummy_fct.f
dummy_fct.f

Cheers,

Olivier

> On 20 Feb 2019, at 09:52, Rui Zhang <email address hidden> wrote:
>
> Question #678667 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/678667
>
> Status: Answered => Open
>
> Rui Zhang is still having a problem:
> Thank you, Olivier!
>
>
> I'm using 2.6.3.2. These cuts are related to all external particles so I think I can't use the basic cut in run_card.dat.
>
> Actually, I tried to change P(0:3,3) to zero in cuts.f but the momentum
> in lhe file is not zero. So I just hesitate think if I change the wrong
> file.
>
> Let me try to use dummy_functions.f... But where is it? I use "find" but
> shell tell me nothing.
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thanks Olivier Mattelaer, that solved my question.

Hi, Olivier!

Sorry to bother you again. But I don't think I use it correctly.
As a simple example, I use
----------
      do iproc=1,maxsproc
      do i=nincoming+1,nexternal
          if (idup(i, 1, iproc).eq.6)then
              topp(0:3)=p(0:3,i)
          endif
          if (idup(i, 1, iproc).eq.-6)then
              topm(0:3)=p(0:3,i)
          endif
          if (idup(i, 1, iproc).eq.25)then
              higgs(0:3)=p(0:3,i)
          endif
      enddo
      enddo
----------
to extract the particles, and use
----------
      if((angle.lt.Pi*0.0d0).or.(angle.ge.Pi*0.1d0)) then
          dummy_cuts=.false.
      endif
----------
to apply the cut.

Finally, the cross section and event weight don't change. I also tried to change the P in dummy_fct.f, but it still doesn't change the final momentum. I don't know how to avoid this.
I just consider to set ANS=0 in matrixY.f when it doesn't pass the cut. Is it safe?

Could you help me with extracting the 4-momentum of particles? Thank you!

Hi,

angle is here not defined. Note also that such type of observables are not lorentz invariant and therefore it's depend of the frame in which the momenta are provided. It is possible that your angle is ALWAYS zero in the frame in which we do the computation.

Since I do not know what process you are looking for neither which angle you use, it is obviously just a guess.
Also the double loop that you use might be dangerous since you might pick particle in a weird way with such loop (picking up particles from different iproc). But again I have no clue what you are doing.

Cheers,

Olivier

Thanks. Now I know I used the wrong momentum.

Actually I don't know how to extract the correct momentum. I think I should apply the cut in each matrix element, i.e. with every iproc. But I don't know how this work in madgraph.
I mean that the cut seems work on every iproc so I needn't use iproc in cuts.f and dummy_fct.f. But idup(i,1,iproc) is the only way that I can find to use particle id, and iproc is not a global variable. So when I use particle id, I need to define a iproc and I don't know the value. That's the reason why I use the double loop.

Also, I check the angle by setting the interval to [0.0,0.1), [0.1,0.2), [0.2,0.3), [0.3,0.4), and so on. I think if it's wrong, the result will show me. But the cross section of each bin are the same. If needed, I can try some simple cuts to check it.

Hi, now I try the variable b4 in ttxh production process for example.
----------
define p = p b b~
generate p p > t t~ h
launch
0
set nevents 50000
set asrwgtflavor 5
set cut_decays False
0
exit
----------
And I modify dummy_fct.f as
----------
      double precision function P3D(P)
      implicit none

      double precision P(0:3)
      P3D=sqrt(P(1)*P(1)+P(2)*P(2)+P(3)*P(3))

      return
      end
----------
      double precision topp(0:3)
      double precision topm(0:3)
      double precision higgs(0:3)
      double precision b4
      include 'maxamps.inc'
      integer idup(nexternal,maxproc,maxsproc)
      integer iproc,i
      double precision P3D
----------
      do iproc=1,maxsproc
      do i=nincoming+1,nexternal
          if (idup(i, 1, iproc).eq.6)then
              topp(0:3)=p(0:3,i)
          endif
          if (idup(i, 1, iproc).eq.-6)then
              topm(0:3)=p(0:3,i)
          endif
      enddo
      enddo
      b4=topp(3)*topm(3)/P3D(topp)/P3D(topm)
      if(b4.lt.0.d0) then
          dummy_cuts=.false.
      endif
----------
Also, I change " if(b4.lt.0.d0) then" to " if(b4.ge.0.d0) then". And these three (1 without cut and 2 with complementary cuts) have the same cross section 0.4 pb.

I don't know if the problem is that I use the wrong momentum. But I actually don't know how to find the correct momentum.

Hi,

Ok I'm going to test this on my machine.

This being said, you should rather run like this:
import model sm-no_bmass
generate p p > t t~ h # since the b is massless the "p" contains the b automatically
launch
0
set nevents 50000
#set asrwgtflavor 5 # this is default since your b is massless
#set cut_decays False # this is actually not relevant for your process

Additionally, in this case since you do not have multiparticles in the final state, you will always have the top for i=3, the anti-top for i=4 and the higgs for i=6. For ALL iproc. So you could in principle hardcode such value and/or pick a single iproc for the identification.

Cheers,

Olivier

Hi,

I have it working in my case:
1) I changed the line

      include 'maxamps.inc'
      integer idup(nexternal,maxproc,maxsproc)

to

      include 'maxamps.inc'
      integer idup(nexternal,maxproc,maxsproc)
      integer mothup(2,nexternal)
      integer icolup(2,nexternal,maxflow,maxsproc)
      include 'leshouche.inc'

otherwise your idup is never initialised ...

2) I changed

      if(b4.lt.0.d0) then
          dummy_cuts=.false.
      endif

to
      dummy_cuts=.true.
      if(b4.lt.0.d0) then
          dummy_cuts=.false.
      endif

in order to be sure that you have the default value for dummy_cuts to be correctly set to true (maybe you have jsut forget to your modification). After that I have different value when I changed .lt. to .ge.
(Note I run on a single channel of integration for my test, so I do not assure that you will still get different number when summing over all the channel --if a symmetry implies that such symmetry is kept.

Concerning the momenta, all the momenta provided are in the partonic center of mass-frame.
So if you want such quantity in another frame, you have to boost your momenta in the frame of your choice.

Cheers,

Olivier

Thanks a lot.

This works and the cross sections are approximately correct.

I still have a problem:
I use the double loop to extract the momenta, is it correct?

In this case the double loop will not create any issue (except a quite irrelevant cpu charge which can be avoided)

Cheers,

Olivier

> On 26 Feb 2019, at 14:42, Rui Zhang <email address hidden> wrote:
>
> Question #678667 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/678667
>
> Status: Answered => Open
>
> Rui Zhang is still having a problem:
> Thanks a lot.
>
> This works and the cross sections are approximately correct.
>
> I still have a problem:
> I use the double loop to extract the momenta, is it correct?
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thanks ,Olivier!

I just want to confirm as I actually want to let the particle decay.

Do you mean it's better to directly use the particle number 3,4,5,..., and it's exactly the order which I write when "generate"?

> Do you mean it's better to directly use the particle number 3,4,5,...,
> and it's exactly the order which I write when "generate"?

That's the case, when you do not use the decay syntax.

> I just want to confirm as I actually want to let the particle decay.

If you use the decay syntax, then again the double loop is useless and you can always pick one.
If you do not use the decay syntax, then you need to be very careful with using the id.

Now in the decay syntax, you will also need to loop over "i" which is negative, since you will want to get the id of the propagator (that we assign to negative number in term of ordering).

The ordering of final state is actually given by the ordering of the final state if you replace each decaying particle by their equivalent decay and then use the ordering.
I will not enter nore in the details here since I do not think that you need that.

Cheers,

Olivier

> On 26 Feb 2019, at 14:57, Rui Zhang <email address hidden> wrote:
>
> Question #678667 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/678667
>
> Status: Answered => Open
>
> Rui Zhang is still having a problem:
> Thanks ,Olivier!
>
> I just want to confirm as I actually want to let the particle decay.
>
> Do you mean it's better to directly use the particle number 3,4,5,...,
> and it's exactly the order which I write when "generate"?
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thanks Olivier. That solves all my question.