Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25pom.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  CALL getenv('INIDAT',inidat)
605 C OPEN(IUNIT,FILE=INIDAT,STATUS='UNKNOWN',ERR=99)
606  OPEN(unit=iunit,file='pomtab.dat'
607  * ,status='UNKNOWN',err=99)
608 C
609  IF (ipomta.EQ.0) THEN
610  DO 1 ii=1,nestep
611  energy = poen(ii)
612  CALL qrblm2(energy)
613  DO 10 jj=0,mxpa25
614  DO 10 kk=0,mxpu50
615  DO 10 ll=0,mxpa13
616  plmnee(jj,kk,ll,ii)=plmncu(jj,kk,ll)
617  10 CONTINUE
618  1 CONTINUE
619 C
620  WRITE(iunit,7102) nestep
621  DO 31 ii=1,nestep
622  WRITE(iunit,7101) poen(ii), poen1(ii), poen2(ii)
623  WRITE(iunit,101)(((plmnee(jj,kk,ll,ii),ll=0,mxpa13),
624  * kk=0,mxpu50),jj=0,mxpa25)
625  31 CONTINUE
626  7102 FORMAT(i13)
627  7101 FORMAT(3e13.5)
628  101 FORMAT(8e13.5)
629 C$
630  ELSEIF (ipomta.EQ.1)THEN
631  READ(iunit,7102) nestep
632  WRITE(6,7102) nestep
633  DO 11 ii=1,nestep
634  READ(iunit,7101) poen(ii), poen1(ii), poen2(ii)
635  WRITE(6,7101) poen(ii), poen1(ii), poen2(ii)
636  READ(iunit,101)(((plmnee(jj,kk,ll,ii),ll=0,mxpa13),
637  * kk=0,mxpu50),jj=0,mxpa25)
638 C WRITE(6,101)(((PLMNEE(JJ,KK,LL,II),LL=0,MXPA13),
639 C * KK=0,MXPU50),JJ=0,MXPA25)
640  11 CONTINUE
641  ENDIF
642  CLOSE(iunit)
643  RETURN
644  99 CONTINUE
645  WRITE(6,'(A)')'Error in PRBLM2 : file pomtab.dat ERROR'
646  CLOSE(iunit)
647  stop
648  END
649 
650 ************************************************************************
651 *
652  SUBROUTINE qrblm2(ECM)
653 * * input:
654 * ECM
655 * output:
656 * PLMN, PLMNCUmmulative, AVSOFN, AVHRDN,SIGDD/D/QEL/EL, PSOFT
657 *
658 *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
659 *
660 * Probabilities for L soft cut, M hard cut pomerons
661 * and N someplace cut trippel pomerons
662 *
663 * Aurenche Maire's version of PRBLM
664 * modified from A.M. to include calculation of x-sections
665 * modified to get:
666 * version with high mass diffraction (Y's and PHI's)
667 * and on both sides 2 kanal low mass diffraction
668 * *** OPTION 21.2.90 FB
669 *----------------------------------------------------------------------
670 *
671  IMPLICIT DOUBLE PRECISION(a-h,o-z)
672  SAVE
673  parameter( zero=0.d0, one=1.d0)
674  parameter(conv=0.38935d0)
675  parameter(pi=3.141592654d0)
676  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
677 * PARAMETRIZATION FOR PTMIN= 3. GEV
678  parameter(mxpa50=250,mxpa51=mxpa50+1)
679 * PARAMETRIZATION FOR PTMIN= 2. GEV
680 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
681  parameter(mxpa96=96)
682 C PARAMETER (MXPA96=480)
683  LOGICAL lsqrt
684 C PARAMETER (MXLMN=5,LSQRT=.false.)
685  parameter(mxlmn=5,lsqrt=.true.)
686  DOUBLE PRECISION dtiny
687 C PARAMETER (TINY = 1.2D-38,DTINY=1.D-70,TIN=1.D-22,TINEXP = -300.D0)
688 C PARAMETER (TINY = 1.2D-38,DTINY=1.D-300,TIN=1.D-22,TINEXP =
689 C -700.D0)
690  parameter(tiny=1.2d-38,dtiny=1.d-70,tin=1.d-22,tinexp=-700.d0)
691 C PARAMETER (TINY = 1.D-38,DTINY=tiny,TIN=1.D-22,TINEXP = -300.D0)
692 * in older version used:
693  parameter(tinyex = -48.d0)
694 * --- -- - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - -
695 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
696  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
697  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
698  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
699  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
700 * *** /POMPAR/*/SIGMA/*/POMTYP/ used only in SIGMAPOM-routines (->POMDI)
701 * (LMAX,MMAX,NMAX (max number of soft/hard/trippel pomerons)
702  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
703  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
704  * sigloo,zloo
705  common/pompar/alfa,alfap,a,c,ak
706  COMMON /singdi/silmsd,sigdi
707 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
708  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
709  COMMON /alala/alalam
710  CHARACTER*80 title
711  CHARACTER*8 projty,targty
712 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
713 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
714  COMMON /user1/title,projty,targty
715  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
716 * --- -- - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - -
717  DOUBLE PRECISION sig,sigp,sigm,sign,sigo
718  dimension sig(0:mxpa25,0:mxpa50,0:mxpa13),
719  &sigp(0:mxpa25,0:mxpa50,0:mxpa13),sigm(0:mxpa25,0:mxpa50,0:mxpa13),
720  &sign(0:mxpa25,0:mxpa50,0:mxpa13),sigo(0:mxpa25,0:mxpa50,0:mxpa13)
721  dimension xpnt(mxpa96),wght(mxpa96),
722  &ssoft(0:mxpa25),shard(0:mxpa50),strpl(0:mxpa25)
723 C - - required MXPA25 > NMAX - -
724  dimension fak(0:mxpa13),cmbin(0:mxpa13,0:mxpa13)
725  double precision
726  & expsop,expsoh,exmsop,exmsoh,exnsop,exnsoh,exosop,exosoh,
727  & exphap,exphah,exmhap,exmhah,exnhap,exnhah,exohap,exohah,
728  & exptrp,exptrh,exmtrp,exmtrh,exntrp,exntrh,exotrp,exotrh,
729  & explop,exploh,exmlop,exmloh,exnlop,exnloh,exolop,exoloh,
730  & expexh,exmexh,exnexh,exoexh,expexp,exmexp,exnexp,exoexp
731  DOUBLE PRECISION fapsof,famsof,fansof,faosof,
732  & faphar,famhar,fanhar,faohar,
733  & faptrp,famtrp,fantrp,faotrp,
734  & faploo,famloo,fanloo,faoloo
735  DOUBLE PRECISION denom,denomi,xpntk,wghtk,rmxlmn
736  & ,sigsum,siginl,sighri
737 *
738 *---------------------------------------------------------------------------
739 *
740 * externe ICON option to internes NMAX=1,2,free
741  IF(icon/10.EQ.4) nmax=2
742  IF(icon/10.EQ.5) nmax=1
743 *
744 * for safty
745  IF( nmax.GT.mxpa13) THEN
746  WRITE(6,*)' arrays limit NMAX set to' , mxpa13
747  nmax=mxpa13
748  ENDIF
749  IF( mmax.GT.mxpa50) THEN
750  WRITE(6,*)' arrays limit MMAX set to' , mxpa50
751  mmax=mxpa50
752  ENDIF
753  IF( lmax.GT.mxpa25) THEN
754  WRITE(6,*)' arrays limit LMAX set to' , mxpa25
755  lmax=mxpa25
756  ENDIF
757 *
758  lmaxi = lmax
759  mmaxi = mmax
760  IF( nmax.GE.3)THEN
761  nmaxi = nmax
762  nnmaxi=(mxpa13-nmaxi)/(1+nmaxi)
763 C aim: MXPA13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
764  nlmaxi=0
765  ELSEIF( nmax.EQ.2)THEN
766  nmaxi=1
767  nnmaxi=1
768  nlmaxi=1
769  ELSEIF( nmax.EQ.1)THEN
770  nmaxi=1
771  nnmaxi=0
772  nlmaxi=1
773  ELSEIF( nmax.LE.0)THEN
774  nmaxi=0
775  nnmaxi=0
776  nlmaxi=0
777  ENDIF
778 *
779 *
780  lentry=0
781 *
782  goto 111
783 *
784 *----------------------------------------------------------------------
785 *
786  entry sigma2(ecm)
787 *
788 *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
789 *
790  lentry=1
791 * externe ICON option to internes NMAX=1,2,free
792  IF(icon/10.EQ.4) nmax=2
793  IF(icon/10.EQ.5) nmax=1
794 * we drop L,M..dependent quantities, the rest is integrated for L=1,rest=0
795  lmaxi = 1
796  mmaxi = 0
797  nmaxi = 0
798  nnmaxi= 0
799  nlmaxi= 0
800 *
801 *----------------------------------------------------------------------
802 *
803 * *** calculate the B-space integral for L/M soft/hard cut-pom.
804 * *** NPNT is the number of integration points of B-space integ.
805 *
806  111 sigtot=0.0
807  siginl=0.0
808  sigel =0.0
809  sigele=0.0
810  sigd =0.0
811  sigdd =0.0
812  sigdi =0.0
813  sigddi=0.0
814  sigine=0.0
815  sihmdd=0.0
816  sighmd=0.
817  siglmd=0.
818  sighin=0.
819  sigin=0.
820  sigsin=0.0
821  sighri=0.0
822  sigsum=0.0
823  sigsme=0.0
824  silmsd=0.0
825  silmdd=0.
826  slhmdd=0.
827  DO 10 l=0,lmaxi
828  ssoft(l)=0.
829  strpl(l)=0.
830  DO 10 m=0,mmaxi
831  shard(m)=0.
832  DO 10 n=0, mxpa13
833  sig(l,m,n)=0.
834  sigp(l,m,n)=0.
835  sigm(l,m,n)=0.
836  sign(l,m,n)=0.
837  sigo(l,m,n)=0.
838  plmn(l,m,n)=0.
839  plmncu(l,m,n)=0.
840 10 CONTINUE
841 *
842 * get bare X-sections
843  s=ecm**2
844  CALL sigshd(ecm)
845 *
846  IF(alalam.LE.1.d-2) THEN
847  alam=0.60
848  ELSE
849  alam=alalam
850  ENDIF
851 *
852 * prepare integration:
853  pi4 = 4.*pi
854 *
855  IF(ecm.LT.2000.d0)THEN
856  npnt=96
857  CALL gset(zero,one,npnt,xpnt,wght)
858  ELSE
859  npnt=mxpa96
860  CALL gset(zero,one,npnt,xpnt,wght)
861  ENDIF
862 *
863 * as low mass diffraction extra, high mass reduced:
864  redu = 1.0
865  IF(ioutpo.GE.0) WRITE (6,*) ' ALAM,REDU= ',alam,redu
866 *
867 * --here the versions started--
868 * prepare factors to enter sum:
869 * notation: Z(HARd/SOFt/LOOp)(Plus/Minus/N=mixed/O=mixed)
870 *
871  zharp=(1.+alam)**2*zhar
872  zsofp=(1.+alam)**2*zsof
873  zloop=(1.+alam)**2*zloo * redu
874  zharm=(1.-alam)**2*zhar
875  zsofm=(1.-alam)**2*zsof
876  zloom=(1.-alam)**2*zloo * redu
877  zharn=(1.-alam**2)*zhar
878  zsofn=(1.-alam**2)*zsof
879  zloon=(1.-alam**2)*zloo * redu
880  zharo=(1.-alam**2)*zhar
881  zsofo=(1.-alam**2)*zsof
882  zlooo=(1.-alam**2)*zloo * redu
883 * one more factor at the top or at the bottom:
884  ztrpp=(1.+alam)**3*ztrp * redu
885  ztrpm=(1.-alam)**3*ztrp * redu
886  ztrpn=(1.-alam**2)*(1.+alam)*ztrp * redu
887  ztrpo=(1.-alam**2)*(1.-alam)*ztrp * redu
888 *
889 * begin M,N,L,LL loop
890 *
891  DO 720 l=0,lmaxi
892  IF(l.EQ.0) THEN
893  fapsof=1.
894  famsof=1.
895  fansof=1.
896  faosof=1.
897  ELSEIF(lsqrt) THEN
898  fapsof=fapsof* sqrt( zsofp/float(l))
899  famsof=famsof* sqrt( zsofm/float(l))
900  fansof=fansof* sqrt( zsofn/float(l))
901  faosof=faosof* sqrt( zsofo/float(l))
902  IF ( fapsof .LT.dtiny ) fapsof=0.
903  IF ( famsof .LT.dtiny ) famsof=0.
904  IF ( fansof .LT.dtiny ) fansof=0.
905  IF ( faosof .LT.dtiny ) faosof=0.
906  ELSEIF(.NOT.lsqrt) THEN
907  fapsof=fapsof*zsofp/float(l)
908  famsof=famsof*zsofm/float(l)
909  fansof=fansof*zsofn/float(l)
910  faosof=faosof*zsofo/float(l)
911  IF (fapsof.LT.dtiny ) fapsof=0.
912  IF (famsof.LT.dtiny ) famsof=0.
913  IF (fansof.LT.dtiny ) fansof=0.
914  IF (faosof.LT.dtiny ) faosof=0.
915  ENDIF
916  DO 730 m=0,mmaxi
917  IF(m.EQ.0) THEN
918  faphar=1.
919  famhar=1.
920  fanhar=1.
921  faohar=1.
922  ELSEIF(lsqrt) THEN
923 C WRITE(6,*)FAPHAR,ZHARP/FLOAT(M),FAMHAR,ZHARM/FLOAT(M)
924  faphar=faphar* sqrt( zharp/float(m) )
925  famhar=famhar* sqrt( zharm/float(m) )
926  fanhar=fanhar* sqrt( zharn/float(m) )
927  faohar=faohar* sqrt( zharo/float(m) )
928  IF ( fapsof*faphar .LT.dtiny ) faphar=0.
929  IF ( famsof*famhar .LT.dtiny ) famhar=0.
930  IF ( fansof*fanhar .LT.dtiny ) fanhar=0.
931  IF ( faosof*faohar .LT.dtiny ) faohar=0.
932  ELSEIF(.NOT.lsqrt) THEN
933  faphar=faphar*zharp/float(m)
934  famhar=famhar*zharm/float(m)
935  fanhar=fanhar*zharn/float(m)
936  faohar=faohar*zharo/float(m)
937  IF (fapsof*faphar.LT.dtiny ) faphar=0.
938  IF (famsof*famhar.LT.dtiny ) famhar=0.
939  IF (fansof*fanhar.LT.dtiny ) fanhar=0.
940  IF (faosof*faohar.LT.dtiny ) faohar=0.
941  ENDIF
942  DO 740 n=0,nmaxi
943  IF( n.EQ.0) THEN
944  faptrp=1.
945  famtrp=1.
946  fantrp=1.
947  faotrp=1.
948  ELSEIF(lsqrt) THEN
949  faptrp=-faptrp* sqrt( ztrpp/float(n) )
950  famtrp=-famtrp* sqrt( ztrpm/float(n) )
951  fantrp=-fantrp* sqrt( ztrpn/float(n) )
952  faotrp=-faotrp* sqrt( ztrpo/float(n) )
953  IF (abs(faptrp*fapsof*faphar).LT.dtiny ) faptrp=0.
954  IF (abs(famtrp*famsof*famhar).LT.dtiny ) famtrp=0.
955  IF (abs(fantrp*fansof*fanhar).LT.dtiny ) fantrp=0.
956  IF (abs(faotrp*faosof*faohar).LT.dtiny ) faotrp=0.
957  ELSEIF(.NOT.lsqrt) THEN
958  faptrp=-faptrp*ztrpp/float(n)
959  famtrp=-famtrp*ztrpm/float(n)
960  fantrp=-fantrp*ztrpn/float(n)
961  faotrp=-faotrp*ztrpo/float(n)
962  IF (abs(faptrp*fapsof*faphar).LT.dtiny ) faptrp=0.
963  IF (abs(famtrp*famsof*famhar).LT.dtiny ) famtrp=0.
964  IF (abs(fantrp*fansof*fanhar).LT.dtiny ) fantrp=0.
965  IF (abs(faotrp*faosof*faohar).LT.dtiny ) faotrp=0.
966  ENDIF
967  DO 750 nn=0,nnmaxi
968 * for compatibility no new subscript is introduced in some arrays:
969  nnn=n+(nmaxi+1)*nn
970 * if only first order option jump out of second order contr.:
971  IF( nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1 ) go to 750
972  IF(nn.EQ.0) THEN
973  faploo=1.
974  famloo=1.
975  fanloo=1.
976  faoloo=1.
977  ELSEIF(lsqrt) THEN
978  faploo=-faploo* sqrt( zloop/float(nn))
979  famloo=-famloo* sqrt( zloom/float(nn))
980  fanloo=-fanloo* sqrt( zloon/float(nn))
981  faoloo=-faoloo* sqrt( zlooo/float(nn))
982  IF(abs(faploo*faptrp*fapsof*faphar).LT.dtiny )faploo=0.
983  IF(abs(famloo*famtrp*famsof*famhar).LT.dtiny )famloo=0.
984  IF(abs(fanloo*fantrp*fansof*fanhar).LT.dtiny )fanloo=0.
985  IF(abs(faoloo*faotrp*faosof*faohar).LT.dtiny )faoloo=0.
986  ELSEIF(.NOT.lsqrt) THEN
987  faploo=-faploo*zloop/float(nn)
988  famloo=-famloo*zloom/float(nn)
989  fanloo=-fanloo*zloon/float(nn)
990  faoloo=-faoloo*zlooo/float(nn)
991  IF(abs(faploo*faptrp*fapsof*faphar).LT.dtiny )faploo=0.
992  IF(abs(famloo*famtrp*famsof*famhar).LT.dtiny )famloo=0.
993  IF(abs(fanloo*fantrp*fansof*fanhar).LT.dtiny )fanloo=0.
994  IF(abs(faoloo*faotrp*faosof*faohar).LT.dtiny )faoloo=0.
995  ENDIF
996 *
997 * elastic processes are not generated
998  IF(l.EQ.0.AND.m.EQ.0.AND.n.EQ.0.AND.nn.EQ.0) go to 750
999 *
1000  denom=dble(m)/dble(bh)+dble(l)/dble(bs)+dble(n)/dble(bt)
1001  & +dble(nn)/dble(bt)
1002 *
1003  DO 735 k=1,npnt
1004 *
1005 C change intergration for large L+M+N+NN?
1006  IF ( (m+l+n+nn) .LE. mxlmn ) THEN
1007  xpntk=dble(xpnt(k))
1008  wghtk=dble(wght(k))
1009  denomi=denom
1010  ELSE
1011  rmxlmn = dble(m+l+n+nn) /dble(mxlmn)
1012  xpntk=dble(xpnt(k))
1013  wghtk= dble(wght(k)) * xpntk**(rmxlmn-1.)
1014  denomi= denom / rmxlmn
1015  ENDIF
1016 *
1017  exposp=-zsofp*xpntk**(1./(denomi*dble(bs)))
1018  exposm=-zsofm*xpntk**(1./(denomi*dble(bs)))
1019  exposn=-zsofn*xpntk**(1./(denomi*dble(bs)))
1020  exposo=-zsofo*xpntk**(1./(denomi*dble(bs)))
1021 *
1022  expohp=-zharp*xpntk**(1./(denomi*dble(bh)))
1023  expohm=-zharm*xpntk**(1./(denomi*dble(bh)))
1024  expohn=-zharn*xpntk**(1./(denomi*dble(bh)))
1025  expoho=-zharo*xpntk**(1./(denomi*dble(bh)))
1026 *
1027  expotp=+ztrpp*xpntk**(1./(denomi*dble(bt)))
1028  expotm=+ztrpm*xpntk**(1./(denomi*dble(bt)))
1029  expotn=+ztrpn*xpntk**(1./(denomi*dble(bt)))
1030  expoto=+ztrpo*xpntk**(1./(denomi*dble(bt)))
1031 *
1032  expolp=+zloop*xpntk**(1./(denomi*dble(bt)))
1033  expolm=+zloom*xpntk**(1./(denomi*dble(bt)))
1034  expoln=+zloon*xpntk**(1./(denomi*dble(bt)))
1035  expolo=+zlooo*xpntk**(1./(denomi*dble(bt)))
1036 *
1037  IF(ioutpo.GE.7) THEN
1038  WRITE(6,*)
1039  * ' K=',k,' EXPOS/H=',exposp,expohp,' DENOMI/BH=',denomi,bh
1040  WRITE(6,*)
1041  * ' K=',k,' EXPOS/H=',exposm,expohm,' DENOMI/BH=',denomi,bh
1042  WRITE(6,*)
1043  * ' K=',k,' EXPOS/H=',exposn,expohn,' DENOMI/BH=',denomi,bh
1044  WRITE(6,*)
1045  * ' K=',k,'XPNT=',xpntk,'WGHT=',wghtk,'DENO=',denomi
1046  ENDIF
1047 *
1048 * notation:
1049 * EX(P=+/M=-/N=mixed/O=mixed)(EX=all/HArd,TRippel,LOop)(P/Half)
1050 *
1051  IF( exposp .GT. tinexp) THEN
1052  expsoh=exp(0.5d00*exposp)
1053  exmsoh=exp(0.5d00*exposm)
1054  exnsoh=exp(0.5d00*exposn)
1055  exosoh=exp(0.5d00*exposo)
1056  ELSE
1057  expsoh=0.
1058  exmsoh=0.
1059  exnsoh=0.
1060  exosoh=0.
1061  ENDIF
1062  expsop=expsoh**2
1063  exmsop=exmsoh**2
1064  exnsop=exnsoh**2
1065  exosop=exosoh**2
1066 *
1067  IF( expohp .GT. tinexp) THEN
1068  exphah=exp(0.5d00*expohp)
1069  exmhah=exp(0.5d00*expohm)
1070  exnhah=exp(0.5d00*expohn)
1071  exohah=exp(0.5d00*expoho)
1072  ELSE
1073  exphah=0.
1074  exmhah=0.
1075  exnhah=0.
1076  exohah=0.
1077  ENDIF
1078  exphap=exphah**2
1079  exmhap=exmhah**2
1080  exnhap=exnhah**2
1081  exohap=exohah**2
1082 *
1083  IF( nmax.GE.3) THEN
1084  IF( expotp .GT. tinexp) THEN
1085  exptrh=exp(0.5d00*expotp)
1086  exmtrh=exp(0.5d00*expotm)
1087  exntrh=exp(0.5d00*expotn)
1088  exotrh=exp(0.5d00*expoto)
1089  ELSE
1090  exptrh=0.
1091  exmtrh=0.
1092  exntrh=0.
1093  exotrh=0.
1094  ENDIF
1095  exptrp= exptrh**2
1096  exmtrp= exmtrh**2
1097  exntrp= exntrh**2
1098  exotrp= exotrh**2
1099  ELSEIF( nmax.LE.2) THEN
1100  exptrh= 1 + 0.5*expotp
1101  exmtrh= 1 + 0.5*expotm
1102  exntrh= 1 + 0.5*expotn
1103  exotrh= 1 + 0.5*expoto
1104  exptrp= 1 + expotp
1105  exmtrp= 1 + expotm
1106  exntrp= 1 + expotn
1107  exotrp= 1 + expoto
1108  ENDIF
1109 *
1110  IF( nmax.GE.3) THEN
1111  IF( expolp .GT. tinexp) THEN
1112  exploh=exp(0.5d00*expolp)
1113  exmloh=exp(0.5d00*expolm)
1114  exnloh=exp(0.5d00*expoln)
1115  exoloh=exp(0.5d00*expolo)
1116  ELSE
1117  exploh=0.
1118  exmloh=0.
1119  exnloh=0.
1120  exoloh=0.
1121  ENDIF
1122  explop=exploh**2
1123  exmlop=exmloh**2
1124  exnlop=exnloh**2
1125  exolop=exoloh**2
1126  ELSEIF( nmax.EQ.2 ) THEN
1127  exploh= 1 + 0.5*expolp
1128  exmloh= 1 + 0.5*expolm
1129  exnloh= 1 + 0.5*expoln
1130  exoloh= 1 + 0.5*expolo
1131  explop= 1 + expolp
1132  exmlop= 1 + expolm
1133  exnlop= 1 + expoln
1134  exolop= 1 + expolo
1135  ELSEIF( nmax.LE.1 ) THEN
1136  exploh= 1
1137  exmloh= 1
1138  exnloh= 1
1139  exoloh= 1
1140  explop= 1
1141  exmlop= 1
1142  exnlop= 1
1143  exolop= 1
1144  ENDIF
1145 *
1146  expexh = expsoh *exphah *exptrh *exploh
1147  exmexh = exmsoh *exmhah *exmtrh *exmloh
1148  exnexh = exnsoh *exnhah *exntrh *exnloh
1149  exoexh = exosoh *exohah *exotrh *exoloh
1150  expexp = expsop *exphap *exptrp *explop
1151  exmexp = exmsop *exmhap *exmtrp *exmlop
1152  exnexp = exnsop *exnhap *exntrp *exnlop
1153  exoexp = exosop *exohap *exotrp *exolop
1154 *
1155  IF( ( nmax.LE.2 .AND. n.EQ.1 ) .OR.
1156  * ( nmax.EQ.2 .AND. nn.EQ.1 ) .OR.
1157  * nmax.EQ.0 ) THEN
1158  sigp(l,m,nnn)=sigp(l,m,nnn)+expsop *exphap *wghtk
1159  sigm(l,m,nnn)=sigm(l,m,nnn)+exmsop *exmhap *wghtk
1160  sign(l,m,nnn)=sign(l,m,nnn)+exnsop *exnhap *wghtk
1161  sigo(l,m,nnn)=sigo(l,m,nnn)+exosop *exohap *wghtk
1162  ELSE
1163  sigp(l,m,nnn)=sigp(l,m,nnn)+expexp*wghtk
1164  sigm(l,m,nnn)=sigm(l,m,nnn)+exmexp*wghtk
1165  sign(l,m,nnn)=sign(l,m,nnn)+exnexp*wghtk
1166  sigo(l,m,nnn)=sigo(l,m,nnn)+exoexp*wghtk
1167  ENDIF
1168 *
1169 * quantities without L,M,N,NN dependence considered ones
1170 * (when, chosen to get suitable weights)
1171  IF(l.EQ.1.AND.m.EQ.0.AND.n.EQ.0.AND.nn.EQ.0) THEN
1172 *
1173  IF ( (m+l+n+nn) .GT. mxlmn ) THEN
1174  WRITE(6,*)' MXLMN too low ' , mxlmn,m,l,n,nn
1175  RETURN
1176  ENDIF
1177  wghfac = wghtk/xpntk *pi4/denomi
1178  IF ( nmax.GE.3 ) THEN
1179  sigele = sigele + wghfac *
1180  * 0.0625*( 1.-expexh + 1.-exmexh
1181  * +1.-exnexh + 1.-exoexh )**2
1182 * low mass diffraction:
1183  silmsd = silmsd + wghfac *
1184  * 0.125*(expexh -exmexh)**2
1185  silmdd = silmdd + wghfac *
1186  * 0.0625*(expexh+exmexh-exnexh-exoexh)**2
1187  ELSEIF( nmax.LE.2 ) THEN
1188  sigele = sigele + wghfac *
1189  * 0.0625*( ( 1.-expexh + 1.-exmexh
1190  * +1.-exnexh + 1.-exoexh
1191 * subtract second order terms in each factor:
1192 C////CHANGED - TO + BRACKET UNNECESSARY
1193  * +(1.-exptrh)*(1-exploh) *expsoh *exphah
1194  * +(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1195  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1196  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah)**2
1197 * subtract second order terms of product:
1198  * - ( (2.-exptrh-exploh) *expsoh *exphah
1199  * +(2.-exmtrh-exmloh) *exmsoh *exmhah
1200  * +(2.-exntrh-exnloh) *exnsoh *exnhah
1201  * +(2.-exotrh-exoloh) *exosoh *exohah ) **2)
1202 * low mass diffraction:
1203  silmsd = silmsd + wghfac *
1204  * 0.125*( ( expexh -exmexh
1205 * subtract second order terms in each factor:
1206  * -(1.-exptrh)*(1-exploh) *expsoh*exphah
1207  * +(1.-exmtrh)*(1-exmloh) *exmsoh*exmhah )**2
1208 * subtract second order terms of product:
1209  * -( (2.-exptrh-exploh) *expsoh *exphah
1210  * -(2.-exmtrh-exmloh) *exmsoh*exmhah ) **2)
1211  silmdd = silmdd + wghfac *
1212  * 0.0625*( (expexh+exmexh-exnexh-exoexh
1213 * subtract second order terms in each factor:
1214  * -(1.-exptrh)*(1-exploh) *expsoh *exphah
1215  * -(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1216  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1217  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah)**2
1218 * subtract second order terms of product:
1219  * - ( (2.-exptrh-exploh) *expsoh *exphah
1220  * +(2.-exmtrh-exmloh) *exmsoh *exmhah
1221  * -(2.-exntrh-exnloh) *exnsoh *exnhah
1222  * -(2.-exotrh-exoloh) *exosoh *exohah ) **2)
1223  ENDIF
1224  IF( nmax.NE.2 ) THEN
1225  sigtot=sigtot+2.*wghfac*
1226  * 0.25*( 1.-expexh + 1.-exmexh +
1227  * 1.-exnexh + 1.-exoexh )
1228  sigine = sigine + wghfac *
1229  * 0.25*( 1.-expexp + 1.-exmexp +
1230  * 1.-exnexp + 1.-exoexp )
1231 * pure-soft-inelastic (hard scatterring is included as absorbtion)
1232  sigsin=sigsin+ wghfac *
1233  * 0.25*( (exphap-expexp)
1234  * +(exmhap-exmexp)
1235  * +(exnhap-exnexp)
1236  * +(exohap-exoexp) )
1237 * hard-inelastic (soft scatterring disregarded)
1238  sighin=sighin+ wghfac*
1239  * 0.25*( 1.-exphap + 1.-exmhap +
1240  * 1.-exnhap + 1.-exohap )
1241  ELSEIF( nmax.EQ.2 ) THEN
1242  sigtot=sigtot+2.*wghfac*
1243  * 0.25*( 1.-expexh + 1.-exmexh +
1244  * 1.-exnexh + 1.-exoexh
1245 * subtract second order terms in .TRH and LOH:
1246 C/////CHANGED - TO +
1247  * +(1.-exptrh)*(1-exploh) *expsoh *exphah
1248  * +(1.-exmtrh)*(1-exmloh) *exmsoh *exmhah
1249  * +(1.-exntrh)*(1-exnloh) *exnsoh *exnhah
1250  * +(1.-exotrh)*(1-exoloh) *exosoh *exohah )
1251  sigine = sigine + wghfac *
1252  * 0.25*( 1.-expexp + 1.-exmexp +
1253  * 1.-exnexp + 1.-exoexp
1254 * subtract second order terms in .TRP and LOP:
1255 C/////CHANGED - TO +
1256  * +(1.-exptrp)*(1-explop) *expsop *exphap
1257  * +(1.-exmtrp)*(1-exmlop) *exmsop *exmhap
1258  * +(1.-exntrp)*(1-exnlop) *exnsop *exnhap
1259  * +(1.-exotrp)*(1-exolop) *exosop *exohap )
1260 * pure-soft-inelastic (hard scatterring is included as absorbtion)
1261  sigsin=sigsin+ wghfac *
1262  * 0.25*( (exphap-expexp)
1263  * +(exmhap-exmexp)
1264  * +(exnhap-exnexp)
1265  * +(exohap-exoexp)
1266 * subtract second order terms of 2nd column:
1267  * +(1.-exptrp)*(1-explop) *expsop *exphap
1268  * +(1.-exmtrp)*(1-exmlop) *exmsop *exmhap
1269  * +(1.-exntrp)*(1-exnlop) *exnsop *exnhap
1270  * +(1.-exotrp)*(1-exolop) *exosop *exohap)
1271 * hard-inelastic (soft scatterring disregarded)
1272  sighin=sighin+ wghfac*
1273  * 0.25*( 1.-exphap + 1.-exmhap +
1274  * 1.-exnhap + 1.-exohap )
1275  ENDIF
1276 * high mass diffraction (sep.low mass diffr. arbitrary)
1277 * (naive 1/EX?TRP -> EX?TR as selected cut counts -1)
1278  IF( nmax.GE.3 ) THEN
1279  sighmd=sighmd + wghfac *
1280  * 0.25*( (exptrp-1.)*expexp
1281  * +(exmtrp-1.)*exmexp
1282  * +(exntrp-1.)*exnexp
1283  * +(exotrp-1.)*exoexp)
1284  ELSE
1285  sighmd=sighmd + wghfac *
1286  * 0.25*( expotp * expsop*exphap
1287  * +expotm * exmsop*exmhap
1288  * +expotn * exnsop*exnhap
1289  * +expoto * exosop*exohap )
1290  ENDIF
1291  IF( nmax.GE.3 ) THEN
1292  sihmdd=sihmdd + wghfac *
1293  * 0.25*( (explop-1.)*expexp
1294  * +(exmlop-1.)*exmexp
1295  * +(exnlop-1.)*exnexp
1296  * +(exolop-1.)*exoexp)
1297  ELSEIF (nmax.EQ.2 ) THEN
1298  sihmdd=sihmdd + wghfac *
1299  * 0.25*( expolp * expsop*exphap
1300  * +expolm * exmsop*exmhap
1301  * +expoln * exnsop*exnhap
1302  * +expolo * exosop*exohap )
1303 * no action:
1304 * ELSEIF (NMAX.LE.1 ) THEN
1305 * SIHMDD=SIHMDD + WGHFAC *
1306 * * 0.25*( 0. * EXPSOP*EXPHAP
1307 * * + 0. * EXMSOP*EXMHAP
1308 * * + 0. * EXNSOP*EXNHAP
1309 * * + 0. * EXOSOP*EXOHAP )
1310  ENDIF
1311  ENDIF
1312 * ending non L,M,N,NN depending part
1313 *
1314 735 CONTINUE
1315 * ending impact-integral loop
1316 *
1317  IF(abs(faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)).LT.dtiny)
1318  & THEN
1319  sigp(l,m,nnn)=0.
1320  ELSEIF(lsqrt) THEN
1321  sigp(l,m,nnn)=faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)
1322  * * abs(faphar*fapsof*faptrp*faploo)/denomi*pi4
1323  ELSEIF(.NOT.lsqrt) THEN
1324  sigp(l,m,nnn)=faphar*fapsof*faptrp*faploo*sigp(l,m,nnn)
1325  * /denomi*pi4
1326  ENDIF
1327  IF(abs(famhar*famsof*famtrp*famloo*sigm(l,m,nnn)).LT.dtiny)
1328  & THEN
1329  sigm(l,m,nnn)=0.
1330  ELSEIF(lsqrt) THEN
1331  sigm(l,m,nnn)=famhar*famsof*famtrp*famloo*sigm(l,m,nnn)
1332  * * abs( famhar*famsof*famtrp*famloo)/denomi*pi4
1333  ELSEIF(.NOT.lsqrt) THEN
1334  sigm(l,m,nnn)=famhar*famsof*famtrp*famloo*sigm(l,m,nnn)
1335  * /denomi*pi4
1336  ENDIF
1337  IF(abs(fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)).LT.dtiny)
1338  & THEN
1339  sign(l,m,nnn)=0.
1340  ELSEIF(lsqrt) THEN
1341  sign(l,m,nnn)=fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)
1342  * * abs( fanhar*fansof*fantrp*fanloo)/denomi*pi4
1343  ELSEIF(.NOT.lsqrt) THEN
1344  sign(l,m,nnn)=fanhar*fansof*fantrp*fanloo*sign(l,m,nnn)
1345  * /denomi*pi4
1346  ENDIF
1347  IF(abs(faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)).LT.dtiny)
1348  & THEN
1349  sigo(l,m,nnn)=0.
1350  ELSEIF(lsqrt) THEN
1351  sigo(l,m,nnn)=faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)
1352  * * abs( faohar*faosof*faotrp*faoloo/denomi)*pi4
1353  ELSEIF(.NOT.lsqrt) THEN
1354  sigo(l,m,nnn)=faohar*faosof*faotrp*faoloo*sigo(l,m,nnn)
1355  * /denomi*pi4
1356  ENDIF
1357 *
1358 750 CONTINUE
1359 740 CONTINUE
1360 730 CONTINUE
1361 720 CONTINUE
1362 *
1363 * *** summing up the three contributions
1364 *
1365  nnnmax=nmaxi+(nmaxi+1)*nnmaxi
1366  DO 820 l=0,lmaxi
1367  DO 830 m=0,mmaxi
1368  DO 830 nnn=0,nnnmax
1369  sig(l,m,nnn)=(sigp(l,m,nnn)+sigm(l,m,nnn)+
1370  * sign(l,m,nnn)+sigo(l,m,nnn) )/4.
1371 830 CONTINUE
1372 820 CONTINUE
1373 *
1374 *
1375 * *** calculate summed quantities for print out
1376 *
1377  DO 3 l=0,lmaxi
1378  DO 4 m=0,mmaxi
1379  DO 4 n=0,nmaxi
1380  DO 4 nn=0,nnmaxi
1381  IF( nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1 ) go to 4
1382  nnn=n+(nmaxi+1)*nn
1383  sigsum=sigsum + sig(l,m,nnn)
1384 * for options outlawing hard without soft:
1385  IF(m.EQ.0.OR.l.GE.1) sigsme=sigsme + sig(l,m,nnn)
1386  shard(m)=shard(m)+sig(l,m,nnn)
1387  ssoft(l)=ssoft(l)+sig(l,m,nnn)
1388  strpl(n)=strpl(n)+sig(l,m,nnn)
1389  siginl = siginl + sig(l,m,nnn)
1390  IF(m.GE.1) sighri = sighri + sig(l,m,nnn)
1391  IF(l.EQ.0.AND.m.EQ.0.AND.nn.EQ.0.AND.n.GE.1) THEN
1392  sigdi = sigdi + (-1)**n*sig(l,m,nnn)
1393  ELSEIF(l.EQ.0.AND.m.EQ.0.AND.n.EQ.0.AND.nn.GE.1) THEN
1394  sigddi= sigddi + (-1)**nn*sig(l,m,nnn)
1395  ENDIF
1396  4 CONTINUE
1397  3 CONTINUE
1398 * elastic processes were not generated, L,M,N,NN=0 no problem
1399 *
1400  siglmd=silmsd+silmdd
1401  sithmd=sighmd+sihmdd
1402  sigd = siglmd + sithmd
1403  slhmdd = sqrt(abs(silmdd*sihmdd))
1404  sigdd= silmdd + sihmdd + slhmdd
1405  sigin=sigine+siglmd
1406  sigel=sigtot-sigin
1407 *
1408 * *** print out
1409 *
1410  IF(lentry.EQ.1.AND.ioutpo.LE.1) RETURN
1411 *
1412  WRITE(6,*)' '
1413  WRITE(6,*)' --- properties of events ---'
1414  WRITE (6,102)
1415  WRITE(6,*)' Energy=',ecm
1416  WRITE (6,102)
1417  WRITE(6,*)' max.contributing soft/hard/diffr./doubl.diffr. cuts'
1418  WRITE(6,*)' LMAXI= MMAXI= NMAXI= NNMAXI='
1419  WRITE(6,'(15X,4I9)') lmaxi,mmaxi,nmaxi,nnmaxi
1420  WRITE(6,*)' methode used: '
1421  WRITE(6,*)' ISIG= ICON= IPIM= '
1422  WRITE(6,'(15X,3I9)') isig,icon,ipim
1423  WRITE (6,102)
1424  WRITE(6,*)' --- bare cross section and eikonal constants ---'
1425 C COMMON/POMPAR/ALFA,ALFAP,A,C,AK
1426 C COMMON /SIGMA/SIGSOF,BS,ZSOF,SIGHAR,BH,ZHAR,SIGTRP,BT,ZTRP,
1427 C * SIGLOO,ZLOO
1428  WRITE(6,*)' ALFA =',alfa,' ALFAP =',alfap,' A =',a
1429  WRITE(6,*)' C =',c,' AK =',ak
1430  WRITE(6,*)' ALALAM =',alalam
1431  WRITE (6,102)
1432  WRITE(6,*)' SIGSOF=',sigsof,' BS=',bs,' ZSOF=',zsof
1433  WRITE(6,*)' SIGHAR=',sighar,' BH=',bh,' ZHAR=',zhar
1434  WRITE(6,*)' SIGTRP=',sigtrp,' BT=',bt,' ZTRP=',ztrp
1435  WRITE(6,*)' SIGLOO=',sigloo,' BT=',bt,' ZLOO=',zloo
1436  WRITE (6,102)
1437  WRITE(6,*)' --- observable cross sections ---'
1438  WRITE (6,102)
1439  WRITE(6,*)' TOTAL X-SECTION = ',sigtot
1440  WRITE(6,*)' ELASTIC X-SECTION = ',sigele
1441  WRITE(6,*)' INELASTIC X-SECTION-LMD = ',sigine
1442  WRITE(6,*)' INELASTIC X-SECTION = ',sigin
1443  WRITE(6,*)' HARD INEL. X-SECTION = ',sighin
1444  WRITE (6,102)
1445  WRITE(6,*)' LOW MASS SING./DOUB.DIFFR.X-SECTION= ',silmsd,silmdd
1446  WRITE(6,*)' => LOW MASS TOTAL DIFFRACTIV.X-SECTION= ',siglmd
1447  WRITE(6,*)' HIGH MASS SING./DOUB.DIFFR.X-SECTION= ',sigdi,sigddi
1448  WRITE(6,*)' => HIGH MASS TOTAL DIFFRACTIV.X-SECTION= ',sithmd
1449  WRITE(6,*)' ESTIMAT.MIXED (LM+HM) DOUBL.DIFFRAC.X.SEC.= ',slhmdd
1450  WRITE(6,*)' => '
1451  WRITE(6,*)' DIFFRACTIVE X-SECTION = ',sigd
1452  WRITE(6,*)' DOUBLY DIFFRACTIVE X-SECT. =',sigdd
1453  WRITE (6,102)
1454 *
1455  IF(ioutpo.GE.0) THEN
1456  WRITE(6,*)' --- observ. x-sections, altern. calculated ---'
1457  WRITE(6,*)' ELASTIC X-SECTION = ',sigel
1458  WRITE(6,*)' INELASTIC X-SECTION-LMD = ',siginl
1459  WRITE(6,*)' HARD INEL. X-SECTION= ',sighri
1460  WRITE(6,*)' HIGH MASS SING./DOUB.DIFFR.X-SECT.=',sighmd,sihmdd
1461  WRITE(6,*)' X-SECTION FOR (L,M,N,NN)= 1000 0100 0010 0001'
1462  WRITE(6,*)' ',sig(1,0,0),sig(0,1,0)
1463  * ,sig(0,0,1),sig(0,0,2)
1464  WRITE (6,102)
1465  ENDIF
1466 *
1467  IF(ioutpo.GE.2) THEN
1468  WRITE (6,102)
1469  nnmaxp=nmaxi/2
1470  IF( nmaxi.LT.2)nnmaxp=1
1471  DO 52 n=0,nnmaxp
1472 * printout loops:
1473  DO 48 l=0,lmaxi
1474  48 WRITE(6,101)(sig(l,m,n),m=0,7)
1475  WRITE (6,102)
1476  DO 50 l=0,lmaxi
1477  50 WRITE(6,101)(sig(l,m,n),m=8,15)
1478  WRITE (6,102)
1479  WRITE(6,*)
1480  & ' # CUT-POMERON SSOFT X-SECT. SHARD X-SECT.'
1481  DO 58 l=0,lmaxi
1482  58 WRITE (6,103)l,ssoft(l),shard(l)
1483  WRITE (6,102)
1484 * printoutloop ends
1485 52 CONTINUE
1486  ENDIF
1487 *
1488 * *** attribute x-sections (SIG) for CUT objects ('s and PHI's)
1489 * to string configurations (PLMN) s: **********
1490 *
1491 C CHANGED 10.1.90 BY J.R.
1492 * preparations:
1493 C just for Y and (Phi) - cuts:
1494  fak(0)=1
1495  DO 500 i=1,nmaxi
1496  fak(i)=fak(i-1)*i
1497  500 CONTINUE
1498  DO 501 i=0,nmaxi
1499  DO 501 j=0,i
1500  cmbin(i,j)=fak(i)/(fak(j)*fak(i-j))
1501  501 CONTINUE
1502 *
1503  tmmp=0.
1504  DO 5 l=0,lmaxi
1505  DO 5 m=0,mmaxi
1506  IF(icon.EQ.44.OR.icon.EQ.46.OR.icon.EQ.48.
1507  * or.icon.EQ.54) THEN
1508 C///test:
1509 * no Y or PHI cut
1510  plmntm=sig(l,m,0)/(sigsum+tin)
1511  plmn(l,m,0) = plmntm + plmn(l,m,0)
1512  tmmp=tmmp+plmntm
1513 * Y but no PHI cut
1514  plmntm=sig(l,m,1)/(sigsum+tin)
1515  tmmp=tmmp+plmntm
1516  IF(l+2.LE.lmaxi) THEN
1517  plmn(l+2,m,0) = (-2.)* plmntm + plmn(l+2,m,0)
1518  plmn(l+1,m,0) = 4. * plmntm + plmn(l+1,m,0)
1519  ELSE
1520  plmn(lmaxi,m,0) = (-2.)* plmntm + plmn(lmaxi,m,0)
1521  plmn(lmaxi,m,0) = 4. * plmntm + plmn(lmaxi,m,0)
1522  ENDIF
1523  IF(l.EQ.0 .AND. m.EQ.0) THEN
1524  plmn(l ,m,1) = (-1.)* plmntm + plmn(l ,m,1)
1525  ELSE
1526  plmn(l ,m,0) = (-1.)* plmntm + plmn(l ,m,0)
1527  ENDIF
1528 * no Y but PHI cut
1529  plmntm=sig(l,m,2)/(sigsum+tin)
1530  tmmp=tmmp+plmntm
1531  IF(l+2.LE.lmaxi) THEN
1532  plmn(l+2,m,0) = (-2.)* plmntm + plmn(l+2,m,0)
1533  plmn(l+1,m,0) = 4. * plmntm + plmn(l+1,m,0)
1534  ELSE
1535  plmn(lmaxi,m,0) = (-2.)* plmntm + plmn(lmaxi,m,0)
1536  plmn(lmaxi,m,0) = 4. * plmntm + plmn(lmaxi,m,0)
1537  ENDIF
1538  IF(l.EQ.0 .AND. m.EQ.0) THEN
1539  plmn(l ,m,2) = (-1.)* plmntm + plmn(l ,m,2)
1540  ELSE
1541  plmn(l ,m,0) = (-1.)* plmntm + plmn(l ,m,0)
1542  ENDIF
1543 C/// test end
1544  ELSE
1545  DO 51 n=0,nmaxi
1546  DO 51 nn=0,nnmaxi
1547  IF(nmax.LE.2 .AND. n.EQ.1 .AND. nn.EQ.1) go to 51
1548  nnn=n+(nmaxi+1)*nn
1549 *
1550 * to be attributed::
1551  plmntm=sig(l,m,nnn)/(sigsum+tin)
1552  tmmp=tmmp+plmntm
1553 *
1554 * attribution loop for Y-cuts
1555  DO 511 n0cut=0,n
1556  DO 511 n1cut=0,n-n0cut
1557  n2cut=n-n0cut-n1cut
1558 * combinatoric weight:
1559  cmb0=cmbin(n,n2cut)
1560  cmb1=cmbin(n-n2cut,n1cut)
1561 *
1562 * attribution loop for PHI-cuts
1563  DO 511 nn0cut=0,nn
1564  DO 511 nn1cut=0,nn-nn0cut
1565  nn2cut=nn-nn0cut-nn1cut
1566 * combinatoric weight:
1567  cmbn0=cmbin(nn,nn2cut)
1568  cmbn1=cmbin(nn-nn2cut,nn1cut)
1569 *
1570 * attributions matrix:
1571 * ("L"soft,"M"hard,"N"1-diffr.,"NN"2-dif.,"NL"diffr.long in.part ):
1572 * obviously:
1573  l2str = l
1574  m2str = m
1575  n2str = n0cut
1576  nn2str= nn0cut
1577  nl2str= 0
1578 * specialties for Y's and PHI's:
1579  l2str=l2str + n1cut + nn1cut + n2cut + nn2cut
1580  IF(nmax.LE.2)THEN
1581 * room to have NL2STR's:
1582  nl2str= n2cut + nn2cut
1583  ELSEIF(nmax.GE.3)THEN
1584 * the extra inner piece counts here like long strings:
1585  l2str=l2str+n2cut+nn2cut
1586  ENDIF
1587  IF((icon.EQ.26.OR.icon.EQ.36.OR.icon.EQ.46.OR.icon.EQ.56)
1588  & .AND. (l2str.GE.1.OR.m2str.GE.1))THEN
1589  l2str=l2str + nl2str
1590  n2str = 0
1591  nn2str = 0
1592  nl2str = 0
1593  ENDIF
1594 *
1595 * getting and checking parameter for storing:
1596  IF(l2str.GT.lmaxi) l2str=lmaxi
1597  IF(m2str.GT.lmaxi) m2str=lmaxi
1598  nnnstr =n2str +(nmaxi+1)*nn2str
1599  * +(nnmaxi+1)*(nmaxi+1)*nl2str
1600  IF(nnnstr.GT.mxpa13) nnnstr=mxpa13
1601 *
1602 * summing contributions
1603  plmn(l2str,m2str,nnnstr) = plmntm
1604  * *cmb0*cmb1 * (-2)**n2cut * (4)**n1cut * (-1)**n0cut
1605  * *cmbn0*cmbn1*(-2)**nn2cut* (4)**nn1cut* (-1)**nn0cut
1606  & + plmn(l2str,m2str,nnnstr)
1607 *
1608  511 CONTINUE
1609  51 CONTINUE
1610  ENDIF
1611 C///// initial general methode ends
1612  5 CONTINUE
1613  IF(abs(tmmp-1.d0).GT..03d0)THEN
1614  WRITE(6,*)
1615  & ' NORMALISATION ERROR SUM PLM before LMD reatribution=',tmmp
1616  ENDIF
1617 *
1618 * *** built in low mass diffraction, get averages and check normalisation:
1619 *
1620 * low mass diffraction was SIGLMD=SILMSD+SILMDD
1621 * and mixed LM/HM diffraction SLHMDD=SQRT(SILMDD*SIHMDD)
1622  plmfac= (sigsum+tin) / (sigsum+tin +siglmd)
1623  plmn(0,0,1)= plmn(0,0,1) +
1624  & ( silmsd - slhmdd ) / (sigsum+tin)
1625  plmn(0,0,2)= plmn(0,0,2) +
1626  & ( silmdd + slhmdd ) / (sigsum+tin)
1627  661 CONTINUE
1628 * AVerage_SOft_N,AVerage_HaRD_N,SUM_over_Pl
1629  avsofn=0.
1630  avharn=0.
1631  avdifn=0.
1632  avddfn=0.
1633  avdlfn=0.
1634  psoft=0.
1635  temp=0.
1636  tmp=0.
1637  tmmp=0.
1638  tmmp1=0.
1639 *
1640 * (L,M,N,NN,NL repl. L2STR,M2... for s.,h.,Y-dif.,PHI-dif.,i.Y-dif.str.#)
1641  DO 6 nl=0,nlmaxi
1642  DO 6 nn=0,nnmaxi
1643  DO 6 n=0,nmaxi
1644  IF(nmax.LE.2 .AND. n+nn+nl.GE.2) go to 6
1645  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1646  DO 63 m=0,mmaxi
1647  DO 63 l=0,lmaxi
1648  IF(nl.EQ.0)tmmp1 = tmmp1 + sig(l,m,nnn)
1649  tmmp = tmmp + sig(l,m,nnn)
1650  plmn(l,m,nnn)=plmn(l,m,nnn) * plmfac
1651  tmp = tmp + plmn(l,m,nnn)
1652 C IF(PLMN(L,M,NNN).LT.-.000001D0)
1653  IF(plmn(l,m,nnn).LT.-.000005d0)
1654  & WRITE(6,*)' 0>PLMN',plmn(l,m,nnn),l,m,n,nn,nl
1655  avsofn=avsofn+plmn(l,m,nnn)*l
1656  avharn=avharn+plmn(l,m,nnn)*m
1657  avdifn=avdifn+plmn(l,m,nnn)*n
1658  avddfn=avddfn+plmn(l,m,nnn)*nn
1659  avdlfn=avdlfn+plmn(l,m,nnn)*nl
1660  IF (m.EQ.0)psoft=psoft+plmn(l,m,nnn)
1661  63 CONTINUE
1662  6 CONTINUE
1663  IF(abs(tmp-1.d0).GT..01d0)THEN
1664  WRITE(6,*)
1665  & ' NORMALISATION ERROR SUM PLM before M reatribution=',tmp
1666  ENDIF
1667  tmmp=tmmp/sigsum
1668  tmmp1=tmmp1/sigsum
1669  IF(abs(tmmp-1.d0).GT..01d0 .OR.abs(tmmp1-1.d0).GT..01d0)THEN
1670  WRITE(6,*)
1671  & ' NORMALISATION ERROR TMMP,TMMP1=',tmmp,tmmp1
1672  ENDIF
1673 *
1674 * *** reattribute purely hard scattering and get cummulant distribution
1675 *
1676 * (L,M,N,NN,NL repl. L2STR,M2...
1677 * for soft,hard,Y-diffr.,PHI-diffr.,inner diffr. string #)
1678  DO 61 nl=0,nlmaxi
1679  DO 61 nn=0,nnmaxi
1680  DO 61 n=0,nmaxi
1681  IF(nmax.LE.2 .AND. n+nn+nl.GE.2) go to 61
1682  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1683  DO 612 m=0,mmaxi
1684  DO 611 l=0,lmaxi
1685 * -- hard: a pure hard scattering gets an extra soft chain, it is considered
1686 * a specialty of fragmentation and therefor implemented only here
1687 * there are are number of other options which are dropped as they
1688 * are hard to implement at this point
1689  IF (l.EQ.0.AND.m.GE.1)THEN
1690  plmn(1,m,nnn)=plmn(1,m,nnn)+plmn(0,m,nnn)
1691  plmn(0,m,nnn)=0.
1692  ENDIF
1693 * -- cummulant of distribution
1694  temp = temp + plmn(l,m,nnn)
1695  plmncu(l,m,nnn)=temp
1696  611 CONTINUE
1697 *
1698  IF(ioutpo.GE.3)WRITE (6,*)' M,(L,PLMN(L,M,N),L=0,LMAX)'
1699  IF(ioutpo.GE.3)WRITE (6,106) m,(l,plmn(l,m,n),l=0,lmaxi)
1700  IF(ioutpo.GE.2)WRITE (6,*)' M,(L,PLMNCU(L,M,N),L=0,LMAX/2)'
1701  IF(ioutpo.GE.2)WRITE (6,106) m,(l,plmncu(l,m,n),l=0,lmaxi/2)
1702  106 FORMAT (i3,9(i3,e11.2))
1703 *
1704  612 CONTINUE
1705  61 CONTINUE
1706 *
1707  IF(abs(temp-1.d0).GT..01d0)THEN
1708  WRITE(6,*)' NORMALISATION ERROR SUM PLM=',temp
1709  plmfac=1./(temp+tin)
1710  go to 661
1711  ENDIF
1712 *
1713  IF(ioutpo.GE.1)WRITE (6,*)
1714  & '(((L,M,N,PLMN(L,M,N),N=0,2),M=0,5),L=0,7)'
1715  IF(ioutpo.GE.1)WRITE (6,1106)
1716  & (((l,m,n,plmn(l,m,n),n=0,2),m=0,5),l=0,7)
1717  IF(ioutpo.GE.1)WRITE (6,*)
1718  & '(((L,M,N,SIG(L,M,N),N=0,2),M=0,5),L=0,7)'
1719  IF(ioutpo.GE.1)WRITE (6,1106)
1720  & (((l,m,n,sig(l,m,n),n=0,2),m=0,5),l=0,7)
1721  1106 FORMAT (1x,3(i5,i5,i5,g12.5))
1722 *
1723  phard=1.-psoft
1724  alfah=sighin/(sigine+0.00001)
1725  betah=1.-alfah
1726  WRITE(6,116)avsofn,avharn,avdifn,avddfn,avdlfn,
1727  & phard,psoft,alfah,betah
1728  116 FORMAT(/'--- various averages:'/
1729  & /' AVSOFN= AVHARN= AVDIFN= AVDDFN= AVDLFN='
1730  & /' ',5f11.3
1731  & /' PHARD= PSOFT= ALFAH= BETAH= '
1732  & /' ',4f11.3)
1733  IF(ioutpo.GE.1)WRITE(6,*)'SIGSUM=SIGINL-LMD',sigsum
1734 *
1735  IF(ioutpo.GE.1)WRITE(6,610) sigtot,sigine,sigd,sigdd,sighin
1736  610 FORMAT (' SIGTOT,SIGINE,SIGD,SIGDD,SIGHIN= '/' ',5e18.6)
1737 *
1738 101 FORMAT(' ',10e10.3)
1739 102 FORMAT(' ')
1740 103 FORMAT(' ',5x,i4,5x,2e15.3)
1741 *
1742  RETURN
1743  END
1744 * ende problm
1745 *
1746 ************************************************************************
1747 *
1748 *
1749  SUBROUTINE samplx(L2STR,M2STR,N2STR,NN2STR,NL2STR)
1750 *
1751 * input:
1752 * PLMNCU
1753 * output:
1754 * samples number of soft (L) and hard (M) cut pomerons from PLMNC
1755 * and of (N=0/1/2) diffractive excitations (for L=M=0
1756 *
1757 *----------------------------------------------------------------------
1758  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1759  SAVE
1760  COMMON /nncms/ gamcm,bgcm,umo,pcm,eproj,pproj
1761 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20),NESTEP
1762  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
1763  parameter(mxpu50=100,mxpu51=mxpu50+1)
1764  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
1765 * PARAMETRIZATION FOR PTMIN= 3. GEV
1766  parameter(mxpa50=250,mxpa51=mxpa50+1)
1767 * PARAMETRIZATION FOR PTMIN= 2. GEV
1768 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
1769 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
1770  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
1771  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
1772 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
1773  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
1774  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
1775 C COMMON /POLMN1/ PLMNEE(0:MXPA25,0:MXPU50,0:MXPA13,20)
1776  COMMON /polmn1/ plmnee(0:mxpa25,0:mxpu50,0:mxpa13,28)
1777  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
1778  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
1779 *
1780  parameter(pi=3.141592654d0)
1781  DATA nprint/0/
1782 C Determine the energy index
1783  ipoen=1
1784  DO 20 ii=1,nestep
1785  IF(umo.GE.poen1(ii).AND.umo.LT.poen2(ii))THEN
1786  ipoen=ii
1787  go to 22
1788  ENDIF
1789  20 CONTINUE
1790  22 CONTINUE
1791 * "L"SOFT,"M"HARD
1792  lmaxi = lmax
1793  mmaxi = mmax
1794  IF(ipim.NE.2) THEN
1795 * "N"diffr.,"NN"dou.dif.,("NL"long inner contr. for 2*cut Y or PHI)
1796  nmaxi = nmax
1797  nnmaxi=0
1798  nlmaxi=0
1799  ELSEIF(ipim.EQ.2) THEN
1800  IF( nmax.GE.3)THEN
1801  nmaxi = nmax
1802  nnmaxi=(13-nmaxi)/(1+nmaxi)
1803 * 13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
1804  nlmaxi=0
1805  ELSEIF( nmax.EQ.2)THEN
1806  nmaxi=1
1807  nnmaxi=1
1808  nlmaxi=1
1809  ELSEIF( nmax.EQ.1)THEN
1810  nmaxi=1
1811  nnmaxi=0
1812  nlmaxi=1
1813  ENDIF
1814  ENDIF
1815  111 CONTINUE
1816 *
1817  x=rndm(v)
1818 *
1819  IF (x.LE.plmncu(0,0,0) .AND. nprint.LT.100)THEN
1820  WRITE(6,*) ' No generator of elastic events '
1821  WRITE(6,*) ' PLMNCU (0,0,0) =!= 0 = ',plmncu(0,0,0)
1822  nprint=nprint+1
1823  goto 111
1824  ENDIF
1825 *
1826  DO 5 nl=0,nlmaxi
1827  DO 5 nn=0,nnmaxi
1828  DO 5 n=0,nmaxi
1829  nnn =n +(nmaxi+1)*nn +(nnmaxi+1)*(nmaxi+1)* nl
1830  DO 6 m=0,mmaxi
1831  DO 7 l=0,lmaxi
1832 *
1833 C IF (X.LE.PLMNCU(L,M,NNN)) THEN
1834  IF (x.LE.plmnee(l,m,nnn,ipoen)) THEN
1835  l2str=l
1836  m2str=m
1837  n2str=n
1838  nn2str=nn
1839  nl2str=nl
1840  RETURN
1841 *
1842  ENDIF
1843  7 CONTINUE
1844  6 CONTINUE
1845  5 CONTINUE
1846 *
1847  nprint=nprint+1
1848  IF(nprint.LT.100) WRITE(6,*)' RAR.IN SAMPLM,PLMNCU,RND=',
1849  & plmncu(lmax, mmax,nnn),x,nprint
1850  IF( plmncu(lmax,mmax,nnn) .GT. 0.1d0 ) RETURN
1851  IF( plmncu(lmax,0,0) .GT. 0.1d0 ) RETURN
1852  WRITE(6,*)' RAR.IN SAMPLM- PROBLEM SEEMS BAD, DECIDE TO STOP'
1853  stop
1854  END
1855 *
1856 *
1857 ************************************************************************
1858 *
1859  SUBROUTINE samplm(L2STR,M2STR,N2STR)
1860 *
1861 * input:
1862 * PLMNCU
1863 * output:
1864 * samples number of soft (L) and hard (M) cut pomerons from PLMNC
1865 * and of (N=0/1/2) diffractive excitations (for L=M=0
1866 *
1867 *----------------------------------------------------------------------
1868  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1869  SAVE
1870  COMMON /nncms/ gamcm,bgcm,umo,pcm,eproj,pproj
1871 C COMMON /POMENE/POEN(20),POEN1(20),POEN2(20),NESTEP
1872  COMMON /pomene/poen(28),poen1(28),poen2(28),nestep
1873  parameter(mxpu50=100,mxpu51=mxpu50+1)
1874 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
1875  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
1876  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
1877 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
1878  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
1879 * PARAMETRIZATION FOR PTMIN= 3. GEV
1880  parameter(mxpa50=250,mxpa51=mxpa50+1)
1881 * PARAMETRIZATION FOR PTMIN= 2. GEV
1882 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
1883  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
1884  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
1885 C COMMON /POLMN1/ PLMNEE(0:MXPA25,0:MXPU50,0:MXPA13,20)
1886  COMMON /polmn1/ plmnee(0:mxpa25,0:mxpu50,0:mxpa13,28)
1887  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
1888  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
1889 *
1890  parameter(pi=3.141592654d0)
1891  ipoen=1
1892  DO 20 ii=1,nestep
1893  IF(umo.GE.poen1(ii).AND.umo.LT.poen2(ii))THEN
1894  ipoen=ii
1895  go to 22
1896  ENDIF
1897  20 CONTINUE
1898  22 CONTINUE
1899  111 CONTINUE
1900 *
1901  x=rndm(v)
1902 *
1903  IF (x.LE.plmncu(0,0,0))THEN
1904  WRITE(6,*) ' No generator of elastic events '
1905  WRITE(6,*) ' PLMNCU (0,0,0) =!= 0 = ',plmncu(0,0,0)
1906  goto 111
1907  ENDIF
1908 *
1909  DO 5 n=0,nmax
1910 C IF (N.GT.1) GO TO 111
1911  DO 6 m=0,mmax
1912  DO 7 l=0,lmax
1913 *
1914 C IF (X.LE.PLMNCU(L,M,N)) THEN
1915  IF (x.LE.plmnee(l,m,n,ipoen)) THEN
1916  l2str=l
1917  m2str=m
1918  n2str=n
1919  RETURN
1920 *
1921  ENDIF
1922  7 CONTINUE
1923  6 CONTINUE
1924  5 CONTINUE
1925 *
1926  WRITE(6,*)' RAR.IN SAMPLM,PLMNCU,RND=',plmncu(lmax,mmax,nmax),x
1927  IF( plmncu(lmax,mmax,nmax) .GT. 0.1d0 ) RETURN
1928  IF( plmncu(lmax,0,0) .GT. 0.1d0 ) RETURN
1929  WRITE(6,*)' RAR.IN SAMPLM- PROBLEM SEEMS BAD, DECIDE TO STOP'
1930  stop
1931  END
1932 *
1933 *
1934 *
1935 *
1936 ************************************************************************
1937 C--------------------------------------------------------------------
1938 C
1939 C dtUpom9h.for
1940 C
1941 C---------------------------------------------------------------------
1942 ************************************************************************
1943 *
1944 * POMDI,SIGMAS,SIGSHD,PRBLM0..9,SAMPLM,SIGMA1
1945 *
1946 * routines called by code word
1947 * SIGMAPOM :
1948 *
1949 
1950 C--------------------------------------------------------------------
1951 C
1952 C dtUpom9h.for
1953 C
1954 C---------------------------------------------------------------------
1955 ************************************************************************
1956 *
1957 * POMDI,SIGMAS,SIGSHD,PRBLM0..9,SAMPLX,SIGMA1
1958 *
1959 * routines called by code word
1960 * SIGMAPOM :
1961 *
1962 ************************************************************************
1963 *
1964 * J.RANFT September 1987
1965  SUBROUTINE pomdi
1966 *
1967 * to calculate the s-dependent X-sections
1968 *
1969 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1970 *
1971 * input:
1972 * ISIG characterizing X sections, transmitted to SIGSHD
1973 * IPIM characterizes the method to calculate SIGMA(LSOFT,mhard)
1974 * IPIM=1 : integral method with low mass dif.matr
1975 * IPIM=2 : int.meth.with low mass dif.matr.(2*2)
1976 * and reduced high mass diffraction
1977 * IPIM=3 : integral methode
1978 * IPIM=4 : integral methode with Y -cuts
1979 * IPIM=5 : integral methode with Y-cuts + 2 CHANNEL EIK.
1980 * rest : not implemented or
1981 * special cases for checking
1982 * LMAX < MXPA25 maximal number of considered soft pomerons
1983 * MMAX < MXPA50 maximal number of considered hard pomerons
1984 * NMAX < MXPA13 maximal number of considered trippel pomerons
1985 * (not used in low masss diffraction formalism)
1986 * output: printet and plottet in SIGMAS and
1987 * printed on the end of this subroutine
1988 *
1989 * card XSECTION causes call to POMDI
1990 * card SIGMAPOM with ITEST=1 causes call to POMDI
1991 * POMDI calling
1992 * SIGMAS ( SIGMA1..3( SIGSD, (1:/3:)GSET ), PLOT)
1993 * PRBLM.. ( SIGSHD, (1:)GSET )
1994 * (2:)SAMPLX, (1:/3:)SAMPLM
1995 *
1996 *----------------------------------------------------------------------
1997  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1998  SAVE
1999 *
2000 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
2001  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
2002 * PARAMETRIZATION FOR PTMIN= 3. GEV
2003  parameter(mxpa50=250,mxpa51=mxpa50+1)
2004 * PARAMETRIZATION FOR PTMIN= 2. GEV
2005 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
2006  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
2007  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
2008  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
2009  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
2010 * *** /HISTOO/
2011  INTEGER*2 ndislm
2012  COMMON /histoo/as(50,9),aecm(50,9),asig(50,9),alos(50,9),
2013  * aloecm(50,9),ndislm(0:mxpa25,0:mxpa50,0:mxpa13)
2014 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
2015  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
2016 * --- used only in SIGMAPOM-routines
2017 * *** /POMPAR/ contains pomeron parameters used in some options
2018 * ALFA,ALFAP,A,BH,C,BS are soft Pomeron parameters chosen in SIGMAS
2019 * appearing in SIGMAS, PRBLM.., AK is K-factor
2020  common/pompar/alfa,alfap,a,c,ak
2021 * *** /SIGMA/ contains variables actually used in iteration
2022 * ZSOF, ZHAR, BS, BH, SIGSOF,SIGHAR ( for soft, hard, trippel-Pom.)
2023 * are input for X-section calculated in SIGSHD
2024 C (/sigma/ is out of P.A.'s CM88, containing variables used in iter.
2025  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
2026  * sigloo,zloo
2027 * *** /POMTYP/ contains parameters determining X-sections
2028 * IPIM,ISIG,LMAX,MMAX,NMAX as described at "CODEWD = SIGMAPOM"
2029 C LMAX,MMAX,NMAX kept open to enable avoiding of numerical problem
2030  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
2031  COMMON /alala/alalam
2032 C COMMON/COLLIS/SS,IJPROJ,IJTAR,PTTHR,PTTHR2,IOPHRD,IJPRLU,IJTALU
2033 C
2034  COMMON /collpo/s,ptthr,ptthr2
2035 C COMMON /COLLIS/SS,IJPROJ,IJTAR,PTTPO,IOPHRD,IJPRLU,IJTALU,PTTPO2
2036  COMMON /collis/ss,ijproj,ijtar,pttpo,pttpo2,iophrd,ijprlu,ijtalu
2037 C
2038 C ECM calculated as SQRT
2039 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2040  parameter(pi=3.141592654d0)
2041  lmaxi = lmax
2042  mmaxi = mmax
2043  nmaxi = nmax
2044  skeep=ss
2045 *
2046 *---------------------------------------------------------------------
2047 * CALL SIGMAS
2048 *
2049 * calculates and plotts out X-section at various energies
2050 * runs thru energies independent of the rest
2051 * and contains no preperation for other work
2052 *
2053  WRITE(6,'(1X/1X)')
2054  WRITE(6,*)' '
2055  WRITE(6,*)
2056  *' ------ testing the energy dependence of x-sections ----------'
2057  WRITE(6,*)' '
2058  IF(ioutpo.GT.-1) WRITE(6,*)
2059  *' (as function of ALAM i.e.a low mass diffr.parameter)'
2060  WRITE(6,*)' -----------------------------------------------'
2061  WRITE(6,'(1X)')
2062 *
2063  DO 1007 iijj=1,10
2064  IF(ioutpo.GT.-1 .OR. iijj.EQ.6)THEN
2065 C WITHOUT SIMGMAS FOR NORMAL USERS
2066 C IF(IOUTPO.GT.-1 )THEN
2067  alalam=iijj*0.1
2068  IF(ioutpo.GT.-1) WRITE(6,1008)alalam
2069  1008 FORMAT (' ALAM= ',f10.3)
2070  CALL sigmas
2071  ENDIF
2072  1007 CONTINUE
2073 *
2074 *---------------------------------------------------------------------
2075 * test sampling L,M
2076 *
2077  1111 CONTINUE
2078 * looping thru the energies
2079 *
2080 C DO 100 III=3,9
2081 C S=10.**III
2082 C ECM=SQRT(S)
2083 *
2084  s=skeep
2085  ecm=sqrt(s)
2086 *
2087 * chosing one options with a standard call
2088 *
2089  IF(ipim.EQ.2) THEN
2090  IF( nmax.GE.3)THEN
2091  nnmaxi=(13-nmaxi)/(1+nmaxi)
2092 * 13 =!= NNNMAX = NMAXI+(NMAXI+1)*NNMAXI
2093  nlmaxi=0
2094  ELSEIF( nmax.EQ.2)THEN
2095  nmaxi=1
2096  nnmaxi=1
2097  nlmaxi=1
2098  ELSEIF( nmax.EQ.1)THEN
2099  nmaxi=1
2100  nnmaxi=0
2101  nlmaxi=1
2102  ENDIF
2103  CALL prblm2(ecm)
2104  ENDIF
2105  IF(ipim.LT.1.AND.ipim.GT.9)THEN
2106  WRITE(6,*) 'RETURN caused by IPIM=',ipim
2107  RETURN
2108  ENDIF
2109 *
2110 * randomly sample events and printout
2111 *
2112  WRITE (6,'(1X)' )
2113  WRITE (6,102)ecm,s
2114  102 FORMAT
2115  * ('--- sample distribution for L soft and M hard inelastic'
2116  * , ' pomerons (string pairs)--- '
2117  * / 20x,'at ECM = ',f10.2,' S = ',f12.1)
2118  DO 31 l=0,lmaxi
2119  DO 32 m=0,mmaxi
2120  DO 32 n=0,13
2121  ndislm(l,m,n)=0
2122  32 CONTINUE
2123  31 CONTINUE
2124 *
2125  IF(icon.EQ.12)go to 100
2126  DO 320 ii=1,10000
2127  IF(ipim.EQ.2) THEN
2128  CALL samplx(l2str,m2str,n2str,nn2str,nl2str)
2129  nnnstr =n2str +(nmaxi+1)*nn2str
2130  * +(nnmaxi+1)*(nmaxi+1)*nl2str
2131  ndislm(l2str,m2str,nnnstr)=ndislm(l2str,m2str,nnnstr)+1
2132  ELSE
2133  CALL samplm(l2str,m2str,n2str)
2134  ndislm(l2str,m2str,n2str)=ndislm(l2str,m2str,n2str)+1
2135  ENDIF
2136  320 CONTINUE
2137 *
2138  WRITE(6,*)
2139  * ' with no diffractive contribution'
2140  WRITE(6,*) ' '
2141  WRITE(6,*)
2142  * ' ....... vertical: NSTR, horizontal MSTR .........'
2143  DO 3344 l=0,min(20,lmaxi)
2144  3344 WRITE(6,34)l,(ndislm(l,m,0),m=0,20)
2145  WRITE(6,*) ' '
2146  IF(ioutpo.GE.0)THEN
2147  DO 333 n=0,5
2148  IF(n.NE.0)THEN
2149  WRITE(6,*)' WITH NSTR=',n
2150  DO 334 l=0,min(20,lmaxi)
2151  WRITE(6,34)l,(ndislm(l,m,n),m=0,20)
2152  334 CONTINUE
2153  WRITE(6,*) ' '
2154  ENDIF
2155  jmpa50 = int(mxpa50/25)
2156 C WRITE(6,*) 'WIDE PLOT 0<L<25, 0<M<'
2157  WRITE(6,*) 'WIDE PLOT 0<L<',mxpa25,' 0<M<'
2158  & ,mxpa50,' IN STEPS OF ',jmpa50
2159 C DO 335 L=0,MIN(25,LMAXI)
2160  DO 335 l=0,mxpa25
2161  WRITE(6,35)l,(ndislm(l,m,n),m=0,mxpa50,jmpa50)
2162  335 CONTINUE
2163  WRITE(6,*) ' '
2164  333 CONTINUE
2165  ENDIF
2166  34 FORMAT (i5,':',21i4)
2167  35 FORMAT (i5,26i4)
2168  100 CONTINUE
2169 * energy loop has ended
2170  RETURN
2171  END
2172 *
2173 *
2174 ************************************************************************
2175 *
2176  SUBROUTINE sigmas
2177 *
2178 * output:
2179 * energy dependence of
2180 * SIGMa-TOT, SIGma-INel, SIGma-Diffractive, SIGma-Hard-INelastic
2181 * print- and plott-out
2182 * using:
2183 * called routines SIGSHD , SIGMA1..3 , PLOT
2184 *
2185 * runs thru energies independent of the rest
2186 * and contains no preperation for other parts
2187 *
2188 *---------------------------------------------------------------------
2189  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2190  SAVE
2191 *
2192 * *** /OUTLEV/ controls output level for POMDI and parton X distribution
2193  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
2194 * *** /POLMN/ arrays having to do with cut soft and hard Pomerons
2195  parameter(mxpa25=30,mxpa26=mxpa25+1,mxpa13=13)
2196  parameter( zero=0.d0, one=1.d0)
2197 * PARAMETRIZATION FOR PTMIN= 3. GEV
2198  parameter(mxpa50=250,mxpa51=mxpa50+1)
2199 * PARAMETRIZATION FOR PTMIN= 2. GEV
2200 C PARAMETER (MXPA50=350,MXPA51=MXPA50+1)
2201  COMMON /polmn/plmn(0:mxpa25,0:mxpa50,0:mxpa13),
2202  * plmncu(0:mxpa25,0:mxpa50,0:mxpa13)
2203  COMMON /polmn0/pdifr,phard,psoft,alfah,betah,
2204  * sigtot,sigqel,sigel,sigine,sighin,sigd,sigdd
2205 *
2206 * *** /POMPAR/*/SIGMA/*/POMTYP/ used only in SIGMAPOM-routines (->POMDI)
2207  COMMON /pomtyp/ipim,icon,isig,lmax,mmax,nmax,difel,difnu
2208  common/pompar/alfa,alfap,a,c,ak
2209  COMMON /sigma/sigsof,bs,zsof,sighar,bh,zhar,sigtrp,bt,ztrp,
2210  * sigloo,zloo
2211 * ***
2212  COMMON /topdr/itopd,idumtp
2213 * ***
2214  INTEGER*2 ndislm
2215  COMMON /histoo/as(50,9),aecm(50,9),asig(50,9),alos(50,9),
2216  * aloecm(50,9),ndislm(0:mxpa25,0:mxpa50,0:mxpa13)
2217 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2218  parameter(pi=3.141592654d0)
2219 *
2220 *---------------------------------------------------------------------
2221 * run thru energies
2222 * ------------------------------------------------------------------
2223 *
2224  istep=2
2225  IF(ioutpo.GT.-1)istep=7
2226  DO 100 i=1,50,istep
2227  s=1.6**i
2228  ecm=sqrt(s)+3.4
2229  s=ecm**2
2230 *
2231 * *** running thru energies initializing AS,AECM,ASIG,ALOS,ALOECM *****
2232 *
2233  DO 111 iii=1,9
2234  as(i,iii)=s
2235  aecm(i,iii)=ecm
2236  alos(i,iii)=log10(s)
2237  aloecm(i,iii)=log10(ecm)
2238  asig(i,iii)=0.
2239  111 CONTINUE
2240 *
2241 * *** calling calculation of x- section
2242 *
2243  IF(ipim.EQ.2 )THEN
2244  CALL sigma2(ecm)
2245  IF(i.EQ.1 .AND. ioutpo.GE.0 ) WRITE(6,*)
2246  & ' s-dep. by integr.with Y,PHI,LMD'
2247  ELSE
2248  CALL sigma2(ecm)
2249  IF(i.EQ.1 .AND. ioutpo.GE.0 ) WRITE(6,*)
2250  & ' s-dep. by integr.with Y,PHI,LMD (DEFAULT)'
2251  ENDIF
2252 *
2253 *
2254 C* the preceeding coresponds to a call to SIGMA2 except extra prin
2255 * *** getting plot array ASIG(I,J) :
2256  asig(i,1)=sigtot
2257  asig(i,2)=sigine
2258  asig(i,3)=sighin
2259  asig(i,4)=sigsof
2260  asig(i,5)=sighar
2261  asig(i,6)=sigtrp
2262  asig(i,7)=sigtot-sigine
2263  asig(i,8)=sigine-sighin
2264  asig(i,9)=sigd
2265  WRITE (6,1007)ecm,sigtot,sigine,sigel,sigd
2266  1007 FORMAT (' ECM,SIGTOT,SIGINE,SIGEL,SIGD',f10.1,4e14.3)
2267  100 CONTINUE
2268 *
2269 * energy loop ends
2270 *---------------------------------------------------------------------
2271 * print out results
2272 * ------------------------------------------------------------------
2273  WRITE (6,991)
2274  991 FORMAT (//' shown as line printer plott'/' with'/
2275 * J: drawn quantities:
2276  1 ' (*) SIGTOT total x-section',
2277  2 ' (2) SIGINE inelastic x-section'/
2278  3 ' (3) SIGHIN hard inelastic cross section, one or more jets',
2279  4 ' (4) SIGSOF input soft x-section'/
2280  5 ' (5) SIGHAR input hard x-sections',
2281  6 ' (6) SIGTRP input diffractive x-section (triple pomeron)'/
2282  7 ' (7) SIGTOT-SIGINE elastic x-section',
2283  8 ' (8) SIGINE-SIGHIN non-hard inelastic x-section, (no jets)'/
2284  9 ' (9) SIGD diffractive xross section '/
2285  * ' are plotted against LOG(10)of(CMENERGY)' //)
2286 *
2287  CALL plot(aloecm,asig,450,9,50,zero, 0.1*one,zero, 2.0*one)
2288 *
2289 * special output on unit number 7
2290 C I kept it as it was
2291  IF (itopd.EQ.1) THEN
2292  WRITE(7,95)
2293  95 FORMAT(' NEW FRAME'/' SET FONT DUPLEX'/' SET SCALE X LOG'/
2294  * ' SET LIMITS X FROM 1.0 TO 1E5 Y FROM 0. TO 200'/
2295  * ' TITLE TOP < TOTAL,INEL. AND HARD (MINIJET) CROSS SECT.<'/
2296  * ' TITLE BOTTOM <C.M.ENERGY [GEV]<'/
2297  * ' TITLE < DUAL UNITARIZATION OF SOFT AND HARD CROSS SECTIONS<'/
2298  * ' TITLE LEFT LINES=-1 <CROSS SECTION [MB]<'/
2299  * ' TITLE 3 8.5 < SOLID = TOTAL X.S. <'/
2300  * ' TITLE < DASHED= INELASTIC X.S. <'/
2301  * ' TITLE < DOTTED= HARD X.S.<'/
2302  * ' TITLE < DOT-DASH= HARD INPUT X.S. <'/
2303  * ' TITLE < DOT-DASH= ELASTIC X.S. <')
2304  92 FORMAT (5f15.5)
2305  DO 94 iuu=1,7
2306  IF (iuu.EQ.4)go to 94
2307  IF (iuu.EQ.6)go to 94
2308  IF (iuu.EQ.1) WRITE(7,97)
2309  97 FORMAT (' SET TEXTURE SOLID')
2310  IF (iuu.EQ.2) WRITE(7,98)
2311  98 FORMAT (' SET TEXTURE DASHES')
2312  IF (iuu.EQ.3) WRITE(7,99)
2313  99 FORMAT (' SET TEXTURE DOTS')
2314  IF (iuu.EQ.5) WRITE(7,197)
2315  197 FORMAT (' SET TEXTURE DOTDASH')
2316  DO 93 iu=2,46
2317  WRITE(7,92)aecm(iu,iuu),asig(iu,iuu)
2318  93 CONTINUE
2319  WRITE(7,96)
2320  96 FORMAT (' JOIN')
2321  94 CONTINUE
2322  ENDIF
2323 * ending IF(ITOP=1) special output
2324  RETURN
2325  END
2326 *
2327 ******************************************************************
2328 * end dtupom90
2329 ******************************************************************
2330 
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