Commit c83b077d authored by Ewelina Maria Lobodzinska's avatar Ewelina Maria Lobodzinska
Browse files

prepare the interface to run with Superchic 4.02 - obo Junichi

parent fa45fb07
......@@ -62,6 +62,7 @@ protected:
std::string m_beam;
std::string m_outtg;
std::string m_sfaci;
std::string m_diff;
std::string m_an;
std::string m_az;
std::string m_rz;
......@@ -76,7 +77,7 @@ protected:
unsigned int m_inccall;
unsigned int m_itend;
int m_iseed;
int m_s2int;
//int m_s2int;
std::string m_genunw;
int m_nev;
......
evgenConfig.description = "SuperChic4 MC gamma + gamma pp collisions at 13000 GeV to 2 gamma ALP mediated"
evgenConfig.keywords = ["2photon","2photon"]
#evgenConfig.weighting = 0
evgenConfig.contact = ["gen.tateno@cern.ch"]
if not os.path.exists('inputs'):
os.makedirs('inputs')
if not os.path.exists('evrecs'):
os.makedirs('evrecs')
from Superchic_i.Superchic_iConf import Superchic_i
genSeq += Superchic_i("Superchic")
genSeq.Superchic.McEventKey = "GEN_EVENT"
evgenConfig.generators += ["Superchic"]
from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
_evgenstream = AthenaPoolOutputStream("StreamEVGEN")
_evgenstream.ItemList = ["2101#*","133273#GEN_EVENT"]
del _evgenstream
# TODO: Sort out proper param setting based on runArgs.ecmEnergy
if int(runArgs.ecmEnergy) != 13000:
evgenLog.error(" Set beam energy in JO initialization with parameter rts ")
sys.exit(1)
genSeq.Superchic.Initialize = [
"rts 13000d0", # set the COM collision energy (in fortran syntax)
"isurv 4", # Model of soft survival
"intag 'in5'", # for input files
"PDFname 'MMHT2014lo68cl'", # PDF set name
"PDFmember 0", # PDF member
"proc 68", # Process number (59 = gg->gg, 56 = gg->ee, 68 = gg->a->gg )
"beam 'prot'", # Beam type ('prot', 'ion')
"outtg 'out'", # for output file name
"sfaci .true.", # Include soft survival effects
"diff 'el'", # interaction: pure elastic (el), single prodton dissociation (ds), double dissociation (dd)
"an 208d0",
"az 82d0",
"rz 6.68d0",
"dz 0.447d0",
"rn 6.7d0",
"dn 0.55d0",
"ionqcd 'coh'",
"ncall 10000", # Number of calls for preconditioning
"itmx 10", # Number of iterations for preconditioning
"prec 0.5d0", # precision
"ncall1 10000",
"inccall 10000",
"itend 1000",
"iseed 34",
"genunw .true.",
"nev 1", # 1 or a few
"erec 'hepevt'",
"readwt .false.",
"wtmax 0d0",
"ymin -2.4d0", # Minimum object rapidity Y_X
"ymax 2.4d0", # Maximum object rapidity Y_X
"mmin 6d0", # Minimum object mass M_X
"mmax 500d0", # Maximum object mass M_X
"gencuts .true.", # Generate cuts below
"scorr .true.",
"fwidth .true.",
"ptxmax 100d0",
"ptamin 3.0d0", # Minimum pT of outgoing object a (gamma)
"ptbmin 3.0d0", # Minimum pT of outgoing object b (gamma)
"etaamin -2.4d0", # Minimum eta of outgoing object a
"etaamax 2.4d0", # Maximum eta of outgoing object a
"etabmin -2.4d0", # Minimum eta of outgoing object b
"etabmax 2.4d0", # Maximum eta of outgoing object b
"acoabmax 100d0",
"rjet 0.6d0",
"jalg 'antikt'",
"malp 1000d0", # ALP mass (GeV)
"gax 0.001d0", # ALP coupling (GeV^-1)
"alpt 'ps'" # ALP type (scalar - 'sc', pseudoscalar - 'ps')
]
This diff is collapsed.
ccc reads input parameters and calculated needed parameters
subroutine calcparam
implicit double precision (a-z)
integer i,j,k,nt,ns
integer i1,i2,nb,ib
integer nhistmax
integer outl
integer iinc,ncallu
logical histol
integer i,j,k
character*100 dum
integer idum
COMMON /ranno/ idum
......@@ -70,17 +65,15 @@ ccc reads input parameters and calculated needed parameters
include 'mn.f'
include 'rech.f'
include 'ptXcuts.f'
include 'nchan.f'
include 'opac.f'
include 'mcharg.f'
include 'varsi.f'
include 'spA.f'
include 'sAA.f'
c include 'mphi.f'
include 'diff.f'
include 'diss.f'
include 'enew.f'
include 'gamma.f'
c include 'init.f'
c include 'initpars.f'
c include 'initsud.f'
ccccccc
open(unit=1,file='input.DAT',status='old', action='read')
......@@ -100,6 +93,11 @@ ccccccc
read(1,*)beam
read(1,*)outtag
read(1,*)sfaci
read(1,*)dum
read(1,*)dum
read(1,*)dum
read(1,*)diff
c read(1,*)elcoll
read(1,*)dum
read(1,*)an
read(1,*)az
......@@ -119,8 +117,6 @@ ccccccc
read(1,*)itend
read(1,*)iseed
read(1,*)dum
read(1,*)s2int
read(1,*)dum
read(1,*)dum
read(1,*)dum
read(1,*)genunw
......@@ -230,14 +226,24 @@ ccccccc
read(1,*)dum
read(1,*)mcharg
read(1,*)mneut
cccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccc
c main/superchic.f
diffsd=diff
if(diff.eq.'sda'.or.diff.eq.'sdb')then
diff='sd'
else
diffsd='n'
endif
offshell=.false.
forward=.false.
call length(outtag,outl)
C call length(outtag,outl)
open(45,file='evrecs/evrec'//outtag(1:outl)//'.dat')
C open(45,file='evrecs/evrec'//outtag(1:outl)//'.dat')
wmax=0d0
evnum=0
......@@ -272,8 +278,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccc
mtau=1.77682d0
mpip=0.13957018d0
mkp=0.493677d0
c mphi=1.019461d0
alpha=7.2974d-3
pi=dacos(-1d0)
conv=389379d3
zi=(0d0,1d0)
......@@ -322,18 +328,61 @@ cccccccccccc
jdahep(2,i)=0
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccc Init LHAPDF
cccccccccccccccccccccccccccccccccccccccccccccccccccc
pi=dacos(-1d0)
cccccccccccccccccccccccccc
c main/superchic.f
call inpdf
print*,'Read PDF successfully...'
cccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccc
call supinit
s2int=8
if(beam.eq.'ionp')s2int=16
if(diff.eq.'el'.and.gamma.eqv..true.)s2int=16
call header
call gaminit
call gaminit_comb
call readcoh
call gaussinit
cccccccccccccccccccccccccc
c main/superchic.f
if(diff.eq.'sd'.or.diff.eq.'dd')then
if(offshell.eqv..false.)then
print*,'Dissociation not currently supported'//
& ' for this process/beam - STOP'
stop
endif
if(erec.eq.'hepevt'.or.erec.eq.'hepmc')then
print*,'Dissociation currently only supported with LHE'
endif
endif
elcoll=.false.
difftot=.false.
if(diff.eq.'el')then
diss1=.false.
diss2=.false.
elseif(diff.eq.'sd')then
ndim=ndim+1
if(genunw)elcoll=.true.
elseif(diff.eq.'dd')then
diss1=.true.
diss2=.true.
ndim=ndim+2
elseif(diff.eq.'tot')then
ndim=ndim+3
difftot=.true.
endif
if(diff.eq.'tot'.or.diff.eq.'sd'.or.diff.eq.'dd')then
call apfelinit
call calcs2diss
endif
cccccccccccccccccccccccccc
c main/superchic.f
if(mes)then
call calcmes
......@@ -341,6 +390,7 @@ cccccccccccccccccccccccccc
endif
ccccccccccccccccccccccccc
c main/superchic.f
if(beam.eq.'el')then
if(sfaci)then
......@@ -350,35 +400,32 @@ ccccccccccccccccccccccccc
endif
endif
ccccccccccccccccccccccccccccccccccccccccccccccccccc
call initparsr(isurv) ! Initialise soft survival parameters
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
call calcop ! proton opacity
cccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc Initialises grids for skewed PDFs and survival factors
c init/init.f
call initpars(isurv) ! Initialise soft survival parameters
call calcop ! proton opacity
call calcscreen ! screening amplitude
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
call initsud ! init sudakov factor
call calcsud ! sudakov factor
call inithg ! init Skewd PDF
call calchg ! skewed PDF
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
print*,'Calculated opacity,screening ampl.,sudakov, & skewed PDF'
call readcoh
call dd
call readscreen
call sdcoh
call apfelinit
call sdincoh
c call inpdf
c call supinit
c call header
c stop
surv=1d0
if(beam.eq.'prot'.or.ionqcd.eq.'coh'.or.ionqcd.eq.'incoh')then
call initparsr(isurv)
call readscreen
if(beam.eq.'prot'.or.ionqcd.eq.'incoh')surv=1d0/norm**2
endif
call calcsud ! sudakov factor
call calchg ! skewed PDF
cccccccccccccccccccccccccccccccccccccccccccccccccccc
print*,'Calculated opacity,screening ampl.,sudakov, & skewed PDF'
if(qcd)then
call calcsud
call calchg
endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c main/superchic.f
C rts ... CMS collision energy (GeV)
s=rts**2
......@@ -431,17 +478,7 @@ c call header
pdgid(4)=pdgid(2)
endif
if(beam.eq.'ion'.or.beam.eq.'ionp')call ioninit
if(beam.eq.'ionp')then
rts=rtspa
s=spa
elseif(beam.eq.'ion')then
rts=rtsaa
s=saa
endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
return
end
......@@ -76,422 +76,87 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
include 'varsi.f'
include 'spA.f'
include 'sAA.f'
c include 'mphi.f'
include 'diff.f'
include 'diss.f'
include 'enew.f'
include 'gamma.f'
cccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccc
forward=.false.
ccccccccccccccccccccccc
c main/superchic.f
call length(outtag,outl)
open(45,file='evrecs/evrec'//outtag(1:outl)//'.dat')
wmax=0d0
evnum=0
if(genunw)then
else
readwt=.false.
endif
if(readwt)wmax=wtmax
if(erec.eq.'hepmc')then
erech=.true.
erec='lhe'
endif
iw=0
gf=1.16639d-5
v=dsqrt(1d0/dsqrt(2d0)/gf)
mt=173d0
mb=4.75d0
mc=1.4d0
mmu=0.10566d0
mpsi=3.096916d0
mpsip=3.686109d0
mups=9.46030d0
mchic0=3.41475d0
mchib0=9.85944d0
mp=0.938272046d0
mn=0.939565413d0
mw=80.318d0
me=0.511d-3
mtau=1.77682d0
mpip=0.13957018d0
mkp=0.493677d0
c mphi=1.019461d0
alpha=7.2974d-3
pi=dacos(-1d0)
conv=389379d3
zi=(0d0,1d0)
mup=0.062d0
md=0.083d0
ms=0.215d0
rmf1( 1) = 1d-10
rmf1( 2) = me
rmf1( 3) = 1d-10
rmf1( 4) = mmu
rmf1( 5) = 1d-10
rmf1( 6) = mtau
rmf1( 7) = 0.062d0
rmf1( 8) = 0.083d0
rmf1( 9) = mc
rmf1(10) = 0.215d0
rmf1(11) = mt
rmf1(12) = mb
rmf1( 1) = me
rmf1( 2) = mmu
rmf1( 3) = mtau
rmf1( 4) = md
rmf1( 5) = mup
rmf1( 6) = ms
rmf1( 7) = mc
rmf1( 8) = mb
rmf1( 9) = mt
mq=0d0
hel=1
mes=.false.
mfact='mx'
forward=.false.
decay2=.false.
decay3=.false.
decay4=.false.
decay6=.false.
ccccccccccccccccccccccccccccccccccccccccccccccc
cccc HEPEVT
ccccccccccccccccccccccccccccccccccccccccccccccc
do k=1,2
do j=1,4
phep(j,k)=q(j,k)
enddo
phep(5,k)=mass(k)
if(beam.eq.'el')then
phep(5,k)=me
elseif(beam.eq.'prot')then
phep(5,k)=mp
elseif(beam.eq.'ion'.or.beam.eq.'ionp')then
phep(5,k)=mion
endif
enddo
do k=1,20
do j=1,4
vhep(j,k)=0d0
enddo
enddo
isthep(1)=2
isthep(2)=2
isthep(3)=1
isthep(4)=1
jmohep(2,5)=2
jdahep(1,1)=0
jdahep(2,1)=0
jdahep(1,2)=0
jdahep(2,2)=0
jdahep(1,3)=0
jdahep(2,3)=0
jdahep(1,4)=0
jdahep(2,4)=0
ccccccccccccccccccccccccccccccccccccccccccccccc
ccccc Les Houches
ccccccccccccccccccccccccccccccccccccccccccccccc
nprup=1
idwtup=3
pdfsup(1)=0
pdfsup(2)=0
pdfgup(1)=0
pdfgup(2)=0
idprup=0
xwgtup=1d0
aqedup=alpha
ccc NEW LHE init
pdfsup(1)=247000
pdfsup(2)=247000
pdfgup(1)=0
pdfgup(2)=0
do k=1,2
do j=1,4
pup(j,k)=q(j,k)
enddo
pup(5,k)=mass(k)
if(beam.eq.'el')then
pup(5,k)=me
elseif(beam.eq.'prot')then
pup(5,k)=mp
elseif(beam.eq.'ion'.or.beam.eq.'ionp')then
pup(5,k)=mion
endif
enddo
istup(1)=-1
istup(2)=-1
istup(3)=1
istup(4)=1
mothup(1,1)=0
mothup(2,1)=0
mothup(1,2)=0
mothup(2,2)=0
mothup(1,3)=0
mothup(2,3)=0
mothup(1,4)=0
mothup(2,4)=0
mothup(1,5)=0
mothup(2,5)=0
icolup(1,1)=0
icolup(2,1)=0
icolup(1,2)=0
icolup(2,2)=0
icolup(1,3)=0
icolup(2,3)=0
icolup(1,4)=0
icolup(2,4)=0
icolup(1,5)=0
icolup(2,5)=0
do i=1,20
vtimup(i)=0
spinup(i)=9
enddo
do i=1,2
do j=1,5
jmohep(i,j)=mothup(i,j)
enddo
enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if(proc.eq.18.or.proc.eq.19.or.proc.eq.20)then
nup=11
elseif(proc.eq.54.or.proc.eq.55)then
nup=11
elseif(proc.eq.73.or.proc.eq.74.or.proc.eq.75)then
nup=11
elseif(proc.eq.76)then
nup=9
elseif(decay4)then
if(proc.eq.53)then
nup=10
else
nup=9
endif
elseif(decay6)then
nup=11
elseif(dps.eq.2.or.decay2)then
nup=7
elseif(dps.eq.3)then
nup=8
elseif(decay3)then
nup=9
elseif(dps.eq.1)then
nup=5
endif
do j=3,nup
do i=1,4
evrec(evnum,j,i)=q(i,j)
enddo
enddo
nhep=nup
nup=nup-2
if(proc.eq.68)then
nhep=7
nup=7
endif
cccccccccccccccccccccccccccccccccccccccccccccc
nhist=0
nhistmax=20
neff=0
neff0=0
do i=1,10
xu(i)=1d0
xl(i)=0d0
enddo