Tag Archive | Fortran

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