Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25pomg4.f
Go to the documentation of this file.
1 C---------------------------------------------------------------------
2 C dtcpomqj.f (July 1993)
3 C-----------------------------------------------------------------------
4 ************************************************************************
5 *
6  SUBROUTINE sigshd(ECM)
7 * May 1991
8 * input:
9 * hard cross sections (see DATA statements)
10 * ECM (independent of CMENER in /USER/)
11 * output:
12 * bar cross sections for soft (SIGSOF) hard (SIGHAR) and diffractive
13 * (SIGTRP) SCATTERING
14 *
15 *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16 *
17 * version determined by ISIG
18 C*********************************************************************
19 C ISIG=1-9 dropped since dpmjet-II.4.2
20 C*********************************************************************
21 * ISIG=1 Capella,Tran Thanh Van,Kwiecinski,PRL 58(1987)2015
22 * ISIG=2 Capella,Tran Thanh Van,Kwiecinski,PRL 58(1987)2015 CHG.SIG
23 * ISIG=3 ALL SUBPROCESSES HARD CROSS SECTIONS PTHR CHANGING WITH ECM
24 * ISIG=4 ALL SUBPROCESSES HARD CROSS SECTIONS PTHR=3 GEV/C MRS1
25 * ISIG=5 ALL SUBPROCESSES HARD CROSS SECTIONS PTHR=2 GEV/C MRS1
26 * ISIG=6 ALL SUBPROCESSES HARD CROSS SECTIONS PTHR=1.3 GEV/C MRS1
27 * ISIG=7 all subprocesses hard cross sections PTHR=2GEV/C MRS1
28 * ISIG=8 program written by Patrick and Maire
29 * ISIG=9 ALL SUBPROCESSES HARD CROSS SECTIONS PTHR=1.0 GEV/C MRS1
30 C*********************************************************************
31 C ISIG=1-9 dropped since dpmjet-II.4.2
32 C*********************************************************************
33 * ISIG=10 Version ISIG=4 including different sets
34 * of structure functions
35 *
36 *-----------------------------------------------------------------------
37  IMPLICIT DOUBLE PRECISION(a-h,o-z)
38  SAVE
39 * CONVERSION FACTOR GEV**-2 TO MILLIBARNS
40  parameter(conv=.38935d0)
41  parameter(pi=3.141592654d0,
42  & three=3.d0,
43  & two =2.d0,
44  & eps=1.d-3)
45  parameter(thousa = 1000.d0)
46 *
47 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
48  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
49 *
50 * *** /POMPAR/*/SIGMA/*/POMTYP/ used only in SIGMAPOM-routines (->POMDI)
51  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
52  common/pompar/alfa,alfap,a,c,ak
53  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
54  * sigloo,zloo
55 *
56  CHARACTER*80 title
57  CHARACTER*8 projty,targty
58 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
59 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
60  COMMON /user1/title,projty,targty
61  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
62 *
63  COMMON /strufu/istrum,istrut
64  COMMON /alala/alalam
65 C COMMON/COLLIS/S,IJPROJ,IJTAR,PTTHR,PTTHR2,IOPHRD,IJPRLU,IJTALU
66 C Prevent initialization energy
67 C to pollute s in /COLLIS/
68  COMMON /collpo/s,ptthr,ptthr2
69 C COMMON /COLLIS/SPO,IJPROJ,IJTAR,PTTPO,IOPHRD,IJPRLU,IJTALU,
70 C * PTTPO2
71  common/collis/spo,ijproj,ijtar,pttpo,pttpo2,iophrd,ijprlu,ijtalu
72  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
73 *
74  dimension sqs(13)
75  dimension sqsj(17)
76  dimension xsqsj(21),xxhhj4(21)
77 *
78 *
79 *
80 * used in ISIG=3,5,6,7,9
81  DATA xsqsj/0.005,0.01,0.02,0.035,0.053,
82  * 0.1,0.2,0.35,0.54,1.,2.,5.,
83  *10.,20.,40.,100.,200.,400.,1000.,2000.,4000./
84 *
85 * used in ISIG 10,40
86  DATA sqs/1.,2.,3.,4.,5.,10.,20.,30.,40.,100.,200.,500.,1000./
87 *
88  pi4=4.*pi
89  s=ecm**2
90 ***********************************************************************
91 *----------------------------------------------------------------------
92 *
93 * ** select option used
94 *
95  go to(10,20,30,40,50,60,70,80,90,100),isig
96 *
97  10 CONTINUE
98  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
99 *----------------------------------------------------------------------
100 * nach: Capella, Tran Thanh Van, Kwiecinski, PRL 58(1987)2015
101 * as used in Ranft et al. SSC 149 Eq. 4,5
102 *
103  20 CONTINUE
104  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
105 *----------------------------------------------------------------------
106 * nach: Capella,Tran Thanh Van,Kwiecinski,PRL 58(1987)2015 CHG.SIG
107 *
108  30 CONTINUE
109  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
110 *----------------------------------------------------------------------
111 * nach: all subprocesses hard cross sections LIKE 70
112 C PTHR RISING WITH ECM
113 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
114 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
115  40 CONTINUE
116  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
117 *----------------------------------------------------------------------
118 * nach: all subprocesses hard cross sections LIKE 70
119 C PTHR RISING WITH ECM ptthr=3gev
120 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
121 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
122  50 CONTINUE
123  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
124 *----------------------------------------------------------------------
125 * nach: all subprocesses hard cross sections LIKE 70
126 C PTHR RISING WITH ECM ptthr=2 gev
127 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
128 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
129  60 CONTINUE
130  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
131 *----------------------------------------------------------------------
132 * nach: all subprocesses hard cross sections LIKE 70
133 C PTHR RISING WITH ECM ptthr=1.3 GEV
134 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
135 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
136  70 CONTINUE
137  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
138 *----------------------------------------------------------------------
139 * nach: all subprocesses hard cross sections
140 *
141  80 CONTINUE
142  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
143 *----------------------------------------------------------------------
144 * nach: Patrick Maires program
145  90 CONTINUE
146  WRITE(6,*)' This value of ISIG no longer available ISIG=',isig
147 *----------------------------------------------------------------------
148 * nach: all subprocesses hard cross sections LIKE 70
149 C PTHR RISING WITH ECM ptthr=1.0 GEV
150 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
151 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
152  RETURN
153  100 CONTINUE
154 *----------------------------------------------------------------------
155 * nach: all subprocesses hard cross sections LIKE 70
156 C PTHR RISING WITH ECM
157 * NEW SIGTRP RISING LIKE LOG(S) A.CAPELLA 30.3.90
158 C ONLY FOR USE WITH PRBLM2!!!!!!!!!!!!!!!!
159 C INTRODUCED 18.12.90
160 C BY DIETER PERTERMANN
161 C modified 11-06-92 (R.Engel)
162 C
163 C default parameter set
164  alfa=1.076
165  alfap=0.24
166  a=40.8
167  bh=3.51
168  bhoo=bh
169  bsoo=bh
170  ak=1.5
171 C ALALAM --> ALAM (see PRBLM2 and POMDI)
172  alalam=0.0
173 C begin fixed ptthr=3GeV
174 C--------------------------------------------------------------------
175  IF(abs(ptthr-three).LT.eps) THEN
176  WRITE(6,*)' PTTHR=3. not available in dpmjet25'
177  WRITE(6,*) ' WARNING: no model parameter set available'
178  WRITE(6,*) ' for this combination of PTCUT and ISTRUF'
179  WRITE(6,*) ' (initialization using default values)'
180  alfa = 1.078
181  a = 42.6
182  alalam=0.740
183  aqqal=1.d0
184  aqqpd=1.d0
185  alfap= 0.24
186  ak=2.0
187  ENDIF
188 C end fixed ptthr=3GeV
189 C begin fixed ptthr=2GeV
190  IF(abs(ptthr-two).LT.eps) THEN
191  WRITE(6,*)' PTTHR=2. not available in dpmjet25'
192  WRITE(6,*) ' WARNING: no model parameter set available'
193  WRITE(6,*) ' for this combination of PTCUT and ISTRUF'
194  WRITE(6,*) ' (initialization using default values)'
195  alfa = 1.042
196  a = 64.54
197  alalam=0.6402
198  aqqal=1.d0
199  aqqpd=1.d0
200  alfap= 0.24
201  ak=2.0
202  ENDIF
203 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
204 C end fixed ptthr=2GeV
205 C-----------------------------------------------------------
206 C-----------------------------------------------------------
207 C-----------------------------------------------------------
208 C begin ptthr= PTTHR=2.1+0.15*(LOG10(ECM/50.))**3
209 C-----------------------------------------------------------
210  IF(istrut.EQ.1) THEN
211  WRITE(6,*)' ISTRUT=1 (PTTHR=2.1+0.15*(LOG10(ECM/50.))**3)',
212  * 'not available in dpmjet25'
213  ptthr=2.1+0.15*(log10(ecm/50.))**3
214  ptthr2=ptthr
215  WRITE(6,*) ' WARNING: no model parameter set available'
216  WRITE(6,*) ' for this combination of PTCUT and ISTRUF'
217  WRITE(6,*) ' (initialization using default values)'
218  alfa = 1.042
219  a = 64.54
220  alalam=0.6402
221  aqqal=1.d0
222  aqqpd=1.d0
223  alfap= 0.24
224  ak=2.0
225  ENDIF
226 C end PTTHR=2.1+0.15*(LOG10(ECM/50.))**3
227 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
228  bhoo=bh
229  bsoo=bh
230 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
231 C-----------------------------------------------------------
232 C-----------------------------------------------------------
233 C-----------------------------------------------------------
234 C begin PTTHR=2.5+0.12*(LOG10(ECM/50.))**3
235 C-----------------------------------------------------------
236  IF(istrut.EQ.2) THEN
237  ptthr=2.5+0.12*(log10(ecm/50.))**3
238  ptthr2=ptthr
239  IF( istruf.EQ.9 ) THEN
240  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
241  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
242  go to 778
243  ELSEIF( istruf.EQ.10 ) THEN
244  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
245  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
246  go to 778
247  ELSEIF( istruf.EQ.11 ) THEN
248  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
249  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
250  go to 778
251  ELSEIF( istruf.EQ.12 ) THEN
252  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
253  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
254  go to 778
255  ELSEIF( istruf.EQ.13 ) THEN
256  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
257  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
258  go to 778
259  ELSEIF( istruf.EQ.14 ) THEN
260  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
261  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
262  go to 778
263  ELSEIF( istruf.EQ.15 ) THEN
264  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
265  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
266 C CETQ PDFs with other scale
267  go to 778
268  ELSEIF( istruf.EQ.16 ) THEN
269  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
270  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
271  ak=2.0
272  go to 778
273  ELSEIF( istruf.EQ.17 ) THEN
274  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
275  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
276  go to 778
277  ELSEIF( istruf.EQ.18 ) THEN
278  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
279  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
280  go to 778
281  ELSEIF( istruf.EQ.19 ) THEN
282  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
283  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
284  go to 778
285  ELSEIF( istruf.EQ.20 ) THEN
286  WRITE(6,*)' ISTRUT=2 (PTTHR=2.5+0.12*(LOG10(ECM/50.))**3)',
287  * 'and ISTRUF= ',istruf ,' not available in dpmjet25'
288 C GRV94LO AK=1. only for ISTRUT = 2
289  go to 778
290  ELSEIF( istruf.EQ.21 ) THEN
291  alfa = 1.0733
292  alfap= 0.171
293  a = 47.84
294  alalam=0.621
295  bsoo=1.58
296  bhoo=3.54
297  ak=1.000
298 C GRV94LO AK=2. only for ISTRUT = 2
299  ELSEIF( istruf.EQ.22 ) THEN
300  alfa = 1.0513
301  alfap= 0.3246
302  a = 55.16
303  alalam=0.5846
304  bsoo=1.114
305  bhoo=1.703
306  ak=2.000
307 C CTEQ96 AK=2. only for ISTRUT = 2
308  ELSEIF( istruf.EQ.23 ) THEN
309  alfa = 1.0448
310  alfap= 0.372
311  a = 57.51
312  alalam=0.566
313  bsoo=0.97
314  bhoo=1.47
315  ak=2.000
316  ELSE
317  778 CONTINUE
318  WRITE(6,*) ' WARNING: no model parameter set available'
319  WRITE(6,*) ' for this combination of PTCUT and ISTRUF'
320  WRITE(6,*) ' (initialization using default values)'
321  alfa = 1.042
322  a = 64.54
323  alalam=0.6402
324  aqqal=1.d0
325  aqqpd=1.d0
326  alfap= 0.24
327  ak=2.0
328  ENDIF
329  bhoo = bhoo/conv
330  bsoo = bsoo/conv
331  ENDIF
332 C end PTTHR=2.5+0.12*(LOG10(ECM/50.))**3
333 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
334 C slopes in GeV**-2
335  bs=bsoo+alfap*log(s)
336  bh=bhoo
337 * BS=BSOO+ALFAP*LOG(S)
338 C change units to mb
339  bh=bh*conv
340  bs=bs*conv
341  bt=bs
342 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
343 C BT=BS/2.
344  c=40.
345 C CHANGED 13.1.90 BY J.R.
346 C C=1.8
347 C\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
348 C C=1.E-8
349 * parametrizations of input cross sections
350  sigsof=a*s**(alfa-1.)
351 * *** hard X-section
352 * read interpolation data for different
353 * sets of structure functions:
354 *
355  CALL rdxsec(xxhhj4)
356  IF(istruf.EQ.21)ak=2.
357  sighar=1.e-8
358  IF (s.GT.2450.d0)
359  * sighar=ak*0.1*(s-2450.)**0.35
360  IF(ecm.GE.thousa*xsqsj(2)) THEN
361  DO 1031 i=1,20
362  iii=i+1
363  IF(ecm.LT.xsqsj(iii)*thousa.AND.
364  * ecm.GE.thousa*xsqsj(i))THEN
365  dsq=ecm-thousa*xsqsj(i)
366  ddsq=thousa*(xsqsj(iii)-xsqsj(i))
367  dhs=(xxhhj4(iii)-xxhhj4(i))
368  sighar=ak*(xxhhj4(i)+dhs*dsq/ddsq)*0.5
369  ENDIF
370  1031 CONTINUE
371  ENDIF
372 C
373 C *** trippel pomeron X-section
374 C
375 C VERSION A.CAPELLA 30.3.90
376  gca=sqrt(a)
377  g3ca=gca**3
378  gaca=0.42
379 C BSDOCA=1.372
380  bsdoca=bsoo*conv
381 C ALSCA=.0925
382  alsca=alfap*conv
383  alns=log(s)
384  bsdca=bsdoca+2.*alsca*alns
385  sigtrp=g3ca*gaca*log(s/10.)/(8.*3.14*bsdca)
386  IF (sigtrp.LT.0.d0)sigtrp=0.01
387 C
388  bddca=2.*alsca*alns
389  alo1sq=(log(s/400.))**2
390  alo2sq=(log(25./s))**2
391  alo3sq=(log(5./20.))**2
392  sigloo=a*gaca**2*(alo1sq+alo2sq-2.*alo3sq)/(32.*3.14*bddca)
393 C
394  zsof=sigsof/(pi4*bs)
395  zhar=sighar/(pi4*bh)
396  ztrp=sigtrp/(pi4*bt)
397  zloo=sigloo/(pi4*bt)
398 C
399 C WRITE(6,'(2(/1X,A))') 'SELECTED PARAMETERS:',
400 C & '===================='
401 C WRITE(6,'(1X,A,E12.3)') ' ALFA ',ALFA
402 C WRITE(6,'(1X,A,E12.3)') ' ALFAP ',ALFAP
403 C WRITE(6,'(1X,A,E12.3)') ' A ',A
404 C WRITE(6,'(1X,A,2E12.3)') ' BS,BSOO',BS,BSOO*CONV
405 C WRITE(6,'(1X,A,2E12.3)') ' BH,BHOO',BH,BHOO*CONV
406 C WRITE(6,'(1X,A,E12.3)') ' GACA ',GACA
407 C WRITE(6,'(1X,A,E12.3,/)') ' AK ',AK
408 C
409  RETURN
410  END
411 *
412 *
413 ************************************************************************
414 ************************************************************************
415 *
416  SUBROUTINE rdxsec(XSEC)
417 C
418 C 18.12.90, Dieter Pertermann
419 C 15.03.93, 27.05.93 modified (R. Engel)
420 C
421 C RDXSEC READS DATA FOR INTERPOLATION
422 C OF THE TOTAL CROSS SECTIONS FOR DIFFERENT
423 C SETS OF STRUCTURE FUNCTIONS. THE CHOICE
424 C OF THE CORRESPONDIG DATA SET IS CONTROLED
425 C BY THE OVERALL STRUCTURE FUNCTION PARAMETER
426 C ISTRUF.
427 C
428 C CMENER taken out but no action necessary as not in SUB
429 C
430 C modified 11-05-92 (R.Engel)
431 *
432  IMPLICIT DOUBLE PRECISION(a-h,o-z)
433  SAVE
434  CHARACTER*80 title
435  CHARACTER*8 projty,targty
436 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
437 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
438  COMMON /user1/title,projty,targty
439  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
440 *
441 C COMMON/COLLIS/S,IJPROJ,IJTAR,PTTHR,PTTHR2,IOPHRD,IJPRLU,IJTALU
442 *
443  COMMON /collpo/s,ptthr,ptthr2
444 C COMMON /COLLIS/SPO,IJPROJ,IJTAR,PTTPO,IOPHRD,IJPRLU,IJTALU,
445 C * PTTPO2
446  common/collis/spo,ijproj,ijtar,pttpo,pttpo2,iophrd,ijprlu,ijtalu
447  COMMON /strufu/istrum,istrut
448 C
449  dimension xsec(21)
450  dimension xs21(189)
451 C
452  parameter(epsil=1.d-4,
453  & three=3.d0,
454  & two=2.d0)
455 C
456  DATA xs21 /
457  & 0.000000e+00,0.137854e-04, .02, .13, .37, 1.32,
458  & 3.88, 8.02, 13.15, 24.32, 43.43, 79.69, 113.13,
459  & 147.5, 180.47, 221.01, 250.37,
460  & 279.4, 320.1, 349.6, 381.6,
461 * total X-section , cteq1M PTTHR=x GEV/C
462  & .000000e+00, .494767e-05, .02, .14, .41,
463  & 1.48, 4.17, 7.92, 11.90, 19.03, 28.59, 42.36,
464  & 52.78, 62.86, 72.65, 85.61, 95.97,
465  & 96., 96., 96., 96.,
466 * total X-section , cteq1MS PTTHR=x GEV/C
467  & 0.000000e+00,
468  & 0.517461e-05, .02, .14, .42, 1.49, 4.14,
469  & 7.87, 11.93, 19.58, 30.67, 48.39, 63.08,
470  & 78.1, 93.28, 114.33, 132.24,
471  & 133., 133., 133., 133.,
472 * total X-section , cteq1ML PTTHR=x GEV/C
473  & 0.000000e+00,
474  & 0.717097e-05, .03, .19, .54, 1.91, 5.33, 10.11,
475  & 16.16, 24.21, 36.41, 54.21, 67.92, 81.44,
476  & 94.81,112.9, 127.63,
477  & 128., 128., 128., 128.,
478 * total X-section , cteq1D PTTHR=x GEV/C
479  & 0.000000e+00,
480  & 0.761464e-05, .02, .17, .47, 1.56, 4.19,
481  & 7.76, 11.48, 18.11, 26.97, 39.82, 49.86, 59.35,
482  & 68.88, 81.65, 91.94,
483  & 92., 92., 92., 92.,
484 * total X-section , cteq1L PTTHR=x GEV/C
485  & .000000e+00,
486  & .620779e-05, .02, .12, .34, 1.19, 3.27,
487  & 6.16, 9.27, 14.99, 23.2, 36.85, 49.45,
488  & 64.43, 82.38, 112.06, 140.36,
489  & 141., 141., 141., 141.,
490 * total X-section , GRV94LO AK=1. PTTHR=x GEV/C
491  & .000000e+00,
492  & .620779e-05, .01, .05, .14, 0.55, 1.87,
493  & 4.29, 7.49, 14.81, 27.8, 55.99, 77.49,
494  & 105.98,138.48, 189.33, 236.37,
495  & 294., 395., 496., 629.,
496 * total X-section , GRV94LO AK=2. PTTHR=x GEV/C
497  & .000000e+00,
498  & .620779e-05, .01, .10, .31, 1.16, 3.76,
499  & 8.31, 14.16, 27.11, 49.3, 90.93,129.77,
500  & 174.16,223.83, 300.20, 370.00,
501  & 455., 600., 746., 936.,
502 * total X-section , CTEC96 AK=2. PTTHR=x GEV/C
503  & .000000e+00,
504  & .620779e-05, .01, .08, .27, 1.17, 4.15,
505  & 9.60, 16.75, 32.88, 61.1,125.98,169.87,
506  & 233.75,308.22, 426.95, 537.90,
507  & 673., 898., 1112., 1379./
508 *******************************************************************
509 *
510  IF( abs(ptthr-three).LT.epsil ) THEN
511  WRITE(6,*) ' ERROR RDXSEC: invalid pdf No. ',istruf
512  stop
513  ELSEIF( abs(ptthr-two).LT.epsil ) THEN
514  WRITE(6,*) ' ERROR RDXSEC: invalid pdf No. ',istruf
515  stop
516  ELSEIF( istrut.EQ.1 ) THEN
517  WRITE(6,*) ' ERROR RDXSEC: invalid pdf No. ',istruf
518  stop
519  ELSEIF( istrut.EQ.2 ) THEN
520  IF( (istruf.GE.9).AND.(istruf.LE.20) ) THEN
521  WRITE(6,*) ' ERROR RDXSEC: invalid pdf No. ',istruf
522  stop
523  ELSEIF( (istruf.GE.21).AND.(istruf.LE.23) ) THEN
524  DO 311 i=1,21
525  nxs = 21*(istruf-15)+i
526  xsec(i)=xs21(nxs)
527  311 CONTINUE
528  ELSE
529  WRITE(6,*) ' ERROR RDXSEC: invalid pdf No. ',istruf
530  stop
531  ENDIF
532  ELSE
533  WRITE(6,*) ' ERROR RDXSEC: PTCUT ',ptthr,' not supported ***'
534  stop
535  ENDIF
536 C
537  RETURN
538  END
539 *
540 *
541 ************************************************************************
542 *
543  BLOCK DATA pomen
544  IMPLICIT DOUBLE PRECISION(a-h,o-z)
545 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20)
546  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
547  DATA poen/20.d0,50.d0,100.d0,200.d0,500.d0,
548  * 1000.d0,1500.d0,
549  * 2000.d0,3000.d0,4000.d0,6000.d0,8000.d0,10000.d0,
550  *15000.d0,20000.d0,30000.d0,40000.d0,60000.d0,
551  *80000.d0,100000.d0,150000.d0,200000.d0,300000.d0
552  *,400000.d0,600000.d0,800000.d0,1000000.d0,2000000.d0/
553  DATA poen1/5.d0,30.d0,70.d0,150.d0,300.d0,
554  * 700.d0,1200.d0,1700.d0,
555  * 2500.d0,3500.d0,5000.d0,7000.d0,9000.d0,
556  *12000.d0,17000.d0,25000.d0,35000.d0,50000.d0,
557  *70000.d0,90000.d0,120000.d0,170000.d0,250000.d0,
558  *250000.d0,500000.d0,700000.d0,900000.d0,1500000.d0/
559  DATA poen2/30.d0,70.d0,150.d0,300.d0,
560  * 700.d0,1200.d0,1700.d0,2500.d0,
561  * 3500.d0,5000.d0,7000.d0,9000.d0,12000.d0,
562  *17000.d0,25000.d0,35000.d0,50000.d0,70000.d0,
563  *90000.d0,120000.d0,170000.d0,250000.d0,350000.d0,
564  *500000.d0,700000.d0,900000.d0,1500000.d0,3000000.d0/
565  DATA nestep/28/
566  END
567 ************************************************************************
568 ************************************************************************
569 ************************************************************************
570 *
571  SUBROUTINE prblm2(ECM)
572  IMPLICIT DOUBLE PRECISION(a-h,o-z)
573 C
574 C Routine to call QRBLM2 for NESTEP energies
575 * j.r.3/94
576 C-----------------------------------------------------------------------
577 C$
578 C$ If IPOMTA=1 file INIDAT already exists : the content of commons
579 C$ POMENE and POLMN1 will be read in (Logic unit is IUNIT=37)
580 C$ If IPOMTA=0 file INIDAT does not exist : the content of COMMON
581 C$ POMENE will be taken from BLOCKDATA POMEN or, when interfaced to
582 C$ HEMAS, directly from hemas input file INPFIL. Then, the content of
583 C$ POMENE and POLMN1 will be written in INIDAT (Logic unit is IUNIT=37)
584 C$ C.Forti 18-nov-94
585 C$
586 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20),NESTEP
587  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
588  COMMON /pomtab/ipomta
589 C$
590  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
591  parameter(mxpa50=250,mxpa51=mxpa50+1)
592  parameter(mxpu50=100,mxpu51=mxpu50+1)
593 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
594  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
595  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
596  COMMON /polmn1/ plmnee(0:mxpa25,0:mxpu50,0:mxpa13,28)
597  DATA iunit/37/
598 C
599  CHARACTER*80 inidat
600 C
601 C$ Statement : setenv INIDAT /nfs/hpmac1/macro02/dpmjet/POMTAB_06.DAT
602 C$ in the job file is required
603 C
604 C CALL GETENV('INIDAT',INIDAT)
605 CC OPEN(IUNIT,FILE=INIDAT,STATUS='UNKNOWN',ERR=99)
606 C OPEN(UNIT=IUNIT,FILE='pomtab.dat'
607 C * ,STATUS='UNKNOWN',ERR=99)
608 Cvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
609  CHARACTER dpmdirname*400, pomfilename*400
610  CALL getenv('G4DPMJET2_5DATA',dpmdirname)
611  pomfilename = dpmdirname(1:len_trim(dpmdirname))//
612  &'/pomtab.dat'
613  OPEN(unit=iunit,file=pomfilename,status='OLD',err=99)
614 C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
615 C
616  IF (ipomta.EQ.0) THEN
617  DO 1 ii=1,nestep
618  energy = poen(ii)
619  CALL qrblm2(energy)
620  DO 10 jj=0,mxpa25
621  DO 10 kk=0,mxpu50
622  DO 10 ll=0,mxpa13
623  plmnee(jj,kk,ll,ii)=plmncu(jj,kk,ll)
624  10 CONTINUE
625  1 CONTINUE
626 C
627  WRITE(iunit,7102) nestep
628  DO 31 ii=1,nestep
629  WRITE(iunit,7101) poen(ii), poen1(ii), poen2(ii)
630  WRITE(iunit,101)(((plmnee(jj,kk,ll,ii),ll=0,mxpa13),
631  * kk=0,mxpu50),jj=0,mxpa25)
632  31 CONTINUE
633  7102 FORMAT(i13)
634  7101 FORMAT(3e13.5)
635  101 FORMAT(8e13.5)
636 C$
637  ELSEIF (ipomta.EQ.1)THEN
638  READ(iunit,7102) nestep
639  WRITE(6,7102) nestep
640  DO 11 ii=1,nestep
641  READ(iunit,7101) poen(ii), poen1(ii), poen2(ii)
642  WRITE(6,7101) poen(ii), poen1(ii), poen2(ii)
643  READ(iunit,101)(((plmnee(jj,kk,ll,ii),ll=0,mxpa13),
644  * kk=0,mxpu50),jj=0,mxpa25)
645 C WRITE(6,101)(((PLMNEE(JJ,KK,LL,II),LL=0,MXPA13),
646 C * KK=0,MXPU50),JJ=0,MXPA25)
647  11 CONTINUE
648  ENDIF
649  CLOSE(iunit)
650  RETURN
651  99 CONTINUE
652  WRITE(6,'(A)')'Error in PRBLM2 : file pomtab.dat ERROR'
653  CLOSE(iunit)
654  stop
655  END
656 
657 ************************************************************************
658 *
659  SUBROUTINE qrblm2(ECM)
660 * * input:
661 * ECM
662 * output:
663 * PLMN, PLMNCUmmulative, AVSOFN, AVHRDN,SIGDD/D/QEL/EL, PSOFT
664 *
665 *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
666 *
667 * Probabilities for L soft cut, M hard cut pomerons
668 * and N someplace cut trippel pomerons
669 *
670 * Aurenche Maire's version of PRBLM
671 * modified from A.M. to include calculation of x-sections
672 * modified to get:
673 * version with high mass diffraction (Y's and PHI's)
674 * and on both sides 2 kanal low mass diffraction
675 * *** OPTION 21.2.90 FB
676 *----------------------------------------------------------------------
677 *
678  IMPLICIT DOUBLE PRECISION(a-h,o-z)
679  SAVE
680  parameter( zero=0.d0, one=1.d0)
681  parameter(conv=0.38935d0)
682  parameter(pi=3.141592654d0)
683  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
684 * PARAMETRIZATION FOR PTMIN= 3. GEV
685  parameter(mxpa50=250,mxpa51=mxpa50+1)
686 * PARAMETRIZATION FOR PTMIN= 2. GEV
687 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
688  parameter(mxpa96=96)
689 C PARAMETER (MXPA96=480)
690  LOGICAL lsqrt
691 C PARAMETER (MXLMN=5,LSQRT=.false.)
692  parameter(mxlmn=5,lsqrt=.true.)
693  DOUBLE PRECISION dtiny
694 C PARAMETER (TINY = 1.2D-38,DTINY=1.D-70,TIN=1.D-22,TINEXP = -300.D0)
695 C PARAMETER (TINY = 1.2D-38,DTINY=1.D-300,TIN=1.D-22,TINEXP =
696 C -700.D0)
697  parameter(tiny=1.2d-38,dtiny=1.d-70,tin=1.d-22,tinexp=-700.d0)
698 C PARAMETER (TINY = 1.D-38,DTINY=tiny,TIN=1.D-22,TINEXP = -300.D0)
699 * in older version used:
700  parameter(tinyex = -48.d0)
701 * --- -- - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - -
702 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
703  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
704  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
705  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
706  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
707 * *** /POMPAR/*/SIGMA/*/POMTYP/ used only in SIGMAPOM-routines (->POMDI)
708 * (LMAX,MMAX,NMAX (max number of soft/hard/trippel pomerons)
709  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
710  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
711  * sigloo,zloo
712  common/pompar/alfa,alfap,a,c,ak
713  COMMON /singdi/silmsd,sigdi
714 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
715  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
716  COMMON /alala/alalam
717  CHARACTER*80 title
718  CHARACTER*8 projty,targty
719 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
720 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
721  COMMON /user1/title,projty,targty
722  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
723 * --- -- - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - -
724  DOUBLE PRECISION sig,sigp,sigm,sign,sigo
725  dimension sig(0:mxpa25,0:mxpa50,0:mxpa13),
726  &sigp(0:mxpa25,0:mxpa50,0:mxpa13),sigm(0:mxpa25,0:mxpa50,0:mxpa13),
727  &sign(0:mxpa25,0:mxpa50,0:mxpa13),sigo(0:mxpa25,0:mxpa50,0:mxpa13)
728  dimension xpnt(mxpa96),wght(mxpa96),
729  &ssoft(0:mxpa25),shard(0:mxpa50),strpl(0:mxpa25)
730 C - - required MXPA25 > NMAX - -
731  dimension fak(0:mxpa13),cmbin(0:mxpa13,0:mxpa13)
732  double precision
733  & expsop,expsoh,exmsop,exmsoh,exnsop,exnsoh,exosop,exosoh,
734  & exphap,exphah,exmhap,exmhah,exnhap,exnhah,exohap,exohah,
735  & exptrp,exptrh,exmtrp,exmtrh,exntrp,exntrh,exotrp,exotrh,
736  & explop,exploh,exmlop,exmloh,exnlop,exnloh,exolop,exoloh,
737  & expexh,exmexh,exnexh,exoexh,expexp,exmexp,exnexp,exoexp
738  DOUBLE PRECISION fapsof,famsof,fansof,faosof,
739  & faphar,famhar,fanhar,faohar,
740  & faptrp,famtrp,fantrp,faotrp,
741  & faploo,famloo,fanloo,faoloo
742  DOUBLE PRECISION denom,denomi,xpntk,wghtk,rmxlmn
743  & ,sigsum,siginl,sighri
744 *
745 *---------------------------------------------------------------------------
746 *
747 * externe ICON option to internes NMAX=1,2,free
748  IF(icon/10.EQ.4) nmax=2
749  IF(icon/10.EQ.5) nmax=1
750 *
751 * for safty
752  IF( nmax.GT.mxpa13) THEN
753  WRITE(6,*)' arrays limit NMAX set to' , mxpa13
754  nmax=mxpa13
755  ENDIF
756  IF( mmax.GT.mxpa50) THEN
757  WRITE(6,*)' arrays limit MMAX set to' , mxpa50
758  mmax=mxpa50
759  ENDIF
760  IF( lmax.GT.mxpa25) THEN
761  WRITE(6,*)' arrays limit LMAX set to' , mxpa25
762  lmax=mxpa25
763  ENDIF
764 *
765  lmaxi = lmax
766  mmaxi = mmax
767  IF( nmax.GE.3)THEN
768  nmaxi = nmax
769  nnmaxi=(mxpa13-nmaxi)/(1+nmaxi)
770 C aim: MXPA13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
771  nlmaxi=0
772  ELSEIF( nmax.EQ.2)THEN
773  nmaxi=1
774  nnmaxi=1
775  nlmaxi=1
776  ELSEIF( nmax.EQ.1)THEN
777  nmaxi=1
778  nnmaxi=0
779  nlmaxi=1
780  ELSEIF( nmax.LE.0)THEN
781  nmaxi=0
782  nnmaxi=0
783  nlmaxi=0
784  ENDIF
785 *
786 *
787  lentry=0
788 *
789  goto 111
790 *
791 *----------------------------------------------------------------------
792 *
793  entry sigma2(ecm)
794 *
795 *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
796 *
797  lentry=1
798 * externe ICON option to internes NMAX=1,2,free
799  IF(icon/10.EQ.4) nmax=2
800  IF(icon/10.EQ.5) nmax=1
801 * we drop L,M..dependent quantities, the rest is integrated for L=1,rest=0
802  lmaxi = 1
803  mmaxi = 0
804  nmaxi = 0
805  nnmaxi= 0
806  nlmaxi= 0
807 *
808 *----------------------------------------------------------------------
809 *
810 * *** calculate the B-space integral for L/M soft/hard cut-pom.
811 * *** NPNT is the number of integration points of B-space integ.
812 *
813  111 sigtot=0.0
814  siginl=0.0
815  sigel =0.0
816  sigele=0.0
817  sigd =0.0
818  sigdd =0.0
819  sigdi =0.0
820  sigddi=0.0
821  sigine=0.0
822  sihmdd=0.0
823  sighmd=0.
824  siglmd=0.
825  sighin=0.
826  sigin=0.
827  sigsin=0.0
828  sighri=0.0
829  sigsum=0.0
830  sigsme=0.0
831  silmsd=0.0
832  silmdd=0.
833  slhmdd=0.
834  DO 10 l=0,lmaxi
835  ssoft(l)=0.
836  strpl(l)=0.
837  DO 10 m=0,mmaxi
838  shard(m)=0.
839  DO 10 n=0, mxpa13
840  sig(l,m,n)=0.
841  sigp(l,m,n)=0.
842  sigm(l,m,n)=0.
843  sign(l,m,n)=0.
844  sigo(l,m,n)=0.
845  plmn(l,m,n)=0.
846  plmncu(l,m,n)=0.
847 10 CONTINUE
848 *
849 * get bare X-sections
850  s=ecm**2
851  CALL sigshd(ecm)
852 *
853  IF(alalam.LE.1.d-2) THEN
854  alam=0.60
855  ELSE
856  alam=alalam
857  ENDIF
858 *
859 * prepare integration:
860  pi4 = 4.*pi
861 *
862  IF(ecm.LT.2000.d0)THEN
863  npnt=96
864  CALL gset(zero,one,npnt,xpnt,wght)
865  ELSE
866  npnt=mxpa96
867  CALL gset(zero,one,npnt,xpnt,wght)
868  ENDIF
869 *
870 * as low mass diffraction extra, high mass reduced:
871  redu = 1.0
872  IF(ioutpo.GE.0) WRITE (6,*) ' ALAM,REDU= ',alam,redu
873 *
874 * --here the versions started--
875 * prepare factors to enter sum:
876 * notation: Z(HARd/SOFt/LOOp)(Plus/Minus/N=mixed/O=mixed)
877 *
878  zharp=(1.+alam)**2*zhar
879  zsofp=(1.+alam)**2*zsof
880  zloop=(1.+alam)**2*zloo * redu
881  zharm=(1.-alam)**2*zhar
882  zsofm=(1.-alam)**2*zsof
883  zloom=(1.-alam)**2*zloo * redu
884  zharn=(1.-alam**2)*zhar
885  zsofn=(1.-alam**2)*zsof
886  zloon=(1.-alam**2)*zloo * redu
887  zharo=(1.-alam**2)*zhar
888  zsofo=(1.-alam**2)*zsof
889  zlooo=(1.-alam**2)*zloo * redu
890 * one more factor at the top or at the bottom:
891  ztrpp=(1.+alam)**3*ztrp * redu
892  ztrpm=(1.-alam)**3*ztrp * redu
893  ztrpn=(1.-alam**2)*(1.+alam)*ztrp * redu
894  ztrpo=(1.-alam**2)*(1.-alam)*ztrp * redu
895 *
896 * begin M,N,L,LL loop
897 *
898  DO 720 l=0,lmaxi
899  IF(l.EQ.0) THEN
900  fapsof=1.
901  famsof=1.
902  fansof=1.
903  faosof=1.
904  ELSEIF(lsqrt) THEN
905  fapsof=fapsof* sqrt( zsofp/float(l))
906  famsof=famsof* sqrt( zsofm/float(l))
907  fansof=fansof* sqrt( zsofn/float(l))
908  faosof=faosof* sqrt( zsofo/float(l))
909  IF ( fapsof .LT.dtiny ) fapsof=0.
910  IF ( famsof .LT.dtiny ) famsof=0.
911  IF ( fansof .LT.dtiny ) fansof=0.
912  IF ( faosof .LT.dtiny ) faosof=0.
913  ELSEIF(.NOT.lsqrt) THEN
914  fapsof=fapsof*zsofp/float(l)
915  famsof=famsof*zsofm/float(l)
916  fansof=fansof*zsofn/float(l)
917  faosof=faosof*zsofo/float(l)
918  IF (fapsof.LT.dtiny ) fapsof=0.
919  IF (famsof.LT.dtiny ) famsof=0.
920  IF (fansof.LT.dtiny ) fansof=0.
921  IF (faosof.LT.dtiny ) faosof=0.
922  ENDIF
923  DO 730 m=0,mmaxi
924  IF(m.EQ.0) THEN
925  faphar=1.
926  famhar=1.
927  fanhar=1.
928  faohar=1.
929  ELSEIF(lsqrt) THEN
930 C WRITE(6,*)FAPHAR,ZHARP/FLOAT(M),FAMHAR,ZHARM/FLOAT(M)
931  faphar=faphar* sqrt( zharp/float(m) )
932  famhar=famhar* sqrt( zharm/float(m) )
933  fanhar=fanhar* sqrt( zharn/float(m) )
934  faohar=faohar* sqrt( zharo/float(m) )
935  IF ( fapsof*faphar .LT.dtiny ) faphar=0.
936  IF ( famsof*famhar .LT.dtiny ) famhar=0.
937  IF ( fansof*fanhar .LT.dtiny ) fanhar=0.
938  IF ( faosof*faohar .LT.dtiny ) faohar=0.
939  ELSEIF(.NOT.lsqrt) THEN
940  faphar=faphar*zharp/float(m)
941  famhar=famhar*zharm/float(m)
942  fanhar=fanhar*zharn/float(m)
943  faohar=faohar*zharo/float(m)
944  IF (fapsof*faphar.LT.dtiny ) faphar=0.
945  IF (famsof*famhar.LT.dtiny ) famhar=0.
946  IF (fansof*fanhar.LT.dtiny ) fanhar=0.
947  IF (faosof*faohar.LT.dtiny ) faohar=0.
948  ENDIF
949  DO 740 n=0,nmaxi
950  IF( n.EQ.0) THEN
951  faptrp=1.
952  famtrp=1.
953  fantrp=1.
954  faotrp=1.
955  ELSEIF(lsqrt) THEN
956  faptrp=-faptrp* sqrt( ztrpp/float(n) )
957  famtrp=-famtrp* sqrt( ztrpm/float(n) )
958  fantrp=-fantrp* sqrt( ztrpn/float(n) )
959  faotrp=-faotrp* sqrt( ztrpo/float(n) )
960  IF (abs(faptrp*fapsof*faphar).LT.dtiny ) faptrp=0.
961  IF (abs(famtrp*famsof*famhar).LT.dtiny ) famtrp=0.
962  IF (abs(fantrp*fansof*fanhar).LT.dtiny ) fantrp=0.
963  IF (abs(faotrp*faosof*faohar).LT.dtiny ) faotrp=0.
964  ELSEIF(.NOT.lsqrt) THEN
965  faptrp=-faptrp*ztrpp/float(n)
966  famtrp=-famtrp*ztrpm/float(n)
967  fantrp=-fantrp*ztrpn/float(n)
968  faotrp=-faotrp*ztrpo/float(n)
969  IF (abs(faptrp*fapsof*faphar).LT.dtiny ) faptrp=0.
970  IF (abs(famtrp*famsof*famhar).LT.dtiny ) famtrp=0.
971  IF (abs(fantrp*fansof*fanhar).LT.dtiny ) fantrp=0.
972  IF (abs(faotrp*faosof*faohar).LT.dtiny ) faotrp=0.
973  ENDIF
974  DO 750 nn=0,nnmaxi
975 * for compatibility no new subscript is introduced in some arrays:
976  nnn=n+(nmaxi+1)*nn
977 * if only first order option jump out of second order contr.:
978  IF( nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1 ) go to 750
979  IF(nn.EQ.0) THEN
980  faploo=1.
981  famloo=1.
982  fanloo=1.
983  faoloo=1.
984  ELSEIF(lsqrt) THEN
985  faploo=-faploo* sqrt( zloop/float(nn))
986  famloo=-famloo* sqrt( zloom/float(nn))
987  fanloo=-fanloo* sqrt( zloon/float(nn))
988  faoloo=-faoloo* sqrt( zlooo/float(nn))
989  IF(abs(faploo*faptrp*fapsof*faphar).LT.dtiny )faploo=0.
990  IF(abs(famloo*famtrp*famsof*famhar).LT.dtiny )famloo=0.
991  IF(abs(fanloo*fantrp*fansof*fanhar).LT.dtiny )fanloo=0.
992  IF(abs(faoloo*faotrp*faosof*faohar).LT.dtiny )faoloo=0.
993  ELSEIF(.NOT.lsqrt) THEN
994  faploo=-faploo*zloop/float(nn)
995  famloo=-famloo*zloom/float(nn)
996  fanloo=-fanloo*zloon/float(nn)
997  faoloo=-faoloo*zlooo/float(nn)
998  IF(abs(faploo*faptrp*fapsof*faphar).LT.dtiny )faploo=0.
999  IF(abs(famloo*famtrp*famsof*famhar).LT.dtiny )famloo=0.
1000  IF(abs(fanloo*fantrp*fansof*fanhar).LT.dtiny )fanloo=0.
1001  IF(abs(faoloo*faotrp*faosof*faohar).LT.dtiny )faoloo=0.
1002  ENDIF
1003 *
1004 * elastic processes are not generated
1005  IF(l.EQ.0.AND.m.EQ.0.AND.n.EQ.0.AND.nn.EQ.0) go to 750
1006 *
1007  denom=dble(m)/dble(bh)+dble(l)/dble(bs)+dble(n)/dble(bt)
1008  & +dble(nn)/dble(bt)
1009 *
1010  DO 735 k=1,npnt
1011 *
1012 C change intergration for large L+M+N+NN?
1013  IF ( (m+l+n+nn) .LE. mxlmn ) THEN
1014  xpntk=dble(xpnt(k))
1015  wghtk=dble(wght(k))
1016  denomi=denom
1017  ELSE
1018  rmxlmn = dble(m+l+n+nn) /dble(mxlmn)
1019  xpntk=dble(xpnt(k))
1020  wghtk= dble(wght(k)) * xpntk**(rmxlmn-1.)
1021  denomi= denom / rmxlmn
1022  ENDIF
1023 *
1024  exposp=-zsofp*xpntk**(1./(denomi*dble(bs)))
1025  exposm=-zsofm*xpntk**(1./(denomi*dble(bs)))
1026  exposn=-zsofn*xpntk**(1./(denomi*dble(bs)))
1027  exposo=-zsofo*xpntk**(1./(denomi*dble(bs)))
1028 *
1029  expohp=-zharp*xpntk**(1./(denomi*dble(bh)))
1030  expohm=-zharm*xpntk**(1./(denomi*dble(bh)))
1031  expohn=-zharn*xpntk**(1./(denomi*dble(bh)))
1032  expoho=-zharo*xpntk**(1./(denomi*dble(bh)))
1033 *
1034  expotp=+ztrpp*xpntk**(1./(denomi*dble(bt)))
1035  expotm=+ztrpm*xpntk**(1./(denomi*dble(bt)))
1036  expotn=+ztrpn*xpntk**(1./(denomi*dble(bt)))
1037  expoto=+ztrpo*xpntk**(1./(denomi*dble(bt)))
1038 *
1039  expolp=+zloop*xpntk**(1./(denomi*dble(bt)))
1040  expolm=+zloom*xpntk**(1./(denomi*dble(bt)))
1041  expoln=+zloon*xpntk**(1./(denomi*dble(bt)))
1042  expolo=+zlooo*xpntk**(1./(denomi*dble(bt)))
1043 *
1044  IF(ioutpo.GE.7) THEN
1045  WRITE(6,*)
1046  * ' K=',k,' EXPOS/H=',exposp,expohp,' DENOMI/BH=',denomi,bh
1047  WRITE(6,*)
1048  * ' K=',k,' EXPOS/H=',exposm,expohm,' DENOMI/BH=',denomi,bh
1049  WRITE(6,*)
1050  * ' K=',k,' EXPOS/H=',exposn,expohn,' DENOMI/BH=',denomi,bh
1051  WRITE(6,*)
1052  * ' K=',k,'XPNT=',xpntk,'WGHT=',wghtk,'DENO=',denomi
1053  ENDIF
1054 *
1055 * notation:
1056 * EX(P=+/M=-/N=mixed/O=mixed)(EX=all/HArd,TRippel,LOop)(P/Half)
1057 *
1058  IF( exposp .GT. tinexp) THEN
1059  expsoh=exp(0.5d00*exposp)
1060  exmsoh=exp(0.5d00*exposm)
1061  exnsoh=exp(0.5d00*exposn)
1062  exosoh=exp(0.5d00*exposo)
1063  ELSE
1064  expsoh=0.
1065  exmsoh=0.
1066  exnsoh=0.
1067  exosoh=0.
1068  ENDIF
1069  expsop=expsoh**2
1070  exmsop=exmsoh**2
1071  exnsop=exnsoh**2
1072  exosop=exosoh**2
1073 *
1074  IF( expohp .GT. tinexp) THEN
1075  exphah=exp(0.5d00*expohp)
1076  exmhah=exp(0.5d00*expohm)
1077  exnhah=exp(0.5d00*expohn)
1078  exohah=exp(0.5d00*expoho)
1079  ELSE
1080  exphah=0.
1081  exmhah=0.
1082  exnhah=0.
1083  exohah=0.
1084  ENDIF
1085  exphap=exphah**2
1086  exmhap=exmhah**2
1087  exnhap=exnhah**2
1088  exohap=exohah**2
1089 *
1090  IF( nmax.GE.3) THEN
1091  IF( expotp .GT. tinexp) THEN
1092  exptrh=exp(0.5d00*expotp)
1093  exmtrh=exp(0.5d00*expotm)
1094  exntrh=exp(0.5d00*expotn)
1095  exotrh=exp(0.5d00*expoto)
1096  ELSE
1097  exptrh=0.
1098  exmtrh=0.
1099  exntrh=0.
1100  exotrh=0.
1101  ENDIF
1102  exptrp= exptrh**2
1103  exmtrp= exmtrh**2
1104  exntrp= exntrh**2
1105  exotrp= exotrh**2
1106  ELSEIF( nmax.LE.2) THEN
1107  exptrh= 1 + 0.5*expotp
1108  exmtrh= 1 + 0.5*expotm
1109  exntrh= 1 + 0.5*expotn
1110  exotrh= 1 + 0.5*expoto
1111  exptrp= 1 + expotp
1112  exmtrp= 1 + expotm
1113  exntrp= 1 + expotn
1114  exotrp= 1 + expoto
1115  ENDIF
1116 *
1117  IF( nmax.GE.3) THEN
1118  IF( expolp .GT. tinexp) THEN
1119  exploh=exp(0.5d00*expolp)
1120  exmloh=exp(0.5d00*expolm)
1121  exnloh=exp(0.5d00*expoln)
1122  exoloh=exp(0.5d00*expolo)
1123  ELSE
1124  exploh=0.
1125  exmloh=0.
1126  exnloh=0.
1127  exoloh=0.
1128  ENDIF
1129  explop=exploh**2
1130  exmlop=exmloh**2
1131  exnlop=exnloh**2
1132  exolop=exoloh**2
1133  ELSEIF( nmax.EQ.2 ) THEN
1134  exploh= 1 + 0.5*expolp
1135  exmloh= 1 + 0.5*expolm
1136  exnloh= 1 + 0.5*expoln
1137  exoloh= 1 + 0.5*expolo
1138  explop= 1 + expolp
1139  exmlop= 1 + expolm
1140  exnlop= 1 + expoln
1141  exolop= 1 + expolo
1142  ELSEIF( nmax.LE.1 ) THEN
1143  exploh= 1
1144  exmloh= 1
1145  exnloh= 1
1146  exoloh= 1
1147  explop= 1
1148  exmlop= 1
1149  exnlop= 1
1150  exolop= 1
1151  ENDIF
1152 *
1153  expexh = expsoh *exphah *exptrh *exploh
1154  exmexh = exmsoh *exmhah *exmtrh *exmloh
1155  exnexh = exnsoh *exnhah *exntrh *exnloh
1156  exoexh = exosoh *exohah *exotrh *exoloh
1157  expexp = expsop *exphap *exptrp *explop
1158  exmexp = exmsop *exmhap *exmtrp *exmlop
1159  exnexp = exnsop *exnhap *exntrp *exnlop
1160  exoexp = exosop *exohap *exotrp *exolop
1161 *
1162  IF( ( nmax.LE.2 .AND. n.EQ.1 ) .OR.
1163  * ( nmax.EQ.2 .AND. nn.EQ.1 ) .OR.
1164  * nmax.EQ.0 ) THEN
1165  sigp(l,m,nnn)=sigp(l,m,nnn)+expsop *exphap *wghtk
1166  sigm(l,m,nnn)=sigm(l,m,nnn)+exmsop *exmhap *wghtk
1167  sign(l,m,nnn)=sign(l,m,nnn)+exnsop *exnhap *wghtk
1168  sigo(l,m,nnn)=sigo(l,m,nnn)+exosop *exohap *wghtk
1169  ELSE
1170  sigp(l,m,nnn)=sigp(l,m,nnn)+expexp*wghtk
1171  sigm(l,m,nnn)=sigm(l,m,nnn)+exmexp*wghtk
1172  sign(l,m,nnn)=sign(l,m,nnn)+exnexp*wghtk
1173  sigo(l,m,nnn)=sigo(l,m,nnn)+exoexp*wghtk
1174  ENDIF
1175 *
1176 * quantities without L,M,N,NN dependence considered ones
1177 * (when, chosen to get suitable weights)
1178  IF(l.EQ.1.AND.m.EQ.0.AND.n.EQ.0.AND.nn.EQ.0) THEN
1179 *
1180  IF ( (m+l+n+nn) .GT. mxlmn ) THEN
1181  WRITE(6,*)' MXLMN too low ' , mxlmn,m,l,n,nn
1182  RETURN
1183  ENDIF
1184  wghfac = wghtk/xpntk *pi4/denomi
1185  IF ( nmax.GE.3 ) THEN
1186  sigele = sigele + wghfac *
1187  * 0.0625*( 1.-expexh + 1.-exmexh
1188  * +1.-exnexh + 1.-exoexh )**2
1189 * low mass diffraction:
1190  silmsd = silmsd + wghfac *
1191  * 0.125*(expexh -exmexh)**2
1192  silmdd = silmdd + wghfac *
1193  * 0.0625*(expexh+exmexh-exnexh-exoexh)**2
1194  ELSEIF( nmax.LE.2 ) THEN
1195  sigele = sigele + wghfac *
1196  * 0.0625*( ( 1.-expexh + 1.-exmexh
1197  * +1.-exnexh + 1.-exoexh
1198 * subtract second order terms in each factor:
1199 C////CHANGED - TO + BRACKET UNNECESSARY
1200  * +(1.-exptrh)*(1-exploh) *expsoh *exphah
1201  * +(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1202  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1203  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah)**2
1204 * subtract second order terms of product:
1205  * - ( (2.-exptrh-exploh) *expsoh *exphah
1206  * +(2.-exmtrh-exmloh) *exmsoh *exmhah
1207  * +(2.-exntrh-exnloh) *exnsoh *exnhah
1208  * +(2.-exotrh-exoloh) *exosoh *exohah ) **2)
1209 * low mass diffraction:
1210  silmsd = silmsd + wghfac *
1211  * 0.125*( ( expexh -exmexh
1212 * subtract second order terms in each factor:
1213  * -(1.-exptrh)*(1-exploh) *expsoh*exphah
1214  * +(1.-exmtrh)*(1-exmloh) *exmsoh*exmhah )**2
1215 * subtract second order terms of product:
1216  * -( (2.-exptrh-exploh) *expsoh *exphah
1217  * -(2.-exmtrh-exmloh) *exmsoh*exmhah ) **2)
1218  silmdd = silmdd + wghfac *
1219  * 0.0625*( (expexh+exmexh-exnexh-exoexh
1220 * subtract second order terms in each factor:
1221  * -(1.-exptrh)*(1-exploh) *expsoh *exphah
1222  * -(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1223  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1224  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah)**2
1225 * subtract second order terms of product:
1226  * - ( (2.-exptrh-exploh) *expsoh *exphah
1227  * +(2.-exmtrh-exmloh) *exmsoh *exmhah
1228  * -(2.-exntrh-exnloh) *exnsoh *exnhah
1229  * -(2.-exotrh-exoloh) *exosoh *exohah ) **2)
1230  ENDIF
1231  IF( nmax.NE.2 ) THEN
1232  sigtot=sigtot+2.*wghfac*
1233  * 0.25*( 1.-expexh + 1.-exmexh +
1234  * 1.-exnexh + 1.-exoexh )
1235  sigine = sigine + wghfac *
1236  * 0.25*( 1.-expexp + 1.-exmexp +
1237  * 1.-exnexp + 1.-exoexp )
1238 * pure-soft-inelastic (hard scatterring is included as absorbtion)
1239  sigsin=sigsin+ wghfac *
1240  * 0.25*( (exphap-expexp)
1241  * +(exmhap-exmexp)
1242  * +(exnhap-exnexp)
1243  * +(exohap-exoexp) )
1244 * hard-inelastic (soft scatterring disregarded)
1245  sighin=sighin+ wghfac*
1246  * 0.25*( 1.-exphap + 1.-exmhap +
1247  * 1.-exnhap + 1.-exohap )
1248  ELSEIF( nmax.EQ.2 ) THEN
1249  sigtot=sigtot+2.*wghfac*
1250  * 0.25*( 1.-expexh + 1.-exmexh +
1251  * 1.-exnexh + 1.-exoexh
1252 * subtract second order terms in .TRH and LOH:
1253 C/////CHANGED - TO +
1254  * +(1.-exptrh)*(1-exploh) *expsoh *exphah
1255  * +(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1256  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1257  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah )
1258  sigine = sigine + wghfac *
1259  * 0.25*( 1.-expexp + 1.-exmexp +
1260  * 1.-exnexp + 1.-exoexp
1261 * subtract second order terms in .TRP and LOP:
1262 C/////CHANGED - TO +
1263  * +(1.-exptrp)*(1-explop) *expsop *exphap
1264  * +(1.-exmtrp)*(1-exmlop) *exmsop *exmhap
1265  * +(1.-exntrp)*(1-exnlop) *exnsop *exnhap
1266  * +(1.-exotrp)*(1-exolop) *exosop *exohap )
1267 * pure-soft-inelastic (hard scatterring is included as absorbtion)
1268  sigsin=sigsin+ wghfac *
1269  * 0.25*( (exphap-expexp)
1270  * +(exmhap-exmexp)
1271  * +(exnhap-exnexp)
1272  * +(exohap-exoexp)
1273 * subtract second order terms of 2nd column:
1274  * +(1.-exptrp)*(1-explop) *expsop *exphap
1275  * +(1.-exmtrp)*(1-exmlop) *exmsop *exmhap
1276  * +(1.-exntrp)*(1-exnlop) *exnsop *exnhap
1277  * +(1.-exotrp)*(1-exolop) *exosop *exohap)
1278 * hard-inelastic (soft scatterring disregarded)
1279  sighin=sighin+ wghfac*
1280  * 0.25*( 1.-exphap + 1.-exmhap +
1281  * 1.-exnhap + 1.-exohap )
1282  ENDIF
1283 * high mass diffraction (sep.low mass diffr. arbitrary)
1284 * (naive 1/EX?TRP -> EX?TR as selected cut counts -1)
1285  IF( nmax.GE.3 ) THEN
1286  sighmd=sighmd + wghfac *
1287  * 0.25*( (exptrp-1.)*expexp
1288  * +(exmtrp-1.)*exmexp
1289  * +(exntrp-1.)*exnexp
1290  * +(exotrp-1.)*exoexp)
1291  ELSE
1292  sighmd=sighmd + wghfac *
1293  * 0.25*( expotp * expsop*exphap
1294  * +expotm * exmsop*exmhap
1295  * +expotn * exnsop*exnhap
1296  * +expoto * exosop*exohap )
1297  ENDIF
1298  IF( nmax.GE.3 ) THEN
1299  sihmdd=sihmdd + wghfac *
1300  * 0.25*( (explop-1.)*expexp
1301  * +(exmlop-1.)*exmexp
1302  * +(exnlop-1.)*exnexp
1303  * +(exolop-1.)*exoexp)
1304  ELSEIF (nmax.EQ.2 ) THEN
1305  sihmdd=sihmdd + wghfac *
1306  * 0.25*( expolp * expsop*exphap
1307  * +expolm * exmsop*exmhap
1308  * +expoln * exnsop*exnhap
1309  * +expolo * exosop*exohap )
1310 * no action:
1311 * ELSEIF (NMAX.LE.1 ) THEN
1312 * SIHMDD=SIHMDD + WGHFAC *
1313 * * 0.25*( 0. * EXPSOP*EXPHAP
1314 * * + 0. * EXMSOP*EXMHAP
1315 * * + 0. * EXNSOP*EXNHAP
1316 * * + 0. * EXOSOP*EXOHAP )
1317  ENDIF
1318  ENDIF
1319 * ending non L,M,N,NN depending part
1320 *
1321 735 CONTINUE
1322 * ending impact-integral loop
1323 *
1324  IF(abs(faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)).LT.dtiny)
1325  & THEN
1326  sigp(l,m,nnn)=0.
1327  ELSEIF(lsqrt) THEN
1328  sigp(l,m,nnn)=faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)
1329  * * abs(faphar*fapsof*faptrp*faploo)/denomi*pi4
1330  ELSEIF(.NOT.lsqrt) THEN
1331  sigp(l,m,nnn)=faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)
1332  * /denomi*pi4
1333  ENDIF
1334  IF(abs(famhar*famsof*famtrp*famloo*sigm(l,m,nnn)).LT.dtiny)
1335  & THEN
1336  sigm(l,m,nnn)=0.
1337  ELSEIF(lsqrt) THEN
1338  sigm(l,m,nnn)=famhar*famsof*famtrp*famloo*sigm(l,m,nnn)
1339  * * abs( famhar*famsof*famtrp*famloo)/denomi*pi4
1340  ELSEIF(.NOT.lsqrt) THEN
1341  sigm(l,m,nnn)=famhar*famsof*famtrp*famloo*sigm(l,m,nnn)
1342  * /denomi*pi4
1343  ENDIF
1344  IF(abs(fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)).LT.dtiny)
1345  & THEN
1346  sign(l,m,nnn)=0.
1347  ELSEIF(lsqrt) THEN
1348  sign(l,m,nnn)=fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)
1349  * * abs( fanhar*fansof*fantrp*fanloo)/denomi*pi4
1350  ELSEIF(.NOT.lsqrt) THEN
1351  sign(l,m,nnn)=fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)
1352  * /denomi*pi4
1353  ENDIF
1354  IF(abs(faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)).LT.dtiny)
1355  & THEN
1356  sigo(l,m,nnn)=0.
1357  ELSEIF(lsqrt) THEN
1358  sigo(l,m,nnn)=faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)
1359  * * abs( faohar*faosof*faotrp*faoloo/denomi)*pi4
1360  ELSEIF(.NOT.lsqrt) THEN
1361  sigo(l,m,nnn)=faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)
1362  * /denomi*pi4
1363  ENDIF
1364 *
1365 750 CONTINUE
1366 740 CONTINUE
1367 730 CONTINUE
1368 720 CONTINUE
1369 *
1370 * *** summing up the three contributions
1371 *
1372  nnnmax=nmaxi+(nmaxi+1)*nnmaxi
1373  DO 820 l=0,lmaxi
1374  DO 830 m=0,mmaxi
1375  DO 830 nnn=0,nnnmax
1376  sig(l,m,nnn)=(sigp(l,m,nnn)+sigm(l,m,nnn)+
1377  * sign(l,m,nnn)+sigo(l,m,nnn) )/4.
1378 830 CONTINUE
1379 820 CONTINUE
1380 *
1381 *
1382 * *** calculate summed quantities for print out
1383 *
1384  DO 3 l=0,lmaxi
1385  DO 4 m=0,mmaxi
1386  DO 4 n=0,nmaxi
1387  DO 4 nn=0,nnmaxi
1388  IF( nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1 ) go to 4
1389  nnn=n+(nmaxi+1)*nn
1390  sigsum=sigsum + sig(l,m,nnn)
1391 * for options outlawing hard without soft:
1392  IF(m.EQ.0.OR.l.GE.1) sigsme=sigsme + sig(l,m,nnn)
1393  shard(m)=shard(m)+sig(l,m,nnn)
1394  ssoft(l)=ssoft(l)+sig(l,m,nnn)
1395  strpl(n)=strpl(n)+sig(l,m,nnn)
1396  siginl = siginl + sig(l,m,nnn)
1397  IF(m.GE.1) sighri = sighri + sig(l,m,nnn)
1398  IF(l.EQ.0.AND.m.EQ.0.AND.nn.EQ.0.AND.n.GE.1) THEN
1399  sigdi = sigdi + (-1)**n*sig(l,m,nnn)
1400  ELSEIF(l.EQ.0.AND.m.EQ.0.AND.n.EQ.0.AND.nn.GE.1) THEN
1401  sigddi= sigddi + (-1)**nn*sig(l,m,nnn)
1402  ENDIF
1403  4 CONTINUE
1404  3 CONTINUE
1405 * elastic processes were not generated, L,M,N,NN=0 no problem
1406 *
1407  siglmd=silmsd+silmdd
1408  sithmd=sighmd+sihmdd
1409  sigd = siglmd + sithmd
1410  slhmdd = sqrt(abs(silmdd*sihmdd))
1411  sigdd= silmdd + sihmdd + slhmdd
1412  sigin=sigine+siglmd
1413  sigel=sigtot-sigin
1414 *
1415 * *** print out
1416 *
1417  IF(lentry.EQ.1.AND.ioutpo.LE.1) RETURN
1418 *
1419  WRITE(6,*)' '
1420  WRITE(6,*)' --- properties of events ---'
1421  WRITE (6,102)
1422  WRITE(6,*)' Energy=',ecm
1423  WRITE (6,102)
1424  WRITE(6,*)' max.contributing soft/hard/diffr./doubl.diffr. cuts'
1425  WRITE(6,*)' LMAXI= MMAXI= NMAXI= NNMAXI='
1426  WRITE(6,'(15X,4I9)') lmaxi,mmaxi,nmaxi,nnmaxi
1427  WRITE(6,*)' methode used: '
1428  WRITE(6,*)' ISIG= ICON= IPIM= '
1429  WRITE(6,'(15X,3I9)') isig,icon,ipim
1430  WRITE (6,102)
1431  WRITE(6,*)' --- bare cross section and eikonal constants ---'
1432 C COMMON/POMPAR/ALFA,ALFAP,A,C,AK
1433 C COMMON /SIGMA/SIGSOF,BS,ZSOF,SIGHAR,BH,ZHAR,SIGTRP,BT,ZTRP,
1434 C * SIGLOO,ZLOO
1435  WRITE(6,*)' ALFA =',alfa,' ALFAP =',alfap,' A =',a
1436  WRITE(6,*)' C =',c,' AK =',ak
1437  WRITE(6,*)' ALALAM =',alalam
1438  WRITE (6,102)
1439  WRITE(6,*)' SIGSOF=',sigsof,' BS=',bs,' ZSOF=',zsof
1440  WRITE(6,*)' SIGHAR=',sighar,' BH=',bh,' ZHAR=',zhar
1441  WRITE(6,*)' SIGTRP=',sigtrp,' BT=',bt,' ZTRP=',ztrp
1442  WRITE(6,*)' SIGLOO=',sigloo,' BT=',bt,' ZLOO=',zloo
1443  WRITE (6,102)
1444  WRITE(6,*)' --- observable cross sections ---'
1445  WRITE (6,102)
1446  WRITE(6,*)' TOTAL X-SECTION = ',sigtot
1447  WRITE(6,*)' ELASTIC X-SECTION = ',sigele
1448  WRITE(6,*)' INELASTIC X-SECTION-LMD = ',sigine
1449  WRITE(6,*)' INELASTIC X-SECTION = ',sigin
1450  WRITE(6,*)' HARD INEL. X-SECTION = ',sighin
1451  WRITE (6,102)
1452  WRITE(6,*)' LOW MASS SING./DOUB.DIFFR.X-SECTION= ',silmsd,silmdd
1453  WRITE(6,*)' => LOW MASS TOTAL DIFFRACTIV.X-SECTION= ',siglmd
1454  WRITE(6,*)' HIGH MASS SING./DOUB.DIFFR.X-SECTION= ',sigdi,sigddi
1455  WRITE(6,*)' => HIGH MASS TOTAL DIFFRACTIV.X-SECTION= ',sithmd
1456  WRITE(6,*)' ESTIMAT.MIXED (LM+HM) DOUBL.DIFFRAC.X.SEC.= ',slhmdd
1457  WRITE(6,*)' => '
1458  WRITE(6,*)' DIFFRACTIVE X-SECTION = ',sigd
1459  WRITE(6,*)' DOUBLY DIFFRACTIVE X-SECT. =',sigdd
1460  WRITE (6,102)
1461 *
1462  IF(ioutpo.GE.0) THEN
1463  WRITE(6,*)' --- observ. x-sections, altern. calculated ---'
1464  WRITE(6,*)' ELASTIC X-SECTION = ',sigel
1465  WRITE(6,*)' INELASTIC X-SECTION-LMD = ',siginl
1466  WRITE(6,*)' HARD INEL. X-SECTION= ',sighri
1467  WRITE(6,*)' HIGH MASS SING./DOUB.DIFFR.X-SECT.=',sighmd,sihmdd
1468  WRITE(6,*)' X-SECTION FOR (L,M,N,NN)= 1000 0100 0010 0001'
1469  WRITE(6,*)' ',sig(1,0,0),sig(0,1,0)
1470  * ,sig(0,0,1),sig(0,0,2)
1471  WRITE (6,102)
1472  ENDIF
1473 *
1474  IF(ioutpo.GE.2) THEN
1475  WRITE (6,102)
1476  nnmaxp=nmaxi/2
1477  IF( nmaxi.LT.2)nnmaxp=1
1478  DO 52 n=0,nnmaxp
1479 * printout loops:
1480  DO 48 l=0,lmaxi
1481  48 WRITE(6,101)(sig(l,m,n),m=0,7)
1482  WRITE (6,102)
1483  DO 50 l=0,lmaxi
1484  50 WRITE(6,101)(sig(l,m,n),m=8,15)
1485  WRITE (6,102)
1486  WRITE(6,*)
1487  & ' # CUT-POMERON SSOFT X-SECT. SHARD X-SECT.'
1488  DO 58 l=0,lmaxi
1489  58 WRITE (6,103)l,ssoft(l),shard(l)
1490  WRITE (6,102)
1491 * printoutloop ends
1492 52 CONTINUE
1493  ENDIF
1494 *
1495 * *** attribute x-sections (SIG) for CUT objects ('s and PHI's)
1496 * to string configurations (PLMN) s: **********
1497 *
1498 C CHANGED 10.1.90 BY J.R.
1499 * preparations:
1500 C just for Y and (Phi) - cuts:
1501  fak(0)=1
1502  DO 500 i=1,nmaxi
1503  fak(i)=fak(i-1)*i
1504  500 CONTINUE
1505  DO 501 i=0,nmaxi
1506  DO 501 j=0,i
1507  cmbin(i,j)=fak(i)/(fak(j)*fak(i-j))
1508  501 CONTINUE
1509 *
1510  tmmp=0.
1511  DO 5 l=0,lmaxi
1512  DO 5 m=0,mmaxi
1513  IF(icon.EQ.44.OR.icon.EQ.46.OR.icon.EQ.48.
1514  * or.icon.EQ.54) THEN
1515 C///test:
1516 * no Y or PHI cut
1517  plmntm=sig(l,m,0)/(sigsum+tin)
1518  plmn(l,m,0) = plmntm + plmn(l,m,0)
1519  tmmp=tmmp+plmntm
1520 * Y but no PHI cut
1521  plmntm=sig(l,m,1)/(sigsum+tin)
1522  tmmp=tmmp+plmntm
1523  IF(l+2.LE.lmaxi) THEN
1524  plmn(l+2,m,0) = (-2.)* plmntm + plmn(l+2,m,0)
1525  plmn(l+1,m,0) = 4. * plmntm + plmn(l+1,m,0)
1526  ELSE
1527  plmn(lmaxi,m,0) = (-2.)* plmntm + plmn(lmaxi,m,0)
1528  plmn(lmaxi,m,0) = 4. * plmntm + plmn(lmaxi,m,0)
1529  ENDIF
1530  IF(l.EQ.0 .AND. m.EQ.0) THEN
1531  plmn(l ,m,1) = (-1.)* plmntm + plmn(l ,m,1)
1532  ELSE
1533  plmn(l ,m,0) = (-1.)* plmntm + plmn(l ,m,0)
1534  ENDIF
1535 * no Y but PHI cut
1536  plmntm=sig(l,m,2)/(sigsum+tin)
1537  tmmp=tmmp+plmntm
1538  IF(l+2.LE.lmaxi) THEN
1539  plmn(l+2,m,0) = (-2.)* plmntm + plmn(l+2,m,0)
1540  plmn(l+1,m,0) = 4. * plmntm + plmn(l+1,m,0)
1541  ELSE
1542  plmn(lmaxi,m,0) = (-2.)* plmntm + plmn(lmaxi,m,0)
1543  plmn(lmaxi,m,0) = 4. * plmntm + plmn(lmaxi,m,0)
1544  ENDIF
1545  IF(l.EQ.0 .AND. m.EQ.0) THEN
1546  plmn(l ,m,2) = (-1.)* plmntm + plmn(l ,m,2)
1547  ELSE
1548  plmn(l ,m,0) = (-1.)* plmntm + plmn(l ,m,0)
1549  ENDIF
1550 C/// test end
1551  ELSE
1552  DO 51 n=0,nmaxi
1553  DO 51 nn=0,nnmaxi
1554  IF(nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1) go to 51
1555  nnn=n+(nmaxi+1)*nn
1556 *
1557 * to be attributed::
1558  plmntm=sig(l,m,nnn)/(sigsum+tin)
1559  tmmp=tmmp+plmntm
1560 *
1561 * attribution loop for Y-cuts
1562  DO 511 n0cut=0,n
1563  DO 511 n1cut=0,n-n0cut
1564  n2cut=n-n0cut-n1cut
1565 * combinatoric weight:
1566  cmb0=cmbin(n,n2cut)
1567  cmb1=cmbin(n-n2cut,n1cut)
1568 *
1569 * attribution loop for PHI-cuts
1570  DO 511 nn0cut=0,nn
1571  DO 511 nn1cut=0,nn-nn0cut
1572  nn2cut=nn-nn0cut-nn1cut
1573 * combinatoric weight:
1574  cmbn0=cmbin(nn,nn2cut)
1575  cmbn1=cmbin(nn-nn2cut,nn1cut)
1576 *
1577 * attributions matrix:
1578 * ("L"soft,"M"hard,"N"1-diffr.,"NN"2-dif.,"NL"diffr.long in.part ):
1579 * obviously:
1580  l2str = l
1581  m2str = m
1582  n2str = n0cut
1583  nn2str= nn0cut
1584  nl2str= 0
1585 * specialties for Y's and PHI's:
1586  l2str=l2str + n1cut + nn1cut + n2cut + nn2cut
1587  IF(nmax.LE.2)THEN
1588 * room to have NL2STR's:
1589  nl2str= n2cut + nn2cut
1590  ELSEIF(nmax.GE.3)THEN
1591 * the extra inner piece counts here like long strings:
1592  l2str=l2str+n2cut+nn2cut
1593  ENDIF
1594  IF((icon.EQ.26.OR.icon.EQ.36.OR.icon.EQ.46.OR.icon.EQ.56)
1595  & .AND. (l2str.GE.1.OR.m2str.GE.1))THEN
1596  l2str=l2str + nl2str
1597  n2str = 0
1598  nn2str = 0
1599  nl2str = 0
1600  ENDIF
1601 *
1602 * getting and checking parameter for storing:
1603  IF(l2str.GT.lmaxi) l2str=lmaxi
1604  IF(m2str.GT.lmaxi) m2str=lmaxi
1605  nnnstr =n2str +(nmaxi+1)*nn2str
1606  * +(nnmaxi+1)*(nmaxi+1)*nl2str
1607  IF(nnnstr.GT.mxpa13) nnnstr=mxpa13
1608 *
1609 * summing contributions
1610  plmn(l2str,m2str,nnnstr) = plmntm
1611  * *cmb0*cmb1 * (-2)**n2cut * (4)**n1cut * (-1)**n0cut
1612  * *cmbn0*cmbn1*(-2)**nn2cut* (4)**nn1cut* (-1)**nn0cut
1613  & + plmn(l2str,m2str,nnnstr)
1614 *
1615  511 CONTINUE
1616  51 CONTINUE
1617  ENDIF
1618 C///// initial general methode ends
1619  5 CONTINUE
1620  IF(abs(tmmp-1.d0).GT..03d0)THEN
1621  WRITE(6,*)
1622  & ' NORMALISATION ERROR SUM PLM before LMD reatribution=',tmmp
1623  ENDIF
1624 *
1625 * *** built in low mass diffraction, get averages and check normalisation:
1626 *
1627 * low mass diffraction was SIGLMD=SILMSD+SILMDD
1628 * and mixed LM/HM diffraction SLHMDD=SQRT(SILMDD*SIHMDD)
1629  plmfac= (sigsum+tin) / (sigsum+tin +siglmd)
1630  plmn(0,0,1)= plmn(0,0,1) +
1631  & ( silmsd - slhmdd ) / (sigsum+tin)
1632  plmn(0,0,2)= plmn(0,0,2) +
1633  & ( silmdd + slhmdd ) / (sigsum+tin)
1634  661 CONTINUE
1635 * AVerage_SOft_N,AVerage_HaRD_N,SUM_over_Pl
1636  avsofn=0.
1637  avharn=0.
1638  avdifn=0.
1639  avddfn=0.
1640  avdlfn=0.
1641  psoft=0.
1642  temp=0.
1643  tmp=0.
1644  tmmp=0.
1645  tmmp1=0.
1646 *
1647 * (L,M,N,NN,NL repl. L2STR,M2... for s.,h.,Y-dif.,PHI-dif.,i.Y-dif.str.#)
1648  DO 6 nl=0,nlmaxi
1649  DO 6 nn=0,nnmaxi
1650  DO 6 n=0,nmaxi
1651  IF(nmax.LE.2 .AND. n+nn+nl.GE.2) go to 6
1652  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1653  DO 63 m=0,mmaxi
1654  DO 63 l=0,lmaxi
1655  IF(nl.EQ.0)tmmp1 = tmmp1 + sig(l,m,nnn)
1656  tmmp = tmmp + sig(l,m,nnn)
1657  plmn(l,m,nnn)=plmn(l,m,nnn) * plmfac
1658  tmp = tmp + plmn(l,m,nnn)
1659 C IF(PLMN(L,M,NNN).LT.-.000001D0)
1660  IF(plmn(l,m,nnn).LT.-.000005d0)
1661  & WRITE(6,*)' 0>PLMN',plmn(l,m,nnn),l,m,n,nn,nl
1662  avsofn=avsofn+plmn(l,m,nnn)*l
1663  avharn=avharn+plmn(l,m,nnn)*m
1664  avdifn=avdifn+plmn(l,m,nnn)*n
1665  avddfn=avddfn+plmn(l,m,nnn)*nn
1666  avdlfn=avdlfn+plmn(l,m,nnn)*nl
1667  IF (m.EQ.0)psoft=psoft+plmn(l,m,nnn)
1668  63 CONTINUE
1669  6 CONTINUE
1670  IF(abs(tmp-1.d0).GT..01d0)THEN
1671  WRITE(6,*)
1672  & ' NORMALISATION ERROR SUM PLM before M reatribution=',tmp
1673  ENDIF
1674  tmmp=tmmp/sigsum
1675  tmmp1=tmmp1/sigsum
1676  IF(abs(tmmp-1.d0).GT..01d0 .OR.abs(tmmp1-1.d0).GT..01d0)THEN
1677  WRITE(6,*)
1678  & ' NORMALISATION ERROR TMMP,TMMP1=',tmmp,tmmp1
1679  ENDIF
1680 *
1681 * *** reattribute purely hard scattering and get cummulant distribution
1682 *
1683 * (L,M,N,NN,NL repl. L2STR,M2...
1684 * for soft,hard,Y-diffr.,PHI-diffr.,inner diffr. string #)
1685  DO 61 nl=0,nlmaxi
1686  DO 61 nn=0,nnmaxi
1687  DO 61 n=0,nmaxi
1688  IF(nmax.LE.2 .AND. n+nn+nl.GE.2) go to 61
1689  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1690  DO 612 m=0,mmaxi
1691  DO 611 l=0,lmaxi
1692 * -- hard: a pure hard scattering gets an extra soft chain, it is considered
1693 * a specialty of fragmentation and therefor implemented only here
1694 * there are are number of other options which are dropped as they
1695 * are hard to implement at this point
1696  IF (l.EQ.0.AND.m.GE.1)THEN
1697  plmn(1,m,nnn)=plmn(1,m,nnn)+plmn(0,m,nnn)
1698  plmn(0,m,nnn)=0.
1699  ENDIF
1700 * -- cummulant of distribution
1701  temp = temp + plmn(l,m,nnn)
1702  plmncu(l,m,nnn)=temp
1703  611 CONTINUE
1704 *
1705  IF(ioutpo.GE.3)WRITE (6,*)' M,(L,PLMN(L,M,N),L=0,LMAX)'
1706  IF(ioutpo.GE.3)WRITE (6,106) m,(l,plmn(l,m,n),l=0,lmaxi)
1707  IF(ioutpo.GE.2)WRITE (6,*)' M,(L,PLMNCU(L,M,N),L=0,LMAX/2)'
1708  IF(ioutpo.GE.2)WRITE (6,106) m,(l,plmncu(l,m,n),l=0,lmaxi/2)
1709  106 FORMAT (i3,9(i3,e11.2))
1710 *
1711  612 CONTINUE
1712  61 CONTINUE
1713 *
1714  IF(abs(temp-1.d0).GT..01d0)THEN
1715  WRITE(6,*)' NORMALISATION ERROR SUM PLM=',temp
1716  plmfac=1./(temp+tin)
1717  go to 661
1718  ENDIF
1719 *
1720  IF(ioutpo.GE.1)WRITE (6,*)
1721  & '(((L,M,N,PLMN(L,M,N),N=0,2),M=0,5),L=0,7)'
1722  IF(ioutpo.GE.1)WRITE (6,1106)
1723  & (((l,m,n,plmn(l,m,n),n=0,2),m=0,5),l=0,7)
1724  IF(ioutpo.GE.1)WRITE (6,*)
1725  & '(((L,M,N,SIG(L,M,N),N=0,2),M=0,5),L=0,7)'
1726  IF(ioutpo.GE.1)WRITE (6,1106)
1727  & (((l,m,n,sig(l,m,n),n=0,2),m=0,5),l=0,7)
1728  1106 FORMAT (1x,3(i5,i5,i5,g12.5))
1729 *
1730  phard=1.-psoft
1731  alfah=sighin/(sigine+0.00001)
1732  betah=1.-alfah
1733  WRITE(6,116)avsofn,avharn,avdifn,avddfn,avdlfn,
1734  & phard,psoft,alfah,betah
1735  116 FORMAT(/'--- various averages:'/
1736  & /' AVSOFN= AVHARN= AVDIFN= AVDDFN= AVDLFN='
1737  & /' ',5f11.3
1738  & /' PHARD= PSOFT= ALFAH= BETAH= '
1739  & /' ',4f11.3)
1740  IF(ioutpo.GE.1)WRITE(6,*)'SIGSUM=SIGINL-LMD',sigsum
1741 *
1742  IF(ioutpo.GE.1)WRITE(6,610) sigtot,sigine,sigd,sigdd,sighin
1743  610 FORMAT (' SIGTOT,SIGINE,SIGD,SIGDD,SIGHIN= '/' ',5e18.6)
1744 *
1745 101 FORMAT(' ',10e10.3)
1746 102 FORMAT(' ')
1747 103 FORMAT(' ',5x,i4,5x,2e15.3)
1748 *
1749  RETURN
1750  END
1751 * ende problm
1752 *
1753 ************************************************************************
1754 *
1755 *
1756  SUBROUTINE samplx(L2STR,M2STR,N2STR,NN2STR,NL2STR)
1757 *
1758 * input:
1759 * PLMNCU
1760 * output:
1761 * samples number of soft (L) and hard (M) cut pomerons from PLMNC
1762 * and of (N=0/1/2) diffractive excitations (for L=M=0
1763 *
1764 *----------------------------------------------------------------------
1765  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1766  SAVE
1767  COMMON /nncms/ gamcm,bgcm,umo,pcm,eproj,pproj
1768 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20),NESTEP
1769  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
1770  parameter(mxpu50=100,mxpu51=mxpu50+1)
1771  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
1772 * PARAMETRIZATION FOR PTMIN= 3. GEV
1773  parameter(mxpa50=250,mxpa51=mxpa50+1)
1774 * PARAMETRIZATION FOR PTMIN= 2. GEV
1775 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
1776 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
1777  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
1778  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
1779 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
1780  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
1781  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
1782 C COMMON /POLMN1/ PLMNEE(0:MXPA25,0:MXPU50,0:MXPA13,20)
1783  COMMON /polmn1/ plmnee(0:mxpa25,0:mxpu50,0:mxpa13,28)
1784  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
1785  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
1786 *
1787  parameter(pi=3.141592654d0)
1788  DATA nprint/0/
1789 C Determine the energy index
1790  ipoen=1
1791  DO 20 ii=1,nestep
1792  IF(umo.GE.poen1(ii).AND.umo.LT.poen2(ii))THEN
1793  ipoen=ii
1794  go to 22
1795  ENDIF
1796  20 CONTINUE
1797  22 CONTINUE
1798 * "L"SOFT,"M"HARD
1799  lmaxi = lmax
1800  mmaxi = mmax
1801  IF(ipim.NE.2) THEN
1802 * "N"diffr.,"NN"dou.dif.,("NL"long inner contr. for 2*cut Y or PHI)
1803  nmaxi = nmax
1804  nnmaxi=0
1805  nlmaxi=0
1806  ELSEIF(ipim.EQ.2) THEN
1807  IF( nmax.GE.3)THEN
1808  nmaxi = nmax
1809  nnmaxi=(13-nmaxi)/(1+nmaxi)
1810 * 13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
1811  nlmaxi=0
1812  ELSEIF( nmax.EQ.2)THEN
1813  nmaxi=1
1814  nnmaxi=1
1815  nlmaxi=1
1816  ELSEIF( nmax.EQ.1)THEN
1817  nmaxi=1
1818  nnmaxi=0
1819  nlmaxi=1
1820  ENDIF
1821  ENDIF
1822  111 CONTINUE
1823 *
1824  x=rndm(v)
1825 *
1826  IF (x.LE.plmncu(0,0,0) .AND. nprint.LT.100)THEN
1827  WRITE(6,*) ' No generator of elastic events '
1828  WRITE(6,*) ' PLMNCU (0,0,0) =!= 0 = ',plmncu(0,0,0)
1829  nprint=nprint+1
1830  goto 111
1831  ENDIF
1832 *
1833  DO 5 nl=0,nlmaxi
1834  DO 5 nn=0,nnmaxi
1835  DO 5 n=0,nmaxi
1836  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1837  DO 6 m=0,mmaxi
1838  DO 7 l=0,lmaxi
1839 *
1840 C IF (X.LE.PLMNCU(L,M,NNN)) THEN
1841  IF (x.LE.plmnee(l,m,nnn,ipoen)) THEN
1842  l2str=l
1843  m2str=m
1844  n2str=n
1845  nn2str=nn
1846  nl2str=nl
1847  RETURN
1848 *
1849  ENDIF
1850  7 CONTINUE
1851  6 CONTINUE
1852  5 CONTINUE
1853 *
1854  nprint=nprint+1
1855  IF(nprint.LT.100) WRITE(6,*)' RAR.IN SAMPLM,PLMNCU,RND=',
1856  & plmncu(lmax, mmax,nnn),x,nprint
1857  IF( plmncu(lmax,mmax,nnn) .GT. 0.1d0 ) RETURN
1858  IF( plmncu(lmax,0,0) .GT. 0.1d0 ) RETURN
1859  WRITE(6,*)' RAR.IN SAMPLM- PROBLEM SEEMS BAD, DECIDE TO STOP'
1860  stop
1861  END
1862 *
1863 *
1864 ************************************************************************
1865 *
1866  SUBROUTINE samplm(L2STR,M2STR,N2STR)
1867 *
1868 * input:
1869 * PLMNCU
1870 * output:
1871 * samples number of soft (L) and hard (M) cut pomerons from PLMNC
1872 * and of (N=0/1/2) diffractive excitations (for L=M=0
1873 *
1874 *----------------------------------------------------------------------
1875  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1876  SAVE
1877  COMMON /nncms/ gamcm,bgcm,umo,pcm,eproj,pproj
1878 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20),NESTEP
1879  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
1880  parameter(mxpu50=100,mxpu51=mxpu50+1)
1881 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
1882  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
1883  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
1884 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
1885  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
1886 * PARAMETRIZATION FOR PTMIN= 3. GEV
1887  parameter(mxpa50=250,mxpa51=mxpa50+1)
1888 * PARAMETRIZATION FOR PTMIN= 2. GEV
1889 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
1890  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
1891  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
1892 C COMMON /POLMN1/ PLMNEE(0:MXPA25,0:MXPU50,0:MXPA13,20)
1893  COMMON /polmn1/ plmnee(0:mxpa25,0:mxpu50,0:mxpa13,28)
1894  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
1895  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
1896 *
1897  parameter(pi=3.141592654d0)
1898  ipoen=1
1899  DO 20 ii=1,nestep
1900  IF(umo.GE.poen1(ii).AND.umo.LT.poen2(ii))THEN
1901  ipoen=ii
1902  go to 22
1903  ENDIF
1904  20 CONTINUE
1905  22 CONTINUE
1906  111 CONTINUE
1907 *
1908  x=rndm(v)
1909 *
1910  IF (x.LE.plmncu(0,0,0))THEN
1911  WRITE(6,*) ' No generator of elastic events '
1912  WRITE(6,*) ' PLMNCU (0,0,0) =!= 0 = ',plmncu(0,0,0)
1913  goto 111
1914  ENDIF
1915 *
1916  DO 5 n=0,nmax
1917 C IF (N.GT.1) GO TO 111
1918  DO 6 m=0,mmax
1919  DO 7 l=0,lmax
1920 *
1921 C IF (X.LE.PLMNCU(L,M,N)) THEN
1922  IF (x.LE.plmnee(l,m,n,ipoen)) THEN
1923  l2str=l
1924  m2str=m
1925  n2str=n
1926  RETURN
1927 *
1928  ENDIF
1929  7 CONTINUE
1930  6 CONTINUE
1931  5 CONTINUE
1932 *
1933  WRITE(6,*)' RAR.IN SAMPLM,PLMNCU,RND=',plmncu(lmax,mmax,nmax),x
1934  IF( plmncu(lmax,mmax,nmax) .GT. 0.1d0 ) RETURN
1935  IF( plmncu(lmax,0,0) .GT. 0.1d0 ) RETURN
1936  WRITE(6,*)' RAR.IN SAMPLM- PROBLEM SEEMS BAD, DECIDE TO STOP'
1937  stop
1938  END
1939 *
1940 *
1941 *
1942 *
1943 ************************************************************************
1944 C--------------------------------------------------------------------
1945 C
1946 C dtUpom9h.for
1947 C
1948 C---------------------------------------------------------------------
1949 ************************************************************************
1950 *
1951 * POMDI,SIGMAS,SIGSHD,PRBLM0..9,SAMPLM,SIGMA1
1952 *
1953 * routines called by code word
1954 * SIGMAPOM :
1955 *
1956 
1957 C--------------------------------------------------------------------
1958 C
1959 C dtUpom9h.for
1960 C
1961 C---------------------------------------------------------------------
1962 ************************************************************************
1963 *
1964 * POMDI,SIGMAS,SIGSHD,PRBLM0..9,SAMPLX,SIGMA1
1965 *
1966 * routines called by code word
1967 * SIGMAPOM :
1968 *
1969 ************************************************************************
1970 *
1971 * J.RANFT September 1987
1972  SUBROUTINE pomdi
1973 *
1974 * to calculate the s-dependent X-sections
1975 *
1976 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1977 *
1978 * input:
1979 * ISIG characterizing X sections, transmitted to SIGSHD
1980 * IPIM characterizes the method to calculate SIGMA(LSOFT,mhard)
1981 * IPIM=1 : integral method with low mass dif.matr
1982 * IPIM=2 : int.meth.with low mass dif.matr.(2*2)
1983 * and reduced high mass diffraction
1984 * IPIM=3 : integral methode
1985 * IPIM=4 : integral methode with Y -cuts
1986 * IPIM=5 : integral methode with Y-cuts + 2 CHANNEL EIK.
1987 * rest : not implemented or
1988 * special cases for checking
1989 * LMAX < MXPA25 maximal number of considered soft pomerons
1990 * MMAX < MXPA50 maximal number of considered hard pomerons
1991 * NMAX < MXPA13 maximal number of considered trippel pomerons
1992 * (not used in low masss diffraction formalism)
1993 * output: printet and plottet in SIGMAS and
1994 * printed on the end of this subroutine
1995 *
1996 * card XSECTION causes call to POMDI
1997 * card SIGMAPOM with ITEST=1 causes call to POMDI
1998 * POMDI calling
1999 * SIGMAS ( SIGMA1..3( SIGSD, (1:/3:)GSET ), PLOT)
2000 * PRBLM.. ( SIGSHD, (1:)GSET )
2001 * (2:)SAMPLX, (1:/3:)SAMPLM
2002 *
2003 *----------------------------------------------------------------------
2004  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2005  SAVE
2006 *
2007 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
2008  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
2009 * PARAMETRIZATION FOR PTMIN= 3. GEV
2010  parameter(mxpa50=250,mxpa51=mxpa50+1)
2011 * PARAMETRIZATION FOR PTMIN= 2. GEV
2012 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
2013  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
2014  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
2015  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
2016  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
2017 * *** /HISTOO/
2018  INTEGER*2 ndislm
2019  COMMON /histoo/as(50,9),aecm(50,9),asig(50,9),alos(50,9),
2020  * aloecm(50,9),ndislm(0:mxpa25,0:mxpa50,0:mxpa13)
2021 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
2022  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
2023 * --- used only in SIGMAPOM-routines
2024 * *** /POMPAR/ contains pomeron parameters used in some options
2025 * ALFA,ALFAP,A,BH,C,BS are soft Pomeron parameters chosen in SIGMAS
2026 * appearing in SIGMAS, PRBLM.., AK is K-factor
2027  common/pompar/alfa,alfap,a,c,ak
2028 * *** /SIGMA/ contains variables actually used in iteration
2029 * ZSOF, ZHAR, BS, BH, SIGSOF,SIGHAR ( for soft, hard, trippel-Pom.)
2030 * are input for X-section calculated in SIGSHD
2031 C (/sigma/ is out of P.A.'s CM88, containing variables used in iter.
2032  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
2033  * sigloo,zloo
2034 * *** /POMTYP/ contains parameters determining X-sections
2035 * IPIM,ISIG,LMAX,MMAX,NMAX as described at "CODEWD = SIGMAPOM"
2036 C LMAX,MMAX,NMAX kept open to enable avoiding of numerical problem
2037  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
2038  COMMON /alala/alalam
2039 C COMMON/COLLIS/SS,IJPROJ,IJTAR,PTTHR,PTTHR2,IOPHRD,IJPRLU,IJTALU
2040 C
2041  COMMON /collpo/s,ptthr,ptthr2
2042 C COMMON /COLLIS/SS,IJPROJ,IJTAR,PTTPO,IOPHRD,IJPRLU,IJTALU,PTTPO2
2043  COMMON /collis/ss,ijproj,ijtar,pttpo,pttpo2,iophrd,ijprlu,ijtalu
2044 C
2045 C ECM calculated as SQRT
2046 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2047  parameter(pi=3.141592654d0)
2048  lmaxi = lmax
2049  mmaxi = mmax
2050  nmaxi = nmax
2051  skeep=ss
2052 *
2053 *---------------------------------------------------------------------
2054 * CALL SIGMAS
2055 *
2056 * calculates and plotts out X-section at various energies
2057 * runs thru energies independent of the rest
2058 * and contains no preperation for other work
2059 *
2060  WRITE(6,'(1X/1X)')
2061  WRITE(6,*)' '
2062  WRITE(6,*)
2063  *' ------ testing the energy dependence of x-sections ----------'
2064  WRITE(6,*)' '
2065  IF(ioutpo.GT.-1) WRITE(6,*)
2066  *' (as function of ALAM i.e.a low mass diffr.parameter)'
2067  WRITE(6,*)' -----------------------------------------------'
2068  WRITE(6,'(1X)')
2069 *
2070  DO 1007 iijj=1,10
2071  IF(ioutpo.GT.-1 .OR. iijj.EQ.6)THEN
2072 C WITHOUT SIMGMAS FOR NORMAL USERS
2073 C IF(IOUTPO.GT.-1 )THEN
2074  alalam=iijj*0.1
2075  IF(ioutpo.GT.-1) WRITE(6,1008)alalam
2076  1008 FORMAT (' ALAM= ',f10.3)
2077  CALL sigmas
2078  ENDIF
2079  1007 CONTINUE
2080 *
2081 *---------------------------------------------------------------------
2082 * test sampling L,M
2083 *
2084  1111 CONTINUE
2085 * looping thru the energies
2086 *
2087 C DO 100 III=3,9
2088 C S=10.**III
2089 C ECM=SQRT(S)
2090 *
2091  s=skeep
2092  ecm=sqrt(s)
2093 *
2094 * chosing one options with a standard call
2095 *
2096  IF(ipim.EQ.2) THEN
2097  IF( nmax.GE.3)THEN
2098  nnmaxi=(13-nmaxi)/(1+nmaxi)
2099 * 13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
2100  nlmaxi=0
2101  ELSEIF( nmax.EQ.2)THEN
2102  nmaxi=1
2103  nnmaxi=1
2104  nlmaxi=1
2105  ELSEIF( nmax.EQ.1)THEN
2106  nmaxi=1
2107  nnmaxi=0
2108  nlmaxi=1
2109  ENDIF
2110  CALL prblm2(ecm)
2111  ENDIF
2112  IF(ipim.LT.1.AND.ipim.GT.9)THEN
2113  WRITE(6,*) 'RETURN caused by IPIM=',ipim
2114  RETURN
2115  ENDIF
2116 *
2117 * randomly sample events and printout
2118 *
2119  WRITE (6,'(1X)' )
2120  WRITE (6,102)ecm,s
2121  102 FORMAT
2122  * ('--- sample distribution for L soft and M hard inelastic'
2123  * , ' pomerons (string pairs)--- '
2124  * / 20x,'at ECM = ',f10.2,' S = ',f12.1)
2125  DO 31 l=0,lmaxi
2126  DO 32 m=0,mmaxi
2127  DO 32 n=0,13
2128  ndislm(l,m,n)=0
2129  32 CONTINUE
2130  31 CONTINUE
2131 *
2132  IF(icon.EQ.12)go to 100
2133  DO 320 ii=1,10000
2134  IF(ipim.EQ.2) THEN
2135  CALL samplx(l2str,m2str,n2str,nn2str,nl2str)
2136  nnnstr =n2str +(nmaxi+1)*nn2str
2137  * +(nnmaxi+1)*(nmaxi+1)*nl2str
2138  ndislm(l2str,m2str,nnnstr)=ndislm(l2str,m2str,nnnstr)+1
2139  ELSE
2140  CALL samplm(l2str,m2str,n2str)
2141  ndislm(l2str,m2str,n2str)=ndislm(l2str,m2str,n2str)+1
2142  ENDIF
2143  320 CONTINUE
2144 *
2145  WRITE(6,*)
2146  * ' with no diffractive contribution'
2147  WRITE(6,*) ' '
2148  WRITE(6,*)
2149  * ' ....... vertical: NSTR, horizontal MSTR .........'
2150  DO 3344 l=0,min(20,lmaxi)
2151  3344 WRITE(6,34)l,(ndislm(l,m,0),m=0,20)
2152  WRITE(6,*) ' '
2153  IF(ioutpo.GE.0)THEN
2154  DO 333 n=0,5
2155  IF(n.NE.0)THEN
2156  WRITE(6,*)' WITH NSTR=',n
2157  DO 334 l=0,min(20,lmaxi)
2158  WRITE(6,34)l,(ndislm(l,m,n),m=0,20)
2159  334 CONTINUE
2160  WRITE(6,*) ' '
2161  ENDIF
2162  jmpa50 = int(mxpa50/25)
2163 C WRITE(6,*) 'WIDE PLOT 0<L<25, 0<M<'
2164  WRITE(6,*) 'WIDE PLOT 0<L<',mxpa25,' 0<M<'
2165  & ,mxpa50,' IN STEPS OF ',jmpa50
2166 C DO 335 L=0,MIN(25,LMAXI)
2167  DO 335 l=0,mxpa25
2168  WRITE(6,35)l,(ndislm(l,m,n),m=0,mxpa50,jmpa50)
2169  335 CONTINUE
2170  WRITE(6,*) ' '
2171  333 CONTINUE
2172  ENDIF
2173  34 FORMAT (i5,':',21i4)
2174  35 FORMAT (i5,26i4)
2175  100 CONTINUE
2176 * energy loop has ended
2177  RETURN
2178  END
2179 *
2180 *
2181 ************************************************************************
2182 *
2183  SUBROUTINE sigmas
2184 *
2185 * output:
2186 * energy dependence of
2187 * SIGMa-TOT, SIGma-INel, SIGma-Diffractive, SIGma-Hard-INelastic
2188 * print- and plott-out
2189 * using:
2190 * called routines SIGSHD , SIGMA1..3 , PLOT
2191 *
2192 * runs thru energies independent of the rest
2193 * and contains no preperation for other parts
2194 *
2195 *---------------------------------------------------------------------
2196  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2197  SAVE
2198 *
2199 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
2200  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
2201 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
2202  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
2203  parameter( zero=0.d0, one=1.d0)
2204 * PARAMETRIZATION FOR PTMIN= 3. GEV
2205  parameter(mxpa50=250,mxpa51=mxpa50+1)
2206 * PARAMETRIZATION FOR PTMIN= 2. GEV
2207 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
2208  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
2209  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
2210  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
2211  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
2212 *
2213 * *** /POMPAR/*/SIGMA/*/POMTYP/ used only in SIGMAPOM-routines (->POMDI)
2214  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
2215  common/pompar/alfa,alfap,a,c,ak
2216  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
2217  * sigloo,zloo
2218 * ***
2219  COMMON /topdr/itopd,idumtp
2220 * ***
2221  INTEGER*2 ndislm
2222  COMMON /histoo/as(50,9),aecm(50,9),asig(50,9),alos(50,9),
2223  * aloecm(50,9),ndislm(0:mxpa25,0:mxpa50,0:mxpa13)
2224 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2225  parameter(pi=3.141592654d0)
2226 *
2227 *---------------------------------------------------------------------
2228 * run thru energies
2229 * ------------------------------------------------------------------
2230 *
2231  istep=2
2232  IF(ioutpo.GT.-1)istep=7
2233  DO 100 i=1,50,istep
2234  s=1.6**i
2235  ecm=sqrt(s)+3.4
2236  s=ecm**2
2237 *
2238 * *** running thru energies initializing AS,AECM,ASIG,ALOS,ALOECM *****
2239 *
2240  DO 111 iii=1,9
2241  as(i,iii)=s
2242  aecm(i,iii)=ecm
2243  alos(i,iii)=log10(s)
2244  aloecm(i,iii)=log10(ecm)
2245  asig(i,iii)=0.
2246  111 CONTINUE
2247 *
2248 * *** calling calculation of x- section
2249 *
2250  IF(ipim.EQ.2 )THEN
2251  CALL sigma2(ecm)
2252  IF(i.EQ.1 .AND. ioutpo.GE.0 ) WRITE(6,*)
2253  & ' s-dep. by integr.with Y,PHI,LMD'
2254  ELSE
2255  CALL sigma2(ecm)
2256  IF(i.EQ.1 .AND. ioutpo.GE.0 ) WRITE(6,*)
2257  & ' s-dep. by integr.with Y,PHI,LMD (DEFAULT)'
2258  ENDIF
2259 *
2260 *
2261 C* the preceeding coresponds to a call to SIGMA2 except extra prin
2262 * *** getting plot array ASIG(I,J) :
2263  asig(i,1)=sigtot
2264  asig(i,2)=sigine
2265  asig(i,3)=sighin
2266  asig(i,4)=sigsof
2267  asig(i,5)=sighar
2268  asig(i,6)=sigtrp
2269  asig(i,7)=sigtot-sigine
2270  asig(i,8)=sigine-sighin
2271  asig(i,9)=sigd
2272  WRITE (6,1007)ecm,sigtot,sigine,sigel,sigd
2273  1007 FORMAT (' ECM,SIGTOT,SIGINE,SIGEL,SIGD',f10.1,4e14.3)
2274  100 CONTINUE
2275 *
2276 * energy loop ends
2277 *---------------------------------------------------------------------
2278 * print out results
2279 * ------------------------------------------------------------------
2280  WRITE (6,991)
2281  991 FORMAT (//' shown as line printer plott'/' with'/
2282 * J: drawn quantities:
2283  1 ' (*) SIGTOT total x-section',
2284  2 ' (2) SIGINE inelastic x-section'/
2285  3 ' (3) SIGHIN hard inelastic cross section, one or more jets',
2286  4 ' (4) SIGSOF input soft x-section'/
2287  5 ' (5) SIGHAR input hard x-sections',
2288  6 ' (6) SIGTRP input diffractive x-section (triple pomeron)'/
2289  7 ' (7) SIGTOT-SIGINE elastic x-section',
2290  8 ' (8) SIGINE-SIGHIN non-hard inelastic x-section, (no jets)'/
2291  9 ' (9) SIGD diffractive xross section '/
2292  * ' are plotted against LOG(10)of(CMENERGY)' //)
2293 *
2294  CALL plot(aloecm,asig,450,9,50,zero, 0.1*one,zero, 2.0*one)
2295 *
2296 * special output on unit number 7
2297 C I kept it as it was
2298  IF (itopd.EQ.1) THEN
2299  WRITE(7,95)
2300  95 FORMAT(' NEW FRAME'/' SET FONT DUPLEX'/' SET SCALE X LOG'/
2301  * ' SET LIMITS X FROM 1.0 TO 1E5 Y FROM 0. TO 200'/
2302  * ' TITLE TOP < TOTAL,INEL. AND HARD (MINIJET) CROSS SECT.<'/
2303  * ' TITLE BOTTOM <C.M.ENERGY [GEV]<'/
2304  * ' TITLE < DUAL UNITARIZATION OF SOFT AND HARD CROSS SECTIONS<'/
2305  * ' TITLE LEFT LINES=-1 <CROSS SECTION [MB]<'/
2306  * ' TITLE 3 8.5 < SOLID = TOTAL X.S. <'/
2307  * ' TITLE < DASHED= INELASTIC X.S. <'/
2308  * ' TITLE < DOTTED= HARD X.S.<'/
2309  * ' TITLE < DOT-DASH= HARD INPUT X.S. <'/
2310  * ' TITLE < DOT-DASH= ELASTIC X.S. <')
2311  92 FORMAT (5f15.5)
2312  DO 94 iuu=1,7
2313  IF (iuu.EQ.4)go to 94
2314  IF (iuu.EQ.6)go to 94
2315  IF (iuu.EQ.1) WRITE(7,97)
2316  97 FORMAT (' SET TEXTURE SOLID')
2317  IF (iuu.EQ.2) WRITE(7,98)
2318  98 FORMAT (' SET TEXTURE DASHES')
2319  IF (iuu.EQ.3) WRITE(7,99)
2320  99 FORMAT (' SET TEXTURE DOTS')
2321  IF (iuu.EQ.5) WRITE(7,197)
2322  197 FORMAT (' SET TEXTURE DOTDASH')
2323  DO 93 iu=2,46
2324  WRITE(7,92)aecm(iu,iuu),asig(iu,iuu)
2325  93 CONTINUE
2326  WRITE(7,96)
2327  96 FORMAT (' JOIN')
2328  94 CONTINUE
2329  ENDIF
2330 * ending IF(ITOP=1) special output
2331  RETURN
2332  END
2333 *
2334 ******************************************************************
2335 * end dtupom90
2336 ******************************************************************
2337 
c2_transformation::zero
static float_type zero(float_type)
utility function f(x)=0 useful in axis transforms
Definition: c2_function.hh:1203
rdxsec
subroutine rdxsec(XSEC)
Definition: dpm25pom.f:416
qrblm2
subroutine qrblm2(ECM)
Definition: dpm25pom.f:652
sigmas
subroutine sigmas
Definition: dpm25pom.f:2176
gset
subroutine gset(AX, BX, NX, Z, W)
Definition: dpm25nulib.f:689
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
s
const XML_Char * s
Definition: include/expat.h:262
z
G4double z
Definition: TRTMaterials.hh:39
HepGeom::BasicVector3D::unit
BasicVector3D< T > unit() const
Definition: BasicVector3D.h:304
a
G4double a
Definition: TRTMaterials.hh:39
HepGeom::Plane3D::d
T d() const
Definition: Plane3D.h:86
samplm
subroutine samplm(L2STR, M2STR, N2STR)
Definition: dpm25pom.f:1859
c2_transformation::one
static float_type one(float_type)
utility function f(x)=1 useful in axis transforms
Definition: c2_function.hh:1201
samplx
subroutine samplx(L2STR, M2STR, N2STR, NN2STR, NL2STR)
Definition: dpm25pom.f:1749
plot
subroutine plot(X, Y, N, M, MM, XO, DX, YO, DY)
Definition: dpm25nulib.f:176
energy
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
nmax
const G4int nmax
Definition: G4UrQMD1_3Interface.hh:69
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
rndm
double precision function rndm(RDUMMY)
Definition: dpm25nulib.f:1460
sigshd
subroutine sigshd(ECM)
Definition: dpm25pom.f:6
status
int status
Definition: tracer.cxx:24
title
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
c2_factory::log
static c2_log_p< float_type > & log()
make a *new object
Definition: c2_factory.hh:138
c2_factory::sqrt
static c2_sqrt_p< float_type > & sqrt()
make a *new object
Definition: c2_factory.hh:142
prblm2
subroutine prblm2(ECM)
Definition: dpm25pom.f:571
pomdi
subroutine pomdi
Definition: dpm25pom.f:1965
c2_factory::exp
static c2_exp_p< float_type > & exp()
make a *new object
Definition: c2_factory.hh:140
Generated on Wed Apr 30 2014 15:54:32 for Geant4.10 by   doxygen 1.8.7