Page MenuHomeHEPForge
Files F11221632

No OneTemporary
Actions

View Options

Index: branches/DMNLO/run_dmnlo/stst2QQ/stst2QQ_MG5.F
===================================================================
--- branches/DMNLO/run_dmnlo/stst2QQ/stst2QQ_MG5.F (revision 0)
+++ branches/DMNLO/run_dmnlo/stst2QQ/stst2QQ_MG5.F (revision 1093)
@@ -0,0 +1,498 @@
+* Template for a new process @ NLO
+* DM@NLO Project - B.Herrmann, K.Kovarik, M.Klasen, Q.Le Boulc`h, J.Harz, P.Steppeler, M. Meinecke, S. Schmiemann
+
+#include "stst2QQ_Couplings.F"
+#include "stst2QQ_Kinematics.F"
+
+********************************
+ subroutine stst2qq_cs(renscale,scalevarpara,Pcm,iflag,myresult)
+ implicit none
+********************************
+
+
+#include "../util/DMNLO_LoopIntegrals.h"
+#include "stst2QQ_Model.h"
+#include "stst2QQ_Kinematics.h"
+#include "stst2QQ_GenCouplings.h"
+#include "stst2QQ_GenCounterterms.h"
+
+ integer ndim, ncomp, iflag(15), inptype, inpgen, iflux, i, j, k, platzhalter, u
+c parameter (ndim = 2)
+ parameter (ncomp = 1)
+
+ double precision Pcm,myresult(9), sinth, costh, avgfac, massth
+ double precision renscale, scalevarpara(5)
+ double precision result(9), error(9)
+ double precision result23(ncomp),error23(ncomp)
+ double precision cut
+
+ double precision theta,thetastep
+
+ double precision AS,zetaDC(2),ASDR
+
+ character buff(4)
+ double precision p(0:3,4)
+
+* subroutines
+ external DMNLO_ModelPara, DMNLO_ModelIni, DMNLO_ModelDigest, Init_RenScheme
+ external DMNLO_CalcPifth, DMNLO_CalcDeltaMb, DMNLO_CalcPiftest
+ external ASzetaDecCoe
+
+ QScale = renscale
+ Qscalealphas = Qscale
+
+ ischeme = 1
+
+** flags assignment
+c number of Sfermion1
+ isf1 = iflag(1)
+
+c type of initial state partical
+ itt1 = iflag(2)
+
+! number of family of initial state partical
+ igen1 = iflag(3)
+
+! number of Sfermion2
+ isf2 = iflag(4)
+
+! type of initial state anti-partical
+ itt2 = iflag(5)
+
+! number of family of initial state anti-partical
+ igen2 = iflag(6)
+
+! type of final state partical
+ ftt1 = iflag(7)
+
+! number of family of final state partical
+ fgen1 = iflag(8)
+
+! type of final state anti-partical
+ ftt2 = iflag(9)
+
+! number of family of final state anti-partical
+ fgen2 = iflag(10)
+
+
+! imass = 0- take masses from MicrOmegas, 1-diagonalize here
+ imass = iflag(11)
+
+! iflux = 0 - calculate v.sigma, 1- standard cross-section sigma
+ iflux = iflag(12)
+
+! itree = 0- full one loop, 1- only tree level (quicker)
+ itree = iflag(13)
+
+
+ if ((itt1.eq.itt2).AND.(ftt1.eq.ftt2))then
+ !Higgs
+ imin=1
+ imax=4
+ jmin=1
+ jmax=4
+ !Vektorboson
+ bmin=2
+ bmax=3
+ cmin=2
+ cmax=3
+! write(*,*)'This do work by now'
+!
+ if(itt1.eq.ftt1) then
+ ! write(*,*)'S-channel and t-channel with neutral particles are activated'
+ !Gaugino
+ pmin=1
+ pmax=4
+ qmin=1
+ qmax=4
+ else
+ !write(*,*)'S-channel with neutral particles is activated and and t-channel with Chargino is activated'
+ !Gaugino
+ pmin=5
+ pmax=6
+ qmin=5
+ qmax=6
+ endif
+ elseif ((itt1.eq.ftt1).AND.(itt2.eq.ftt2)) then
+ if(itt1.eq.itt2-1) then
+ !write(*,*)'S-channel with charged particles is activated and and t-channel with neutral particles is activated'
+ elseif (itt2.eq.itt1-1) then
+ !write(*,*)'Changing up and downtype quark as particle and antiparticle in both state'
+ platzhalter = itt1
+ itt1 = itt2
+ itt2 = platzhalter
+ platzhalter = ftt1
+ ftt1 = ftt2
+ ftt2=platzhalter
+ endif
+ !setzten der Propagator-Indizes
+ !Higgs
+ imin=5
+ imax=6
+ jmin=5
+ jmax=6
+ !Vekotrboson
+ bmin=4
+ bmax=4
+ cmin=4
+ cmax=4
+ !Gaugino
+ pmin=1
+ pmax=4
+ qmin=1
+ qmax=4
+
+ elseif ((itt2.eq.ftt1).AND.(itt2.eq.ftt1)) then
+ if(itt2.eq.itt1-1) then
+ !write(*,*)'Changing up and downtype quark as particle and antiparticle in initial state'
+ platzhalter = itt1
+ itt1 = itt2
+ itt2 = platzhalter
+ elseif (ftt2.eq.ftt1-1) then
+ !write(*,*)'Changing up and downtype quark as particle and antiparticle in final state'
+ platzhalter = ftt1
+ ftt1 = ftt2
+ ftt2=platzhalter
+ endif
+ !setzten der Propagator-Indizes
+ !Higgs
+ imin=5
+ imax=6
+ jmin=5
+ jmax=6
+ !Vekotrboson
+ bmin=4
+ bmax=4
+ cmin=4
+ cmax=4
+ !Gaugino
+ pmin=1
+ pmax=4
+ qmin=1
+ qmax=4
+
+ else
+ write(*,*)'This do not work - please enter first the uptype-quark and then the downtype quark'
+ endif
+
+
+! ===== Set model parameters =====
+! reading the MSSM parameters from micrOmegas and setting the masses & mixings & all shortforms e.g. TB,SA
+
+ call DMNLO_ModelPara()
+ call DMNLO_ModelIni()
+! call DMNLO_ModelDigest
+
+ if (scalevarpara(1).eq.1) then
+! print*,"scalevariation"
+ QScale = scalevarpara(2)
+! switch for alpha_s scale
+ Qscalealphas = scalevarpara(2)
+! Qscalealphas = renscale
+ Af(3,3) = scalevarpara(3)
+ Af(4,3) = scalevarpara(4)
+ endif
+
+! write(*,*)'TEST2'
+! ===== set sfermion-widths to zero/10 ===========
+ do i = 1,2
+ do j = 3,4
+ do k = 1,3
+
+ if(WSf(i,j,k).lt.10.0) then
+ !see Warning below for purpose of this
+ !print*,'WARNING: MicrOMEGAs sQuark-width small (=',WSf(i,j,k),'GeV )=> set width WSf(',i,',',j,',',k,') to 10 GeV'
+ WSf(i,j,k)=10.0
+ endif
+
+ enddo
+ enddo
+ enddo
+
+
+
+
+!-----------------Comparing with MadGraph:-------------------------------------------------------------------------
+
+!Pass in the check_sa.f file after the declaration of the 4-Momenta and declare u as an Integer
+! Open the data file
+ open (u, FILE='points.dat', STATUS='OLD')
+ do i=1,4
+ read (u,'(i2,1x,5e15.7)') k, p(0,i),p(1,i),p(2,i),p(3,i)
+ enddo
+! Close the file
+ close (u)
+
+ do i=1,4
+ write (*,'(i2,1x,5e15.7)') i, p(0,i),p(1,i),p(2,i),p(3,i)
+ enddo
+
+ theta= dacos((p(1,1)*p(1,3)+p(2,1)*p(2,3)+p(3,1)*p(3,3))/(dsqrt(p(1,1)**2+p(2,1)**2+p(3,1)**2)*dsqrt(p(1,3)**2+p(2,3)**2+p(3,3)**2)))
+ sqrtS= p(0,1)+p(0,2)
+
+!-----------------------------------------------Ende MadGraph Vergleich--------------------------------------------------------------
+
+! ===== Definition of Renormalization Scheme, Scale & Input =====
+ call Init_RenScheme(1)
+ ! -------------------Ausgabe der parameter----------------------
+ write(*,*) 'theta',theta
+ write(*,*)'sqrtS',sqrtS
+ write(*,*)'pcm',Pcm
+ write(*,*)'msf1',MSf(isf1,itt1,igen1)
+ write(*,*)'msf2',MSf(isf2,itt2,igen2)
+ write(*,*)'mf1',Mf(ftt1,fgen1)
+ write(*,*)'mf2',Mf(ftt2,fgen2)
+
+! ===== Renormalization Scales for Loop integrals =====
+! Divergent UV & IR poles for Loop integrals
+ UVdiv = 0d0
+ IRdiv = 0d0
+
+! Renormalization scale \mu & switch between dim.reg and mass regularization for IR divergence - xiIR
+! xiIR enables a check with mass regularization if set to 0d0 -> DEFAULT is xiIR = 1d0 !!!!
+ xiIR = 1d0
+ muSc = QScale
+
+! EpsPoles = -2,-1,0,1,... steers LoopFunctions
+! e.g. EpsPoles = -2/-1 lets all LoopFunctions return the IR diveregent coefficient corresponding to double/single pole
+! checking UV divergence requires setting EpsPole=1 and UVdiv=1 (loop amplitude should be zero or extremely small)
+ EpsPole = 0d0
+
+! flux factor (in mykinematics.h)
+ fluxflag = iflux
+! write(*,*)'TEST5'
+! write(*,*)itt3, igen3, itt4,igen4
+! Generic couplings (here the type and generation is set !!! )
+ call stst2qq_SetCouplings
+! write(*,*)'TEST6'
+! write(*,*)itt3, igen3, itt4,igen4
+! FA kinematics & reduced masses
+ call stst2qq_SetKinematics
+
+! self-energies & counterterms
+! call Template_SetCounterTerms()
+
+*******************************************
+
+! Mass threshold
+ massth = sqrtS*(muu3+muu4)
+
+
+
+
+ if (sqrtS.le.massth) then
+ do i=1,9
+ result(i) = 0d0
+ enddo
+ result23(1) = 0d0
+ else
+
+
+
+! Compute differential cross-section
+
+
+ call stst2qq_DiffCS(result, dcos(theta))
+
+
+ result23(1) = 0d0
+
+
+! Copy to output vectors
+
+ if (itree.eq.0) then
+
+ do i=1,6
+ myresult(i) = result(i)
+ enddo
+
+ myresult(7) = result23(1)
+ myresult(8) = error23(1)
+ myresult(9) = result(6) + result23(1)
+
+ else
+
+ do i=1,9
+ myresult(i) = result(i)
+ enddo
+
+ endif
+
+ endif
+
+ end
+
+**************************************************************
+
+ subroutine stst2qq_IntCS(result,error,mytree)
+
+ implicit none
+
+#include "stst2QQ_Kinematics.h"
+#include "stst2QQ_Model.h"
+
+ double precision result(9), error(9)
+
+ double precision costhmin, costhmax, relaccuracy, absaccuracy
+ double precision kappa
+ integer c, neval, fail, ncomp, mytree
+ external stst2qq_TreeDiffCS, stst2qq_DiffCS, Patterson
+
+ ncomp = 7
+
+! Integration limits for cos(theta)
+ costhmin = -1d0
+ costhmax = 1d0
+
+! Set requested accuracy
+ absaccuracy = 1d-12
+ relaccuracy = 1d-5
+
+! Full calculation
+ if (mytree.eq.0) then
+
+ call Patterson(ncomp, costhmin, costhmax, stst2qq_DiffCS,
+ & relaccuracy, absaccuracy, neval, fail, result, error)
+
+! Only tree-level calculation
+ else
+
+ call Patterson(ncomp, costhmin, costhmax, stst2qq_TreeDiffCS,
+ & relaccuracy, absaccuracy, neval, fail, result, error)
+
+ endif
+
+
+c if( fail .ne. 0 )
+c & print *, "Failed to reach the desired accuracy."
+
+ end
+
+*************************************
+
+ subroutine stst2qq_DiffCS(result, costh)
+
+ implicit none
+
+#include "stst2QQ_Kinematics.h"
+#include "stst2QQ_Model.h"
+#include "stst2QQ_GenCouplings.h"
+
+ double precision costh, intfac, intfacMO, avgfac, testmin, testmax
+ double precision result(9)
+
+ double precision kappa
+ double complex stst2qq_M2Tree
+
+ double precision test
+
+
+
+! write(*,*) 'cos(theta)', costh
+! Mandelstam variables t=(p1-k1)^2 , u=(p1-k2)^2
+ tman = sqrtS**2*(muu1**2 + muu3**2 - 0.5d0*(1d0+muu1**2-muu2**2)*(1d0+muu3**2-muu4**2) +
+ & 0.5d0*kappa(1d0,muu1**2,muu2**2)*kappa(1d0,muu3**2,muu4**2)*costh)
+ tred = (muu1**2 + muu3**2 - 0.5d0*(1d0+muu1**2-muu2**2)*(1d0+muu3**2-muu4**2) +
+ & 0.5d0*kappa(1d0,muu1**2,muu2**2)*kappa(1d0,muu3**2,muu4**2)*costh)
+
+ uman = sqrtS**2*(muu1**2 + muu2**2 + muu3**2 + muu4**2 - 1d0 - tred)
+ ured = muu1**2 + muu2**2 + muu3**2 + muu4**2 - 1d0 - tred
+
+
+! Integration factor for 2->2 phase-space integration in CMS
+! (the azimuthal integration constant 2pi is inlcuded in intfac)
+
+ intfac = 2*pi*(kappa(1d0,muu3**2,muu4**2)/2d0)/(4*(2*pi)**2)
+ intfacMO = 2*pi*(kappa(sqrtS**2,Mfm(4)**2,Mfm(4)**2)/(2d0*sqrtS))/(4*(2*pi)**2*sqrtS)
+
+! Flux factors for different units of cross-section
+ if (fluxflag.eq.1) then
+! in units of pb
+cc flux = hbar_c2/(4d0*(kappa(sqrtS**2,muu1**2*sqrtS**2,muu2**2*sqrtS**2)/(2d0*sqrtS))*sqrtS)
+! in units of GeV^-2
+ flux = 1d0/(4d0*(kappa(1d0,muu1**2,muu2**2)/2d0)*sqrtS**2)
+ else
+
+! in units of cm^3/sec
+c flux = 2.9979d-26*hbar_c2/sqrtS**2
+! in units of GeV^-2
+ flux = 1d0/sqrtS**2
+ endif
+!!distinguish between identical and non identical particels
+! if(indk.ge.indl) then
+! if(indk.lt.5)then
+!! Average factor+Symmetryfactor from equal final particels
+! avgfac = 1d0/4d0*1d0/2d0
+! end if
+! else
+! Average factor
+ avgfac = 1d0
+! end if
+! ***************************************
+! contributions to the matrix element
+! ***************************************
+! (1-tree, 2-full one loop virtual, 3-vertex, 4-propagator, 5-box, 6-dipole)
+
+ result(1) = intfac*flux*avgfac*stst2qq_M2Tree()
+
+
+ end
+
+*************************************
+
+ subroutine stst2qq_TreeDiffCS(result, costh)
+
+ implicit none
+
+#include "stst2QQ_Kinematics.h"
+#include "stst2QQ_Model.h"
+#include "stst2QQ_GenCouplings.h"
+
+ double precision costh, intfac, intfacMO, avgfac
+ double precision result(9)
+
+ double precision kappa, stst2qq_M2TreeM0, stst2qq_M2TreeHO, stst2qq_M2TreePerc
+ double complex stst2qq_M2Tree
+
+
+
+! Mandelstam variables t=(p2-k1)**2 , u=(p1-k1)**2
+! Scattering angle theta between p1 and k1
+
+ tman = sqrtS**2*(muu2**2 + muu3**2 - 0.5d0*(1d0-muu1**2+muu2**2)*(1d0-muu4**2+muu3**2) -
+ & 0.5d0*kappa(1d0,muu1**2,muu2**2)*kappa(1d0,muu3**2,muu4**2)*costh)
+ tred = (muu2**2 + muu3**2 - 0.5d0*(1d0-muu1**2+muu2**2)*(1d0-muu4**2+muu3**2) -
+ & 0.5d0*kappa(1d0,muu1**2,muu2**2)*kappa(1d0,muu3**2,muu4**2)*costh)
+
+ uman = sqrtS**2*(muu1**2 + muu2**2 + muu3**2 + muu4**2 - 1d0 - tred)
+ ured = muu1**2 + muu2**2 + muu3**2 + muu4**2 - 1d0 - tred
+
+
+
+! Integration factor for 2->2 phase-space integration in CMS
+! the azimuthal integration constant 2pi is inlcuded in intfac
+ intfac = 2*pi*(kappa(1d0,muu3**2,muu4**2)/2d0)/(4*(2*pi)**2)
+
+! Flux factor
+ if (fluxflag.eq.1) then
+! in units of pb
+cc flux = hbar_c2/(4d0*(kappa(sqrtS**2,muu1**2*sqrtS**2,muu2**2*sqrtS**2)/(2d0*sqrtS))*sqrtS)
+! in units of GeV^-2
+ flux = 1d0/(4d0*(kappa(1d0,muu1**2,muu2**2)/2d0)*sqrtS**2)
+ else
+
+! in units of cm^3/sec
+c flux = 2.9979d-26*hbar_c2/sqrtS**2
+! in units of GeV^-2
+ flux = 1d0/sqrtS**2
+ endif
+ avgfac = 1d0
+
+! no flux and averaging factor!
+
+ result(1) = intfac*flux*avgfac*stst2qq_M2Tree()
+
+
+ end
+
+#include "stst2QQ_M2Tree.F"

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:40 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111224
Default Alt Text
(13 KB)

Event Timeline