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 c 2/24/99 Revised version for muon arm ---Yajun Mao c 12/6/00 slightly revised -- DMLee c c Notes: taken from routine by Naohito Saito c c********************************************************** c jmm 9/26/01 put only muon momentum and highest energy hadron momentum into c pythia.dat. pythia.hist contains muon and highest energy hadron momenta as in c 9/24/01 version. c 12/18/2001 jmm. modification of gen_pythia_muon.f_charmtomuhad to do bbar->J/psi implicit none include 'event.inc' 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' integer egzinit integer write_code_bank real*8 bimpact /100./ !minimun bias impact par. 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, isignal real py_varstore(8) integer iwrite, ntries, interm1, interm2 real zmu, ptp, thp, pdot, px1st, py1st, pz1st, p41st c input values integer imsel, imsub(2,5) !msel, mstp's to change integer imstp(2,5),imstj(2,5) real ickin(2,4),iparp(2,5) !ckin's to change integer imrlu(2,2) !MRLU's to change real ntstore(13) !ntuple storage real sqrts /200./ logical verbose /.true./ logical xverbose /.false./ logical zebrastore /.false./ ! logical for new zebra store function logical muon_found /.false./ logical all_keep /.false./ character*50 hist_file /'pythia.hist'/ integer nsize /13/ character*8 chform (13) data chform /'Event', 'Pid1', 'Px1','Py1','Pz1','Energy1', + 'vertx','verty','vertz', + 'x1','x2','xf','tau'/ real nstore(13) character*60 file_out ! output file name, added 30 apr 04 HvH integer iseq ! sequence number in the file name c c Par file section c CHARACTER*50 PAR_FILE NAMELIST /PYTH_PAR/ tnumevt, imsel, iMSUB, imstp, ickin, + imrlu,imstj,iparp, sqrts,verbose, xverbose, + zebrastore, hist_file EXTERNAL LUDATA, PYDATA, IACCEPT integer NWPAWC, H ! initialize cern stuff parameter (NWPAWC=80000) COMMON /PAWC/H(NWPAWC) call hlimit(NWPAWC) ! end declarations *---------------------------------------------------------------------------------------------* 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,5 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 c c tell me if the process is signal or background c isignal = 0 if(iMSUB(1,1).eq.1.and.iMSUB(2,1).eq.1)isignal=1 ! Drell-Yan if(iMSUB(1,1).eq.86)isignal=1 ! J/Psi if(msel.eq.4)isignal=1 ! Charm if(msel.eq.5)isignal=1 ! Bottom write(*,*)'Filling MSTP array.' do i=1,5 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,2 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 write(*,*)'Filling MSTJ array.' do i=1,5 if (iMSTJ(1,i).ne.0) then MSTJ(iMSTJ(1,i))=iMSTJ(2,i) if (verbose) then write(*,*) 'MSTJ(', iMSTJ(1,i),')',iMSTJ(2,i) endif endif enddo write(*,*) 'Filling PARP array ' do i=1,5 if (iparp(1,i).ne.0)then PARP(iparp(1,i))=iparp(2,i) if (verbose) then write(*,*) 'PARP(', iparp(1,i),')',iparp(2,i) endif endif enddo ********************************************************* c Initialization Section c Initialize ZEBRA bank c call MZEBRA(0) c c Initialize zebra writeout c file = pythia.dat is the file for the subsequent pisa run c if (.not.zebrastore) then open(unit=15, status='UNKNOWN', file='sequence.dat') iseq = 0 ! Pick up the next sequence read (15,'(i4)',err=40) iseq ! number for the output file. 40 iseq = iseq+1 ! Where to write the output file? rewind(unit=15) ! Run script disk_space.pl to see write (15,'(i4)') iseq ! which of the /phenix/data??/ disks close (unit=15) ! have free space. Is this legal? write (file_out,50) iseq 50 format('/phenix/data07/hubert/pythia/pythia_',i3.3,'.dat') write (6,'(''Output file name:'',a60)') file_out open(unit=15,file=file_out, + status='NEW', access='sequential', form='unformatted') write (hist_file,60) iseq ! same for the ntuple file 60 format('pythia_',i3.3,'.hist') ! call hropen(16,'LUN1',hist_file,'N',1024,istat) endif ! end output files opened ! Now declare the ntuple: call hbookn(1000,'Pythia Particles',nsize,'LUN1',8192,chform) write (15) sqrts,tnumevt ! write file header to pythia_xx.dat file call PYINIT('CMS','p+','p+',sqrts) ! Initialize Pythia if(verbose) Call PYSTAT(2) write(*,*)'AAR::PYTHIA EVENT GENREATOR INITIALIZED' call ugpini !Define standard GEANT particles chevt_name = 'Pythia/Isajet generator' ! store name of generator (in common?) c*************************************************************************************** c c Main even loop c c 1. Call event generator Pythia c 2. Histogram particles c 3. Record events in output file ntries = 0 ! jmm 10/9/01 count number of tries in order to evaluate !acceptance all_keep = .false. ! true : old version with all tracks written out ! false: new version (30 apr 2004), only one muon kept do evtnum=1,tnumevt ! master event loop muon_found = .false. do while (.not.muon_found) ntries = ntries + 1 ! number of tries to get a muon in the event call PYEVNT ! produce a new Pythia event ! jmm 9/18/01 check to see if there is a muon in either arm. ! use angular range of so. arm, 12-34.5 deg. do Ilp=1,N ! loop over particle list ! look for a muon. k(xx,2)=+-13 means mu+- pg.49 if(abs(k(ilp,2)).eq.13) then ptp = sqrt(Ppyth(ilp,1)**2+Ppyth(ilp,2)**2) ! muon pt thp = atan(ptp/abs(Ppyth(ilp,3))) ! muon angle if (thp.gt.0.2126.and.thp.lt.0.6873.and. ! acceptance cut + ppyth(ilp,4).gt.2.6) then ! p>2.6 Gev ?? muon_found = .true. ! good muon zmu = vpyth(ilp,3) ! save vertex endif ! good angle, good momentum endif ! this is a muon enddo ! end search for a muon in the list enddo ! muon found in the list NEVHEP=NEVHEP+1 ! NEVHEP pg.30 npart=0 if(evtnum.eq.1) call LULIST(1) ! ??? not in Pythia manual c c Pythia returns non-trackable particles. Exclude these... c ! count up # of particles nptls=0 ! nptls in pythia common (?) do ilp=1,N If(K(Ilp,1).GE.1.and.K(Ilp,1).LE.10) nptls=nptls+1 enddo if (mod(evtnum,100).eq.0) & write (6,70) evtnum, tnumevt, ntries, nptls 70 format(' event ',I4,'/',I4,' for ',I4, & ' events tried; nptls in this event: ',I4) 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(.not.all_keep) nptls = 2 if(.not.zebrastore) then ! Write out global event properties write (15) evtnum,nptls ! event number, multiplicity 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 interm1 = 5 ! Write to output file interm2 = 8 ! values for intermediate particles if ( N. lt. 8 )then ! gluon,gluon,c and cbar, indices 5-8. interm1 = 1 ! or 1-4 ?? interm2 = 4 endif do i=7,11 ! zero out 2nd muon ntuple for each event nstore(i)=0. ! enddo ! zmu = 1000. ! set zvertex to some unreasonable value do Ilp=1,N ! loop over all particles if(Ilp.ge.interm1.and.ILP.le.interm2) then xyzvert(1,Ilp) = Vpyth(Ilp,1) ! store info on intermediary particles (5-8) xyzvert(2,Ilp) = Vpyth(Ilp,2) ! or 1-4?? 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) endif ! end store gluons, c-quarks etc. If(K(Ilp,1).GE.1.and.K(Ilp,1).LE.10) then ! all particles c this is the 2nd muon c if(vpyth(ilp,3).eq.zmu.and.abs(k(ilp,2)).eq.13) then c nstore(7)=real(K(ilp,2)) c nstore(8)=Ppyth(ilp,1) c nstore(9)=Ppyth(ilp,2) c nstore(10)=Ppyth(ilp,3) c nstore(11)=Ppyth(ilp,4) c pdot = Ppyth(ilp,1)*px1st+Ppyth(ilp,2)*py1st+ c + Ppyth(ilp,3)*pz1st c nstore(19)=sqrt(2.*(Ppyth(ilp,4)*p41st-pdot)+0.02233) c print *, ilp,zmu,nstore(7),nstore(11) c endif c ! 1st muon (=13), z=1000 if(abs(k(ilp,2)).eq.13.and.zmu.eq.1000.) then zmu = vpyth(ilp,3) ! this is the first muon, nstore(1)=real(evtnum) ! store ntuple info. nstore(2)=real(K(ilp,2)) nstore(3)=Ppyth(ilp,1) nstore(4)=Ppyth(ilp,2) nstore(5)=Ppyth(ilp,3) nstore(6)=Ppyth(ilp,4) nstore(7)=vpyth(ilp,1) nstore(8)=vpyth(ilp,2) nstore(9)=vpyth(ilp,3) nstore(10)=PARI(33) ! x1 nstore(11)=PARI(34) ! x2 nstore(12)=PARI(35) ! xf=x1-x2 nstore(13)=PARI(36) ! tau=x1*x2 px1st = Ppyth(ilp,1) ! save 1st muon momenta for calculating mass py1st = Ppyth(ilp,2) pz1st = Ppyth(ilp,3) p41st = Ppyth(ilp,4) endif mxtot = mxtot + 1 ! mxtot must be in common c jmm 9/26/01, added to simplify pythia.par to correspond to ntuple c storage -- only muon and highest energy hadron. c if( all_keep ) then write (6,*)'>>>>>>> ERROR: all_keep turned up .true.' 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(.not.zebrastore) then write(15) idtot(ILP) write(15) p4vec(1,Ilp) ,p4vec(2,Ilp),p4vec(3,Ilp), + p4vec(4,Ilp) endif endif ! jmm, end of all_keep = true loop Endif ! end if all particles 20 enddo ! loop over all particles call hfn(1000,nstore) c print *, 'nstore 1, 6, 14', nstore(1),nstore(6),nstore(14) if(.not.all_keep) then p4vec(1,1) = nstore(3) p4vec(2,1) = nstore(4) p4vec(3,1) = nstore(5) p4vec(4,1) = nstore(6) idtot(1) = 5 if(nstore(2).lt.0) idtot(1) = 6 write(15) idtot(1) write(15) p4vec(1,1) ,p4vec(2,1),p4vec(3,1), + p4vec(4,1) p4vec(1,2) = nstore(8) p4vec(2,2) = nstore(9) p4vec(3,2) = nstore(10) p4vec(4,2) = nstore(11) idtot(2) = 5 if(nstore(7).lt.0) idtot(2) = 6 write(15) idtot(2) write(15) p4vec(1,2) ,p4vec(2,2),p4vec(3,2), + p4vec(4,2) endif ! keep-only-1 case enddo ! end master event loop write(6,*)'##### PYTHIA EVENT GENERATOR END #####' goto 1001 ! successfull end, go to exit 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 (.not.zebrastore) then close (15) call hrout(0,ICYCLE,' ') call hrend('LUN1') endif *----- LUND Statistics Call PYSTAT(1) print *, 'number of events requested =', tnumevt print *, 'number of events attempted =', ntries c 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 END *=======================================================================================* FUNCTION IACCEPT(numevt,isignal) Common/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) INTEGER*2 PI_SOUTH_P,PI_NORTH_P,PI_SOUTH_N,PI_NORTH_N PMCUT=2.0 ! set threshold 2.0GeV/c for muon PPCUT=2.5 ! set threshold 2.5GeV/c for pion IACCEPT=0 ! set the event unaccepted primarily MU_SOUTH_P=0 MU_NORTH_P=0 PI_SOUTH_P=0 PI_NORTH_P=0 MU_SOUTH_N=0 MU_NORTH_N=0 PI_SOUTH_N=0 PI_NORTH_N=0 DO I=1,N PTOTAL=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) THETA=ACOS(P(I,3)/PTOTAL)*57.29578 ID=K(I,2) c c ID = +-13, +-mu; ID=+-211, +-pi; ID=+-321,+-K c IF(THETA.GE.9.0.AND.THETA.LE.37.0)THEN ! particle in north arm IF(ID.EQ.-13.AND.PTOTAL.GE.PMCUT)MU_NORTH_P=MU_NORTH_P+1 IF(ID.EQ.13.AND.PTOTAL.GE.PMCUT)MU_NORTH_N=MU_NORTH_N+1 IF(ID.EQ.211.AND.PTOTAL.GE.PPCUT)PI_NORTH_P=PI_NORTH_P+1 IF(ID.EQ.321.AND.PTOTAL.GE.PPCUT)PI_NORTH_P=PI_NORTH_P+1 IF(ID.EQ.-211.AND.PTOTAL.GE.PPCUT)PI_NORTH_N=PI_NORTH_N+1 IF(ID.EQ.-321.AND.PTOTAL.GE.PPCUT)PI_NORTH_N=PI_NORTH_N+1 ELSEIF(THETA.GE.143.0.AND.THETA.LE.171.0)THEN ! particle in south arm IF(ID.EQ.-13.AND.PTOTAL.GE.PMCUT)MU_SOUTH_P=MU_SOUTH_P+1 IF(ID.EQ.13.AND.PTOTAL.GE.PMCUT)MU_SOUTH_N=MU_SOUTH_N+1 IF(ID.EQ.211.AND.PTOTAL.GE.PPCUT)PI_SOUTH_P=PI_SOUTH_P+1 IF(ID.EQ.321.AND.PTOTAL.GE.PPCUT)PI_SOUTH_P=PI_SOUTH_P+1 IF(ID.EQ.-211.AND.PTOTAL.GE.PPCUT)PI_SOUTH_N=PI_SOUTH_N+1 IF(ID.EQ.-321.AND.PTOTAL.GE.PPCUT)PI_SOUTH_N=PI_SOUTH_N+1 ENDIF ENDDO if(numevt.le.20)then C print*,'mu_n, mu_s, pi_n, pi_s, mu+n, mu+s, pi+n, pi+s: ', C + mu_south_n, mu_north_n, pi_south_n, pi_north_n, C + mu_south_p, mu_north_p, pi_south_p, pi_north_p endif if(isignal.eq.1)then C-----Check accepted muons IF(MU_SOUTH_N.GE.1.AND.MU_SOUTH_P.GE.1)GOTO 1 IF(MU_NORTH_N.GE.1.AND.MU_NORTH_P.GE.1)GOTO 1 else C-----Check accepted pions IF(PI_SOUTH_N.GE.1.OR.PI_SOUTH_P.GE.1)GOTO 1 IF(PI_NORTH_N.GE.1.OR.PI_NORTH_P.GE.1)GOTO 1 endif RETURN 1 IACCEPT=1 RETURN END *=====================================================================================*