program pythia_event c********************************************************** c c Pythia_event c c program to call the pythia event generator and store c the event in a PISA standard event format. c c 6/10/97 - 2/11/99 Andrew Rose c c 2/17/99 Corrected version from Yuji Goto c 6/19/99 Corrected version from Yuji Goto c c Notes: taken from routine by Naohito Saito c c********************************************************** * implicit none include 'event.inc' c c include 'param_utg.inc' c include 'param_uts.inc' c include 'hepevt.inc' include 'gcflag.inc' include 'gen_pythia.inc' include 'guevgen.inc' include 'evntcode.inc' c Externals integer egzinit integer write_code_bank real*8 bimpact /100./ !minimun bias impact par. *KEND. c c utsgkine portion of code c integer i,ii,id_geant, irr real p(4),v(3),dumy,d logical lok character text*132,name*20 Integer ILP,Jlp,Iacc,NEVHEP,NPART, icycle real py_varstore(8) c c counters integer Iwrite c input values integer imsel, imsub(2,20) !msel, mstp's to change integer imstp(2,20) real ickin(2,4) !ckin's to change integer imrlu(2,6) !MRLU's to change real ntstore(10) !ntuple storage c integer istat real sqrts /500./ logical verbose /.true./ logical xverbose /.false./ logical zebrastore /.false./ ! logical for new zebra store function character*20 hist_file /'pythia.hist'/ integer nsize /10/ character*8 chform (10) data chform /'Event','Line','Status','PID','Parent', + 'Px','Py','Pz','Energy','Mass'/ real nstore(10) c c Par file section c CHARACTER*50 PAR_FILE NAMELIST /PYTH_PAR/ tnumevt, imsel, iMSUB, imstp, ickin, + imrlu, sqrts,verbose, xverbose, zebrastore, hist_file EXTERNAL LUDATA,PYDATA * initialize cernstuff parameter (NWPAWC=80000) COMMON /PAWC/H(NWPAWC) call hlimit(NWPAWC) c c Read in input name list c c par_file='pythia.par' open(unit=15,file=par_file,status='OLD', + err=997) read(15,nml=pyth_par,err=999) close(unit=15) c c Resolve input into Pythia control arrays c msel=imsel write (6,*)'MSEL value set to ', msel write(*,*)'Filling MSUB array.' do i=1,20 if (iMSUB(1,i).ne.0) then MSUB(iMSUB(1,i))=iMSUB(2,i) if (verbose) then write(*,*) 'MSUB(', iMSUB(1,i),')',iMSUB(2,i) endif endif enddo write(*,*)'Filling MSTP array.' do i=1,20 if (iMSTP(1,i).ne.0) then MSTP(iMSTP(1,i))=iMSTP(2,i) if (verbose) then write(*,*) 'MSTP(', iMSTP(1,i),')',iMSTP(2,i) endif endif enddo write(*,*) 'Filling CKIN array.' do i=1,4 ii=ickin(1,i) if (ickin(1,i).ne.0) then ckin(ii)=ickin(2,i) if (verbose) then write(*,*) 'CKIN(',ii,')',ickin(2,i) endif endif enddo write(*,*) 'Filling MRLU array.' do i=1,6 if (imrlu(1,i).ne.0) then mrlu(imrlu(1,i))=imrlu(2,i) if (verbose) then write(*,*) 'MRLU(',imrlu(1,i),')',imrlu(2,i) endif endif enddo c********************************************************* c c Initialization Section c Initialize ZEBRA bank c call MZEBRA(0) c c Initialize zebra writeout c if (zebrastore) then else open(unit=15,file='pythia.dat',status='NEW', + access='sequential', form='unformatted') call hropen(16,'LUN1',hist_file,'N',1024,istat) endif * declare ntuple call hbookn(1000,'Pythia Particles',nsize,'LUN1',8192,chform) c write file header write(15) sqrts,tnumevt c c Initialize Pythia c Call PYINIT('CMS','p+','p+',sqrts) if(verbose) Call PYSTAT(2) write(*,*)'AAR::PYTHIA EVENT GENREATOR INITIALIZED' call ugpini !Define standard GEANT particles c c Store title of user event generator c chevt_name = 'Pythia/Isajet generator' ! c********************************************************** c c Main even loop c c 1. Call event generator Pythia c 2. Histogram particles c 3. Record events in output file c c loop over event number Do evtnum=1,tnumevt mxtot=0 write(6,*)'Event number ', evtnum Call PYEVNT If(NEVHEP.LT.10) Call LULIST(1) NEVHEP=NEVHEP+1 NPART=0 c write event header information c nptls=N !set common nptls c c Pythia returns non-trackable particles. Exclude these... c c count up # of particles nptls=0 do ilp=1,N If(K(Ilp,1).GE.1.and.K(Ilp,1).LE.10) nptls=nptls+1 enddo c write event header to output file. Header will contain c 1. event number c 2. number of particles in event c 3.Process ID c 4. Bjorken's x1, x2 (2 values) c 5. partonic s,t,u (3 values) c 6. Q squared c 7. transverse momentum if(zebrastore) then else write (15) evtnum,nptls write (6,*) 'Number of particles: ',nptls write(15)MSTI(1) !Process ID write(15)(PARI(Iwrite),Iwrite=33,34) !Bjorken's x1, x2 (2 values) write(15)(PARI(Iwrite),Iwrite=14,16)!partonic s,t,u (3 values) write(15)PARI(22) !Q squared write(15)PARI(17) !transverse momentum endif DO Ilp=5,8 c store info on intermediary particles (5-8) xyzvert(1,Ilp) = Vpyth(Ilp,1) xyzvert(2,Ilp) = Vpyth(Ilp,2) xyzvert(3,Ilp) = Vpyth(Ilp,3) p4vec(1,Ilp) = Ppyth(Ilp,1) p4vec(2,Ilp) = Ppyth(Ilp,2) p4vec(3,Ilp) = Ppyth(Ilp,3) p4vec(4,Ilp) = Ppyth(Ilp,4) write(15) K(ILP,2) write(15) p4vec(1,Ilp) ,p4vec(2,Ilp),p4vec(3,Ilp), + p4vec(4,Ilp) END DO Do Ilp=1,N c store ntuple info nstore(1)=real(evtnum) nstore(2)=real(Ilp) nstore(3)=real(K(ilp,1)) nstore(4)=real(K(ilp,2)) nstore(5)=real(K(ilp,3)) nstore(6)=Ppyth(ilp,1) nstore(7)=Ppyth(ilp,2) nstore(8)=Ppyth(ilp,3) nstore(9)=Ppyth(ilp,4) nstore(10)=Ppyth(ilp,5) call hfn(1000,nstore) If(K(Ilp,1).GE.1.and.K(Ilp,1).LE.10) Then mxtot = mxtot + 1 Call GETPAR(K(Ilp,2),id_geant) if (xverbose) Then write(6,*)'Particle #', ilp,' of ', N, 'ID: ',id_geant write(6,*)Vpyth(Ilp,1),Vpyth(Ilp,2),Vpyth(Ilp,3) write(6,*)Ppyth(Ilp,1),Ppyth(Ilp,2),Ppyth(Ilp,3) endif c Store vertex position xyzvert(1,Ilp) = Vpyth(Ilp,1) xyzvert(2,Ilp) = Vpyth(Ilp,2) xyzvert(3,Ilp) = Vpyth(Ilp,3) c Store momentum of particle p4vec(1,Ilp) = Ppyth(Ilp,1) p4vec(2,Ilp) = Ppyth(Ilp,2) p4vec(3,Ilp) = Ppyth(Ilp,3) p4vec(4,Ilp) = Ppyth(Ilp,4) idtot(ILP) = id_geant id_parent(ILP) = 0 NPART=NPART+1 if(zebrastore) then else write(15) idtot(ILP) write(15) p4vec(1,Ilp) ,p4vec(2,Ilp),p4vec(3,Ilp), + p4vec(4,Ilp) endif Else Endif 20 enddo enddo CALL PYSTAT(1) write(6,*)'##### PYTHIA EVENT GENERATOR END #####' c successfull end, go to exit goto 1001 c************************************************ c c Error handling and Utility closeout c 996 continue write(*,*) 'Error opening output file. Move old data file' goto 1001 997 continue write(6,998) 998 format(3x,'Unable to open pythia.par') stop ' Cannot find user input file' 999 continue write(6,1000) 1000 format(3x,'Read error in pyth_par segment of pythia.par',3x, + ' Namelist mis-match in pyth_par segment of pythia.par ?',3x, + 'The pythia.par file will be re-read to pinpoint the erroneous', + ' line',3x,'****This will cause the program to crash.****') close(unit=15) open(unit=15,file=par_file,status='OLD',err=997) read(15,nml=pyth_par) close(unit=15) stop ' PC23GEM: stop because of pythia.par file error.' 1001 continue if (zebrastore) then else close (15) call hrout(0,ICYCLE,' ') call hrend('LUN1') endif end ************************************ SUBROUTINE GETPAR(IPLUND,IPRG) ************************************ C C Translate LUND Particle Code To GEANT3 Code C (PYTHIA 5.4, JETSET 7.3) C COMMON/LPCOM/LPBARY,LPMESN INTEGER*4 LPBARY(6,6,6,2,2),LPMESN(6,6,4,2) C INTEGER IMES,JMES,IBARY,JBARY,KBARY,ISPIN,ISPINT INTEGER IANT,IPLHM,IPLTMP,IMESBR,ILEPGB C IPRG=0 C * check anti-particle (1:particle 2:anti-particle) IANT=1 IF (IPLUND.LT.0) IANT=2 * check higher multiplet IPLHM=ABS(IPLUND/10000) * particle ID without indication of higher multiplet IPLTMP=ABS(IPLUND-IPLHM*10000) * check leptons or gauge bosons (0:leptons or gauge bosons 1:others) ILEPGB=0 IF (IPLTMP/100.NE.0) ILEPGB=1 IF (ILEPGB.EQ.0) THEN IPRG=0 IF (IPLTMP.EQ.22) IPRG=1 IF (IPLTMP.EQ.11) THEN IF (IANT.EQ.1) IPRG=3 IF (IANT.EQ.2) IPRG=2 ENDIF IF (IPLTMP.EQ.13) THEN IF (IANT.EQ.1) IPRG=6 IF (IANT.EQ.2) IPRG=5 ENDIF IF (IPLTMP.EQ.15) THEN IF (IANT.EQ.1) IPRG=34 IF (IANT.EQ.2) IPRG=33 ENDIF IF ((IPLTMP.EQ.12).OR.(IPLTMP.EQ.14).OR.(IPLTMP.EQ.16)) + IPRG=4 IF (IPLTMP.EQ.24) THEN IF (IANT.EQ.1) IPRG=42 IF (IANT.EQ.2) IPRG=43 ENDIF IF (IPLTMP.EQ.23) IPRG=44 RETURN ENDIF * check diffractive states IF (IPLTMP.EQ.210) THEN IPRG=8 RETURN ELSEIF (IPLTMP.EQ.2110) THEN IPRG=13 RETURN ELSEIF (IPLTMP.EQ.2210) THEN IPRG=14 RETURN ENDIF * check meson or baryon (1:baryon 2:meson) IMESBR=1 IF (IPLTMP/1000.EQ.0) IMESBR=2 IF (IPLHM.NE.0) THEN IPRG=0 RETURN ENDIF IF (IMESBR.EQ.1) THEN IBARY=IPLTMP/1000 JBARY=(IPLTMP-IBARY*1000)/100 KBARY=(IPLTMP-IBARY*1000-JBARY*100)/10 ISPIN=(IPLTMP-IBARY*1000-JBARY*100-KBARY*10)/2 IPRG=LPBARY(IBARY,JBARY,KBARY,ISPIN,IANT) ELSEIF (IMESBR.EQ.2) THEN IMES=IPLTMP/100 JMES=(IPLTMP-IMES*100)/10 ISPIN=(IPLTMP-IMES*100-JMES*10) IF (ISPIN.EQ.0) THEN ISPINT=1 ELSE ISPINT=(ISPIN-1)/2+2 ENDIF IPRG=LPMESN(IMES,JMES,ISPINT,IANT) ENDIF C RETURN END *********************** SUBROUTINE ugpini *********************** COMMON/LPCOM/LPBARY,LPMESN INTEGER*4 LPBARY(6,6,6,2,2),LPMESN(6,6,4,2) LPBARY(2,1,1,1,1)=13 LPBARY(3,1,1,1,1)=21 LPBARY(2,2,1,1,1)=14 LPBARY(3,2,1,1,1)=20 LPBARY(3,3,1,1,1)=23 LPBARY(2,1,2,1,1)=18 LPBARY(3,1,2,1,1)=41 LPBARY(3,2,2,1,1)=19 LPBARY(3,3,2,1,1)=22 LPBARY(3,3,3,2,1)=24 LPBARY(2,1,1,1,2)=25 LPBARY(3,1,1,1,2)=29 LPBARY(2,2,1,1,2)=15 LPBARY(3,2,1,1,2)=28 LPBARY(3,3,1,1,2)=31 LPBARY(3,1,2,1,2)=26 LPBARY(3,2,2,1,2)=27 LPBARY(3,3,2,1,2)=30 LPBARY(3,3,3,2,2)=32 LPMESN(3,1,1,1)=16 LPMESN(1,3,1,1)=10 LPMESN(1,1,2,1)=7 LPMESN(2,1,2,1)=8 LPMESN(4,1,2,1)=35 LPMESN(2,2,2,1)=17 LPMESN(3,2,2,1)=11 LPMESN(4,2,2,1)=37 LPMESN(4,3,2,1)=39 LPMESN(2,1,2,2)=9 LPMESN(4,1,2,2)=36 LPMESN(3,2,2,2)=12 LPMESN(4,2,2,2)=38 LPMESN(4,3,2,2)=40 RETURN END