diff --git a/nnpdf31_proc/NNPDF40_JETS_BBAR/analysis.f b/nnpdf31_proc/NNPDF40_JETS_BBAR/analysis.f new file mode 100644 index 00000000..1f8ed559 --- /dev/null +++ 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 new file mode 100644 index 00000000..f7bdc5d4 --- /dev/null +++ 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 new file mode 100644 index 00000000..9dcd7eab --- /dev/null +++ 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 new file mode 100644 index 00000000..6096997d --- /dev/null +++ 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 new file mode 100644 index 00000000..c31f8b2c --- /dev/null +++ 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 new file mode 100644 index 00000000..cfe66461 --- /dev/null +++ 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 new file mode 100644 index 00000000..a33c66dc --- /dev/null +++ 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