[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / Pythia.F

C--------------------------------------------------------------------------
C
C Environment:
C      This software is part of the EvtGen package developed jointly
C      for the BaBar and CLEO collaborations.  If you use all or part
C      of it, please give an appropriate acknowledgement.
C
C Copyright Information: See BelEvtGen/COPYRIGHT
C      Copyright (C) 1998      Caltech, UCSB
C
C Module: EvtJetSetInit.F
C
C Description:
C
C Modification history:
C
C    DJL/RYD     August 11, 1998         Module created
C    RS          October 28, 2002        copied and modified from Jetset
C
C------------------------------------------------------------------------
      SUBROUTINE EVTPYTHIAINIT(FNAME)
      IMPLICIT NONE
      EXTERNAL PYDATA
      
      CHARACTER*(*) FNAME
      
      PRINT *,'PYUPDA : : ',FNAME
      OPEN(54,STATUS='OLD',FILE=FNAME)
      CALL PYUPDA(3,54)
      CLOSE(54)      
      WRITE (*,*) 'PYUPDA DONE'
      RETURN
      END
      SUBROUTINE INITPYTHIA(GAGA)
      IMPLICIT NONE
      INTEGER PYCOMP
      EXTERNAL PYDATA,PYCOMP
      LOGICAL FIRST
      SAVE FIRST
      DATA FIRST / .TRUE. /
      DOUBLE PRECISION MAXIMUM
      COMMON/CBBEAM/MAXIMUM
      INTEGER IDC,GAGA
C...Decay information.
      INTEGER MDCY,MDME,KFDP
      DOUBLE PRECISION BRAT
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      
      CHARACTER CH1*16,CH2*16
      
      INTEGER DC,CDC
C     order of DC is d u s c b t b' t' ...
C     ... e nu_e mu nu_mu tau nu_tau tau' nu_tau'
      COMMON/DECAYCH/DC(18)
      
      IF(FIRST) THEN
         MAXIMUM=10.6
         WRITE (*,*) "Decay of Z0/gamma to"
C...  Only allow Z0 decay to quarks (i.e. no leptonic final states).
C     RS veto also b quarks
         DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
C     first switch channel off
            MDME(IDC,1)=MIN(0,MDME(IDC,1))
            CDC=IABS(KFDP(IDC,1))
C     then switch on again if requested
            CALL PYNAME(CDC,CH1)
            CALL PYNAME(-CDC,CH2)
            IF(DC(CDC).EQ.1) THEN
               MDME(IDC,1)=1
               WRITE (*,*) "  ",ch1,"+ ",ch2," ALLOWED"
            ELSE
               WRITE (*,*) "  ",ch1,"+ ",ch2," DISABLED"
            ENDIF
         ENDDO
C     CALL PYSTAT(2)
         
         IF(GAGA.EQ.0) THEN
            CALL PYINIT('CMS','E+','E-',MAXIMUM)
         ELSE
            CALL PYINIT('CMS','GAMMA/E+','GAMMA/E-',MAXIMUM)
         ENDIF
C     CALL PYLIST(0)
C         WRITE (*,*) 'done PYTHIA initialization'//
C     >        ' with varying beam energy'
C         WRITE (*,*) 'maximum allowed energy is ',MAXIMUM,' GeV'
         FIRST=.FALSE.
      ENDIF
      RETURN
      END
      
C--------------------------------------------------------------------------
C
C Environment:
C      This software is part of the EvtGen package developed jointly
C      for the BaBar and CLEO collaborations.  If you use all or part
C      of it, please give an appropriate acknowledgement.
C
C Copyright Information: See BelEvtGen/COPYRIGHT
C      Copyright (C) 1998      Caltech, UCSB
C
C Module: continuum.F
C
C Description:
C
C Modification history:
C
C    DJL/RYD     August 11, 1998         Module created
C     26-Sep-2002 - RS : changed to PYTHIA
C
C------------------------------------------------------------------------
      SUBROUTINE PYCONTINUUM(ENERGY,NDAUG,KF,PX,PY,PZ,E)
      IMPLICIT NONE
      DOUBLE PRECISION MAXIMUM
      COMMON/CBBEAM/MAXIMUM
      
      DOUBLE PRECISION P,V
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      
      INTEGER MSTU,MSTJ
      DOUBLE PRECISION PARU,PARJ
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      
      INTEGER MSTP,MSTI
      DOUBLE PRECISION PARP,PARI
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
C...Decay information.
      INTEGER MDCY,MDME,KFDP
      DOUBLE PRECISION BRAT
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      
      INTEGER N,NPAD,K,IMSTJ
      
      REAL*8 PXSUM,PYSUM,PZSUM
      
      REAL*8 ENERGY
      
      INTEGER PYCOMP
      EXTERNAL PYCOMP
      
      INTEGER KF(100),I,NDAUG
      REAL*8  PX(100),PY(100),PZ(100),E(100)

C     RS Particles should not decay at this stage
C     but remember the previous setting
      IMSTJ=MSTJ(21)
      MSTJ(21)=0
      
C      IF(FLAVOR.NE.LastProdFLav) THEN
C...  Only allow Z0 decay to decay to certain flavour
C         DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
C            IF(IABS(KFDP(IDC,1)).NE.Flavor)
C     >           MDME(IDC,1)=MIN(0,MDME(IDC,1))
C         ENDDO
C      ENDIF
      
 4    PARP(171)=Energy/Maximum
      CALL PYEVNT
      
C--- only primary quarks and stable particles remain in the event record
      CALL PYEDIT(5)
      
C     RS and allow decays
      MSTJ(21)=IMSTJ
      
      NDAUG=0
C     sum to check that we preserve momentum
      PXSUM=0.0
      PYSUM=0.0
      PZSUM=0.0
      
      DO 1,I=1,N
         NDAUG=NDAUG+1
         KF(NDAUG)=K(I,2)
         PX(NDAUG)=P(I,1)
         PY(NDAUG)=P(I,2)
         PZ(NDAUG)=P(I,3)
         E(NDAUG) =P(I,4)
         IF(K(I,1).LT.10)THEN
            PXSUM=PXSUM+PX(NDAUG)
            PYSUM=PYSUM+PY(NDAUG)
            PZSUM=PZSUM+PZ(NDAUG)
         ENDIF
1     CONTINUE
      
      IF (ABS(PXSUM).GT.0.001.OR.
     +       ABS(PYSUM).GT.0.001.OR.
     +       ABS(PZSUM).GT.0.001) THEN


         PRINT *, 'Momentum not conserved in jetset fragmentation:'
         PRINT *,'dPx:',pxsum,' dPy:',pysum,' dPz:',pzsum

         CALL PYLIST(2)

         GOTO 4
         
      ENDIF
      RETURN
      END
C--------------------------------------------------------------------------
C
C Environment:
C      This software is part of the EvtGen package developed jointly
C      for the BaBar and CLEO collaborations.  If you use all or part
C      of it, please give an appropriate acknowledgement.
C
C Copyright Information: See BelEvtGen/COPYRIGHT
C      Copyright (C) 1998      Caltech, UCSB
C
C Module: jetset1.F
C
C Description:
C
C Modification history:
C
C    DJL/RYD     August 11, 1998         Module created
C     26-Sep-2002 - RS : changed to PYTHIA
C
C------------------------------------------------------------------------
      SUBROUTINE PYTHIADEC(IP,M,NDAUG,KF,KM,PX,PY,PZ,E)

C
C interface to JETSET 7.4 to have one particle decayed
C including possibly fragmentation, if the decay products include
C partons.
C     
      IMPLICIT NONE
      
      INTEGER MSTU,MSTJ
      DOUBLE PRECISION PARU,PARJ
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      SAVE /PYDAT1/
      DOUBLE PRECISION P,V
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      
      INTEGER N,NPAD,K
      
      INTEGER NDMAX
      PARAMETER (NDMAX=100)
      
      INTEGER PYCOMP
      EXTERNAL PYCOMP
      
      INTEGER IP,KF(NDMAX),I,NDAUG,KM(NDMAX)
      INTEGER KP,KID,IPART1,IPART
      REAL*8  M,PX(NDMAX),PY(NDMAX),PZ(NDMAX),E(NDMAX)
      
      INTEGER IP1,IMSTJ
      DOUBLE PRECISION QMAX
      
C      CALL PY1ENT(1,IP,M,0D0,0D0)
 
      K(1,1)=1
      K(1,2)=ip
      P(1,5)=m
      P(1,4)=m
      P(1,1)=0.0
      P(1,2)=0.0
      P(1,3)=0.0
      n=1
     
C     RS Particles should not decay at this stage
C     but remember the previous setting
      IMSTJ=MSTJ(21)
      MSTJ(21)=0
      
      CALL PYDECY(1)
      
C     RS and allow decays
      MSTJ(21)=IMSTJ


C code copied from LUEXEC to avoid error with shower
C switched on
C...Decay products may develop a shower.
        IF(MSTJ(92).GT.0) THEN
          IP1=MSTJ(92)
          QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
     &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
          CALL PYSHOW(IP1,IP1+1,QMAX)
          CALL PYPREP(IP1)
          MSTJ(92)=0
        ELSEIF(MSTJ(92).LT.0) THEN
          IP1=-MSTJ(92)
          CALL PYSHOW(IP1,-3,P(IP,5))
          CALL PYPREP(IP1)
          MSTJ(92)=0
        ENDIF

c
c for debugging:
c      call lulist(1)
c
      mstj(21)=0
      call pyexec
      mstj(21)=2
 
c find partons, delete secondary partons, set mother pointers

      ndaug = 0
      ipart1 = 1
      ipart = 1

      do 10 i=2,n
         kp = k(i,3)
         kid = k(i,2)
         if (abs(kid) .ge. 1 .and. abs(kid) .le. 8
     1        .or. kid .eq. 21
     2        .or. kid .ge. 91 .and. kid .le. 94)  then
            if (ipart1 .eq. 1) ipart1 = i
            ipart = i
            if (kp .ne. 1)  goto 10
            kp = 0
         else
            if (kp .gt. ipart)  then
               goto 10
            elseif (kp .ge. ipart1)  then
               kp = ipart1-1
            else
               kp = 0
            endif
         endif
         ndaug = ndaug + 1
         km(ndaug)=kp
         kf(ndaug)=kid
         px(ndaug)=p(i,1)
         py(ndaug)=p(i,2)
         pz(ndaug)=p(i,3)
         e(ndaug)=p(i,4)

C        print '( 2I5,I12,4F12.4 )',ndaug,km(ndaug),kf(ndaug),
C     1    px(ndaug),py(ndaug),pz(ndaug),e(ndaug)

 10   CONTINUE

C      NDAUG=0
C      DO I=2,N
C        NDAUG = NDAUG + 1
C        KM(NDAUG)=K(I,3)
C        KF(NDAUG)=K(I,2)
C        PX(NDAUG)=P(I,1)
C        PY(NDAUG)=P(I,2)
C        PZ(NDAUG)=P(I,3)
C        E(NDAUG)=P(I,4)
C      ENDDO


      RETURN
      END
Commit	Line	Data
da0e9ce3	1	C--------------------------------------------------------------------------
	2	C
	3	C Environment:
	4	C This software is part of the EvtGen package developed jointly
	5	C for the BaBar and CLEO collaborations. If you use all or part
	6	C of it, please give an appropriate acknowledgement.
	7	C
	8	C Copyright Information: See BelEvtGen/COPYRIGHT
	9	C Copyright (C) 1998 Caltech, UCSB
	10	C
	11	C Module: EvtJetSetInit.F
	12	C
	13	C Description:
	14	C
	15	C Modification history:
	16	C
	17	C DJL/RYD August 11, 1998 Module created
	18	C RS October 28, 2002 copied and modified from Jetset
	19	C
	20	C------------------------------------------------------------------------
	21	SUBROUTINE EVTPYTHIAINIT(FNAME)
	22	IMPLICIT NONE
	23	EXTERNAL PYDATA
	24
	25	CHARACTER() FNAME
	26
	27	PRINT *,'PYUPDA : : ',FNAME
	28	OPEN(54,STATUS='OLD',FILE=FNAME)
	29	CALL PYUPDA(3,54)
	30	CLOSE(54)
	31	WRITE (,) 'PYUPDA DONE'
	32	RETURN
	33	END
	34	SUBROUTINE INITPYTHIA(GAGA)
	35	IMPLICIT NONE
	36	INTEGER PYCOMP
	37	EXTERNAL PYDATA,PYCOMP
	38	LOGICAL FIRST
	39	SAVE FIRST
	40	DATA FIRST / .TRUE. /
	41	DOUBLE PRECISION MAXIMUM
	42	COMMON/CBBEAM/MAXIMUM
	43	INTEGER IDC,GAGA
	44	C...Decay information.
	45	INTEGER MDCY,MDME,KFDP
	46	DOUBLE PRECISION BRAT
	47	COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
	48
	49	CHARACTER CH116,CH216
	50
	51	INTEGER DC,CDC
	52	C order of DC is d u s c b t b' t' ...
	53	C ... e nu_e mu nu_mu tau nu_tau tau' nu_tau'
	54	COMMON/DECAYCH/DC(18)
	55
	56	IF(FIRST) THEN
	57	MAXIMUM=10.6
	58	WRITE (,) "Decay of Z0/gamma to"
	59	C... Only allow Z0 decay to quarks (i.e. no leptonic final states).
	60	C RS veto also b quarks
	61	DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
	62	C first switch channel off
	63	MDME(IDC,1)=MIN(0,MDME(IDC,1))
	64	CDC=IABS(KFDP(IDC,1))
65	C then switch on again if requested
66	CALL PYNAME(CDC,CH1)
67	CALL PYNAME(-CDC,CH2)
68	IF(DC(CDC).EQ.1) THEN
69	MDME(IDC,1)=1
70	WRITE (,) " ",ch1,"+ ",ch2," ALLOWED"
71	ELSE
72	WRITE (,) " ",ch1,"+ ",ch2," DISABLED"
73	ENDIF
74	ENDDO
75	C CALL PYSTAT(2)
76
77	IF(GAGA.EQ.0) THEN
78	CALL PYINIT('CMS','E+','E-',MAXIMUM)
79	ELSE
80	CALL PYINIT('CMS','GAMMA/E+','GAMMA/E-',MAXIMUM)
81	ENDIF
82	C CALL PYLIST(0)
83	C WRITE (,) 'done PYTHIA initialization'//
84	C > ' with varying beam energy'
85	C WRITE (,) 'maximum allowed energy is ',MAXIMUM,' GeV'
86	FIRST=.FALSE.
87	ENDIF
88	RETURN
89	END
90
91	C--------------------------------------------------------------------------
92	C
93	C Environment:
94	C This software is part of the EvtGen package developed jointly
95	C for the BaBar and CLEO collaborations. If you use all or part
96	C of it, please give an appropriate acknowledgement.
97	C
98	C Copyright Information: See BelEvtGen/COPYRIGHT
99	C Copyright (C) 1998 Caltech, UCSB
100	C
101	C Module: continuum.F
102	C
103	C Description:
104	C
105	C Modification history:
106	C
107	C DJL/RYD August 11, 1998 Module created
108	C 26-Sep-2002 - RS : changed to PYTHIA
109	C
110	C------------------------------------------------------------------------
111	SUBROUTINE PYCONTINUUM(ENERGY,NDAUG,KF,PX,PY,PZ,E)
112	IMPLICIT NONE
113	DOUBLE PRECISION MAXIMUM
114	COMMON/CBBEAM/MAXIMUM
115
116	DOUBLE PRECISION P,V
117	COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
118
119	INTEGER MSTU,MSTJ
120	DOUBLE PRECISION PARU,PARJ
121	COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
122
123	INTEGER MSTP,MSTI
124	DOUBLE PRECISION PARP,PARI
125	COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
126	C...Decay information.
127	INTEGER MDCY,MDME,KFDP
128	DOUBLE PRECISION BRAT
129	COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
130
131	INTEGER N,NPAD,K,IMSTJ
132
133	REAL*8 PXSUM,PYSUM,PZSUM
134
135	REAL*8 ENERGY
136
137	INTEGER PYCOMP
138	EXTERNAL PYCOMP
139
140	INTEGER KF(100),I,NDAUG
141	REAL*8 PX(100),PY(100),PZ(100),E(100)
142
143	C RS Particles should not decay at this stage
144	C but remember the previous setting
145	IMSTJ=MSTJ(21)
146	MSTJ(21)=0
147
148	C IF(FLAVOR.NE.LastProdFLav) THEN
149	C... Only allow Z0 decay to decay to certain flavour
150	C DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
151	C IF(IABS(KFDP(IDC,1)).NE.Flavor)
152	C > MDME(IDC,1)=MIN(0,MDME(IDC,1))
153	C ENDDO
154	C ENDIF
155
156	4 PARP(171)=Energy/Maximum
157	CALL PYEVNT
158
159	C--- only primary quarks and stable particles remain in the event record
160	CALL PYEDIT(5)
161
162	C RS and allow decays
163	MSTJ(21)=IMSTJ
164
165	NDAUG=0
166	C sum to check that we preserve momentum
167	PXSUM=0.0
168	PYSUM=0.0
169	PZSUM=0.0
170
171	DO 1,I=1,N
172	NDAUG=NDAUG+1
173	KF(NDAUG)=K(I,2)
174	PX(NDAUG)=P(I,1)
175	PY(NDAUG)=P(I,2)
176	PZ(NDAUG)=P(I,3)
177	E(NDAUG) =P(I,4)
178	IF(K(I,1).LT.10)THEN
179	PXSUM=PXSUM+PX(NDAUG)
180	PYSUM=PYSUM+PY(NDAUG)
181	PZSUM=PZSUM+PZ(NDAUG)
182	ENDIF
183	1 CONTINUE
184
185	IF (ABS(PXSUM).GT.0.001.OR.
186	+ ABS(PYSUM).GT.0.001.OR.
187	+ ABS(PZSUM).GT.0.001) THEN
188
189
190	PRINT *, 'Momentum not conserved in jetset fragmentation:'
191	PRINT *,'dPx:',pxsum,' dPy:',pysum,' dPz:',pzsum
192
193	CALL PYLIST(2)
194
195	GOTO 4
196
197	ENDIF
198	RETURN
199	END
200	C--------------------------------------------------------------------------
201	C
202	C Environment:
203	C This software is part of the EvtGen package developed jointly
204	C for the BaBar and CLEO collaborations. If you use all or part
205	C of it, please give an appropriate acknowledgement.
206	C
207	C Copyright Information: See BelEvtGen/COPYRIGHT
208	C Copyright (C) 1998 Caltech, UCSB
209	C
210	C Module: jetset1.F
211	C
212	C Description:
213	C
214	C Modification history:
215	C
216	C DJL/RYD August 11, 1998 Module created
217	C 26-Sep-2002 - RS : changed to PYTHIA
218	C
219	C------------------------------------------------------------------------
220	SUBROUTINE PYTHIADEC(IP,M,NDAUG,KF,KM,PX,PY,PZ,E)
221
222	C
223	C interface to JETSET 7.4 to have one particle decayed
224	C including possibly fragmentation, if the decay products include
225	C partons.
226	C
227	IMPLICIT NONE
228
229	INTEGER MSTU,MSTJ
230	DOUBLE PRECISION PARU,PARJ
231	COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
232	SAVE /PYDAT1/
233	DOUBLE PRECISION P,V
234	COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
235
236	INTEGER N,NPAD,K
237
238	INTEGER NDMAX
239	PARAMETER (NDMAX=100)
240
241	INTEGER PYCOMP
242	EXTERNAL PYCOMP
243
244	INTEGER IP,KF(NDMAX),I,NDAUG,KM(NDMAX)
245	INTEGER KP,KID,IPART1,IPART
246	REAL*8 M,PX(NDMAX),PY(NDMAX),PZ(NDMAX),E(NDMAX)
247
248	INTEGER IP1,IMSTJ
249	DOUBLE PRECISION QMAX
250
251	C CALL PY1ENT(1,IP,M,0D0,0D0)
252
253	K(1,1)=1
254	K(1,2)=ip
255	P(1,5)=m
256	P(1,4)=m
257	P(1,1)=0.0
258	P(1,2)=0.0
259	P(1,3)=0.0
260	n=1
261
262	C RS Particles should not decay at this stage
263	C but remember the previous setting
264	IMSTJ=MSTJ(21)
265	MSTJ(21)=0
266
267	CALL PYDECY(1)
268
269	C RS and allow decays
270	MSTJ(21)=IMSTJ
271
272
273
274
275	C code copied from LUEXEC to avoid error with shower
276	C switched on
277	C...Decay products may develop a shower.
278	IF(MSTJ(92).GT.0) THEN
279	IP1=MSTJ(92)
280	QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
281	& 1))2-(P(IP1,2)+P(IP1+1,2))2-(P(IP1,3)+P(IP1+1,3))**2))
282	CALL PYSHOW(IP1,IP1+1,QMAX)
283	CALL PYPREP(IP1)
284	MSTJ(92)=0
285	ELSEIF(MSTJ(92).LT.0) THEN
286	IP1=-MSTJ(92)
287	CALL PYSHOW(IP1,-3,P(IP,5))
288	CALL PYPREP(IP1)
289	MSTJ(92)=0
290	ENDIF
291
292	c
293	c for debugging:
294	c call lulist(1)
295	c
296	mstj(21)=0
297	call pyexec
298	mstj(21)=2
299
300	c find partons, delete secondary partons, set mother pointers
301
302	ndaug = 0
303	ipart1 = 1
304	ipart = 1
305
306	do 10 i=2,n
307	kp = k(i,3)
308	kid = k(i,2)
309	if (abs(kid) .ge. 1 .and. abs(kid) .le. 8
310	1 .or. kid .eq. 21
311	2 .or. kid .ge. 91 .and. kid .le. 94) then
312	if (ipart1 .eq. 1) ipart1 = i
313	ipart = i
314	if (kp .ne. 1) goto 10
315	kp = 0
316	else
317	if (kp .gt. ipart) then
318	goto 10
319	elseif (kp .ge. ipart1) then
320	kp = ipart1-1
321	else
322	kp = 0
323	endif
324	endif
325	ndaug = ndaug + 1
326	km(ndaug)=kp
327	kf(ndaug)=kid
328	px(ndaug)=p(i,1)
329	py(ndaug)=p(i,2)
330	pz(ndaug)=p(i,3)
331	e(ndaug)=p(i,4)
332
333	C print '( 2I5,I12,4F12.4 )',ndaug,km(ndaug),kf(ndaug),
334	C 1 px(ndaug),py(ndaug),pz(ndaug),e(ndaug)
335
336	10 CONTINUE
337
338	C NDAUG=0
339	C DO I=2,N
340	C NDAUG = NDAUG + 1
341	C KM(NDAUG)=K(I,3)
342	C KF(NDAUG)=K(I,2)
343	C PX(NDAUG)=P(I,1)
344	C PY(NDAUG)=P(I,2)
345	C PZ(NDAUG)=P(I,3)
346	C E(NDAUG)=P(I,4)
347	C ENDDO
348
349
350	RETURN
351	END