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 2/17/99 Corrected version from Yuji Goto c 2/24/99 Revised version for muon arm ---Yajun Mao c Notes: taken from routine by Naohito Saito c 3/25/02 Add control flags in input parameter file - Frédéric Fleuret c Add ntuple before kinematic cuts c c 01/09/03 Add momentum cuts on output particles - FF c 01/13/03 make it for pythia6205 c 08/30/03 Add the dump of the seed state to file - SG c c========================================= *=== All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) *=== Pythia event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) *=== Pythia Parameters. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) *=== global variables INTEGER I,ID_GEANT, IWRITE, ILP, ICYCLE LOGICAL EVT_IDS/.FALSE./,EVT_GEO/.FALSE./ *=== event and particle information INTEGER IDTOT, NPTLS REAL*4 P4VEC(4) COMMON/EVT/NPTLS *=== counters ! NUMTOT = total number of generated events ! NUMSIG = total number of signal events ! EVTNUM = total number of events after selection INTEGER NUMTOT/0/,NUMSIG/0/,EVTNUM/0/ *=== par file input variables DATA SQRTS/200./ LOGICAL SEED/.FALSE./ *=============================================================================== *=== Geometrical description REAL THMI_NO,THMA_NO,THMI_SU,THMA_SU LOGICAL NO_ON,SU_ON COMMON/COMGLO/THMI_NO,THMA_NO,NO_ON,THMI_SU,THMA_SU,SU_ON NAMELIST/NMLGLO/THMI_NO,THMA_NO,NO_ON,THMI_SU,THMA_SU,SU_ON DATA THMI_NO/9./,THMA_NO/37./,NO_ON/.TRUE./, + THMI_SU/9./,THMA_SU/37./,SU_ON/.TRUE./ *=== Signal description INTEGER ISIG1,IPROD1,ISIG2,IPROD2,NPROD1(2),NPROD2(2) REAL MI(2),P1MIMA(2),P2MIMA(2) COMMON/COMSIG/ISIG1,IPROD1,NPROD1,P1MIMA, + ISIG2,IPROD2,NPROD2,P2MIMA,MI NAMELIST/NMLSIG/ISIG1,IPROD1,NPROD1,P1MIMA, + ISIG2,IPROD2,NPROD2,P2MIMA,MI *=== Output parameters LOGICAL VERB, WRIFLAG, DO1, DO2 CHARACTER*20 H_FILE COMMON/COMWRI/VERB,WRIFLAG,H_FILE,DO1,DO2 NAMELIST/NMLWRI/VERB,WRIFLAG,H_FILE,DO1,DO2 DATA VERB/.FALSE./, WRIFLAG/.FALSE./, H_FILE/'pythia.hist'/, + DO1/.FALSE./,DO2/.FALSE./ *=============================================================================== NAMELIST/PYTH_PAR/ ! Pythia parameters + TNUMEVT,SQRTS,SEED *=== to read selection file CHARACTER*25 PARAM CHARACTER TEST *===== (before selection) ntuple definition *= This ntuple is intended to store J/psi and muon information *= coming from all generated events before any kinematic selection. CHARACTER*20 HISTF/'pythia.hist'/ INTEGER NSIZE1/19/ CHARACTER*8 CHFORM1(19) DATA CHFORM1 / + 'Event','ID','Px','Py','Pz','E','Vx','Vy','Vz', + 'pID','pPx','pPy','pPz','pE', + 'gpID','gpPx','gpPy','gpPz','gpE'/ REAL NSTORE1(19) *===== output (after selection) ntuple *= This ntuple is intended to store particle information *= coming from events which are stored in the output data file. *= This is a subsample of the previous ntuple. INTEGER NSIZE2/19/ CHARACTER*8 CHFORM2(19) DATA CHFORM2 / + 'Event','ID','Px','Py','Pz','E','Vx','Vy','Vz', + 'pID','pPx','pPy','pPz','pE', + 'gpID','gpPx','gpPy','gpPz','gpE'/ REAL NSTORE2(19) COMMON/NTVAR/NSTORE1,NSTORE2 *=== Externals LOGICAL PART_OK EXTERNAL PYDATA EXTERNAL PART_OK *=== Initialize cernstuff parameter (NWPAWC=100000000) REAL H COMMON /PAWC/H(NWPAWC) call hlimit(NWPAWC) *================================================== initialization * 1. Read input parameters * 2. Initalize output files and ntuples * 3. Initialize Pythia *================================================== initialization *====================================== Reading input name list OPEN(UNIT=15,FILE='pythia.par',STATUS='OLD',ERR=997) 10 READ(15,*) PARAM TEST=PARAM IF (TEST.NE.'!'.AND.PARAM.NE.'$pyth_par') CALL PYGIVE(PARAM) IF (PARAM.NE.'$pyth_par') GOTO 10 BACKSPACE(15) READ(15,NML=PYTH_PAR,ERR=999) READ(15,NML=NMLGLO,ERR=999) READ(15,NML=NMLSIG,ERR=999) READ(15,NML=NMLWRI,ERR=999) CLOSE(UNIT=15) *===================================== Initialize zebra writeout OPEN(UNIT=15, FILE='pythia.dat', STATUS='NEW', + ACCESS='sequential', FORM='unformatted') WRITE(15) SQRTS,TNUMEVT ! Write output file header *========================================== Declare ntuples CALL HROPEN(16,'LUN1',H_FILE,'N',1024,ISTAT) IF(DO1) CALL HBOOKN(100,'Before sel.',NSIZE1,'LUN1',8192,CHFORM1) IF(DO2) CALL HBOOKN(1000,'After sel.',NSIZE2,'LUN1',8192,CHFORM2) *========================================== Open seed file OPEN(UNIT=20, FILE='seed.state', STATUS='OLD', + ACCESS='sequential', FORM='unformatted') *========================================== Initialize Pythia Call PYINIT('CMS','p+','p+',SQRTS) IF(SEED) CALL PYRSET(20,0) *=== Pythia print out IF (VERB) THEN CALL PYSTAT(2) CALL PYSTAT(5) CALL PYSTAT(6) CALL PYLIST(12) ENDIF *================================ Define standard GEANT particles CALL UGPINI *============================================================ event loop * 1. Call event generator Pythia * 2. Histogram particles * 3. Record events in output file *============================================================ event loop 1 CONTINUE *=== Call event generator Pythia Call PYEVNT NUMTOT=NUMTOT+1 ! total number of generated events NPTLS = 0 ! number of stable particles to be written in the output CALL CHK_EVT(EVT_IDS,EVT_GEO) *==============Starting here : treat events passing EVT_IDS selection only IF (EVT_IDS) NUMSIG = NUMSIG+1 ! total number of signal events IF (DO1.AND.EVT_IDS) THEN DO ILP=1,N IF (K(ILP,1).GE.1.AND.K(ILP,1).LE.10) THEN IF (.NOT.WRIFLAG.OR. + (WRIFLAG.AND.PART_OK(ILP,'IDS'))) + CALL FILL_NT(100,NUMSIG,ILP) ! fill ntuple 100 (before selection) ENDIF ENDDO ENDIF *==============Starting here : treat events passing EVT_GEO selection only NWRIT = 0 ! number of indeed written particle IF (.NOT.EVT_IDS.OR..NOT.EVT_GEO) GOTO 1 EVTNUM=EVTNUM+1 ! number of events after selections *==============Event has been accepted, time to write the seed state CALL PYRGET(20,-1) *=== Print full event story (10 first events only) * IF (VERB) THEN IF (EVTNUM.LT.10) Call PYLIST(1) WRITE (6,*) 'Number of particles=',NPTLS,' in event ',EVTNUM * ENDIF *=== Write output file header for this event WRITE(15) EVTNUM,NPTLS ! event number, number of particles in event 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 *=== Write info on intermediary particles (5-8); needed by PISA (don't know why ???) INTERM1 = 5 INTERM2 = 8 IF (N.LT.8) THEN INTERM1 = 1 INTERM2 = 4 ENDIF DO ILP=1,N IF (ILP.GE.INTERM1.AND.ILP.LE.INTERM2) THEN P4VEC(1) = P(ILP,1) P4VEC(2) = P(ILP,2) P4VEC(3) = P(ILP,3) P4VEC(4) = P(ILP,4) WRITE(15) K(ILP,2) WRITE(15) P4VEC(1), P4VEC(2), P4VEC(3), P4VEC(4) ENDIF *=== Write info on stable particles and fill ntuple 1000 (after selections) IF (K(ILP,1).GE.1.AND.K(ILP,1).LE.10) THEN IF (.NOT.WRIFLAG.OR. + (WRIFLAG.AND.PART_OK(ILP,'IDS'))) THEN NWRIT = NWRIT + 1 CALL GETPAR(K(ILP,2),ID_GEANT) ! change particle ID code for PISA P4VEC(1) = P(ILP,1) P4VEC(2) = P(ILP,2) P4VEC(3) = P(ILP,3) P4VEC(4) = P(ILP,4) IDTOT = ID_GEANT WRITE(15) IDTOT WRITE(15) P4VEC(1) ,P4VEC(2),P4VEC(3),P4VEC(4) IF (DO2) CALL FILL_NT(1000,EVTNUM,ILP) ! fill ntuple 1000 (after selection) ENDIF ENDIF 100 CONTINUE ENDDO IF(NWRIT.NE.NPTLS) + WRITE(*,*)'WARNING : NUMBER OF WRITTEN PARTICLE = ',NWRIT, + ' AND NUMBER OF COUNTED PARTICLE = ',NPTLS,' ARE DIFFERENT' IF(EVTNUM.LT.TNUMEVT) GOTO 1 ! number of event not reached yet, new event generation CLOSE(UNIT=20) CLOSE(UNIT=15) *=== Write ending output OPEN(UNIT=15,FILE='pythia.par',STATUS='OLD',ERR=997) WRITE(*,*) 'PYTHIA input parameters are :' 11 READ(15,*) PARAM TEST=PARAM IF (TEST.NE.'!'.AND.PARAM.NE.'$pyth_par') WRITE(*,*) PARAM IF (PARAM.NE.'$pyth_par') GOTO 11 CLOSE(UNIT=15) WRITE(*,*) WRITE(*,*) 'Selection cuts are :' WRITE(*,*) ' tnumevt = ',tnumevt WRITE(*,*) ' sqrts = ' ,sqrts WRITE(*,*) ' seed = ' ,seed WRITE(*,*) ' thmi_no = ',thmi_no WRITE(*,*) ' thma_no = ',thma_no WRITE(*,*) ' no_on = ',no_on WRITE(*,*) ' thmi_su = ',thmi_su WRITE(*,*) ' thma_su = ',thma_su WRITE(*,*) ' su_on = ',su_on WRITE(*,*) ' isig1 = ',isig1 WRITE(*,*) ' iprod1 = ',iprod1 WRITE(*,*) ' p1mima(1) = ',p1mima(1), + ' p1mima(2) = ',p1mima(2) WRITE(*,*) ' isig2 = ',isig2 WRITE(*,*) ' iprod2 = ',iprod2 WRITE(*,*) ' p2mima(1) = ',p2mima(1), + ' p2mima(2) = ',p2mima(2) WRITE(*,*) ' mi(1) = ',mi(1),' mi(2) = ',mi(2) WRITE(*,*) ' verb = ',verb WRITE(*,*) ' wriflag = ',wriflag WRITE(*,*) ' h_file = ',h_file WRITE(*,*) ' do1 = ',do1 WRITE(*,*) ' do2 = ',do2 WRITE(*,*) '====================================' WRITE(*,*)'Total number of generated events = ',NUMTOT WRITE(*,*)'Total number of signal events = ',NUMSIG WRITE(*,*)'Total number of stored events = ',EVTNUM ACCEPTA = REAL(EVTNUM)/REAL(NUMSIG) WRITE(*,*)'Acceptance after selections = ',EVTNUM,' / ', + NUMSIG,' = ',ACCEPTA *=== successfull end, go to exit goto 1001 *======================================= Error handling and Utility closeout 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='pythia.par',STATUS='OLD',ERR=997) READ(15,NML=PYTH_PAR) CLOSE(UNIT=15) STOP ' PC23GEM: stop because of pythia.par file error.' 1001 CONTINUE *=== close output file CLOSE(15) *=== Fill ntuple files IF (DO1.OR.DO2) THEN CALL HROUT(0,ICYCLE,' ') CALL HREND('LUN1') ENDIF CLOSE(16) *=== LUND Statistics CALL PYSTAT(1) END *============================================================== ROUTINES *= FUNCTION PART_OK(I,TYPE) *============================================================== ROUTINES *= FUNCTION SIG_OK(I,ISIGNAL) *============================================================== ROUTINES *= SUBROUTINE CHK_EVT(EVT_IDS,EVT_GEO) *============================================================== ROUTINES *= SUBROUTINE FILL_NT : *= Fill both ntuples *============================================================== ROUTINES *= SUBROUTINE GETPAR : *= Translate Lund Particle Code To GEANT3 Code *= (PYTHIA 5.4, JETSET) *============================================================== ROUTINES *= SUBROUTINE UGPINI : *= Define standard GEANT particles *============================================================== ROUTINES *===================== PART_OK =========================* *=== Test particle acceptance and ID *=======================================================* *= I = particle line number in pythia *= ID = particle ID (mu- = 13, mu+ = -13, ...) *= IDP = particle's parent ID *= IDGP= particle's grand parent ID (needed for DrellYan process) *=======================================================* LOGICAL FUNCTION PART_OK(I,TYPE) *=== All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) CHARACTER*3 TYPE LOGICAL VERB, WRIFLAG, DO1, DO2 CHARACTER*20 H_FILE COMMON/COMWRI/VERB,WRIFLAG,H_FILE,DO1,DO2 *=== Signal description INTEGER ISIG1,IPROD1,ISIG2,IPROD2,NPROD1(2),NPROD2(2) REAL MI(2), P1MIMA(2), P2MIMA(2) COMMON/COMSIG/ISIG1,IPROD1,NPROD1,P1MIMA, + ISIG2,IPROD2,NPROD2,P2MIMA,MI *=== Pythia event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) REAL THETAP, PTOTAL *=== Geometrical description REAL THMI_NO,THMA_NO,THMI_SU,THMA_SU LOGICAL NO_ON,SU_ON COMMON/COMGLO/THMI_NO,THMA_NO,NO_ON,THMI_SU,THMA_SU,SU_ON LOGICAL CONT, SIG_OK EXTERNAL SIG_OK PART_OK = .FALSE. ID = K(I,2) IDP = K(K(I,3),2) IDGP= K(K(K(I,3),3),2) PTOTAL=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) THETAP=ACOS(P(I,3)/PTOTAL)*57.29578 *==============================Check particle ID IF (TYPE.EQ.'IDS') THEN *=== Must be a stable particle IF (K(I,1).GE.1.AND.K(I,1).LE.10) THEN *= If no selection required, all particle pass IF (ISIG1.EQ.0.AND.IPROD1.EQ.0.AND. + ISIG2.EQ.0.AND.IPROD2.EQ.0) PART_OK=.TRUE. *= Test first particle selection IF (IPROD1.NE.0.AND.ID.EQ.IPROD1) THEN IF (ISIG1.EQ.0) PART_OK = .TRUE. IF (ISIG1.NE.0.AND.SIG_OK(I,ISIG1)) PART_OK=.TRUE. ENDIF *= Test second particle selection IF (IPROD2.NE.0.AND.ID.EQ.IPROD2) THEN IF (ISIG2.EQ.0) PART_OK = .TRUE. IF (ISIG2.NE.0.AND.SIG_OK(I,ISIG2)) PART_OK=.TRUE. ENDIF ENDIF ENDIF *====================== Check particle acceptance IF (TYPE.EQ.'GEO') THEN *=== Particle momentum (only if IPROD1/2 is defined) IF (ID.EQ.IPROD1.OR.ID.EQ.IPROD2) THEN IF ( + ((ID.EQ.IPROD1.AND. + PTOTAL.GE.P1MIMA(1).AND.PTOTAL.LE.P1MIMA(2)).OR. + (ID.EQ.IPROD2.AND. + PTOTAL.GE.P2MIMA(1).AND.PTOTAL.LE.P2MIMA(2))).AND. *=== North selection + ((NO_ON.AND. + THETAP.GE.THMI_NO.AND.THETAP.LE.THMA_NO).OR. *=== South selection + (SU_ON.AND. + THETAP.GE.THMI_SU.AND.THETAP.LE.THMA_SU)) + ) PART_OK=.TRUE. ELSE IF ( *=== North selection + ((NO_ON.AND. + THETAP.GE.THMI_NO.AND.THETAP.LE.THMA_NO).OR. *=== South selection + (SU_ON.AND. + THETAP.GE.THMI_SU.AND.THETAP.LE.THMA_SU)) + ) PART_OK=.TRUE. ENDIF ENDIF * RETURN END *======================= SIG_OK =======================* *= Check particle Signal *======================= SIG_OK =======================* LOGICAL FUNCTION SIG_OK(I,ISIGNAL) *=== All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) *=== Pythia event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) INTEGER ID,I_P,IDP ID = K(I,2) I_P = K(I,3) IDP = K(I_P,2) SIG_OK = .FALSE. *=== for all signals (including J/psi) IF (IDP.EQ.ISIGNAL) SIG_OK = .TRUE. *=== for Drell-Yan; look at grand parents IF (ISIGNAL.EQ.23) THEN I_P = K(I_P,3) IDP = K(I_P,2) IF (IDP.EQ.ISIGNAL) SIG_OK = .TRUE. ENDIF *=== for DDbar (resp. BBbar); look for a c (resp. b) or a cbar (resp. bbar) quark IF (ABS(ISIGNAL).LE.5) THEN 31 CONTINUE I_P = K(I_P,3) IDP = K(I_P,2) IF (I_P.GT.2) THEN ! if parents are incoming protons, stop IF (IDP.EQ.ISIGNAL) SIG_OK = .TRUE. IF (IDP.NE.ISIGNAL) GOTO 31 ENDIF ENDIF RETURN END *======================= CHK_EVT =======================* *= Make event selection *= EVT_IDS = .TRUE. if event pass IDs selections *= EVT_GEO = .TRUE. if event pass IDs and geom selections *======================= CHK_EVT =======================* SUBROUTINE CHK_EVT(EVT_IDS,EVT_GEO) *=== All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) *=== Pythia event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) REAL P1(4), P2(4), Minv LOGICAL VERB, WRIFLAG, DO1, DO2 CHARACTER*20 H_FILE COMMON/COMWRI/VERB,WRIFLAG,H_FILE,DO1,DO2 *=== Signal description INTEGER ISIG1,IPROD1,ISIG2,IPROD2,NPROD1(2),NPROD2(2) REAL MI(2), P1MIMA(2), P2MIMA(2) COMMON/COMSIG/ISIG1,IPROD1,NPROD1,P1MIMA, + ISIG2,IPROD2,NPROD2,P2MIMA,MI COMMON/EVT/NPTLS INTEGER NTOT(2) LOGICAL EVT_IDS, EVT_GEO LOGICAL PART_OK, IDS_OK, GEO_OK EXTERNAL PART_OK,Minv INTEGER N_ID, N_ID1, N_ID2 INTEGER N_GEO, N_GEO0, N_GEO1, N_GEO2 INTEGER PARTIDX(100), IDX, IDX1, IDX2, KIDX1, KIDX2 * IDX=0 DO IPART=1,100 PARTIDX(IPART)=0 ENDDO EVT_IDS = .FALSE. EVT_GEO = .FALSE. N_ID = 0 N_ID1 = 0 N_ID2 = 0 N_GEO = 0 N_GEO0 = 0 N_GEO1 = 0 N_GEO2 = 0 DO ILP=1,N ID = K(ILP,2) *=== Check particle IDs IF (PART_OK(ILP,'IDS')) THEN IF (WRIFLAG) N_GEO0 = N_GEO0 + 1 *=== Fill signal particles for resonance selection IF (ISIG1.EQ.ISIG2.AND.ISIG1.NE.0) THEN IDX=IDX+1 PARTIDX(IDX)=ILP ENDIF *=== Check particle geometrical acceptance IF (IPROD1.EQ.0.AND.IPROD2.EQ.0) THEN N_ID = N_ID+1 IF (PART_OK(ILP,'GEO')) N_GEO = N_GEO+1 ELSE IF (IPROD1.EQ.0.OR.ID.EQ.IPROD1) THEN DO I=1,4 P1(I) = P(ILP,I) ENDDO N_ID1 = N_ID1+1 IF (PART_OK(ILP,'GEO')) N_GEO1 = N_GEO1+1 ENDIF IF (IPROD2.EQ.0.OR.ID.EQ.IPROD2) THEN DO I=1,4 P2(I) = P(ILP,I) ENDDO N_ID2 = N_ID2+1 IF (PART_OK(ILP,'GEO')) N_GEO2 = N_GEO2+1 ENDIF ENDIF ENDIF IF (.NOT.WRIFLAG.AND.K(ILP,1).GE.1.AND.K(ILP,1).LE.10) ! number of part to be written in output + N_GEO0 = N_GEO0 + 1 ENDDO NTOT(1) = NPROD1(1)+NPROD2(1) NTOT(2) = NPROD1(2)+NPROD2(2) *=== IDs selection only IF ((N_ID.GE.NTOT(1).AND.N_ID.LE.NTOT(2)).OR. + (N_ID1.GE.NPROD1(1).AND.N_ID1.LE.NPROD1(2).AND. + N_ID2.GE.NPROD2(1).AND.N_ID2.LE.NPROD2(2).AND. + Minv(P1,P2).GE.MI(1).AND.Minv(P1,P2).LE.MI(2))) + EVT_IDS=.TRUE. ! IDs selection only *=== IDs + acceptance selection IF ((N_GEO.GE.NTOT(1).AND.N_GEO.LE.NTOT(2)).OR. + (N_GEO1.GE.NPROD1(1).AND.N_GEO1.LE.NPROD1(2).AND. + N_GEO2.GE.NPROD2(1).AND.N_GEO2.LE.NPROD2(2))) EVT_GEO=.TRUE. ! IDs + acceptance selection NPTLS = N_GEO0 ! set number of particle to be written in the output *=== Specific selection for resonances, check that both signal particles are ok IF (ISIG1.EQ.ISIG2.AND.ISIG1.NE.0.AND.IPROD1.EQ.-IPROD2) THEN EVT_GEO=.FALSE. IF (IDX.GE.2) THEN IDX1 = IDX-1 DO ILP1=1,IDX1 IDX2=ILP1+1 DO ILP2=IDX2,IDX KIDX1=K(PARTIDX(ILP1),3) KIDX2=K(PARTIDX(ILP2),3) *=== In case of Drell-yan, look at grand-parents IF (ISIG1.EQ.23) THEN KIDX1=K(K(PARTIDX(ILP1),3),3) KIDX2=K(K(PARTIDX(ILP2),3),3) ENDIF IF (PART_OK(PARTIDX(ILP1),'IDS').AND. + PART_OK(PARTIDX(ILP2),'IDS').AND. + PART_OK(PARTIDX(ILP1),'GEO').AND. + PART_OK(PARTIDX(ILP2),'GEO').AND. + (KIDX1.EQ.KIDX2)) THEN EVT_GEO=.TRUE. ENDIF ENDDO ENDDO ENDIF ENDIF *=== RETURN END *========================= Minv ===========================* *= Compute invariante mass *========================= Minv ===========================* REAL FUNCTION Minv(P1,P2) REAL P1(4), P2(4), P(4) Minv = 0 DO I=1,4 P(I) = P1(I) + P2(I) ENDDO Minv = SQRT(P(4)**2-P(1)**2-P(2)**2-P(3)**2) RETURN END *====================== FILL_NT ===========================* *= Fill ntuple variables *= if parameter = 100, fill event before fiducial selection *= if parameter = 1000, fill event after fiducial selection *====================== FILL_NT ===========================* SUBROUTINE FILL_NT(INTUPLE,NUMTOT,ILP) *=== All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) REAL NSTORE1(19), NSTORE2(19) COMMON/NTVAR/NSTORE1,NSTORE2 *=== Pythia event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) INTEGER INTUPLE,NUMTOT,ILP IF (INTUPLE.EQ.100) THEN * Store ntuple info NSTORE1(1) = REAL(NUMTOT) ! event number NSTORE1(2) = REAL(K(ILP,2)) ! Particle ID NSTORE1(3) = P(ILP,1) ! Particle Px NSTORE1(4) = P(ILP,2) ! Particle Py NSTORE1(5) = P(ILP,3) ! Particle Pz NSTORE1(6) = P(ILP,4) ! Particle E NSTORE1(7) = V(ILP,1) ! Particle Vx NSTORE1(8) = V(ILP,2) ! Particle Vy NSTORE1(9) = V(ILP,3) ! Particle Vz NSTORE1(10) = REAL(K(K(ILP,3),2)) ! Parent ID NSTORE1(11) = P(K(ILP,3),1) ! Parent Px NSTORE1(12) = P(K(ILP,3),2) ! Parent Py NSTORE1(13) = P(K(ILP,3),3) ! Parent Pz NSTORE1(14) = P(K(ILP,3),4) ! Parent E NSTORE1(15) = REAL(K(K(K(ILP,3),3),2)) ! Grand Parent ID NSTORE1(16) = P(K(K(ILP,3),3),1) ! Grand Parent Px NSTORE1(17) = P(K(K(ILP,3),3),2) ! Grand Parent Py NSTORE1(18) = P(K(K(ILP,3),3),3) ! Grand Parent Pz NSTORE1(19) = P(K(K(ILP,3),3),4) ! Grand Parent E CALL HFN(100,NSTORE1) ELSEIF (INTUPLE.EQ.1000) THEN NSTORE2(1) = REAL(NUMTOT) ! event number NSTORE2(2) = REAL(K(ILP,2)) ! Particle ID NSTORE2(3) = P(ILP,1) ! Particle Px NSTORE2(4) = P(ILP,2) ! Particle Py NSTORE2(5) = P(ILP,3) ! Particle Pz NSTORE2(6) = P(ILP,4) ! Particle E NSTORE2(7) = V(ILP,1) ! Particle Vx NSTORE2(8) = V(ILP,2) ! Particle Vy NSTORE2(9) = V(ILP,3) ! Particle Vz NSTORE2(10) = REAL(K(K(ILP,3),2)) ! Parent ID NSTORE2(11) = P(K(ILP,3),1) ! Parent Px NSTORE2(12) = P(K(ILP,3),2) ! Parent Py NSTORE2(13) = P(K(ILP,3),3) ! Parent Pz NSTORE2(14) = P(K(ILP,3),4) ! Parent E NSTORE2(15) = REAL(K(K(K(ILP,3),3),2)) ! Grand Parent ID NSTORE2(16) = P(K(K(ILP,3),3),1) ! Grand Parent Px NSTORE2(17) = P(K(K(ILP,3),3),2) ! Grand Parent Py NSTORE2(18) = P(K(K(ILP,3),3),3) ! Grand Parent Pz NSTORE2(19) = P(K(K(ILP,3),3),4) ! Grand Parent E call HFN(1000,NSTORE2) ENDIF END *====================== GETPAR ========================* *= Translate LUND Particle Code To GEANT3 Code *= (PYTHIA 5.4, JETSET 7.3) *====================== GETPAR ========================* SUBROUTINE GETPAR(IPLUND,IPRG) COMMON/LPCOM/LPBARY,LPMESN INTEGER*4 LPBARY(6,6,6,2,2),LPMESN(6,6,4,2) INTEGER IMES,JMES,IBARY,JBARY,KBARY,ISPIN,ISPINT INTEGER IANT,IPLHM,IPLTMP,IMESBR,ILEPGB IPRG=0 * 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 RETURN END *====================== UGPINI ========================* *= Define standard GEANT particles *====================== UGPINI ========================* 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