Charm production

Hi Shahryar,
in order to do things consistently, you need to have a massive charm quark in the model.
to do this, do
> import model loop_sm-c_mass

then, always for consistency, define the multiparticle p and j (light jet) without the c

> define p = g u u~ d d~ s s~
> define j = p

and now
> generate p p > c a [QCD]

should normally work.

Let me know if you need further help.

Cheers,

Marco

Dear Marco

 Thanks for your consideration. I did your comment step by step and unfortunately get an error this time: Born diagrams could not generated for this process. Furthermore, It should be noted that I want to compare the MadGraph5_aMC@NLO results for this process to experimental measurements for example from D0 collaboration (arXiv:1210.5033 [hep-ex]). According to this article, photons in association with a charm (c) quark are produced primarily through the Compton-like scattering process gc→γ c, which dominates at lower photon transverse momenta. Therefor I need having charm in initial state and so it should be mass less. I neglected the warning appeared when I generated this process as
generate p p > a c [QCD]
without any change in the model and p and j definition and I get reasonable results comparing to the SHERPA results and D0 measurements specially at low and medium photon transverse momenta. But there is still a problem. According to D0 article (arXiv:1210.5033 [hep-ex]), photons in association with a charm (c) quark are produced also through quark-antiuqrks annihilation q q~->gamma g->gamma c c~ which dominates at higher photon transverse momenta. When I generated this process as above mentioned way, this important sub-process not included in the cross section calculation and so I get very small result at higher p_T.
 Do you have any comment for me? How I can calculate cross section exactly corresponding D0 measurement for instance?

Dear Shahryar,

There are two ways to make predictions for photon + charm jet:

1. three-flavour scheme
In this scheme, you treat the c-quark as a massive quark, which can be distinguished from the other light quarks in the matrix elements. In that case you should not make the incoming protons are outgoing jets contain c quarks:

import model loop_sm-c_mass
define p = g u u~ d d~ s s~
define j = p
generate p p > a c c~ [QCD]

Note that you have to specify both the charm and the anti-charm in the final state, without any generation cuts on the charm quarks (although you need generation cuts for the photon). Then, at the level of the analysis, one (or even two) of the charm quarks can be very soft or collinear to the beam and therefore goes undetected.

2. four (or five) flavour scheme
The charm is treated as massless and is indistinguishable from the other light quarks at the matrix element level. Only after showering you can try to discriminate the c-jets from the light jets (Note that you cannot make predictions for fNLO for this process, only at NLO+PS level)

import model loop_sm
define p = g u u~ d d~ s s~ c c~
define j = p
generate p p > a j [QCD]

In this case, you need a generation cut on the jet (as well as on the photon), otherwise you'll run into divergencies. Your generation cuts should be loose enough so that there is no dependence on them after showering.

You'll not get exactly the same results for the two scheme. Only if you would include all orders in perturbation theory the results will be identical. In most cases at the LHC solution 2 is the preferred one because large logs are resummed in the PDFs. On the other hand, in particular at low pT for the c-quark, solution 1 might be better as it takes the mass of the quark into account and does not need any generation cut.

Best regards,
Rikkert

Dear Rikkert

Thanks for your consideration. At present, I choose the solution 1 (as you mentioned it is better at low pT) and did following steps:

import model loop_sm-c_mass
define p = g u u~ d d~ s s~
define j = p
generate p p > a c c~ [QCD]

Then, According to the D0 publication (arXiv:1210.5033 [hep-ex]), I put in run_card.dat:

    1 = lpp1 ! beam 1 type (0 = no PDF)
    -1 = lpp2 ! beam 2 type (0 = no PDF)
 980 = ebeam1 ! beam 1 energy in GeV
 980 = ebeam2 ! beam 2 energy in GeV

 T = fixed_ren_scale ! if .true. use fixed ren scale
 T = fixed_fac_scale ! if .true. use fixed fac scale
 34.2 = muR_ref_fixed ! fixed ren reference scale
 34.2 = muF1_ref_fixed ! fixed fact reference scale for pdf1
 34.2 = muF2_ref_fixed ! fixed fact reference scale for pdf2
#***********************************************************************
# Renormalization and factorization scales (advanced and NLO options) *
#***********************************************************************
 F = fixed_QES_scale ! if .true. use fixed Ellis-Sexton scale
 34.2 = QES_ref_fixed ! fixed Ellis-Sexton reference scale

   1 = jetalgo ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
 0.7 = jetradius ! The radius parameter for the jet algorithm
  15 = ptj ! Min jet transverse momentum
  1.5 = etaj ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)

  30 = ptgmin ! Min photon transverse momentum
  1 = etagamma ! Max photon abs(pseudo-rap)
 0.4 = R0gamma ! Radius of isolation code
 1.0 = xn ! n parameter of eq.(3.4) in hep-ph/9801442
 1.0 = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
 .true. = isoEM ! isolate photons from EM energy (photons and leptons)

and I also manually changed the file cuts.f as follows

            if(ptg.lt.ptgmin .or. ptg.gt.40D0)then
               passcuts_user=.false.
               return
            endif

because I want to calculate the cross section for photon transverse momenta between 30-40 GeV.
Finally, I choose

Fixed_order=ON

and get 1.912 +- 0.029 as a result for the cross section. But this value is very small comparing SHERPA (10) and NLO QCD result (10.5). I'll be grateful if you guide me again.

Best
Shahryar

Dear Shahryar,

If you kept the maxjetflavour to 4, it means that in your jet definition you also include the c quark. Hence, you are requiring two jets with

  1 = jetalgo ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
 0.7 = jetradius ! The radius parameter for the jet algorithm
  15 = ptj ! Min jet transverse momentum
  1.5 = etaj ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)

You'll have to change the cuts.f to require only one c-charged jet.

Best regards,
Rikkert

Dear Rikkert

Thanks for your comment. But as I mentioned in previous question, I get a very small value for the cross section comparing to SHERPA and QCD NLO results an also experimental measurement from D0. By doing your new coment we will get even a smaller value for the cross section than last time. I put in run_card.dat

maxjetflavour=3

and changed cuts.f as follows

c find the jets
      do i=1,nexternal
         if (istatus(i).eq.1 .and.
     & (abs(ipdg(i)).le.maxjetflavor .or. ipdg(i).eq.21)
     & .or. ipdg(i).eq.4 ) then
            is_a_j(i)=.true.
         else
            is_a_j(i)=.false.
         endif
      enddo

In this way, we get 0.7371 +- 0.093 for the cross section.

Best regards,
Shahryar

Dear Shahryar,

Your suggestion is not correct. In this way you require at least one jet that can be made up by any quarks/gluons or charm quarks. (but not anti-charm quarks). This does not mean that the jet you get is charm flavoured. It could very well be made from only a single gluon.

You have to include *both* charms in the jet-definition, but then check into which jet the charm went and cut on that jet. This can be done using the 'jet' array as filled by the call to amcatnlo_fastjetppgenkt_etamax_timed().
Here you might (or not) exclude jets that contain both the charm and anti-charm. But this is up to you, because it is not clear if

Best,
Rikkert