high-pT top quark cross section

Asked by Sergei Chekanov

 Hello,

 We are sporting one problem in generation of Madgraph files. We run Madgraph for ttbar, for 100 TeV CM energy, and set the cut on pT of top quarks at 2.5 TeV. For NLO program, we use SubProcesses/cuts.f to set pT cut.

Here are LO Monte Carlo events for ttbar with pT>2500 GeV
 http://atlaswww.hep.anl.gov/hepsim/info.php?item=16
 (sigma = 0.2 pb)

 And this is NLO prediction for ttbar with pT>2500 GeV:
 http://atlaswww.hep.anl.gov/hepsim/info.php?item=57
 (sigma = 1 pb)

As you can see "k-factor" is 5, which looks too large. Do you know, maybe we do not compare the same things since implementations of LO and NLO cuts on top are different?

best, Sergei

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Solved by:
Sergei Chekanov
Solved:
Last query:
Last reply:
Revision history for this message
marco zaro (marco-zaro) said :
#1

Hi,
running amcatnlo (with hand-coded cuts on the tops) i get
LO: 1.208e-01 +- 1.7e-03 pb (so i confirm what you guys get)
and the NLO takes a long long time.

I would like anyway to stress that such hard cuts can lead to unstable results at the NLO (a thing that can be seen by the fact that the NLO integration takes a very long time, and of course by looking at the logs), corresponding to 'almost' hitting the collinear singularity when a hard parton is emitted in the real radiation and the top-antitop pair recoils against it.
Given the strong hierarchy of scales, in particular pt_cut >> mtop, the tops almost behave as massless quarks, and the cross-section will show a log(mt/pt_cut).
I am not sure what the best way of doing such a computation can be, i have two ideas, one which can be implemented right now, the other which needs more time:

the easy one is to require a minimum separation between the tops, in order not to get close to the singularity, and can be implemented easily in cuts.f (or a variant which is a bit more involved, but for sure more clean, would be to cluster the top, antitop and the optional extra partons into jets with some suitable R value, and to ask for 2 jets with pt > 2500 gev and with a top/antitop as constituent, as one would do for b-jets)

the other solution is to use a 6flavour scheme, which is not currently available in mg5_amc.

i have added Rik to the discussion, if he knows more about this...

Cheers,

Marco

Revision history for this message
Rikkert Frederix (frederix) said :
#2

Hi Sergei, Marco,

The cuts you apply are at the level of the Les Houches event generation, or at the level of the analysis. Of course, the parton shower will move the tops around (a bit) and therefore a 2500 GeV top at the LHEF level, is not the same as at the hadron level.

The only thing I can add to improve the generation is that you should definitely use a 5 flavour scheme (massless b quark) and not the default 4 flavour scheme, if you aren't using this yet.

Best regards,
Rik

Revision history for this message
Sergei Chekanov (chekanov) said :
#3

 Hi, Rik

 I apply the pT cut as:

 c pT cut
       passcuts=.false.
       do i=nincoming+1,nexternal
               if (pt(p(0,i)).gt.2500) then
                  passcuts=.true.
                  return
               endif
       enddo
       if (passcuts .eqv. .false.) return

in

 SubProcesses/cuts.f

Is this "generator" level cut?

best, Sergei

Revision history for this message
Rikkert Frederix (frederix) said :
#4

Dear Sergei,

That is indeed at the level of the LHEF event generation.

However, what you are doing is not correct.
First, I strongly advice to update to the latest version of the code. In the latest version, applying cuts in cuts.f has become a lot easier for NLO processes. In fact, there is already an example for applying a cuts on the pT of the top quark.
Your mistakes, are that you are applying the cuts to all final state particles, including the possible extra parton present in the real-emission corrections. Moreover, in your if statement

                if (pt(p(0,i)).gt.2500) then

you should put a double precision number, not an integer. Ie.,

                if (pt(p(0,i)).gt.2500d0) then

Best regards,
Rikkert

Revision history for this message
Sergei Chekanov (chekanov) said :
#5

Hello, Rikkert

 Thanks for the suggesting. I've followed it: I've download MG5_aMC_v2.1.2,
 changed SubProcesses/cuts.f to

 c pT cut on top
       passcuts_user=.false.
       do i=nincoming+1,nexternal
               if (pt(p(0,i)).gt.2500d0) then
                  passcuts_user=.true.
                  return
               endif
       enddo
       if (passcuts_user .eqv. .false.) return

Then I've changed: Template/NLO/Cards/run_card.dat to

 50000 = ebeam1 ! beam 1 energy in GeV
 50000 = ebeam2 ! beam 2 energy in GeV

and

  lhapdf = pdlabel ! PDF set
  21100 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number (MSTW2008nlo68cl.LHgrid)

I ran for 30 min (5000 events) and I've got this:

Cross=9.84906270e-01 pb
Error=1.4878e-02 pb

which is still ~4 larger than at LO (~0.2 pb). I hope this is not ralated to the PDF choice..

best, Sergei

Revision history for this message
marco zaro (marco-zaro) said :
#6

Ciao Sergei,

I think that the large k-factor you get is precisely due to the
quasi-collinear singularity that you encounter with these cuts, i.e. the
top pair recoiling against a hard gluon (see my previous mail).
Can you try imposing an extra distance cut on the two tops?

