MacOSX: Solution to PYTHIA / HBOOK Problem in Cernlib
Last month, I was trying to compile some pythia codes to calculate snutau (scalar tau) cross section, Pt and eta distributions. But unfortunately my trusted g77 compiler insistently rejected to compile my code with the following errors:
Undefined symbols:
"_pylist_", referenced from:
_MAIN__ in ccJ7qu36.o
"_pyevnt_", referenced from:
_MAIN__ in ccJ7qu36.o
"_pyinit_", referenced from:
_MAIN__ in ccJ7qu36.o
"_hropen_", referenced from:
_MAIN__ in ccJ7qu36.o
"_hlimit_", referenced from:
_MAIN__ in ccJ7qu36.o
"_pystat_", referenced from:
_MAIN__ in ccJ7qu36.o
_MAIN__ in ccJ7qu36.o
"_pydata_", referenced from:
___g77_forceload_0.0 in ccJ7qu36.o
"_hfill_", referenced from:
_MAIN__ in ccJ7qu36.o
_MAIN__ in ccJ7qu36.o
_MAIN__ in ccJ7qu36.o
"_hrout_", referenced from:
_MAIN__ in ccJ7qu36.o
"_hbook1_", referenced from:
_MAIN__ in ccJ7qu36.o
_MAIN__ in ccJ7qu36.o
_MAIN__ in ccJ7qu36.o
"_hrend_", referenced from:
_MAIN__ in ccJ7qu36.o
ld: symbol(s) not found
collect2: ld returned 1 exit status
After a little research about HBOOK technology, i discovered that it’s a quite old component even before Fortran 77 and the only way to make it run in this way is to include Cernlib components while compiling. But that didn’t work either for MacOsX because probably MacosX version of Cernlib doesn’t support HBOOK libraries although it’s possible for Linux.
Anyway why do we have to use HBOOK? At the end of the day, we just want to have some data files to draw using Paw or Root and it’s possible using Pythia’s own components. So the principle must be “Use Pythia Commands to get an Output file because that’s the easiest way!” Here is my solution: I changed all HBOOK commands to PY..s in my code below and it worked perfect.
Write “g77 -o mycode.x -w mycode.f libpythia6421.a” to compile.
C...All real arithmetic in double precision.
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
C...Three Pythia functions return integers, so need declaring.
INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers
C...(left- and righthanded SUSY, technicolor, excited fermions,
C...extra dimensions).
PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
&KEXCIT=4000000,KDIMEN=5000000)
C...EXTERNAL statement links PYDATA on most machines.
EXTERNAL PYDATA
DIMENSION IHI(10)
C...Commonblocks.
C...The event record.
COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
C...Parameters.
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
C...Particle properties + some flavour parameters.
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
C...Decay information.
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
C...Selection of hard scattering subprocesses.
COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
C...Parameters.
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
C...Supersymmetry parameters.
COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
C....Real definitions...
REAL CTETM,CTETP,ETAUP,ETAUM,PTAUP,PTAUP2,PTAUM,PTAUM2
REAL PTAUMX, PTAUMY, PTAUMZ
REAL PTAUPX, PTAUPY, PTAUPZ
REAL PTX, ETAX, XMTAU
C ....rapidity: YX
C--------------------------------------------
C...First section: initialization.
NEV=10000
C...Select generic SUSY generation.
C...39:all MSSM processes except Higgs prod.
C...41:stop pair(ISUB=261-265), 42:slepton pair(201-214), 45:sbottom(281-296)
C...ISUB=210,211,212 for slepton+sneutrino
MSEL=0
MSUB(214)=1
C...Set SUSY parameters in SUGRA scenario.
IMSS(1)=2 !mSUGRA parameters given to PYTHIA
RMSS(8)=230D0 !m_0
RMSS(1)=360D0 !m_1/2
RMSS(5)=10D0 !tan(beta)
RMSS(4)=1D0 !sign(mu)
RMSS(16)=0D0 !A_0
C...Channels
do kk=1949,1974
MDME(kk,1)=0
enddo
MDME(1950,1)=1
C...Channels for W boson
do kw=206,208
MDME(kw,1)=0
enddo
C...If interested only in cross sections and resonance decays:
C...switch on/off initial and final state radiation,
C...multiple interactions and hadronization.
MSTP(11)=0 ! 1:QED radiation
MSTP(61)=0 ! 2:ISR
MSTP(71)=0 ! 1:FSR
MSTP(81)=0 ! 1:multiple int.
MSTP(111)=0 ! 1:hadronization
C...Initialization for the LHC.
CALL PYINIT('CMS','e-','e+',3000D0)
C...List resonance data: decay channels, widths etc.
CALL PYSTAT(2)
C...Book Histograms
CALL PYBOOK(10,'PT',100,0D0,1000D0)
CALL PYBOOK(20,'ETA',100,-5D0,5D0)
CALL PYBOOK(30,'MTAUTAU',100,0D0,2000D0)
C--------------------------------------------------
C...Second section: event loop.
C...Loop over the number of events.
DO 200 IEV=1,NEV
IF(MOD(IEV,500).EQ.0) WRITE(6,*)
& 'Now at event number',IEV
C...Event generation.
CALL PYEVNT
C...List first few events.
IF(IEV.LE.5) CALL PYLIST(1)
C...Fill the masses of interesting (s)particles.
C...Fill pt of particle
DO I=1,N
C...Catch tau- lepton
IF((K(I,2).EQ.15).AND.(K(I,1).EQ.1)) THEN
PTX=SQRT(P(I,1)**2+P(I,2)**2)
CTETM=P(I,3)/SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
ETAX=-LOG(TAN(ACOS(CTETM)/2.))
etaum=P(I,4)
ptaum2=P(I,1)**2+P(I,2)**2+P(I,3)**2
ptaum=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
ptaumx=P(I,1)
ptaumy=P(I,2)
ptaumz=P(I,3)
CALL PYFILL(10,DBLE(PTX),1D0)
CALL PYFILL(20,DBLE(ETAX),1D0)
C...... CALL HFILL(10,PTX,0.,1.)
C...... CALL HFILL(20,ETAX,0.,1.)
ENDIF
C...Catch tau+ lepton
IF((K(I,2).EQ.-15).AND.(K(I,1).EQ.1)) THEN
etaup=P(I,4)
ptaup2=P(I,1)**2+P(I,2)**2+P(I,3)**2
ptaup=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
ptaupx=P(I,1)
ptaupy=P(I,2)
ptaupz=P(I,3)
ENDIF
XMTAU=SQRT((ETAUP+ETAUM)**2-(ptaup2+ptaum2
. +2.0*(ptaupx*ptaumx+ptaupy*ptaumy+ptaupz*ptaumz)))
C... CALL HFILL(30,XMTAU,0.,1.)
CALL PYFILL(30,DBLE(XMTAU),1D0)
ENDDO
C...End of documentation and event loops.
200 CONTINUE
C--------------------------------------------
C...Third section: produce output and end.
C...Cross section table.
CALL PYSTAT(1)
C...Histogram close
C... CALL HROUT(0,ICYCLE,' ')
C... CALL HREND('SUSY')
C...Histograms.
OPEN(11,file='pt.dat',STATUS='unknown')
IHI(1)=10
CALL PYDUMP(3,11,1,IHI)
CLOSE(11)
OPEN(22,file='eta.dat',STATUS='unknown')
IHI(1)=20
CALL PYDUMP(3,22,1,IHI)
CLOSE(22)
OPEN(33,file='MTau.dat',STATUS='unknown')
IHI(1)=30
CALL PYDUMP(3,33,1,IHI)
CLOSE(33)
CALL PYHIST
END
Author: Sinan Kuday
Click image to visit my personal site.
Kuday WeBlog 
Recent Comments