Unable to create proper grid using MADGraph+amcfast+applgrid

Asked by Rohit Gupta on 2017-08-18

Dear Experts,

I am trying to create grid for eta distribution. I also downloaded some example grid files (https://www.herafitter.org/HERAFitter/HERAFitter/DownloadPage) and in that grids the Q2 axis (actually tau axis where tau= ln(ln(Q2/Lambda0)) ) has very small range for example from 2.4676169 to 2.4676171 but in the grid which i am creating, this range is 2.42 to 2.54. Also in downloaded example grids have points only at two values of Q2 (it look like a sheet) but in my grid it is evenly spread over the entire grid.
 Below is the run_card and the analysis file.
I also modified setscale.f
   elseif(dynamical_scale_choice.eq.0) then
       tmp = 91.118
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc USER-DEFINED SCALE: ENTER YOUR CODE HERE cc
cc to use this code you must set cc
cc dynamical_scale_choice = 0 cc
cc in the run_card (run_card.dat) cc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

======================================run_card.dat====================================================

#***********************************************************************
# MadGraph5_aMC@NLO *
# *
# run_card.dat aMC@NLO *
# *
# This file is used to set the parameters of the run. *
# *
# Some notation/conventions: *
# *
# Lines starting with a hash (#) are info or comments *
# *
# mind the format: value = variable ! comment *
#***********************************************************************
#
#*******************
# Running parameters
#*******************
#
#***********************************************************************
# Tag name for the run (one word) *
#***********************************************************************
  tag_1 = run_tag ! name of the run
#***********************************************************************
# Number of LHE events (and their normalization) and the required *
# (relative) accuracy on the Xsec. *
# These values are ignored for fixed order runs *
#***********************************************************************
  0 = nevents ! Number of unweighted events requested
 -1.0 = req_acc ! Required accuracy (-1=auto determined from nevents)
 -1 = nevt_job! Max number of events per job in event generation.
                 ! (-1= no split).
#***********************************************************************
# Normalize the weights of LHE events such that they sum or average to *
# the total cross section *
#***********************************************************************
 average = event_norm ! average or sum
#***********************************************************************
# Number of points per itegration channel (ignored for aMC@NLO runs) *
#***********************************************************************
 0.001 = req_acc_FO ! Required accuracy (-1=ignored, and use the
                     ! number of points and iter. below)
# These numbers are ignored except if req_acc_FO is equal to -1
 5000 = npoints_FO_grid ! number of points to setup grids
 4 = niters_FO_grid ! number of iter. to setup grids
 10000 = npoints_FO ! number of points to compute Xsec
 6 = niters_FO ! number of iter. to compute Xsec
#***********************************************************************
# Random number seed *
#***********************************************************************
 0 = iseed ! rnd seed (0=assigned automatically=default))
#***********************************************************************
# Collider type and energy *
#***********************************************************************
 1 = lpp1 ! beam 1 type (0 = no PDF)
 1 = lpp2 ! beam 2 type (0 = no PDF)
 4000.0 = ebeam1 ! beam 1 energy in GeV
 4000.0 = ebeam2 ! beam 2 energy in GeV
#***********************************************************************
# PDF choice: this automatically fixes also alpha_s(MZ) and its evol. *
#***********************************************************************
 lhapdf = pdlabel ! PDF set
 11000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
#***********************************************************************
# Include the NLO Monte Carlo subtr. terms for the following parton *
# shower (HERWIG6 | HERWIGPP | PYTHIA6Q | PYTHIA6PT | PYTHIA8) *
# WARNING: PYTHIA6PT works only for processes without FSR!!!! *
#***********************************************************************
  HERWIG6 = parton_shower
  1.0 = shower_scale_factor ! multiply default shower starting
                                  ! scale by this factor
#***********************************************************************
# Renormalization and factorization scales *
# (Default functional form for the non-fixed scales is the sum of *
# the transverse masses of all final state particles and partons. This *
# can be changed in SubProcesses/set_scales.f) *
#***********************************************************************
 true = fixed_ren_scale ! if .true. use fixed ren scale
 true = fixed_fac_scale ! if .true. use fixed fac scale
 91.118 = muR_ref_fixed ! fixed ren reference scale
 91.118 = muF1_ref_fixed ! fixed fact reference scale for pdf1
 91.118 = muF2_ref_fixed ! fixed fact reference scale for pdf2
  0 = dynamical_scale_choice ! Choose one of the preselected dynamical choices
#***********************************************************************
# Renormalization and factorization scales (advanced and NLO options) *
#***********************************************************************
 true = fixed_QES_scale ! if .true. use fixed Ellis-Sexton scale
 91.118 = QES_ref_fixed ! fixed Ellis-Sexton reference scale
 1.0 = muR_over_ref ! ratio of current muR over reference muR
 1.0 = muF1_over_ref ! ratio of current muF1 over reference muF1
 1.0 = muF2_over_ref ! ratio of current muF2 over reference muF2
 1.0 = QES_over_ref ! ratio of current QES over reference QES