Cheers,

Marco

Marco Zaro

On Tue, Sep 2, 2014 at 2:08 AM, Sergei Chekanov <
<email address hidden>> wrote:

> Question #253711 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/253711
>
> Sergei Chekanov posted a new comment:
> Hello, Rikkert
>
> Thanks for the suggesting. I've followed it: I've download
> MG5_aMC_v2.1.2,
> changed SubProcesses/cuts.f to
>
> c pT cut on top
> passcuts_user=.false.
> do i=nincoming+1,nexternal
> if (pt(p(0,i)).gt.2500d0) then
> passcuts_user=.true.
> return
> endif
> enddo
> if (passcuts_user .eqv. .false.) return
>
> Then I've changed: Template/NLO/Cards/run_card.dat to
>
> 50000 = ebeam1 ! beam 1 energy in GeV
> 50000 = ebeam2 ! beam 2 energy in GeV
>
> and
>
> lhapdf = pdlabel ! PDF set
> 21100 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
> (MSTW2008nlo68cl.LHgrid)
>
> I ran for 30 min (5000 events) and I've got this:
>
> Cross=9.84906270e-01 pb
> Error=1.4878e-02 pb
>
> which is still ~4 larger than at LO (~0.2 pb). I hope this is not
> ralated to the PDF choice..
>
> best, Sergei
>
> --
> You received this question notification because you are a direct
> subscriber of the question.
>

Revision history for this message
Rikkert Frederix (frederix) said :
#7

Dear Sergei,

One more thing. You should NOT apply the cuts on all the final state particles, but only on the top and anti-top. In your current setup in cuts.f, you also apply the cuts on the possible extra parton that is present in the real emission corrections. You can simply uncomment this part of cuts.f:

C***************************************************************
C***************************************************************
C PUT HERE YOUR USER-DEFINED CUTS
C***************************************************************
C***************************************************************
C
c$$$C EXAMPLE: cut on top quark pT
c$$$ do i=1,nexternal ! loop over all external particles
c$$$ if (istatus(i).eq.1 ! final state particle
c$$$ & .and. abs(ipdg(i)).eq.6) then ! top quark
c$$$C apply the pT cut (pT should be large than 200 GeV for the event to
c$$$C pass cuts)
c$$$ if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then
c$$$C momenta do not pass cuts. Set passcuts_user to false and return
c$$$ passcuts_user=.false.
c$$$ return
c$$$ endif
c$$$ endif
c$$$ enddo
c

and change the 200d0**2 to 2500d0**2.

Did you try running the LO with the same PDF set?

Cheers,
Rik

Revision history for this message
Sergei Chekanov (chekanov) said :
#8

 Thanks a lot!

 I've just uncommented these lines:

 c$$$C EXAMPLE: cut on top quark pT
      do i=1,nexternal ! loop over all external particles
         if (istatus(i).eq.1 ! final state particle
     & .and. abs(ipdg(i)).eq.6) then ! top quark
c$$$C apply the pT cut (pT should be large than 200 GeV for the event to
c$$$C pass cuts)
            if ( p(1,i)**2+p(2,i)**2 .lt. 2500d0**2 ) then
c$$$C momenta do not pass cuts. Set passcuts_user to false and return
               passcuts_user=.false.
               return
            endif
         endif
      enddo

and reran the whole thing. Now my cross section looks rather low:

Cross=7.07072233e-02
Error=4.5571e-03

This is smaller than the cross section at LO that was done using CTEQ6L (0.2 pb)

 best, Sergei

Revision history for this message
Sergei Chekanov (chekanov) said :
#9

 Hello,

 I think the problem was solved. I made this code in cuts.f (this only works for MG5_aMC_v2_1_2!):

 c pT cut on top
       passcuts_user=.false.
       do i=nincoming+1,nexternal
                   if (istatus(i).eq.1 ! final state particle
     & .and. abs(ipdg(i)).eq.6) then ! top quark
                        if ( (p(1,i)**2+p(2,i)**2) .gt. 2500d0**2 ) then
                        passcuts_user=.true.
                        return
                   endif
                   endif
       enddo
       if (passcuts_user .eqv. .false.) return

And I've got this result (using MSTW2008NLO):

     Summary:
      Process p p > t t~ [QCD]
      Run at p-p collider (50000 + 50000 GeV)
      Total cross-section: 5.665e-01 +- 9.4e-03 pb
      Number of events generated: 5000
      Parton shower to be used: HERWIG6
      Fraction of negative weights: 0.13
      Total running time : 44m 47s

The k-factor is large (at LO, I had 0.2), but this is what we should see.
 Thanks to Michelangelo Mangano, who checked this cross section. According to his program, he has
 0.66 pb, which is close to the above result (the difference can be due to PDF).
 He also observes quite large k-factor (~3)

 I'm closing this problem report.

best, Sergei

Revision history for this message
Rikkert Frederix (frederix) said :
#10

Dear Sergei,

Ok. Your code should work as it requires that at least one top quark should have a pT of 2500 GeV, while what you tried before was requiring that both top quarks should have a pT of 2500 GeV.

Cheers,
Rik