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