#***********************************************************************
# Reweight flags to get scale dependence and PDF uncertainty *
# For scale dependence: factor rw_scale_up/down around central scale *
# For PDF uncertainty: use LHAPDF with supported set *
#***********************************************************************
 True = reweight_scale ! reweight to get scale dependence
 0.5 = rw_Rscale_down ! lower bound for ren scale variations
 2.0 = rw_Rscale_up ! upper bound for ren scale variations
 0.5 = rw_Fscale_down ! lower bound for fact scale variations
 2.0 = rw_Fscale_up ! upper bound for fact scale variations
 False = reweight_PDF ! reweight to get PDF uncertainty
 244601 = PDF_set_min ! First of the error PDF sets
 244700 = PDF_set_max ! Last of the error PDF sets
#***********************************************************************
# ickkw parameter: *
# 0: No merging *
# 3: FxFx Merging - WARNING! Applies merging only at the hard-event *
# level. After showering an MLM-type merging should be applied as *
# well. See http://amcatnlo.cern.ch/FxFx_merging.htm for details. *
# 4: UNLOPS merging (with pythia8 only). No interface from within *
# MG5_aMC available, but available in Pythia8. *
# -1: NNLL+NLO jet-veto computation. See arxiv:1412.8408 [hep-ph]. *
#***********************************************************************
 0 = ickkw
#***********************************************************************
#
#***********************************************************************
# BW cutoff (M+/-bwcutoff*Gamma) *
#***********************************************************************
 15.0 = bwcutoff
#***********************************************************************
# Cuts on the jets *
# Jet clustering is performed by FastJet.
# When matching to a parton shower, these generation cuts should be *
# considerably softer than the analysis cuts. *
# (more specific cuts can be specified in SubProcesses/cuts.f) *
#***********************************************************************
 1.0 = jetalgo ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
 0.7 = jetradius ! The radius parameter for the jet algorithm
 10.0 = ptj ! Min jet transverse momentum
 -1.0 = etaj ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)
#***********************************************************************
# Cuts on the charged leptons (e+, e-, mu+, mu-, tau+ and tau-) *
# (more specific gen cuts can be specified in SubProcesses/cuts.f) *
#***********************************************************************
  0.0 = ptl ! Min lepton transverse momentum
  -1.0 = etal ! Max lepton abs(pseudo-rap) (a value .lt.0 means no cut)
  0.0 = drll ! Min distance between opposite sign lepton pairs
  0.0 = drll_sf ! Min distance between opp. sign same-flavor lepton pairs
  0.0 = mll ! Min inv. mass of all opposite sign lepton pairs
  30.0 = mll_sf ! Min inv. mass of all opp. sign same-flavor lepton pairs
#***********************************************************************
# Photon-isolation cuts, according to hep-ph/9801442 *
# When ptgmin=0, all the other parameters are ignored *
#***********************************************************************
  20.0 = ptgmin ! Min photon transverse momentum
  -1.0 = 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)
#***********************************************************************
# Maximal PDG code for quark to be considered a jet when applying cuts.*
# At least all massless quarks of the model should be included here. *
#***********************************************************************
 4 = maxjetflavor
#***********************************************************************
# For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
#***********************************************************************
 2 = iappl ! aMCfast switch (0=OFF, 1=prepare APPLgrids, 2=fill grids)
#***********************************************************************

====================================================================================================

=========================================Analysis file===================================================

************************************************************************
*
* Analysis routines for vector-boson production.
*
************************************************************************
      subroutine analysis_begin(nwgt,weights_info)
*
      implicit none
*
      include "dbook.inc"
*
      integer nwgt
      character*(*) weights_info(*)
      integer j,kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      character*5 cc(2)
      data cc/" ","Born "/
c*
c* aMCfast common.
c* Needed to redefine the grid parameters
c*
       include "reweight_appl.inc"
       include "appl_common.inc"
c*
c* Grid parameters
c*
       appl_Q2min = 2500d0
       appl_Q2max = 40000d0
       appl_xmin = 2d-7
       appl_xmax = 1d0
       appl_nQ2 = 10
       appl_Q2order = 3
       appl_nx = 50
       appl_xorder = 3
*
* Initialize histograms
*
      call inihist
*
* Check whether the user is trying to fill to many histograms
*
      nwgt_analysis = nwgt
      if (nwgt_analysis*5.gt.nplots/2) then
         write(*,*) "error in analysis_begin: ",
     1 "too many histograms, increase NPLOTS to",
     2 nwgt_analysis * 5 * 2
         call exit(-10)
      endif
