NNPDF · cschwan · Mar 23, 2021
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/analysis.f b/nnpdf31_proc/NNPDF40_JETS_BBAR/analysis.f
@@ -0,0 +1,144 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+
+      call set_error_estimation(1)
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(1,'int',1, 0d0, 1d0)
+
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+c     subroutine parameters
+      double precision p(0:4,nexternal)
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision wgts(*)
+      integer ibody
+
+c     variables for amcatnlo_fastjetppgenkt
+      double precision pQCD(0:3,nexternal),pjet(0:3,nexternal),
+     $                 etajet(nexternal)
+      integer nQCD,jet(nexternal),njet
+
+c     observables
+      integer xbin
+      double precision xystar,xptavg,ptjet(nexternal),yjet(nexternal),
+     $                 xyboost
+
+c     functions
+      double precision getptv4,getinvm,getrapidityv4
+      external getptv4,getinvm,getrapidityv4
+
+c     miscellaneous
+      integer i,j
+
+      nQCD=0
+      do j=nincoming+1,nexternal
+         if (abs(ipdg(j)).le.5 .or. ipdg(j).eq.21 .or.
+     $                              ipdg(j).eq.22) then
+            nQCD=nQCD+1
+            do i=0,3
+               pQCD(i,nQCD)=p(i,j)
+            enddo
+         endif
+      enddo
+
+      do i=1,nexternal
+         do j=0,3
+            pjet(j,i)=0d0
+         enddo
+         jet(i)=0
+      enddo
+
+c     recombine momenta
+      call amcatnlo_fastjetppgenkt_etamax(pQCD,nQCD,jetradius,
+     $     ptj,etaj,jetalgo,pjet,njet,jet)
+
+      if (njet.lt.2) then
+        write (*,*) "error: event contains less than two jets"
+        return
+      endif
+
+      call HwU_fill(1,0.5d0,wgts)
+
+ 999  return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      function getptv4(p)
+      implicit none
+      real*8 getptv4,p(0:3)
+c
+      getptv4=sqrt(p(1)**2+p(2)**2)
+      return
+      end
+
+      function getinvm(en,ptx,pty,pl)
+      implicit none
+      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
+      parameter (tiny=1.d-5)
+c
+      tmp=en**2-ptx**2-pty**2-pl**2
+      if(tmp.gt.0.d0)then
+        tmp=sqrt(tmp)
+      elseif(tmp.gt.-tiny)then
+        tmp=0.d0
+      else
+        write(*,*)'Attempt to compute a negative mass'
+        stop
+      endif
+      getinvm=tmp
+      return
+      end
+
+      function getrapidity(en,pl)
+      implicit none
+      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
+      parameter (tiny=1.d-8)
+c
+      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(1.d0,pl)*1.d8
+        endif
+      else
+        y=sign(1.d0,pl)*1.d8
+      endif
+      getrapidity=y
+      return
+      end
+
+      function getrapidityv4(p)
+      implicit none
+      real*8 getrapidityv4,p(0:3)
+      real*8 getrapidity
+c
+      getrapidityv4=getrapidity(p(0),p(3))
+      return
+      end
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/change_etaj_to_rapj.patch b/nnpdf31_proc/NNPDF40_JETS_BBAR/change_etaj_to_rapj.patch
@@ -0,0 +1,48 @@
+--- NLO/SubProcesses/fastjetfortran_madfks_core.cc.orig	2020-11-27 18:11:04.027146365 +0100
++++ NLO/SubProcesses/fastjetfortran_madfks_core.cc	2020-11-27 18:13:47.285826142 +0100
+@@ -76,7 +76,7 @@
+   /// and the extraction of the jets
+   void amcatnlo_transfer_cluster_transfer(const double * p, const int & npart, 
+                                  const JetDefinition & jet_def,
+-                                 const double & ptmin, const double & etamax,
++                                 const double & ptmin, const double & rapmax,
+ 				 double * f77jets, int & njets, int * whichjet) {
+
+     // transfer p[4*ipart+0..3] -> input_particles[i]
+@@ -90,9 +90,9 @@
+     jets = sorted_by_pt(cs->inclusive_jets(ptmin));
+
+     //apply the eta selector if etamax >0
+-    Selector select_eta = SelectorAbsEtaMax(etamax);
+-    if (etamax > 0.) {
+-        jets = select_eta(jets);
++    Selector select_rap = SelectorAbsRapMax(rapmax);
++    if (rapmax > 0.) {
++        jets = select_rap(jets);
+     }
+
+     // transfer jets -> f77jets[4*ijet+0..3]
+--- NLO/SubProcesses/fastjetfortran_madfks_full.cc.orig	2020-11-27 18:11:12.687075955 +0100
++++ NLO/SubProcesses/fastjetfortran_madfks_full.cc	2020-11-27 18:13:47.289826110 +0100
+@@ -79,7 +79,7 @@
+   /// and the extraction of the jets
+   void amcatnlo_transfer_cluster_transfer(const double * p, const int & npart, 
+                                  const JetDefinition & jet_def,
+-                                 const double & ptmin, const double & etamax,
++                                 const double & ptmin, const double & rapmax,
+ 				 double * f77jets, int & njets, int * whichjet,
+ 				 const double & ghost_maxrap = 0.0,  
+ 				 const int & nrepeat = 0, const double & ghost_area = 0.0) {
+@@ -101,9 +101,9 @@
+     jets = sorted_by_pt(cs->inclusive_jets(ptmin));
+
+     //apply the eta selector if etamax >0
+-    Selector select_eta = SelectorAbsEtaMax(etamax);
+-    if (etamax > 0.) {
+-        jets = select_eta(jets);
++    Selector select_rap = SelectorAbsRapMax(rapmax);
++    if (rapmax > 0.) {
++        jets = select_rap(jets);
+     }
+
+     // transfer jets -> f77jets[4*ijet+0..3]
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/change_scale_to_mjj.patch b/nnpdf31_proc/NNPDF40_JETS_BBAR/change_scale_to_mjj.patch
@@ -0,0 +1,102 @@
+--- CMS_2JET_7TEV/SubProcesses/setscales.f.orig	2020-11-30 09:45:53.669448000 +0100
++++ CMS_2JET_7TEV/SubProcesses/setscales.f	2020-12-09 15:04:10.610477654 +0100
+@@ -527,6 +527,29 @@
+       integer i,j
+       character*80 temp_scale_id
+       common/ctemp_scale_id/temp_scale_id
++
++      integer iPDG_reco(nexternal)
++      double precision ppl(0:3), pplb(0:3), ppv(0:3), xmll
++      double precision p_reco(0:4,nexternal), p_in(0:4,nexternal)
++c     les houches accord stuff to identify particles
++c
++      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
++     &        icolup(2,nexternal,maxflow),niprocs
++      common /c_leshouche_inc/idup,mothup,icolup,niprocs
++c Masses of external particles
++      double precision pmass(nexternal)
++      common/to_mass/pmass
++      double precision pQCD(0:3,nexternal),
++     $       pjet(0:3,nexternal),etajet(nexternal)
++      integer nQCD,jet(nexternal),njet
++      double precision ptjet(nexternal),yjet(nexternal)
++c For boosts
++      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
++      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
++     #                        sqrtshat,shat
++      double precision chybst,shybst,chybstmo
++      double precision xd(1:3)
++      data (xd(i),i=1,3)/0,0,1/
+ c
+       tmp=0
+       if(ickkw.eq.-1)then
+@@ -562,16 +587,59 @@
+           tmp=muR_ref_fixed
+           temp_scale_id='fixed scale'
+       elseif(dynamical_scale_choice.eq.10) then
+-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+-cc      USER-DEFINED SCALE: ENTER YOUR CODE HERE                                 cc
+-cc      to use this code you must set                                            cc
+-cc                 dynamical_scale_choice = 10                                   cc
+-cc      in the run_card (run_card.dat)                                           cc
+-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+-         write(*,*) "User-defined scale not set"
+-         stop 1
+-         temp_scale_id='User-defined dynamical scale' ! use a meaningful string
+-         tmp = 0d0
++          temp_scale_id='Mjj'
++
++          do i=1,nexternal
++            p_in(0:3,i) = pp(0:3,i)
++            p_in(4,i) = pmass(i)
++          enddo
++          call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
++     $                           p_in, idup(1,1), p_reco, iPDG_reco)
++
++          nQCD=0
++          chybst=cosh(ybst_til_tolab)
++          shybst=sinh(ybst_til_tolab)
++          chybstmo=chybst-1.d0
++          do j=nincoming+1,nexternal
++             if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
++     $                                  ipdg_reco(j).eq.22) then
++                nQCD=nQCD+1
++                call boostwdir2(chybst,shybst,chybstmo,xd,
++     $                          p_reco(0,j),pQCD(0,nQCD))
++             endif
++          enddo
++
++          do i=1,nexternal
++             do j=0,3
++                pjet(j,i)=0d0
++             enddo
++             jet(i)=0
++          enddo
++
++c         recombine momenta
++          call amcatnlo_fastjetppgenkt_etamax(pQCD,nQCD,jetradius,ptj,
++     $         etaj,jetalgo,pjet,njet,jet)
++
++          do i=1,njet
++            ptjet(i)=sqrt(pjet(1,i)**2+pjet(2,i)**2)
++            if(i.gt.1)then
++              if (ptjet(i).gt.ptjet(i-1)) then
++                write (*,*) "Error 3: jets should be ordered in pt"
++                write (*,*) ptjet(i), ptjet(i-1), i
++                stop
++              endif
++            endif
++          enddo
++
++         if (njet .lt. 2) then
++c          something is wrong with the cuts - they should check for two jets
++           tmp=4d0
++           write (*,*) "Error 4: scale couldn't find two jets"
++         else
++           tmp=sqrt((pjet(0,1)+pjet(0,2))**2-(pjet(1,1)+pjet(1,2))**2
++     $                                      -(pjet(2,1)+pjet(2,2))**2
++     $                                      -(pjet(3,1)+pjet(3,2))**2)
++         endif
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ cc      USER-DEFINED SCALE: END OF USER CODE                                     cc
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/enable_nf_5_to_4_ct.patch b/nnpdf31_proc/NNPDF40_JETS_BBAR/enable_nf_5_to_4_ct.patch
@@ -0,0 +1,29 @@
+--- NLO/SubProcesses/fks_singular.f.orig	2021-01-14 09:08:28.216954859 +0100
++++ NLO/SubProcesses/fks_singular.f	2021-01-14 09:06:13.450217667 +0100
+@@ -112,7 +112,7 @@
+       double precision wgtborn, alphas
+       ! switch on/off here
+       logical include_6to5_cnt 
+-      data include_6to5_cnt /.false./ 
++      data include_6to5_cnt /.true./ 
+
+ CMZMZ REMEMBER!!!!
+ c     wgt1 : weight of the contribution not multiplying a scale log
+@@ -126,7 +126,7 @@
+
+       ! skip if we don't want this piece or if the scale is
+       ! below mt
+-      if (.not.include_6to5_cnt.or.scale.lt.mdl_mt) return
++      if (.not.include_6to5_cnt.or.scale.lt.mdl_mb) return
+
+ C the contribution is the following (if mu > mt):
+ C      Add a term -alphas n TF/3pi log (muR^2/mt^2) sigma(0) 
+@@ -162,7 +162,7 @@
+      &    - alphas / 3d0 / pi * TF * dble(alphasbpow) * amp_split(iamp) 
+
+           amp_split_6to5f(orders_to_amp_split_pos(orders)) = 
+-     &    dlog(qes2/mdl_mt**2) * 
++     &    dlog(qes2/mdl_mb**2) * 
+      &     (alphas / 3d0 / pi * TF * dble(niglu)   
+      &    - alphas / 3d0 / pi * TF * dble(alphasbpow)) * amp_split(iamp)
+         endif
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/launch.txt b/nnpdf31_proc/NNPDF40_JETS_BBAR/launch.txt
@@ -0,0 +1,28 @@
+launch @OUTPUT@
+fixed_order = ON
+set maxjetflavor 4
+set gf @GF@
+set mh @MH@
+set mt @MT@
+set mw @MW@
+set mz @MZ@
+set wh @WH@
+set wt @WT@
+set ww @WW@
+set wz @WZ@
+set ebeam1 7000
+set ebeam2 7000
+set pdlabel lhapdf
+set lhaid 324900
+set dynamical_scale_choice 10
+set reweight_scale True
+set jetalgo = 1.0
+set jetradius = 0.4
+set ptj = 50.0
+set maxjetflavor = 5
+# `change_etaj_to_rapj.patch` changes the following to rapidity instead of pseudo-rapidity
+set etaj = 4.5
+set req_acc_FO 0.01
+set pineappl True
+done
+quit
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/metadata.txt b/nnpdf31_proc/NNPDF40_JETS_BBAR/metadata.txt
@@ -0,0 +1,7 @@
+description=cross section at 14 TeV
+x1_label=
+x1_label_tex=$$
+x1_unit=GeV
+y_label=dsig/d
+y_label_tex=$\mathrm{d}\sigma/\mathrm{d} $
+y_unit=pb
diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/output.txt b/nnpdf31_proc/NNPDF40_JETS_BBAR/output.txt
@@ -0,0 +1,7 @@
+set complex_mass_scheme True
+import model loop_qcd_qed_sm_Gmu-with_b_mass
+define j = p b b~
+# TODO: include NLO EW corrections
+generate p p > b b~ [QCD]
+output @OUTPUT@
+quit