*
* Define the observables
*
      do j=1,2
         do kk=1,nwgt_analysis
            l = ( kk - 1 ) * 2 + ( j - 1 ) * 1
            call bookup(l+1,"y_l "//cc(j)//weights_info(kk),
     1 0.2d0,0d0,2.0d0)
         enddo
      enddo
*
      return
      end
*
**********************************************************************
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
*
      implicit none
*
      include "nexternal.inc"
*
      integer i,kk,l
      integer ibody
      integer istatus(nexternal)
      integer iPDG(nexternal)
      double precision p(0:4,nexternal)
      double precision wgts(*)
      integer nwgt_analysis
      common/c_analysis/nwgt_analysis
      double precision www
c** double precision yW,mTW,mT2
      double precision pe(0:3),pn(0:3),pw(0:3)
      double precision ye,pte,etmiss,mtr,cphi,ywen,ptw,var
      double precision ye2530,ye3035,ye3540,ye4045,ye45
      double precision getrapidity
      external getrapidity
*
* All information can be found in analysis_td_template.f
*
* The ibody variable is:
* ibody = 1 : (n+1)-body contribution
* ibody = 2 : n-body contribution (excluding the Born)
* ibody = 3 : Born contribution
*
      if(ibody.ne.1.and.ibody.ne.2.and.ibody.ne.3)then
         write(6,*) "Error: Invalid value for ibody =",ibody
         call exit(-10)
      endif
*
* Check that the number of outcoming particles is correct
*
      if (nexternal.ne.5)then
         write(*,*) "error #1 in analysis_fill: ",
     1 "only for process p p > mu+ vm [QCD]"
         call exit(-10)
      endif
*
* Check that the third particle is either a Z/gamma or a W
*
c** if (abs(ipdg(3)).ne.22.and.abs(ipdg(3)).ne.23.and.
c** 1 abs(ipdg(3)).ne.24) then
c** write(*,*) "error in analysis_fill: ",
c** 1 "The outcoming particle must be a vector boson, ",
c** 2 "ID =",ipdg(3)
c** call exit(-10)
c** endif
*
      if (.not. (abs(ipdg(1)).le.4 .or. ipdg(1).eq.21)) then
         write (*,*) 'error #2 in analysis_fill: '/
     & /'only for process "p p > l vl [QCD]"'
         call exit(-10)
      endif
      if (.not. (abs(ipdg(2)).le.4 .or. ipdg(2).eq.21)) then
         write (*,*) 'error #3 in analysis_fill: '/
     & /'only for process "p p > l vl [QCD]"'
         call exit(-10)
      endif
      if (.not. (abs(ipdg(5)).le.4 .or. ipdg(5).eq.21)) then
         write (*,*) 'error #4 in analysis_fill: '/
     & /'only for process "p p > l vl [QCD]"'
         call exit(-10)
      endif
      if( (abs(abs(ipdg(3))-abs(ipdg(4))).ne.1) .or.
     & (sign(1d0,dble(ipdg(3))).eq.sign(1d0,dble(ipdg(4)))) .or.
     & (abs(ipdg(3)).le.10.or.abs(ipdg(3)).ge.16) .or.
     & (abs(ipdg(4)).le.10.or.abs(ipdg(4)).ge.16) )then
         write(*,*)'analysis not suited for this process',ipdg(3),ipdg(4)
      endif
* Compute observables
*
*
      do i=0,3
         if(abs(ipdg(3)).eq.11.or.abs(ipdg(3)).eq.13.or.abs(ipdg(3)).eq.
     1 15)then
            pe(i)=p(i,3)
            pn(i)=p(i,4)
         else
            pe(i)=p(i,4)
            pn(i)=p(i,3)
         endif
         pw(i)=pe(i)+pn(i)
      enddo
      ye = getrapidity(pe(0), pe(3))
      pte = dsqrt(pe(1)**2 + pe(2)**2)
      etmiss = dsqrt(pn(1)**2 + pn(2)**2)
      ptw = dsqrt(pw(1)**2 + pw(2)**2)
      mtr = dsqrt(2d0*pte*etmiss-2d0*pe(1)*pn(1)-2d0*pe(2)*pn(2))
      if(pte.gt.25d0.and.mtr.gt.50d0.and.etmiss.gt.15d0)then
        ye2530=getrapidity(pe(0), pe(3))
      endif

*
* Fill the histograms
*
      do i=1,2
         do kk=1,nwgt_analysis
            www = wgts(kk)
            l = ( kk - 1 ) * 2 + ( i - 1 ) * 1
* Only fill Born histograms for ibody=3
            if (ibody.ne.3.and.i.eq.2) cycle
*
            call mfill(l+1,ye2530,www)
         enddo
      enddo
*
      return
      end
*
************************************************************************
      subroutine analysis_end(xnorm)
*
      implicit none
*
      include "dbook.inc"
*
      character*14 ytit
      double precision xnorm
      integer i
      integer kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
*
      call open_topdrawer_file
      call mclear
*
      do i=1,NPLOTS
         call mopera(i,"+",i,i,xnorm,0.d0)
         call mfinal(i)
      enddo
      ytit = "sigma per bin "
      do i=1,2
         do kk=1,nwgt_analysis
            l = ( kk - 1 ) * 2 + ( i - 1 ) * 1
            call multitop(l+1,3,2,"y_l ", ytit,"LIN")
         enddo
      enddo
*
      call close_topdrawer_file
*
      return
      end
*
************************************************************************
*
* Auxiliary functions for the kinematics
*
************************************************************************
      function getrapidity(en,pl)
*
      implicit none
*
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
      parameter (tiny=1d-8)
*
      xplus = en + pl
      xminus = en - pl
*
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
         if((xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
            y = 0.5d0 * log( xplus / xminus )
         else
            y = sign(1d0,pl) * tiny
         endif
      else
         y = sign(1d0,pl) * tiny
      endif
*
      getrapidity = y
*
      return
      end

=====================================================================================================

Thanks in advance.

Question information

Language:
English Edit question
Status:
Answered
For:
MadGraph5_aMC@NLO Edit question
Assignee:
Rikkert Frederix Edit question
Last query:
2017-08-26
Last reply:
2017-08-29
Rikkert Frederix (frederix) said : #1

Dear Rohit,

What is your question?

best,
Rikkert

Bajarang (bajjubaba) said : #2

I want to create grids similar to the downloaded grids.
Actually I want to create grid with fixed Q2 values instead of dynamic scale variation. In downloaded example grid also all points are at two Q2 values but the grid which I created has points spread across full Q2 range.

Rohit Gupta (rohithep) said : #3

Dear Rikkert

Please ignore previous message.....

I want to create grids similar to the downloaded grids.
Actually I want to create grid with fixed Q2 values instead of dynamic scale variation. In downloaded example grid also all points are at two Q2 values but the grid which I created has points spread across full Q2 range.

Rikkert Frederix (frederix) said : #4

Dear Rohit,

Which version of MadGraph5_aMC@NLO are you using?

Do you get the right grids with setting 'reweight_scale' to 'False' in the run_card?

Best,
Rikkert

Rohit Gupta (rohithep) said : #5

I am using version 2.3.3

I tried with 'reweight_scale' to False but it is giving this error

Command "launch auto " interrupted with error:
InvalidRunCard : APPLgrid generation only possible with including the reweighting to get scale dependence

Rikkert Frederix (frederix) said : #6

Dear Rohit,

I'm still not sure what your problem exactly is. Do you have problems reproducing the results (plots) with the grids? Why should the grids be exactly identical? I mean, in the end, the grids are something non-physical. The plots and cross sections that can be computed with the grids are physical. Are you saying that these latter are not correct?

Best,
Rikkert

Rohit Gupta (rohithep) said : #7

Hi Rikkert,

Sorry for the late reply and also for the confusion.

In one line the problem is : How to create the grids with fixed Q^2, just like in the example file I found for xfitter : http://xfitter.hepforge.org/data.html (1206.2598)

Because when I create it, the cross section values are all spread over the grid from tau~2.3 to tau~2.7 whereas in the example file I see that, its mainly located at tau~2.4 (which corresponds to Q^2~70GeV). How can I obtain that?

I tried many variations, changing the renormalization and factorization scales, tried providing the dynamical_scale_choice (~70GeV), tried changing the rw_Rscale_down etc, but everytime, I get a grid which has the cross-section values spread all over in the grid.

Forgive me but I am not exactly sure, if this spread over of the cross-section is responsible or is not responsible for me not being able to replicate the results in xfitter, but intuitively I think if I get the grid right, I will get matching results in xfitter.

Thanks for your patience.
Regards,
Rohit

Rikkert Frederix (frederix) said : #8

Dear Rohit,

Still, I don't really understand the problem. In the end, are you able to reproduce plots you made with

1. a single run of MG5_aMC,

2. using the grids generated with MG5_aMC+aMCfast+APPLGrid?

MG5_aMC automatically includes a factor 2 scale variation, and this information is also available in the grids. However, of course, one can also only look at the central value.

I don't know what's inside the xfitter grids. But this is not really important, since the grids themselves are not physical observables. Hence, the comparison should be made at the level of plots created from those grids in a consistent way.

We have made many self-consistency tests within the MG5_aMC+aMCfast+APPLGrid code, see arXiv:1406.7693. All checks have been done at the level of physical observables, which is the only thing that needs to be correct.

Best,
Rikkert

Can you help with this problem?

Provide an answer of your own, or ask Rohit Gupta for more information if necessary.

To post a message you must log in.