Double check if SM is running added. Some redundant output removed from SM
[u/mrichter/AliRoot.git] / DPMJET / dpmjet3.0-5.f
CommitLineData
9aaba0d6 1*$ CREATE DT_INIT.FOR
2*COPY DT_INIT
3*
4* +-------------------------------------------------------------+
5* | |
6* | |
7* | DPMJET 3.0 |
8* | |
9* | |
10* | S. Roesler+), R. Engel#), J. Ranft*) |
11* | |
12* | +) CERN, SC-RP |
13* | CH-1211 Geneva 23, Switzerland |
14* | Email: Stefan.Roesler@cern.ch |
15* | |
16* | #) Institut fuer Kernphysik |
17* | Forschungszentrum Karlsruhe |
18* | D-76021 Karlsruhe, Germany |
19* | |
20* | *) University of Siegen, Dept. of Physics |
21* | D-57068 Siegen, Germany |
22* | |
23* | |
24* | http://home.cern.ch/sroesler/dpmjet3.html |
25* | |
26* | |
27* | Monte Carlo models used for event generation: |
28* | PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1 |
29* | |
30* +-------------------------------------------------------------+
31*
32*
33*===init===============================================================*
34*
35 SUBROUTINE DT_INIT(NCASES,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
36 & IDP,IGLAU)
37
38************************************************************************
39* Initialization of event generation *
40* This version dated 7.4.98 is written by S. Roesler. *
41* *
42* Last change 27.12.2006 by S. Roesler. *
43************************************************************************
44
45 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
46 SAVE
47
48 PARAMETER ( LINP = 10 ,
49 & LOUT = 6 ,
50 & LDAT = 9 )
51 PARAMETER (ZERO=0.0D0,ONE=1.0D0)
52
53* particle properties (BAMJET index convention)
54 CHARACTER*8 ANAME
55 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
56 & IICH(210),IIBAR(210),K1(210),K2(210)
57* names of hadrons used in input-cards
58 CHARACTER*8 BTYPE
59 COMMON /DTPAIN/ BTYPE(30)
60* (original name: PAREVT)
61 LOGICAL LDIFFR, LINCTV, LEVPRT, LHEAVY, LDEEXG, LGDHPR, LPREEX,
62 & LHLFIX, LPRFIX, LPARWV, LPOWER, LSNGCH, LLVMOD, LSCHDF
63 PARAMETER ( NALLWP = 39 )
64 COMMON /FKPARE/ DPOWER, FSPRD0, FSHPFN, RN1GSC, RN2GSC,
65 & LDIFFR (NALLWP),LPOWER, LINCTV, LEVPRT, LHEAVY,
66 & LDEEXG, LGDHPR, LPREEX, LHLFIX, LPRFIX, LPARWV,
67 & ILVMOD, JLVMOD, LLVMOD, LSNGCH, LSCHDF
68* (original name: INPFLG)
69 COMMON /FKINPF/ IANG,IFISS,IB0,IGEOM,ISTRAG,KEYDK
70* (original name: FRBKCM)
71 PARAMETER ( MXFFBK = 6 )
72 PARAMETER ( MXZFBK = 9 )
73 PARAMETER ( MXNFBK = 10 )
74 PARAMETER ( MXAFBK = 16 )
75 PARAMETER ( NXZFBK = MXZFBK + MXFFBK / 3 )
76 PARAMETER ( NXNFBK = MXNFBK + MXFFBK / 3 )
77 PARAMETER ( NXAFBK = MXAFBK + 1 )
78 PARAMETER ( MXPSST = 300 )
79 PARAMETER ( MXPSFB = 41000 )
80 LOGICAL LFRMBK, LNCMSS
81 COMMON /FKFRBK/ AMUFBK, EEXFBK (MXPSST), AMFRBK (MXPSST),
82 & EXFRBK (MXPSFB), SDMFBK (MXPSFB), COUFBK (MXPSFB),
83 & EXMXFB, R0FRBK, R0CFBK, C1CFBK, C2CFBK,
84 & IFRBKN (MXPSST), IFRBKZ (MXPSST),
85 & IFBKSP (MXPSST), IFBKPR (MXPSST), IFBKST (MXPSST),
86 & IPSIND (0:MXNFBK,0:MXZFBK,2), JPSIND (0:MXAFBK),
87 & IFBIND (0:NXNFBK,0:NXZFBK,2), JFBIND (0:NXAFBK),
88 & IFBCHA (5,MXPSFB), IPOSST, IPOSFB, IFBSTF,
89 & IFBFRB, NBUFBK, LFRMBK, LNCMSS
90 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
91* emulsion treatment
92 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
93 & NCOMPO,IEMUL
94* Glauber formalism: parameters
95 COMMON /DTGLAM/ RASH(NCOMPX),RBSH(NCOMPX),
96 & BMAX(NCOMPX),BSTEP(NCOMPX),
97 & SIGSH,ROSH,GSH,BSITE(0:NEB,NQB,NCOMPX,KSITEB),
98 & NSITEB,NSTATB
99* Glauber formalism: cross sections
100 COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
101 & XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
102 & XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
103 & XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
104 & XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
105 & XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
106 & XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
107 & XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
108 & XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
109 & BSLOPE,NEBINI,NQBINI
110* interface HADRIN-DPM
111 COMMON /HNTHRE/ EHADTH,EHADLO,EHADHI,INTHAD,IDXTA
112* central particle production, impact parameter biasing
113 COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
114* parameter for intranuclear cascade
115 LOGICAL LPAULI
116 COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
117* various options for treatment of partons (DTUNUC 1.x)
118* (chain recombination, Cronin,..)
119 LOGICAL LCO2CR,LINTPT
120 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
121 & LCO2CR,LINTPT
122* threshold values for x-sampling (DTUNUC 1.x)
123 COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
124 & SSMIMQ,VVMTHR
125* flags for input different options
126 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
127 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
128 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
129* nuclear potential
130 LOGICAL LFERMI
131 COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
132 & EBINDP(2),EBINDN(2),EPOT(2,210),
133 & ETACOU(2),ICOUL,LFERMI
134* n-n cross section fluctuations
135 PARAMETER (NBINS = 1000)
136 COMMON /DTXSFL/ FLUIXX(NBINS),IFLUCT
137* flags for particle decays
138 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
139 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
140 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
141* diquark-breaking mechanism
142 COMMON /DTDIQB/ DBRKR(3,8),DBRKA(3,8),CHAM1,CHAM3,CHAB1,CHAB3
143* nucleon-nucleon event-generator
144 CHARACTER*8 CMODEL
145 LOGICAL LPHOIN
146 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
147* properties of interacting particles
148 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
149* properties of photon/lepton projectiles
150 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
151* flags for diffractive interactions (DTUNUC 1.x)
152 COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
153* parameters for hA-diffraction
154 COMMON /DTDIHA/ DIBETA,DIALPH
155* Lorentz-parameters of the current interaction
156 COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
157 & UMO,PPCM,EPROJ,PPROJ
158* kinematical cuts for lepton-nucleus interactions
159 COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
160 & Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
161* VDM parameter for photon-nucleus interactions
162 COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
163* Glauber formalism: flags and parameters for statistics
164 LOGICAL LPROD
165 CHARACTER*8 CGLB
166 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
167* cuts for variable energy runs
168 COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
169* flags for activated histograms
170 COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
171 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
bd378884 172 COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
9aaba0d6 173* LEPTO
174**LUND single / double precision
175 REAL CUT,PARL,TMPX,TMPY,TMPW2,TMPQ2,TMPU
176 COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),
177 & TMPX,TMPY,TMPW2,TMPQ2,TMPU
178* LEPTO
179 REAL RPPN
180 COMMON /LEPTOI/ RPPN,LEPIN,INTER
181* steering flags for qel neutrino scattering modules
182 COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
183* event flag
184 COMMON /DTEVNO/ NEVENT,ICASCA
185
186 INTEGER PYCOMP
187
188C DIMENSION XPARA(5)
189 DIMENSION XDUMB(40),IPRANG(5)
190
191 PARAMETER (MXCARD=58)
192 CHARACTER*78 CLINE,CTITLE
193 CHARACTER*60 CWHAT
194 CHARACTER*8 BLANK,SDUM
195 CHARACTER*10 CODE,CODEWD
196 CHARACTER*72 HEADER
197 LOGICAL LSTART,LEINP,LXSTAB
198 DIMENSION WHAT(6),CODE(MXCARD)
199 DATA CODE/
200 & 'TITLE ','PROJPAR ','TARPAR ','ENERGY ',
201 & 'MOMENTUM ','CMENERGY ','EMULSION ','FERMI ',
202 & 'TAUFOR ','PAULI ','COULOMB ','HADRIN ',
203 & 'EVAP ','EMCCHECK ','MODEL ','PHOINPUT ',
204 & 'GLAUBERI ','FLUCTUAT ','CENTRAL ','RECOMBIN ',
205 & 'COMBIJET ','XCUTS ','INTPT ','CRONINPT ',
206 & 'SEADISTR ','SEASU3 ','DIQUARKS ','RESONANC ',
207 & 'DIFFRACT ','SINGLECH ','NOFRAGME ','HADRONIZE ',
208 & 'POPCORN ','PARDECAY ','BEAM ','LUND-MSTU ',
209 & 'LUND-MSTJ ','LUND-MDCY ','LUND-PARJ ','LUND-PARU ',
210 & 'OUTLEVEL ','FRAME ','L-TAG ','L-ETAG ',
211 & 'ECMS-CUT ','VDM-PAR1 ','HISTOGRAM ','XS-TABLE ',
212 & 'GLAUB-PAR ','GLAUB-INI ','VDM-PAR2 ','XS-QELPRO ',
213 & 'RNDMINIT ','LEPTO-CUT ','LEPTO-LST ','LEPTO-PARL',
214 & 'START ','STOP '/
215 DATA BLANK /' '/
216
217 DATA LSTART,LXSTAB,IFIRST /.TRUE.,.FALSE.,1/
218 DATA CMEOLD /0.0D0/
219
220*---------------------------------------------------------------------
221* at the first call of INIT: initialize event generation
222 EPNSAV = EPN
223 IF (LSTART) THEN
224 CALL DT_TITLE
225* initialization and test of the random number generator
226 IF (ITRSPT.NE.1) THEN
227 CALL DT_RNDMST(22,54,76,92)
228 CALL DT_RNDMTE(1)
229 ENDIF
230* initialization of BAMJET, DECAY and HADRIN
231 CALL DT_DDATAR
232 CALL DT_DHADDE
233 CALL DT_DCHANT
234 CALL DT_DCHANH
235* set default values for input variables
236 CALL DT_DEFAUL(EPN,PPN)
237 IGLAU = 0
238 IXSQEL = 0
239* flag for collision energy input
240 LEINP = .FALSE.
241 LSTART = .FALSE.
242 ENDIF
243
244*---------------------------------------------------------------------
245 10 CONTINUE
246
247* bypass reading input cards (e.g. for use with Fluka)
248* in this case Epn is expected to carry the beam momentum
249 IF (NCASES.EQ.-1) THEN
250 IP = NPMASS
251 IPZ = NPCHAR
252 PPN = EPNSAV
253 EPN = ZERO
254 CMENER = ZERO
255 LEINP = .TRUE.
256 MKCRON = 0
257 WHAT(1) = 1
258 WHAT(2) = 0
259 CODEWD = 'START '
260 GOTO 900
261 ENDIF
262
263* read control card from input-unit LINP
264 READ(LINP,'(A78)',END=9999) CLINE
265 IF (CLINE(1:1).EQ.'*') THEN
266* comment-line
267 WRITE(LOUT,'(A78)') CLINE
268 GOTO 10
269 ENDIF
270C READ(CLINE,1000,END=9999) CODEWD,(WHAT(I),I=1,6),SDUM
271C1000 FORMAT(A10,6E10.0,A8)
272 DO 1008 I=1,6
273 WHAT(I) = ZERO
274 1008 CONTINUE
275 READ(CLINE,1006,END=9999) CODEWD,CWHAT,SDUM
276 1006 FORMAT(A10,A60,A8)
277 READ(CWHAT,*,END=1007) (WHAT(I),I=1,6)
278 1007 CONTINUE
279 WRITE(LOUT,1001) CODEWD,(WHAT(I),I=1,6),SDUM
280 1001 FORMAT(A10,6G10.3,A8)
281
282 900 CONTINUE
283
284* check for valid control card and get card index
285 ICW = 0
286 DO 11 I=1,MXCARD
287 IF (CODEWD.EQ.CODE(I)) ICW = I
288 11 CONTINUE
289 IF (ICW.EQ.0) THEN
290 WRITE(LOUT,1002) CODEWD
291 1002 FORMAT(/,1X,'---> ',A10,': invalid control-card !',/)
292 GOTO 10
293 ENDIF
294
295 GOTO(
296*------------------------------------------------------------
297* TITLE , PROJPAR , TARPAR , ENERGY , MOMENTUM,
298 & 100 , 110 , 120 , 130 , 140 ,
299*
300*------------------------------------------------------------
301* CMENERGY, EMULSION, FERMI , TAUFOR , PAULI ,
302 & 150 , 160 , 170 , 180 , 190 ,
303*
304*------------------------------------------------------------
305* COULOMB , HADRIN , EVAP , EMCCHECK, MODEL ,
306 & 200 , 210 , 220 , 230 , 240 ,
307*
308*------------------------------------------------------------
309* PHOINPUT, GLAUBERI, FLUCTUAT, CENTRAL , RECOMBIN,
310 & 250 , 260 , 270 , 280 , 290 ,
311*
312*------------------------------------------------------------
313* COMBIJET, XCUTS , INTPT , CRONINPT, SEADISTR,
314 & 300 , 310 , 320 , 330 , 340 ,
315*
316*------------------------------------------------------------
317* SEASU3 , DIQUARKS, RESONANC, DIFFRACT, SINGLECH,
318 & 350 , 360 , 370 , 380 , 390 ,
319*
320*------------------------------------------------------------
321* NOFRAGME, HADRONIZE, POPCORN , PARDECAY, BEAM ,
322 & 400 , 410 , 420 , 430 , 440 ,
323*
324*------------------------------------------------------------
325* LUND-MSTU, LUND-MSTJ, LUND-MDCY, LUND-PARJ, LUND-PARU,
326 & 450 , 451 , 452 , 460 , 470 ,
327*
328*------------------------------------------------------------
329* OUTLEVEL, FRAME , L-TAG , L-ETAG , ECMS-CUT,
330 & 480 , 490 , 500 , 510 , 520 ,
331*
332*------------------------------------------------------------
333* VDM-PAR1, HISTOGRAM, XS-TABLE , GLAUB-PAR, GLAUB-INI,
334 & 530 , 540 , 550 , 560 , 565 ,
335*
336*------------------------------------------------------------
337* , , VDM-PAR2, XS-QELPRO, RNDMINIT ,
338 & 570 , 580 , 590 ,
339*
340*------------------------------------------------------------
341* LEPTO-CUT, LEPTO-LST,LEPTO-PARL, START , STOP )
342 & 600 , 610 , 620 , 630 , 640 ) , ICW
343*
344*------------------------------------------------------------
345
346 GOTO 10
347
348*********************************************************************
349* *
350* control card: codewd = TITLE *
351* *
352* what (1..6), sdum no meaning *
353* *
354* Note: The control-card following this must consist of *
355* a string of characters usually giving the title of *
356* the run. *
357* *
358*********************************************************************
359
360 100 CONTINUE
361 READ(LINP,'(A78)') CTITLE
362 WRITE(LOUT,'(//,5X,A78,//)') CTITLE
363 GOTO 10
364
365*********************************************************************
366* *
367* control card: codewd = PROJPAR *
368* *
369* what (1) = mass number of projectile nucleus default: 1 *
370* what (2) = charge of projectile nucleus default: 1 *
371* what (3..6) no meaning *
372* sdum projectile particle code word *
373* *
374* Note: If sdum is defined what (1..2) have no meaning. *
375* *
376*********************************************************************
377
378 110 CONTINUE
379 IF (SDUM.EQ.BLANK) THEN
380 IP = INT(WHAT(1))
381 IPZ = INT(WHAT(2))
382 IJPROJ = 1
383 IBPROJ = 1
384 ELSE
385 IJPROJ = 0
386 DO 111 II=1,30
387 IF (SDUM.EQ.BTYPE(II)) THEN
388 IP = 1
389 IPZ = 1
390 IF (II.EQ.26) THEN
391 IJPROJ = 135
392 ELSEIF (II.EQ.27) THEN
393 IJPROJ = 136
394 ELSEIF (II.EQ.28) THEN
395 IJPROJ = 133
396 ELSEIF (II.EQ.29) THEN
397 IJPROJ = 134
398 ELSE
399 IJPROJ = II
400 ENDIF
401 IBPROJ = IIBAR(IJPROJ)
402* photon
403 IF ((IJPROJ.EQ.7).AND.(WHAT(1).GT.ZERO)) VIRT = WHAT(1)
404* lepton
405 IF (((IJPROJ.EQ. 3).OR.(IJPROJ.EQ. 4).OR.
406 & (IJPROJ.EQ.10).OR.(IJPROJ.EQ.11)).AND.
407 & (WHAT(1).GT.ZERO)) Q2HI = WHAT(1)
408 ENDIF
409 111 CONTINUE
410 IF (IJPROJ.EQ.0) THEN
411 WRITE(LOUT,1110)
412 1110 FORMAT(/,1X,'invalid PROJPAR card !',/)
413 GOTO 9999
414 ENDIF
415 ENDIF
416 GOTO 10
417
418*********************************************************************
419* *
420* control card: codewd = TARPAR *
421* *
422* what (1) = mass number of target nucleus default: 1 *
423* what (2) = charge of target nucleus default: 1 *
424* what (3..6) no meaning *
425* sdum target particle code word *
426* *
427* Note: If sdum is defined what (1..2) have no meaning. *
428* *
429*********************************************************************
430
431 120 CONTINUE
432 IF (SDUM.EQ.BLANK) THEN
433 IT = INT(WHAT(1))
434 ITZ = INT(WHAT(2))
435 IJTARG = 1
436 IBTARG = 1
437 ELSE
438 IJTARG = 0
439 DO 121 II=1,30
440 IF (SDUM.EQ.BTYPE(II)) THEN
441 IT = 1
442 ITZ = 1
443 IJTARG = II
444 IBTARG = IIBAR(IJTARG)
445 ENDIF
446 121 CONTINUE
447 IF (IJTARG.EQ.0) THEN
448 WRITE(LOUT,1120)
449 1120 FORMAT(/,1X,'invalid TARPAR card !',/)
450 GOTO 9999
451 ENDIF
452 ENDIF
453 GOTO 10
454
455*********************************************************************
456* *
457* control card: codewd = ENERGY *
458* *
459* what (1) = energy (GeV) of projectile in Lab. *
460* if what(1) < 0: |what(1)| = kinetic energy *
461* default: 200 GeV *
462* if |what(2)| > 0: min. energy for variable *
463* energy runs *
464* what (2) = max. energy for variable energy runs *
465* if what(2) < 0: |what(2)| = kinetic energy *
466* *
467*********************************************************************
468
469 130 CONTINUE
470 EPN = WHAT(1)
471 PPN = ZERO
472 CMENER = ZERO
473 IF ((ABS(WHAT(2)).GT.ZERO).AND.
474 & (ABS(WHAT(2)).GT.ABS(WHAT(1)))) THEN
475 VARELO = WHAT(1)
476 VAREHI = WHAT(2)
477 EPN = VAREHI
478 ENDIF
479 LEINP = .TRUE.
480 GOTO 10
481
482*********************************************************************
483* *
484* control card: codewd = MOMENTUM *
485* *
486* what (1) = momentum (GeV/c) of projectile in Lab. *
487* default: 200 GeV/c *
488* what (2..6), sdum no meaning *
489* *
490*********************************************************************
491
492 140 CONTINUE
493 EPN = ZERO
494 PPN = WHAT(1)
495 CMENER = ZERO
496 LEINP = .TRUE.
497 GOTO 10
498
499*********************************************************************
500* *
501* control card: codewd = CMENERGY *
502* *
503* what (1) = energy in nucleon-nucleon cms. *
504* default: none *
505* what (2..6), sdum no meaning *
506* *
507*********************************************************************
508
509 150 CONTINUE
510 EPN = ZERO
511 PPN = ZERO
512 CMENER = WHAT(1)
513 LEINP = .TRUE.
514 GOTO 10
515
516*********************************************************************
517* *
518* control card: codewd = EMULSION *
519* *
520* definition of nuclear emulsions *
521* *
522* what(1) mass number of emulsion component *
523* what(2) charge of emulsion component *
524* what(3) fraction of events in which a scattering on a *
525* nucleus of this properties is performed *
526* what(4,5,6) as what(1,2,3) but for another component *
527* default: no emulsion *
528* sdum no meaning *
529* *
530* Note: If this input-card is once used with valid parameters *
531* TARPAR is obsolete. *
532* Not the absolute values of the fractions are important *
533* but only the ratios of fractions of different comp. *
534* This control card can be repeatedly used to define *
535* emulsions consisting of up to 10 elements. *
536* *
537*********************************************************************
538
539 160 CONTINUE
540 IF ((WHAT(1).GT.ZERO).AND.(WHAT(2).GT.ZERO)
541 & .AND.(ABS(WHAT(3)).GT.ZERO)) THEN
542 NCOMPO = NCOMPO+1
543 IF (NCOMPO.GT.NCOMPX) THEN
544 WRITE(LOUT,1600)
545 STOP
546 ENDIF
547 IEMUMA(NCOMPO) = INT(WHAT(1))
548 IEMUCH(NCOMPO) = INT(WHAT(2))
549 EMUFRA(NCOMPO) = WHAT(3)
550 IEMUL = 1
551C CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
552 ENDIF
553 IF ((WHAT(4).GT.ZERO).AND.(WHAT(5).GT.ZERO)
554 & .AND.(ABS(WHAT(6)).GT.ZERO)) THEN
555 NCOMPO = NCOMPO+1
556 IF (NCOMPO.GT.NCOMPX) THEN
557 WRITE(LOUT,1001)
558 STOP
559 ENDIF
560 IEMUMA(NCOMPO) = INT(WHAT(4))
561 IEMUCH(NCOMPO) = INT(WHAT(5))
562 EMUFRA(NCOMPO) = WHAT(6)
563C CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
564 ENDIF
565 1600 FORMAT(1X,'too many emulsion components - program stopped')
566 GOTO 10
567
568*********************************************************************
569* *
570* control card: codewd = FERMI *
571* *
572* what (1) = -1 Fermi-motion of nucleons not treated *
573* default: 1 *
574* what (2) = scale factor for Fermi-momentum *
575* default: 0.75 *
576* what (3..6), sdum no meaning *
577* *
578*********************************************************************
579
580 170 CONTINUE
581 IF (WHAT(1).EQ.-1.0D0) THEN
582 LFERMI = .FALSE.
583 ELSE
584 LFERMI = .TRUE.
585 ENDIF
586 XMOD = WHAT(2)
587 IF (XMOD.GE.ZERO) FERMOD = XMOD
588 GOTO 10
589
590*********************************************************************
591* *
592* control card: codewd = TAUFOR *
593* *
594* formation time supressed intranuclear cascade *
595* *
596* what (1) formation time (in fm/c) *
597* note: what(1)=10. corresponds roughly to an *
598* average formation time of 1 fm/c *
599* default: 5. fm/c *
600* what (2) number of generations followed *
601* default: 25 *
602* what (3) = 1. p_t-dependent formation zone *
603* = 2. constant formation zone *
604* default: 1 *
605* what (4) modus of selection of nucleus where the *
606* cascade if followed first *
607* = 1. proj./target-nucleus with probab. 1/2 *
608* = 2. nucleus with highest mass *
609* = 3. proj. nucleus if particle is moving in pos. z *
610* targ. nucleus if particle is moving in neg. z *
611* default: 1 *
612* what (5..6), sdum no meaning *
613* *
614*********************************************************************
615
616 180 CONTINUE
617 TAUFOR = WHAT(1)
618 KTAUGE = INT(WHAT(2))
619 INCMOD = 1
620 IF ((WHAT(3).GE.1.0D0).AND.(WHAT(3).LE.2.0D0))
621 & ITAUVE = INT(WHAT(3))
622 IF ((WHAT(4).GE.1.0D0).AND.(WHAT(4).LE.3.0D0))
623 & INCMOD = INT(WHAT(4))
624 GOTO 10
625
626*********************************************************************
627* *
628* control card: codewd = PAULI *
629* *
630* what (1) = -1 Pauli's principle for secondary *
631* interactions not treated *
632* default: 1 *
633* what (2..6), sdum no meaning *
634* *
635*********************************************************************
636
637 190 CONTINUE
638 IF (WHAT(1).EQ.-1.0D0) THEN
639 LPAULI = .FALSE.
640 ELSE
641 LPAULI = .TRUE.
642 ENDIF
643 GOTO 10
644
645*********************************************************************
646* *
647* control card: codewd = COULOMB *
648* *
649* what (1) = -1. Coulomb-energy treatment switched off *
650* default: 1 *
651* what (2..6), sdum no meaning *
652* *
653*********************************************************************
654
655 200 CONTINUE
656 ICOUL = 1
657 IF (WHAT(1).EQ.-1.0D0) THEN
658 ICOUL = 0
659 ELSE
660 ICOUL = 1
661 ENDIF
662 GOTO 10
663
664*********************************************************************
665* *
666* control card: codewd = HADRIN *
667* *
668* HADRIN module *
669* *
670* what (1) = 0. elastic/inelastic interactions with probab. *
671* as defined by cross-sections *
672* = 1. inelastic interactions forced *
673* = 2. elastic interactions forced *
674* default: 1 *
675* what (2) upper threshold in total energy (GeV) below *
676* which interactions are sampled by HADRIN *
677* default: 5. GeV *
678* what (3..6), sdum no meaning *
679* *
680*********************************************************************
681
682 210 CONTINUE
683 IWHAT = INT(WHAT(1))
684 IF ((IWHAT.GE.0).AND.(IWHAT.LE.2)) INTHAD = IWHAT
685 IF ((WHAT(2).GT.ZERO).AND.(WHAT(2).LT.15.0D0)) EHADTH = WHAT(2)
686 GOTO 10
687
688*********************************************************************
689* *
690* control card: codewd = EVAP *
691* *
692* evaporation module *
693* *
694* what (1) =< -1 ==> evaporation is switched off *
695* >= 1 ==> evaporation is performed *
696* *
697* what (1) = i1 + i2*10 + i3*100 + i4*10000 *
698* (i1, i2, i3, i4 >= 0 ) *
699* *
700* i1 is the flag for selecting the T=0 level density option used *
701* = 1: standard EVAP level densities with Cook pairing *
702* energies *
703* = 2: Z,N-dependent Gilbert & Cameron level densities *
704* (default) *
705* = 3: Julich A-dependent level densities *
706* = 4: Z,N-dependent Brancazio & Cameron level densities *
707* *
708* i2 >= 1: high energy fission activated *
709* (default high energy fission activated) *
710* *
711* i3 = 0: No energy dependence for level densities *
712* = 1: Standard Ignyatuk (1975, 1st) energy dependence *
713* for level densities (default) *
714* = 2: Standard Ignyatuk (1975, 1st) energy dependence *
715* for level densities with NOT used set of parameters *
716* = 3: Standard Ignyatuk (1975, 1st) energy dependence *
717* for level densities with NOT used set of parameters *
718* = 4: Second Ignyatuk (1975, 2nd) energy dependence *
719* for level densities *
720* = 5: Second Ignyatuk (1975, 2nd) energy dependence *
721* for level densities with fit 1 Iljinov & Mebel set of *
722* parameters *
723* = 6: Second Ignyatuk (1975, 2nd) energy dependence *
724* for level densities with fit 2 Iljinov & Mebel set of *
725* parameters *
726* = 7: Second Ignyatuk (1975, 2nd) energy dependence *
727* for level densities with fit 3 Iljinov & Mebel set of *
728* parameters *
729* = 8: Second Ignyatuk (1975, 2nd) energy dependence *
730* for level densities with fit 4 Iljinov & Mebel set of *
731* parameters *
732* *
733* i4 >= 1: Original Gilbert and Cameron pairing energies used *
734* (default Cook's modified pairing energies) *
735* *
736* what (2) = ig + 10 * if (ig and if must have the same sign) *
737* *
738* ig =< -1 ==> deexcitation gammas are not produced *
739* (if the evaporation step is not performed *
740* they are never produced) *
741* if =< -1 ==> Fermi Break Up is not invoked *
742* (if the evaporation step is not performed *
743* it is never invoked) *
744* The default is: deexcitation gamma produced and Fermi break up *
745* activated for the new preequilibrium, not *
746* activated otherwise. *
747* what (3..6), sdum no meaning *
748* *
749*********************************************************************
750
751 220 CONTINUE
752 WRITE(LOUT,1009)
753 1009 FORMAT(1X,/,'Warning! Evaporation request rejected since',
754 & ' evaporation modules not available with this version.')
755 LEVPRT = .FALSE.
756 LDEEXG = .FALSE.
757 LHEAVY = .FALSE.
758 LFRMBK = .FALSE.
759 IFISS = 0
760 IEVFSS = 0
761
762 GOTO 10
763
764*********************************************************************
765* *
766* control card: codewd = EMCCHECK *
767* *
768* extended energy-momentum / quantum-number conservation check *
769* *
770* what (1) = -1 extended check not performed *
771* default: 1. *
772* what (2..6), sdum no meaning *
773* *
774*********************************************************************
775
776 230 CONTINUE
777 IF (WHAT(1).EQ.-1) THEN
778 LEMCCK = .FALSE.
779 ELSE
780 LEMCCK = .TRUE.
781 ENDIF
782 GOTO 10
783
784*********************************************************************
785* *
786* control card: codewd = MODEL *
787* *
788* Model to be used to treat nucleon-nucleon interactions *
789* *
790* sdum = DTUNUC two-chain model *
791* = PHOJET multiple chains including minijets *
792* = LEPTO DIS *
793* = QNEUTRIN quasi-elastic neutrino scattering *
794* default: PHOJET *
795* *
796* if sdum = LEPTO: *
797* what (1) (variable INTER) *
798* = 1 gamma exchange *
799* = 2 W+- exchange *
800* = 3 Z0 exchange *
801* = 4 gamma/Z0 exchange *
802* *
803* if sdum = QNEUTRIN: *
804* what (1) = 0 elastic scattering on nucleon and *
805* tau does not decay (default) *
806* = 1 decay of tau into mu.. *
807* = 2 decay of tau into e.. *
808* = 10 CC events on p and n *
809* = 11 NC events on p and n *
810* *
811* what (2..6) no meaning *
812* *
813*********************************************************************
814
815 240 CONTINUE
816 IF (SDUM.EQ.CMODEL(1)) THEN
817 MCGENE = 1
818 ELSEIF (SDUM.EQ.CMODEL(2)) THEN
819 MCGENE = 2
820 ELSEIF (SDUM.EQ.CMODEL(3)) THEN
821 MCGENE = 3
822 IF ((WHAT(1).GE.1.0D0).AND.(WHAT(1).LE.4.0D0))
823 & INTER = INT(WHAT(1))
824 ELSEIF (SDUM.EQ.CMODEL(4)) THEN
825 MCGENE = 4
826 IWHAT = INT(WHAT(1))
827 IF ((IWHAT.EQ.1 ).OR.(IWHAT.EQ.2 ).OR.
828 & (IWHAT.EQ.10).OR.(IWHAT.EQ.11))
829 & NEUDEC = IWHAT
830 ELSE
831 STOP ' Unknown model !'
832 ENDIF
833 GOTO 10
834
835*********************************************************************
836* *
837* control card: codewd = PHOINPUT *
838* *
839* Start of input-section for PHOJET-specific input-cards *
840* Note: This section will not be finished before giving *
841* ENDINPUT-card *
842* what (1..6), sdum no meaning *
843* *
844*********************************************************************
845
846 250 CONTINUE
847 IF (LPHOIN) THEN
848 CALL PHO_INIT(LINP,LOUT,IREJ1)
849 IF (IREJ1.NE.0) THEN
850 WRITE(LOUT,'(1X,A)')'INIT: reading PHOJET-input failed'
851 STOP
852 ENDIF
853 LPHOIN = .FALSE.
854 ENDIF
855 GOTO 10
856
857*********************************************************************
858* *
859* control card: codewd = GLAUBERI *
860* *
861* Pre-initialization of impact parameter selection *
862* *
863* what (1..6), sdum no meaning *
864* *
865*********************************************************************
866
867 260 CONTINUE
868 IF (IFIRST.NE.99) THEN
869 CALL DT_RNDMST(12,34,56,78)
870 CALL DT_RNDMTE(1)
871 OPEN(40,FILE='outdata0/shm.out',STATUS='UNKNOWN')
872C OPEN(11,FILE='outdata0/shm.dbg',STATUS='UNKNOWN')
873 IFIRST = 99
874 ENDIF
875
876 IPPN = 8
877 PLOW = 10.0D0
878C IPPN = 1
879C PLOW = 100.0D0
880 PHI = 1.0D5
881 APLOW = LOG10(PLOW)
882 APHI = LOG10(PHI)
883 ADP = (APHI-APLOW)/DBLE(IPPN)
884
885 IPLOW = 1
886 IDIP = 1
887 IIP = 5
888C IPLOW = 1
889C IDIP = 1
890C IIP = 1
891 IPRANG(1) = 1
892 IPRANG(2) = 2
893 IPRANG(3) = 5
894 IPRANG(4) = 10
895 IPRANG(5) = 20
896
897 ITLOW = 30
898 IDIT = 3
899 IIT = 60
900C IDIT = 10
901C IIT = 21
902
903 DO 473 NCIT=1,IIT
904 IT = ITLOW+(NCIT-1)*IDIT
905C IPHI = IT
906C IDIP = 10
907C IIP = (IPHI-IPLOW)/IDIP
908C IF (IIP.EQ.0) IIP = 1
909C IF (IT.EQ.IPLOW) IIP = 0
910
911 DO 472 NCIP=1,IIP
912 IP = IPRANG(NCIP)
913CC IF (NCIP.LE.IIP) THEN
914C IP = IPLOW+(NCIP-1)*IDIP
915CC ELSE
916CC IP = IT
917CC ENDIF
918 IF (IP.GT.IT) GOTO 472
919
920 DO 471 NCP=1,IPPN+1
921 APPN = APLOW+DBLE(NCP-1)*ADP
922 PPN = 10**APPN
923
924 OPEN(12,FILE='outdata0/shm.sta',STATUS='UNKNOWN')
925 WRITE(12,'(1X,2I5,E15.3)') IP,IT,PPN
926 CLOSE(12)
927
928 XLIM1 = 0.0D0
929 XLIM2 = 50.0D0
930 XLIM3 = ZERO
931 IBIN = 50
932 CALL DT_NEWHGR(XDUM,XDUM,XDUM,XDUMB,-1,IHDUM)
933 CALL DT_NEWHGR(XLIM1,XLIM2,XLIM3,XDUMB,IBIN,IHSHMA)
934
935 NEVFIT = 5
936C IF ((IP.GT.10).OR.(IT.GT.10)) THEN
937C NEVFIT = 5
938C ELSE
939C NEVFIT = 10
940C ENDIF
941 SIGAV = 0.0D0
942
943 DO 478 I=1,NEVFIT
944 CALL DT_SHMAKI(IP,IDUM1,IT,IDUM1,IJPROJ,PPN,99)
945 SIGAV = SIGAV+XSPRO(1,1,1)
946 DO 479 J=1,50
947 XC = DBLE(J)
948 CALL DT_FILHGR(XC,BSITE(1,1,1,J),IHSHMA,I)
949 479 CONTINUE
950 478 CONTINUE
951
952 CALL DT_EVTHIS(IDUM)
953 HEADER = ' BSITE'
954C CALL OUTGEN(IHSHMA,0,0,0,0,0,HEADER,0,NEVFIT,ONE,0,1,-1)
955
956C CALL GENFIT(XPARA)
957C WRITE(40,'(2I4,E11.3,F6.0,5E11.3)')
958C & IP,IT,PPN,SIGAV/DBLE(NEVFIT),XPARA
959
960 471 CONTINUE
961
962 472 CONTINUE
963
964 473 CONTINUE
965
966 STOP
967
968*********************************************************************
969* *
970* control card: codewd = FLUCTUAT *
971* *
972* Treatment of cross section fluctuations *
973* *
974* what (1) = 1 treat cross section fluctuations *
975* default: 0. *
976* what (1..6), sdum no meaning *
977* *
978*********************************************************************
979
980 270 CONTINUE
981 IFLUCT = 0
982 IF (WHAT(1).EQ.ONE) THEN
983 IFLUCT = 1
984 CALL DT_FLUINI
985 ENDIF
986 GOTO 10
987
988*********************************************************************
989* *
990* control card: codewd = CENTRAL *
991* *
992* what (1) = 1. central production forced default: 0 *
993* if what (1) < 0 and > -100 *
994* what (2) = min. impact parameter default: 0 *
995* what (3) = max. impact parameter default: b_max *
996* if what (1) < -99 *
997* what (2) = fraction of cross section default: 1 *
998* if what (1) = -1 : evaporation/fzc suppressed *
999* if what (1) < -1 : evaporation/fzc allowed *
1000* *
1001* what (4..6), sdum no meaning *
1002* *
1003*********************************************************************
1004
1005 280 CONTINUE
1006 ICENTR = INT(WHAT(1))
1007 IF (ICENTR.LT.0) THEN
1008 IF (ICENTR.GT.-100) THEN
1009 BIMIN = WHAT(2)
1010 BIMAX = WHAT(3)
1011 ELSE
1012 XSFRAC = WHAT(2)
1013 ENDIF
1014 ENDIF
1015 GOTO 10
1016
1017*********************************************************************
1018* *
1019* control card: codewd = RECOMBIN *
1020* *
1021* Chain recombination *
1022* (recombine S-S and V-V chains to V-S chains) *
1023* *
1024* what (1) = -1. recombination switched off default: 1 *
1025* what (2..6), sdum no meaning *
1026* *
1027*********************************************************************
1028
1029 290 CONTINUE
1030 IRECOM = 1
1031 IF (WHAT(1).EQ.-1.0D0) IRECOM = 0
1032 GOTO 10
1033
1034*********************************************************************
1035* *
1036* control card: codewd = COMBIJET *
1037* *
1038* chain fusion (2 q-aq --> qq-aqaq) *
1039* *
1040* what (1) = 1 fusion treated *
1041* default: 0. *
1042* what (2) minimum number of uncombined chains from *
1043* single projectile or target nucleons *
1044* default: 0. *
1045* what (3..6), sdum no meaning *
1046* *
1047*********************************************************************
1048
1049 300 CONTINUE
1050 LCO2CR = .FALSE.
1051 IF (INT(WHAT(1)).EQ.1) LCO2CR = .TRUE.
1052 IF (WHAT(2).GE.ZERO) CUTOF = WHAT(2)
1053 GOTO 10
1054
1055*********************************************************************
1056* *
1057* control card: codewd = XCUTS *
1058* *
1059* thresholds for x-sampling *
1060* *
1061* what (1) defines lower threshold for val.-q x-value (CVQ) *
1062* default: 1. *
1063* what (2) defines lower threshold for val.-qq x-value (CDQ) *
1064* default: 2. *
1065* what (3) defines lower threshold for sea-q x-value (CSEA) *
1066* default: 0.2 *
1067* what (4) sea-q x-values in S-S chains (SSMIMA) *
1068* default: 0.14 *
1069* what (5) not used *
1070* default: 2. *
1071* what (6), sdum no meaning *
1072* *
1073* Note: Lower thresholds (what(1..3)) are def. as x_thr=CXXX/ECM *
1074* *
1075*********************************************************************
1076
1077 310 CONTINUE
1078 IF (WHAT(1).GE.0.5D0) CVQ = WHAT(1)
1079 IF (WHAT(2).GE.ONE) CDQ = WHAT(2)
1080 IF (WHAT(3).GE.0.1D0) CSEA = WHAT(3)
1081 IF (WHAT(4).GE.ZERO) THEN
1082 SSMIMA = WHAT(4)
1083 SSMIMQ = SSMIMA**2
1084 ENDIF
1085 IF (WHAT(5).GT.2.0D0) VVMTHR = WHAT(5)
1086 GOTO 10
1087
1088*********************************************************************
1089* *
1090* control card: codewd = INTPT *
1091* *
1092* what (1) = -1 intrinsic transverse momenta of partons *
1093* not treated default: 1 *
1094* what (2..6), sdum no meaning *
1095* *
1096*********************************************************************
1097
1098 320 CONTINUE
1099 IF (WHAT(1).EQ.-1.0D0) THEN
1100 LINTPT = .FALSE.
1101 ELSE
1102 LINTPT = .TRUE.
1103 ENDIF
1104 GOTO 10
1105
1106*********************************************************************
1107* *
1108* control card: codewd = CRONINPT *
1109* *
1110* Cronin effect (multiple scattering of partons at chain ends) *
1111* *
1112* what (1) = -1 Cronin effect not treated default: 1 *
1113* what (2) = 0 scattering parameter default: 0.64 *
1114* what (3..6), sdum no meaning *
1115* *
1116*********************************************************************
1117
1118 330 CONTINUE
1119 IF (WHAT(1).EQ.-1.0D0) THEN
1120 MKCRON = 0
1121 ELSE
1122 MKCRON = 1
1123 ENDIF
1124 CRONCO = WHAT(2)
1125 GOTO 10
1126
1127*********************************************************************
1128* *
1129* control card: codewd = SEADISTR *
1130* *
1131* what (1) (XSEACO) sea(x) prop. 1/x**what (1) default: 1. *
1132* what (2) (UNON) default: 2. *
1133* what (3) (UNOM) default: 1.5 *
1134* what (4) (UNOSEA) default: 5. *
1135* qdis(x) prop. (1-x)**what (1) etc. *
1136* what (5..6), sdum no meaning *
1137* *
1138*********************************************************************
1139
1140 340 CONTINUE
1141 XSEACO = WHAT(1)
1142 XSEACU = 1.05D0-XSEACO
1143 UNON = WHAT(2)
1144 IF (UNON.LT.0.1D0) UNON = 2.0D0
1145 UNOM = WHAT(3)
1146 IF (UNOM.LT.0.1D0) UNOM = 1.5D0
1147 UNOSEA = WHAT(4)
1148 IF (UNOSEA.LT.0.1D0) UNOSEA = 5.0D0
1149 GOTO 10
1150
1151*********************************************************************
1152* *
1153* control card: codewd = SEASU3 *
1154* *
1155* Treatment of strange-quarks at chain ends *
1156* *
1157* what (1) (SEASQ) strange-quark supression factor *
1158* iflav = 1.+rndm*(2.+SEASQ) *
1159* default: 1. *
1160* what (2..6), sdum no meaning *
1161* *
1162*********************************************************************
1163
1164 350 CONTINUE
1165 SEASQ = WHAT(1)
1166 GOTO 10
1167
1168*********************************************************************
1169* *
1170* control card: codewd = DIQUARKS *
1171* *
1172* what (1) = -1. sea-diquark/antidiquark-pairs not treated *
1173* default: 1. *
1174* what (2..6), sdum no meaning *
1175* *
1176*********************************************************************
1177
1178 360 CONTINUE
1179 IF (WHAT(1).EQ.-1.0D0) THEN
1180 LSEADI = .FALSE.
1181 ELSE
1182 LSEADI = .TRUE.
1183 ENDIF
1184 GOTO 10
1185
1186*********************************************************************
1187* *
1188* control card: codewd = RESONANC *
1189* *
1190* treatment of low mass chains *
1191* *
1192* what (1) = -1 low chain masses are not corrected for resonance *
1193* masses (obsolete for BAMJET-fragmentation) *
1194* default: 1. *
1195* what (2) = -1 massless partons default: 1. (massive) *
1196* default: 1. (massive) *
1197* what (3) = -1 chain-system containing chain of too small *
1198* mass is rejected (note: this does not fully *
1199* apply to S-S chains) default: 0. *
1200* what (4..6), sdum no meaning *
1201* *
1202*********************************************************************
1203
1204 370 CONTINUE
1205 IRESCO = 1
1206 IMSHL = 1
1207 IRESRJ = 0
1208 IF (WHAT(1).EQ.-ONE) IRESCO = 0
1209 IF (WHAT(2).EQ.-ONE) IMSHL = 0
1210 IF (WHAT(3).EQ.-ONE) IRESRJ = 1
1211 GOTO 10
1212
1213*********************************************************************
1214* *
1215* control card: codewd = DIFFRACT *
1216* *
1217* Treatment of diffractive events *
1218* *
1219* what (1) = (ISINGD) 0 no single diffraction *
1220* 1 single diffraction included *
1221* +-2 single diffractive events only *
1222* +-3 projectile single diffraction only *
1223* +-4 target single diffraction only *
1224* -5 double pomeron exchange only *
1225* (neg. sign applies to PHOJET events) *
1226* default: 0. *
1227* *
1228* what (2) = (IDOUBD) 0 no double diffraction *
1229* 1 double diffraction included *
1230* 2 double diffractive events only *
1231* default: 0. *
1232* what (3) = 1 projectile diffraction treated (2-channel form.) *
1233* default: 0. *
1234* what (4) = alpha-parameter in projectile diffraction *
1235* default: 0. *
1236* what (5..6), sdum no meaning *
1237* *
1238*********************************************************************
1239
1240 380 CONTINUE
1241 IF (ABS(WHAT(1)).GT.ZERO) ISINGD = INT(WHAT(1))
1242 IF (ABS(WHAT(2)).GT.ZERO) IDOUBD = INT(WHAT(2))
1243 IF ((ISINGD.GT.1).AND.(IDOUBD.GT.1)) THEN
1244 WRITE(LOUT,1380)
1245 1380 FORMAT(1X,'INIT: inconsistent DIFFRACT - input !',/,
1246 & 11X,'IDOUBD is reset to zero')
1247 IDOUBD = 0
1248 ENDIF
1249 IF (WHAT(3).GT.ZERO) DIBETA = WHAT(3)
1250 IF (WHAT(4).GT.ZERO) DIALPH = WHAT(4)
1251 GOTO 10
1252
1253*********************************************************************
1254* *
1255* control card: codewd = SINGLECH *
1256* *
1257* what (1) = 1. Regge contribution (one chain) included *
1258* default: 0. *
1259* what (2..6), sdum no meaning *
1260* *
1261*********************************************************************
1262
1263 390 CONTINUE
1264 ISICHA = 0
1265 IF (WHAT(1).EQ.ONE) ISICHA = 1
1266 GOTO 10
1267
1268*********************************************************************
1269* *
1270* control card: codewd = NOFRAGME *
1271* *
1272* biased chain hadronization *
1273* *
1274* what (1..6) = -1 no of hadronizsation of S-S chains *
1275* = -2 no of hadronizsation of D-S chains *
1276* = -3 no of hadronizsation of S-D chains *
1277* = -4 no of hadronizsation of S-V chains *
1278* = -5 no of hadronizsation of D-V chains *
1279* = -6 no of hadronizsation of V-S chains *
1280* = -7 no of hadronizsation of V-D chains *
1281* = -8 no of hadronizsation of V-V chains *
1282* = -9 no of hadronizsation of comb. chains *
1283* default: complete hadronization *
1284* sdum no meaning *
1285* *
1286*********************************************************************
1287
1288 400 CONTINUE
1289 DO 401 I=1,6
1290 ICHAIN = INT(WHAT(I))
1291 IF ((ICHAIN.LE.-1).AND.(ICHAIN.GE.-9))
1292 & LHADRO(ABS(ICHAIN)) = .FALSE.
1293 401 CONTINUE
1294 GOTO 10
1295
1296*********************************************************************
1297* *
1298* control card: codewd = HADRONIZE *
1299* *
1300* hadronization model and parameter switch *
1301* *
1302* what (1) = 1 hadronization via BAMJET *
1303* = 2 hadronization via JETSET *
1304* default: 2 *
1305* what (2) = 1..3 parameter set to be used *
1306* JETSET: 3 sets available *
1307* ( = 3 default JETSET-parameters) *
1308* BAMJET: 1 set available *
1309* default: 1 *
1310* what (3..6), sdum no meaning *
1311* *
1312*********************************************************************
1313
1314 410 CONTINUE
1315 IWHAT1 = INT(WHAT(1))
1316 IWHAT2 = INT(WHAT(2))
1317 IF ((IWHAT1.EQ.1).OR.(IWHAT1.EQ.2)) IFRAG(1) = IWHAT1
1318 IF ((IWHAT1.EQ.2).AND.(IWHAT2.GE.1).AND.(IWHAT2.LE.3))
1319 & IFRAG(2) = IWHAT2
1320 GOTO 10
1321
1322*********************************************************************
1323* *
1324* control card: codewd = POPCORN *
1325* *
1326* "Popcorn-effect" in fragmentation and diquark breaking diagrams *
1327* *
1328* what (1) = (PDB) frac. of diquark fragmenting directly into *
1329* baryons (PYTHIA/JETSET fragmentation) *
1330* (JETSET: = 0. Popcorn mechanism switched off) *
1331* default: 0.5 *
1332* what (2) = probability for accepting a diquark breaking *
1333* diagram involving the generation of a u/d quark- *
1334* antiquark pair default: 0.0 *
1335* what (3) = same a what (2), here for s quark-antiquark pair *
1336* default: 0.0 *
1337* what (4..6), sdum no meaning *
1338* *
1339*********************************************************************
1340
1341 420 CONTINUE
1342 IF (WHAT(1).GE.0.0D0) PDB = WHAT(1)
1343 IF (WHAT(2).GE.0.0D0) THEN
1344 PDBSEA(1) = WHAT(2)
1345 PDBSEA(2) = WHAT(2)
1346 ENDIF
1347 IF (WHAT(3).GE.0.0D0) PDBSEA(3) = WHAT(3)
1348 DO 421 I=1,8
1349 DBRKA(1,I) = DBRKR(1,I)*PDBSEA(1)/(1.D0-PDBSEA(1))
1350 DBRKA(2,I) = DBRKR(2,I)*PDBSEA(2)/(1.D0-PDBSEA(2))
1351 DBRKA(3,I) = DBRKR(3,I)*PDBSEA(3)/(1.D0-PDBSEA(3))
1352 421 CONTINUE
1353 GOTO 10
1354
1355*********************************************************************
1356* *
1357* control card: codewd = PARDECAY *
1358* *
1359* what (1) = 1. Sigma0/Asigma0 are decaying within JETSET *
1360* = 2. pion^0 decay after intranucl. cascade *
1361* default: no decay *
1362* what (2..6), sdum no meaning *
1363* *
1364*********************************************************************
1365
1366 430 CONTINUE
1367 IF (WHAT(1).EQ.ONE) ISIG0 = 1
1368 IF (WHAT(1).EQ.2.0D0) IPI0 = 1
1369 GOTO 10
1370
1371*********************************************************************
1372* *
1373* control card: codewd = BEAM *
1374* *
1375* definition of beam parameters *
1376* *
1377* what (1/2) > 0 : energy of beam 1/2 (GeV) *
1378* < 0 : abs(what(1/2)) energy per charge of *
1379* beam 1/2 (GeV) *
1380* (beam 1 is directed into positive z-direction) *
1381* what (3) beam crossing angle, defined as 2x angle between *
1382* one beam and the z-axis (micro rad) *
1383* what (4) angle with x-axis defining the collision plane *
1384* what (5..6), sdum no meaning *
1385* *
1386* Note: this card requires previously defined projectile and *
1387* target identities (PROJPAR, TARPAR) *
1388* *
1389*********************************************************************
1390
1391 440 CONTINUE
1392 CALL DT_BEAMPR(WHAT,PPN,1)
1393 EPN = ZERO
1394 CMENER = ZERO
1395 LEINP = .TRUE.
1396 GOTO 10
1397
1398*********************************************************************
1399* *
1400* control card: codewd = LUND-MSTU *
1401* *
1402* set parameter MSTU in JETSET-common /LUDAT1/ *
1403* *
1404* what (1) = index according to LUND-common block *
1405* what (2) = new value of MSTU( int(what(1)) ) *
1406* what (3), what(4) and what (5), what(6) further *
1407* parameter in the same way as what (1) and *
1408* what (2) *
1409* default: default-Lund or corresponding to *
1410* the set given in HADRONIZE *
1411* *
1412*********************************************************************
1413
1414 450 CONTINUE
1415 IF (WHAT(1).GT.ZERO) THEN
1416 NMSTU = NMSTU+1
1417 IMSTU(NMSTU) = INT(WHAT(1))
1418 MSTUX(NMSTU) = INT(WHAT(2))
1419 ENDIF
1420 IF (WHAT(3).GT.ZERO) THEN
1421 NMSTU = NMSTU+1
1422 IMSTU(NMSTU) = INT(WHAT(3))
1423 MSTUX(NMSTU) = INT(WHAT(4))
1424 ENDIF
1425 IF (WHAT(5).GT.ZERO) THEN
1426 NMSTU = NMSTU+1
1427 IMSTU(NMSTU) = INT(WHAT(5))
1428 MSTUX(NMSTU) = INT(WHAT(6))
1429 ENDIF
1430 GOTO 10
1431
1432*********************************************************************
1433* *
1434* control card: codewd = LUND-MSTJ *
1435* *
1436* set parameter MSTJ in JETSET-common /LUDAT1/ *
1437* *
1438* what (1) = index according to LUND-common block *
1439* what (2) = new value of MSTJ( int(what(1)) ) *
1440* what (3), what(4) and what (5), what(6) further *
1441* parameter in the same way as what (1) and *
1442* what (2) *
1443* default: default-Lund or corresponding to *
1444* the set given in HADRONIZE *
1445* *
1446*********************************************************************
1447
1448 451 CONTINUE
1449 IF (WHAT(1).GT.ZERO) THEN
1450 NMSTJ = NMSTJ+1
1451 IMSTJ(NMSTJ) = INT(WHAT(1))
1452 MSTJX(NMSTJ) = INT(WHAT(2))
1453 ENDIF
1454 IF (WHAT(3).GT.ZERO) THEN
1455 NMSTJ = NMSTJ+1
1456 IMSTJ(NMSTJ) = INT(WHAT(3))
1457 MSTJX(NMSTJ) = INT(WHAT(4))
1458 ENDIF
1459 IF (WHAT(5).GT.ZERO) THEN
1460 NMSTJ = NMSTJ+1
1461 IMSTJ(NMSTJ) = INT(WHAT(5))
1462 MSTJX(NMSTJ) = INT(WHAT(6))
1463 ENDIF
1464 GOTO 10
1465
1466*********************************************************************
1467* *
1468* control card: codewd = LUND-MDCY *
1469* *
1470* set parameter MDCY(I,1) for particle decays in JETSET-common *
1471* /LUDAT3/ *
1472* *
1473* what (1-6) = PDG particle index of particle which should *
1474* not decay *
1475* default: default-Lund or forced in *
1476* DT_INITJS *
1477* *
1478*********************************************************************
1479
1480 452 CONTINUE
1481 DO 4521 I=1,6
1482 IF (WHAT(I).NE.ZERO) THEN
1483 KC = PYCOMP(INT(WHAT(I)))
1484 MDCY(KC,1) = 0
1485 ENDIF
1486 4521 CONTINUE
1487 GOTO 10
1488
1489*********************************************************************
1490* *
1491* control card: codewd = LUND-PARJ *
1492* *
1493* set parameter PARJ in JETSET-common /LUDAT1/ *
1494* *
1495* what (1) = index according to LUND-common block *
1496* what (2) = new value of PARJ( int(what(1)) ) *
1497* what (3), what(4) and what (5), what(6) further *
1498* parameter in the same way as what (1) and *
1499* what (2) *
1500* default: default-Lund or corresponding to *
1501* the set given in HADRONIZE *
1502* *
1503*********************************************************************
1504
1505 460 CONTINUE
1506 IF (WHAT(1).NE.ZERO) THEN
1507 NPARJ = NPARJ+1
1508 IPARJ(NPARJ) = INT(WHAT(1))
1509 PARJX(NPARJ) = WHAT(2)
1510 ENDIF
1511 IF (WHAT(3).NE.ZERO) THEN
1512 NPARJ = NPARJ+1
1513 IPARJ(NPARJ) = INT(WHAT(3))
1514 PARJX(NPARJ) = WHAT(4)
1515 ENDIF
1516 IF (WHAT(5).NE.ZERO) THEN
1517 NPARJ = NPARJ+1
1518 IPARJ(NPARJ) = INT(WHAT(5))
1519 PARJX(NPARJ) = WHAT(6)
1520 ENDIF
1521 GOTO 10
1522
1523*********************************************************************
1524* *
1525* control card: codewd = LUND-PARU *
1526* *
1527* set parameter PARJ in JETSET-common /LUDAT1/ *
1528* *
1529* what (1) = index according to LUND-common block *
1530* what (2) = new value of PARU( int(what(1)) ) *
1531* what (3), what(4) and what (5), what(6) further *
1532* parameter in the same way as what (1) and *
1533* what (2) *
1534* default: default-Lund or corresponding to *
1535* the set given in HADRONIZE *
1536* *
1537*********************************************************************
1538
1539 470 CONTINUE
1540 IF (WHAT(1).GT.ZERO) THEN
1541 NPARU = NPARU+1
1542 IPARU(NPARU) = INT(WHAT(1))
1543 PARUX(NPARU) = WHAT(2)
1544 ENDIF
1545 IF (WHAT(3).GT.ZERO) THEN
1546 NPARU = NPARU+1
1547 IPARU(NPARU) = INT(WHAT(3))
1548 PARUX(NPARU) = WHAT(4)
1549 ENDIF
1550 IF (WHAT(5).GT.ZERO) THEN
1551 NPARU = NPARU+1
1552 IPARU(NPARU) = INT(WHAT(5))
1553 PARUX(NPARU) = WHAT(6)
1554 ENDIF
1555 GOTO 10
1556
1557*********************************************************************
1558* *
1559* control card: codewd = OUTLEVEL *
1560* *
1561* output control switches *
1562* *
1563* what (1) = internal rejection informations default: 0 *
1564* what (2) = energy-momentum conservation check output *
1565* default: 0 *
1566* what (3) = internal warning messages default: 0 *
1567* what (4..6), sdum not yet used *
1568* *
1569*********************************************************************
1570
1571 480 CONTINUE
1572 DO 481 K=1,6
1573 IOULEV(K) = INT(WHAT(K))
1574 481 CONTINUE
1575 GOTO 10
1576
1577*********************************************************************
1578* *
1579* control card: codewd = FRAME *
1580* *
1581* frame in which final state is given in DTEVT1 *
1582* *
1583* what (1) = 1 target rest frame (laboratory) *
1584* = 2 nucleon-nucleon cms *
1585* default: 1 *
1586* *
1587*********************************************************************
1588
1589 490 CONTINUE
1590 KFRAME = INT(WHAT(1))
1591 IF ((KFRAME.GE.1).AND.(KFRAME.LE.2)) IFRAME = KFRAME
1592 GOTO 10
1593
1594*********************************************************************
1595* *
1596* control card: codewd = L-TAG *
1597* *
1598* lepton tagger: *
1599* definition of kinematical cuts for radiated photon and *
1600* outgoing lepton detection in lepton-nucleus interactions *
1601* *
1602* what (1) = y_min *
1603* what (2) = y_max *
1604* what (3) = Q^2_min *
1605* what (4) = Q^2_max *
1606* what (5) = theta_min (Lab) *
1607* what (6) = theta_max (Lab) *
1608* default: no cuts *
1609* sdum no meaning *
1610* *
1611*********************************************************************
1612
1613 500 CONTINUE
1614 YMIN = WHAT(1)
1615 YMAX = WHAT(2)
1616 Q2MIN = WHAT(3)
1617 Q2MAX = WHAT(4)
1618 THMIN = WHAT(5)
1619 THMAX = WHAT(6)
1620 GOTO 10
1621
1622*********************************************************************
1623* *
1624* control card: codewd = L-ETAG *
1625* *
1626* lepton tagger: *
1627* what (1) = min. outgoing lepton energy (in Lab) *
1628* what (2) = min. photon energy (in Lab) *
1629* what (3) = max. photon energy (in Lab) *
1630* default: no cuts *
1631* what (2..6), sdum no meaning *
1632* *
1633*********************************************************************
1634
1635 510 CONTINUE
1636 ELMIN = MAX(WHAT(1),ZERO)
1637 EGMIN = MAX(WHAT(2),ZERO)
1638 EGMAX = MAX(WHAT(3),ZERO)
1639 GOTO 10
1640
1641*********************************************************************
1642* *
1643* control card: codewd = ECMS-CUT *
1644* *
1645* what (1) = min. c.m. energy to be sampled *
1646* what (2) = max. c.m. energy to be sampled *
1647* what (3) = min x_Bj to be sampled *
1648* default: no cuts *
1649* what (3..6), sdum no meaning *
1650* *
1651*********************************************************************
1652
1653 520 CONTINUE
1654 ECMIN = WHAT(1)
1655 ECMAX = WHAT(2)
1656 IF (ECMIN.GT.ECMAX) ECMIN = ECMAX
1657 XBJMIN = MAX(WHAT(3),ZERO)
1658 GOTO 10
1659
1660*********************************************************************
1661* *
1662* control card: codewd = VDM-PAR1 *
1663* *
1664* parameters in gamma-nucleus cross section calculation *
1665* *
1666* what (1) = Lambda^2 default: 2. *
1667* what (2) lower limit in M^2 integration *
1668* = 1 (3m_pi)^2 *
1669* = 2 (m_rho0)^2 *
1670* = 3 (m_phi)^2 default: 1 *
1671* what (3) upper limit in M^2 integration *
1672* = 1 s/2 *
1673* = 2 s/4 *
1674* = 3 s default: 3 *
1675* what (4) CKMT F_2 structure function *
1676* = 2212 proton *
1677* = 100 deuteron default: 2212 *
1678* what (5) calculation of gamma-nucleon xsections *
1679* = 1 according to CKMT-parametrization of F_2 *
1680* = 2 integrating SIGVP over M^2 *
1681* = 3 using SIGGA *
1682* = 4 PHOJET cross sections default: 4 *
1683* *
1684* what (6), sdum no meaning *
1685* *
1686*********************************************************************
1687
1688 530 CONTINUE
1689 IF (WHAT(1).GE.ZERO) RL2 = WHAT(1)
1690 IF ((WHAT(2).GE.1).AND.(WHAT(2).LE.3)) INTRGE(1) = INT(WHAT(2))
1691 IF ((WHAT(3).GE.1).AND.(WHAT(3).LE.3)) INTRGE(2) = INT(WHAT(3))
1692 IF ((WHAT(4).EQ.2212).OR.(WHAT(4).EQ.100)) IDPDF = INT(WHAT(4))
1693 IF ((WHAT(5).GE.1).AND.(WHAT(5).LE.4)) MODEGA = INT(WHAT(5))
1694 GOTO 10
1695
1696*********************************************************************
1697* *
1698* control card: codewd = HISTOGRAM *
1699* *
1700* activate different classes of histograms *
1701* *
1702* default: no histograms *
1703* *
1704*********************************************************************
1705
1706 540 CONTINUE
1707 DO 541 J=1,6
1708 IF ((WHAT(J).GE.100).AND.(WHAT(J).LE.150)) THEN
1709 IHISPP(INT(WHAT(J))-100) = 1
1710 ELSEIF ((ABS(WHAT(J)).GE.200).AND.(ABS(WHAT(J)).LE.250)) THEN
1711 IHISXS(INT(ABS(WHAT(J)))-200) = 1
1712 IF (WHAT(J).LT.ZERO) IXSTBL = 1
1713 ENDIF
1714 541 CONTINUE
1715 GOTO 10
1716
1717*********************************************************************
1718* *
1719* control card: codewd = XS-TABLE *
1720* *
1721* output of cross section table for requested interaction *
1722* - particle production deactivated ! - *
1723* *
1724* what (1) lower energy limit for tabulation *
1725* > 0 Lab. frame *
1726* < 0 nucleon-nucleon cms *
1727* what (2) upper energy limit for tabulation *
1728* > 0 Lab. frame *
1729* < 0 nucleon-nucleon cms *
1730* what (3) > 0 # of equidistant lin. bins in E *
1731* < 0 # of equidistant log. bins in E *
1732* what (4) lower limit of particle virtuality (photons) *
1733* what (5) upper limit of particle virtuality (photons) *
1734* what (6) > 0 # of equidistant lin. bins in Q^2 *
1735* < 0 # of equidistant log. bins in Q^2 *
1736* *
1737*********************************************************************
1738
1739 550 CONTINUE
1740 IF (WHAT(1).EQ.99999.0D0) THEN
1741 IRATIO = INT(WHAT(2))
1742 GOTO 10
1743 ENDIF
1744 CMENER = ABS(WHAT(2))
1745 IF (.NOT.LXSTAB) THEN
1746 CALL DT_BERTTP
1747 CALL DT_INCINI
1748 ENDIF
1749 IF ((.NOT.LXSTAB).OR.(CMENER.NE.CMEOLD)) THEN
1750 CMEOLD = CMENER
1751 IF (WHAT(2).GT.ZERO)
1752 & CMENER = SQRT(2.0D0*AAM(1)**2+2.0D0*WHAT(2)*AAM(1))
1753 EPN = ZERO
1754 PPN = ZERO
1755C WRITE(LOUT,*) 'CMENER = ',CMENER
1756 CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,CMENER,1)
1757 CALL DT_PHOINI
1758 ENDIF
1759 CALL DT_XSTABL(WHAT,IXSQEL,IRATIO)
1760 IXSQEL = 0
1761 LXSTAB = .TRUE.
1762 GOTO 10
1763
1764*********************************************************************
1765* *
1766* control card: codewd = GLAUB-PAR *
1767* *
1768* parameters in Glauber-formalism *
1769* *
1770* what (1) # of nucleon configurations sampled in integration *
1771* over nuclear desity default: 1000 *
1772* what (2) # of bins for integration over impact-parameter and *
1773* for profile-function calculation default: 49 *
1774* what (3) = 1 calculation of tot., el. and qel. cross sections *
1775* default: 0 *
1776* what (4) = 1 read pre-calculated impact-parameter distrib. *
1777* from "sdum".glb *
1778* =-1 dump pre-calculated impact-parameter distrib. *
1779* into "sdum".glb *
1780* = 100 read pre-calculated impact-parameter distrib. *
1781* for variable projectile/target/energy runs *
1782* from "sdum".glb *
1783* default: 0 *
1784* what (5..6) no meaning *
1785* sdum if |what (4)| = 1 name of in/output-file (sdum.glb) *
1786* *
1787*********************************************************************
1788
1789 560 CONTINUE
1790 IF (WHAT(1).GT.ZERO) JSTATB = INT(WHAT(1))
1791 IF (WHAT(2).GT.ZERO) JBINSB = INT(WHAT(2))
1792 IF (WHAT(3).EQ.ONE) LPROD = .FALSE.
1793 IF ((ABS(WHAT(4)).EQ.ONE).OR.(WHAT(4).EQ.100)) THEN
1794 IOGLB = INT(WHAT(4))
1795 CGLB = SDUM
1796 ENDIF
1797 GOTO 10
1798
1799*********************************************************************
1800* *
1801* control card: codewd = GLAUB-INI *
1802* *
1803* pre-initialization of profile function *
1804* *
1805* what (1) lower energy limit for initialization *
1806* > 0 Lab. frame *
1807* < 0 nucleon-nucleon cms *
1808* what (2) upper energy limit for initialization *
1809* > 0 Lab. frame *
1810* < 0 nucleon-nucleon cms *
1811* what (3) > 0 # of equidistant lin. bins in E *
1812* < 0 # of equidistant log. bins in E *
1813* what (4) maximum projectile mass number for which the *
1814* Glauber data are initialized for each *
1815* projectile mass number *
1816* (if <= mass given with the PROJPAR-card) *
1817* default: 18 *
1818* what (5) steps in mass number starting from what (4) *
1819* up to mass number defined with PROJPAR-card *
1820* for which Glauber data are initialized *
1821* default: 5 *
1822* what (6) no meaning *
1823* sdum no meaning *
1824* *
1825*********************************************************************
1826
1827 565 CONTINUE
1828 IOGLB = -100
1829 CALL DT_GLBINI(WHAT)
1830 GOTO 10
1831
1832*********************************************************************
1833* *
1834* control card: codewd = VDM-PAR2 *
1835* *
1836* parameters in gamma-nucleus cross section calculation *
1837* *
1838* what (1) = 0 no suppression of shadowing by direct photon *
1839* processes *
1840* = 1 suppression .. default: 1 *
1841* what (2) = 0 no suppression of shadowing by anomalous *
1842* component if photon-F_2 *
1843* = 1 suppression .. default: 1 *
1844* what (3) = 0 no suppression of shadowing by coherence *
1845* length of the photon *
1846* = 1 suppression .. default: 1 *
1847* what (4) = 1 longitudinal polarized photons are taken into *
1848* account *
1849* eps*R*Q^2/M^2 = what(4)*Q^2/M^2 default: 0 *
1850* what (5..6), sdum no meaning *
1851* *
1852*********************************************************************
1853
1854 570 CONTINUE
1855 IF ((WHAT(1).EQ.ZERO).OR.(WHAT(1).EQ.ONE)) ISHAD(1) = INT(WHAT(1))
1856 IF ((WHAT(2).EQ.ZERO).OR.(WHAT(2).EQ.ONE)) ISHAD(2) = INT(WHAT(2))
1857 IF ((WHAT(3).EQ.ZERO).OR.(WHAT(3).EQ.ONE)) ISHAD(3) = INT(WHAT(3))
1858 EPSPOL = WHAT(4)
1859 GOTO 10
1860
1861*********************************************************************
1862* *
1863* control card: XS-QELPRO *
1864* *
1865* what (1..6), sdum no meaning *
1866* *
1867*********************************************************************
1868
1869 580 CONTINUE
1870 IXSQEL = ABS(WHAT(1))
1871 GOTO 10
1872
1873*********************************************************************
1874* *
1875* control card: RNDMINIT *
1876* *
1877* initialization of random number generator *
1878* *
1879* what (1..4) values for initialization (= 1..168) *
1880* what (5..6), sdum no meaning *
1881* *
1882*********************************************************************
1883
1884 590 CONTINUE
1885 IF ((WHAT(1).LT.1.0D0).OR.(WHAT(1).GT.168.0D0)) THEN
1886 NA1 = 22
1887 ELSE
1888 NA1 = WHAT(1)
1889 ENDIF
1890 IF ((WHAT(2).LT.1.0D0).OR.(WHAT(2).GT.168.0D0)) THEN
1891 NA2 = 54
1892 ELSE
1893 NA2 = WHAT(2)
1894 ENDIF
1895 IF ((WHAT(3).LT.1.0D0).OR.(WHAT(3).GT.168.0D0)) THEN
1896 NA3 = 76
1897 ELSE
1898 NA3 = WHAT(3)
1899 ENDIF
1900 IF ((WHAT(4).LT.1.0D0).OR.(WHAT(4).GT.168.0D0)) THEN
1901 NA4 = 92
1902 ELSE
1903 NA4 = WHAT(4)
1904 ENDIF
1905 CALL DT_RNDMST(NA1,NA2,NA3,NA4)
1906 GOTO 10
1907
1908*********************************************************************
1909* *
1910* control card: codewd = LEPTO-CUT *
1911* *
1912* set parameter CUT in LEPTO-common /LEPTOU/ *
1913* *
1914* what (1) = index in CUT-array *
1915* what (2) = new value of CUT( int(what(1)) ) *
1916* what (3), what(4) and what (5), what(6) further *
1917* parameter in the same way as what (1) and *
1918* what (2) *
1919* default: default-LEPTO parameters *
1920* *
1921*********************************************************************
1922
1923 600 CONTINUE
1924 IF (WHAT(1).GT.ZERO) CUT(INT(WHAT(1))) = WHAT(2)
1925 IF (WHAT(3).GT.ZERO) CUT(INT(WHAT(3))) = WHAT(4)
1926 IF (WHAT(5).GT.ZERO) CUT(INT(WHAT(5))) = WHAT(6)
1927 GOTO 10
1928
1929*********************************************************************
1930* *
1931* control card: codewd = LEPTO-LST *
1932* *
1933* set parameter LST in LEPTO-common /LEPTOU/ *
1934* *
1935* what (1) = index in LST-array *
1936* what (2) = new value of LST( int(what(1)) ) *
1937* what (3), what(4) and what (5), what(6) further *
1938* parameter in the same way as what (1) and *
1939* what (2) *
1940* default: default-LEPTO parameters *
1941* *
1942*********************************************************************
1943
1944 610 CONTINUE
1945 IF (WHAT(1).GT.ZERO) LST(INT(WHAT(1))) = INT(WHAT(2))
1946 IF (WHAT(3).GT.ZERO) LST(INT(WHAT(3))) = INT(WHAT(4))
1947 IF (WHAT(5).GT.ZERO) LST(INT(WHAT(5))) = INT(WHAT(6))
1948 GOTO 10
1949
1950*********************************************************************
1951* *
1952* control card: codewd = LEPTO-PARL *
1953* *
1954* set parameter PARL in LEPTO-common /LEPTOU/ *
1955* *
1956* what (1) = index in PARL-array *
1957* what (2) = new value of PARL( int(what(1)) ) *
1958* what (3), what(4) and what (5), what(6) further *
1959* parameter in the same way as what (1) and *
1960* what (2) *
1961* default: default-LEPTO parameters *
1962* *
1963*********************************************************************
1964
1965 620 CONTINUE
1966 IF (WHAT(1).GT.ZERO) PARL(INT(WHAT(1))) = WHAT(2)
1967 IF (WHAT(3).GT.ZERO) PARL(INT(WHAT(3))) = WHAT(4)
1968 IF (WHAT(5).GT.ZERO) PARL(INT(WHAT(5))) = WHAT(6)
1969 GOTO 10
1970
1971*********************************************************************
1972* *
1973* control card: codewd = START *
1974* *
1975* what (1) = number of events default: 100. *
1976* what (2) = 0 Glauber initialization follows *
1977* = 1 Glauber initialization supressed, fitted *
1978* results are used instead *
1979* (this does not apply if emulsion-treatment *
1980* is requested) *
1981* = 2 Glauber initialization is written to *
1982* output-file shmakov.out *
1983* = 3 Glauber initialization is read from input-file *
1984* shmakov.out default: 0 *
1985* what (3..6) no meaning *
1986* what (3..6) no meaning *
1987* *
1988*********************************************************************
1989
1990 630 CONTINUE
1991
1992* check for cross-section table output only
1993 IF (LXSTAB) STOP
1994
1995 NCASES = INT(WHAT(1))
1996 IF (NCASES.LE.0) NCASES = 100
1997 IGLAU = INT(WHAT(2))
1998 IF ((IGLAU.NE.1).AND.(IGLAU.NE.2).AND.(IGLAU.NE.3))
1999 & IGLAU = 0
2000
2001 NPMASS = IP
2002 NPCHAR = IPZ
2003 NTMASS = IT
2004 NTCHAR = ITZ
2005 IDP = IJPROJ
2006 IDT = IJTARG
2007 IF (IDP.LE.0) IDP = 1
2008* muon neutrinos: temporary (missing index)
2009* (new patch in projpar: therefore the following this is probably not
2010* necessary anymore..)
2011C IF (IDP.EQ.26) IDP = 5
2012C IF (IDP.EQ.27) IDP = 6
2013
2014* redefine collision energy
2015 IF (LEINP) THEN
2016 IF (ABS(VAREHI).GT.ZERO) THEN
2017 PDUM = ZERO
2018 IF (VARELO.LT.EHADLO) VARELO = EHADLO
2019 CALL DT_LTINI(IDP,IDT,VARELO,PDUM,VARCLO,1)
2020 PDUM = ZERO
2021 CALL DT_LTINI(IDP,IDT,VAREHI,PDUM,VARCHI,1)
2022 ENDIF
2023 CALL DT_LTINI(IDP,IDT,EPN,PPN,CMENER,1)
2024 ELSE
2025 WRITE(LOUT,1003)
2026 1003 FORMAT(1X,'INIT: collision energy not defined!',/,
2027 & 1X,' -program stopped- ')
2028 STOP
2029 ENDIF
2030
2031* switch off evaporation (even if requested) if central coll. requ.
2032 IF ((ICENTR.EQ.-1).OR.(ICENTR.GT.0).OR.(XSFRAC.LT.0.5D0)) THEN
2033 IF (LEVPRT) THEN
2034 WRITE(LOUT,1004)
2035 1004 FORMAT(1X,/,'Warning! Evaporation request rejected since',
2036 & ' central collisions forced.')
2037 LEVPRT = .FALSE.
2038 LDEEXG = .FALSE.
2039 LHEAVY = .FALSE.
2040 ENDIF
2041 ENDIF
2042
2043* initialization of evaporation-module
2044
2045 WRITE(LOUT,1010)
2046 1010 FORMAT(1X,/,'Warning! No evaporation performed since',
2047 & ' evaporation modules not available with this version.')
2048 LEVPRT = .FALSE.
2049 LDEEXG = .FALSE.
2050 LHEAVY = .FALSE.
2051 LFRMBK = .FALSE.
2052 IFISS = 0
2053 IEVFSS = 0
2054 CALL DT_BERTTP
2055 CALL DT_INCINI
2056
2057* save the default JETSET-parameter
2058 CALL DT_JSPARA(0)
2059
2060* force use of phojet for g-A
2061 IF ((IDP.EQ.7).AND.(MCGENE.NE.3)) MCGENE = 2
2062* initialization of nucleon-nucleon event generator
2063 IF (MCGENE.EQ.2) CALL DT_PHOINI
2064* initialization of LEPTO event generator
2065 IF (MCGENE.EQ.3) THEN
2066
2067 STOP ' This version does not contain LEPTO !'
2068
2069 ENDIF
2070
2071* initialization of quasi-elastic neutrino scattering
2072 IF (MCGENE.EQ.4) THEN
2073 IF (IJPROJ.EQ.5) THEN
2074 NEUTYP = 1
2075 ELSEIF (IJPROJ.EQ.6) THEN
2076 NEUTYP = 2
2077 ELSEIF (IJPROJ.EQ.135) THEN
2078 NEUTYP = 3
2079 ELSEIF (IJPROJ.EQ.136) THEN
2080 NEUTYP = 4
2081 ELSEIF (IJPROJ.EQ.133) THEN
2082 NEUTYP = 5
2083 ELSEIF (IJPROJ.EQ.134) THEN
2084 NEUTYP = 6
2085 ENDIF
2086 ENDIF
2087
2088* normalize fractions of emulsion components
2089 IF (NCOMPO.GT.0) THEN
2090 SUMFRA = ZERO
2091 DO 491 I=1,NCOMPO
2092 SUMFRA = SUMFRA+EMUFRA(I)
2093 491 CONTINUE
2094 IF (SUMFRA.GT.ZERO) THEN
2095 DO 492 I=1,NCOMPO
2096 EMUFRA(I) = EMUFRA(I)/SUMFRA
2097 492 CONTINUE
2098 ENDIF
2099 ENDIF
2100
2101* disallow Cronin's multiple scattering for nucleus-nucleus interactions
6cf1df4c 2102 IF ((IP.GT.1).AND. (IT.GT.1) .AND. (MKCRON.GT.0)) THEN
9aaba0d6 2103 WRITE(LOUT,1005)
2104 1005 FORMAT(/,1X,'INIT: multiple scattering disallowed',/)
2105 MKCRON = 0
2106 ENDIF
2107
2108* initialization of Glauber-formalism (moved to xAEVT, sr 26.3.96)
2109C IF (NCOMPO.LE.0) THEN
2110C CALL DT_SHMAKI(IP,IPZ,IT,ITZ,IDP,PPN,IGLAU)
2111C ELSE
2112C DO 493 I=1,NCOMPO
2113C CALL DT_SHMAKI(IP,IPZ,IEMUMA(I),IEMUCH(I),IDP,PPN,0)
2114C 493 CONTINUE
2115C ENDIF
2116
2117* pre-tabulation of elastic cross-sections
2118 CALL DT_SIGTBL(JDUM,JDUM,DUM,DUM,-1)
2119
2120 CALL DT_XTIME
2121
2122 RETURN
2123
2124*********************************************************************
2125* *
2126* control card: codewd = STOP *
2127* *
2128* stop of the event generation *
2129* *
2130* what (1..6) no meaning *
2131* *
2132*********************************************************************
2133
2134 9999 CONTINUE
2135 WRITE(LOUT,9000)
2136 9000 FORMAT(1X,'---> unexpected end of input !')
2137
2138 640 CONTINUE
2139 STOP
2140
2141 END
2142
2143*$ CREATE DT_KKINC.FOR
2144*COPY DT_KKINC
2145*
2146*===kkinc==============================================================*
2147*
2148 SUBROUTINE DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,
2149 & IREJ)
2150
2151************************************************************************
2152* Treatment of complete nucleus-nucleus or hadron-nucleus scattering *
2153* This subroutine is an update of the previous version written *
2154* by J. Ranft/ H.-J. Moehring. *
2155* This version dated 19.11.95 is written by S. Roesler *
2156************************************************************************
2157
2158 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2159 SAVE
2160 PARAMETER ( LINP = 10 ,
2161 & LOUT = 6 ,
2162 & LDAT = 9 )
2163 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY5=1.0D-5,
2164 & TINY2=1.0D-2,TINY3=1.0D-3)
2165
2166 LOGICAL LFZC
2167
2168* event history
09b429a4 2169
2170 PARAMETER (NMXHEP=4000)
2171 COMMON /POEVT1/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
2172 & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),
2173 & VHEP(4,NMXHEP), NSD1, NSD2, NDD
2174
9aaba0d6 2175 PARAMETER (NMXHKK=200000)
2176 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
2177 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
2178 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
2179* extended event history
2180 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
2181 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
2182 & IHIST(2,NMXHKK)
2183* particle properties (BAMJET index convention)
2184 CHARACTER*8 ANAME
2185 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2186 & IICH(210),IIBAR(210),K1(210),K2(210)
2187* properties of interacting particles
2188 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2189* Lorentz-parameters of the current interaction
2190 COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
2191 & UMO,PPCM,EPROJ,PPROJ
2192* flags for input different options
2193 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2194 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2195 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2196* flags for particle decays
2197 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
2198 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
2199 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
2200* cuts for variable energy runs
2201 COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
2202* Glauber formalism: flags and parameters for statistics
2203 LOGICAL LPROD
2204 CHARACTER*8 CGLB
2205 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
e3f546f5 2206 COMMON /DTGLCP/ RPROJ,RTARG,BIMPAC,
2207 & NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC,
2208 & NCP,NCT
9aaba0d6 2209
2210 DIMENSION WHAT(6)
2211
2212 IREJ = 0
2213 ILOOP = 0
09b429a4 2214 NSD1 = 0
2215 NSD2 = 0
2216 NDD = 0
9aaba0d6 2217 100 CONTINUE
2218 IF (ILOOP.EQ.4) THEN
2219 WRITE(LOUT,1000) NEVHKK
2220 1000 FORMAT(1X,'KKINC: event ',I8,' rejected!')
2221 GOTO 9999
2222 ENDIF
2223 ILOOP = ILOOP+1
2224
2225* variable energy-runs, recalculate parameters for LT's
2226 IF ((ABS(VAREHI).GT.ZERO).OR.(IOGLB.EQ.100)) THEN
2227 PDUM = ZERO
2228 CDUM = ZERO
2229 CALL DT_LTINI(IDP,1,EPN,PDUM,CDUM,1)
2230 ENDIF
2231 IF (EPN.GT.EPROJ) THEN
2232 WRITE(LOUT,'(A,E9.3,2A,E9.3,A)')
2233 & ' Requested energy (',EPN,'GeV) exceeds',
2234 & ' initialization energy (',EPROJ,'GeV) !'
2235 STOP
2236 ENDIF
2237
2238* re-initialize /DTPRTA/
2239 IP = NPMASS
2240 IPZ = NPCHAR
2241 IT = NTMASS
2242 ITZ = NTCHAR
2243 IJPROJ = IDP
2244 IBPROJ = IIBAR(IJPROJ)
2245
2246* calculate nuclear potentials (common /DTNPOT/)
2247 CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
2248
2249* initialize treatment for residual nuclei
2250 CALL DT_RESNCL(EPN,NLOOP,1)
2251
2252* sample hadron/nucleus-nucleus interaction
2253 CALL DT_KKEVNT(KKMAT,IREJ1)
2254 IF (IREJ1.GT.0) THEN
2255 IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 1 in KKINC'
2256 GOTO 9999
2257 ENDIF
2258
2259 IF ((NPMASS.GT.1).OR.(NTMASS.GT.1)) THEN
2260
2261* intranuclear cascade of final state particles for KTAUGE generations
2262* of secondaries
2263 CALL DT_FOZOCA(LFZC,IREJ1)
2264 IF (IREJ1.GT.0) THEN
2265 IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2 in KKINC'
2266 GOTO 9999
2267 ENDIF
2268
2269* baryons unable to escape the nuclear potential are treated as
2270* excited nucleons (ISTHKK=15,16)
2271 CALL DT_SCN4BA
2272
2273* decay of resonances produced in intranuclear cascade processes
2274**sr 15-11-95 should be obsolete
2275C IF (LFZC) CALL DT_DECAY1
2276
2277 101 CONTINUE
2278* treatment of residual nuclei
2279 CALL DT_RESNCL(EPN,NLOOP,2)
2280
2281* evaporation / fission / fragmentation
2282* (if intranuclear cascade was sampled only)
2283 IF (LFZC) THEN
2284 CALL DT_FICONF(IJPROJ,IP,IPZ,IT,ITZ,NLOOP,IREJ1)
2285 IF (IREJ1.GT.1) GOTO 101
2286 IF (IREJ1.EQ.1) GOTO 100
2287 ENDIF
2288
2289 ENDIF
2290
2291* rejection of unphysical configurations
2292 CALL DT_REJUCO(1,IREJ1)
2293 IF (IREJ1.GT.0) THEN
2294 IF (IOULEV(1).GT.0)
2295 & WRITE(LOUT,*) 'rejected 3 in KKINC: too large x'
2296 GOTO 100
2297 ENDIF
2298
2299* transform finale state into Lab.
2300 IFLAG = 2
2301 CALL DT_BEAMPR(WHAT,DUM,IFLAG)
2302 IF ((IFRAME.EQ.1).AND.(IFLAG.EQ.-1)) CALL DT_LT2LAB
2303
2304 IF (IPI0.EQ.1) CALL DT_DECPI0
2305
2306C IF (NEVHKK.EQ.5) CALL DT_EVTOUT(4)
9aaba0d6 2307 RETURN
e3f546f5 2308
9aaba0d6 2309 9999 CONTINUE
2310 IREJ = 1
09b429a4 2311
9aaba0d6 2312 RETURN
2313 END
2314
2315*$ CREATE DT_DEFAUL.FOR
2316*COPY DT_DEFAUL
2317*
2318*===defaul=============================================================*
2319*
2320 SUBROUTINE DT_DEFAUL(EPN,PPN)
2321
2322************************************************************************
2323* Variables are set to default values. *
2324* This version dated 8.5.95 is written by S. Roesler. *
2325************************************************************************
2326
2327 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2328 SAVE
2329 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
2330 PARAMETER (TWOPI = 6.283185307179586454D+00)
2331
2332* particle properties (BAMJET index convention)
2333 CHARACTER*8 ANAME
2334 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2335 & IICH(210),IIBAR(210),K1(210),K2(210)
2336* nuclear potential
2337 LOGICAL LFERMI
2338 COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
2339 & EBINDP(2),EBINDN(2),EPOT(2,210),
2340 & ETACOU(2),ICOUL,LFERMI
2341* interface HADRIN-DPM
2342 COMMON /HNTHRE/ EHADTH,EHADLO,EHADHI,INTHAD,IDXTA
2343* central particle production, impact parameter biasing
2344 COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
2345* properties of interacting particles
2346 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2347* properties of photon/lepton projectiles
2348 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
2349 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2350* emulsion treatment
2351 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2352 & NCOMPO,IEMUL
2353* parameter for intranuclear cascade
2354 LOGICAL LPAULI
2355 COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
2356* various options for treatment of partons (DTUNUC 1.x)
2357* (chain recombination, Cronin,..)
2358 LOGICAL LCO2CR,LINTPT
2359 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
2360 & LCO2CR,LINTPT
2361* threshold values for x-sampling (DTUNUC 1.x)
2362 COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
2363 & SSMIMQ,VVMTHR
2364* flags for input different options
2365 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2366 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2367 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2368* n-n cross section fluctuations
2369 PARAMETER (NBINS = 1000)
2370 COMMON /DTXSFL/ FLUIXX(NBINS),IFLUCT
2371* flags for particle decays
2372 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
2373 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
2374 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
2375* diquark-breaking mechanism
2376 COMMON /DTDIQB/ DBRKR(3,8),DBRKA(3,8),CHAM1,CHAM3,CHAB1,CHAB3
2377* nucleon-nucleon event-generator
2378 CHARACTER*8 CMODEL
2379 LOGICAL LPHOIN
2380 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
2381* flags for diffractive interactions (DTUNUC 1.x)
2382 COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
2383* VDM parameter for photon-nucleus interactions
2384 COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
2385* Glauber formalism: flags and parameters for statistics
2386 LOGICAL LPROD
2387 CHARACTER*8 CGLB
2388 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
2389* kinematical cuts for lepton-nucleus interactions
2390 COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
2391 & Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
2392* flags for activated histograms
2393 COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
2394* cuts for variable energy runs
2395 COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
2396* parameters for hA-diffraction
2397 COMMON /DTDIHA/ DIBETA,DIALPH
2398* LEPTO
2399 REAL RPPN
2400 COMMON /LEPTOI/ RPPN,LEPIN,INTER
2401* steering flags for qel neutrino scattering modules
2402 COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
2403* event flag
2404 COMMON /DTEVNO/ NEVENT,ICASCA
2405
2406 DATA POTMES /0.002D0/
2407
2408* common /DTNPOT/
2409 DO 10 I=1,2
2410 PFERMP(I) = ZERO
2411 PFERMN(I) = ZERO
2412 EBINDP(I) = ZERO
2413 EBINDN(I) = ZERO
2414 DO 11 J=1,210
2415 EPOT(I,J) = ZERO
2416 11 CONTINUE
2417* nucleus independent meson potential
2418 EPOT(I,13) = POTMES
2419 EPOT(I,14) = POTMES
2420 EPOT(I,15) = POTMES
2421 EPOT(I,16) = POTMES
2422 EPOT(I,23) = POTMES
2423 EPOT(I,24) = POTMES
2424 EPOT(I,25) = POTMES
2425 10 CONTINUE
2426 FERMOD = 0.55D0
2427 ETACOU(1) = ZERO
2428 ETACOU(2) = ZERO
2429 ICOUL = 1
2430 LFERMI = .TRUE.
2431
2432* common /HNTHRE/
2433 EHADTH = -99.0D0
2434 EHADLO = 4.06D0
2435 EHADHI = 6.0D0
2436 INTHAD = 1
2437 IDXTA = 2
2438
2439* common /DTIMPA/
2440 ICENTR = 0
2441 BIMIN = ZERO
2442 BIMAX = 1.0D10
2443 XSFRAC = 1.0D0
2444
2445* common /DTPRTA/
2446 IP = 1
2447 IPZ = 1
2448 IT = 1
2449 ITZ = 1
2450 IJPROJ = 1
2451 IBPROJ = 1
2452 IJTARG = 1
2453 IBTARG = 1
2454* common /DTGPRO/
2455 VIRT = ZERO
2456 DO 14 I=1,4
2457 PGAMM(I) = ZERO
2458 PLEPT0(I) = ZERO
2459 PLEPT1(I) = ZERO
2460 PNUCL(I) = ZERO
2461 14 CONTINUE
2462 IDIREC = 0
2463
2464* common /DTFOTI/
2465**sr 7.4.98: changed after corrected B-sampling
2466C TAUFOR = 4.4D0
2467 TAUFOR = 3.5D0
2468 KTAUGE = 25
2469 ITAUVE = 1
2470 INCMOD = 1
2471 LPAULI = .TRUE.
2472
2473* common /DTCHAI/
2474 SEASQ = ONE
2475 MKCRON = 1
2476 CRONCO = 0.64D0
2477 ISICHA = 0
2478 CUTOF = 100.0D0
2479 LCO2CR = .FALSE.
2480 IRECOM = 1
2481 LINTPT = .TRUE.
2482
2483* common /DTXCUT/
2484* definition of soft quark distributions
2485 XSEACU = 0.05D0
2486 UNON = 2.0D0
2487 UNOM = 1.5D0
2488 UNOSEA = 5.0D0
2489* cutoff parameters for x-sampling
2490 CVQ = 1.0D0
2491 CDQ = 2.0D0
2492C CSEA = 0.3D0
2493 CSEA = 0.1D0
2494 SSMIMA = 1.2D0
2495 SSMIMQ = SSMIMA**2
2496 VVMTHR = 2.0D0
2497
2498* common /DTXSFL/
2499 IFLUCT = 0
2500
2501* common /DTFRPA/
2502 PDB = 0.15D0
2503 PDBSEA(1) = 0.0D0
2504 PDBSEA(2) = 0.0D0
2505 PDBSEA(3) = 0.0D0
2506 ISIG0 = 0
2507 IPI0 = 0
2508 NMSTU = 0
2509 NPARU = 0
2510 NMSTJ = 0
2511 NPARJ = 0
2512
2513* common /DTDIQB/
2514 DO 15 I=1,8
2515 DBRKR(1,I) = 5.0D0
2516 DBRKR(2,I) = 5.0D0
2517 DBRKR(3,I) = 10.0D0
2518 DBRKA(1,I) = ZERO
2519 DBRKA(2,I) = ZERO
2520 DBRKA(3,I) = ZERO
2521 15 CONTINUE
2522 CHAM1 = 0.2D0
2523 CHAM3 = 0.5D0
2524 CHAB1 = 0.7D0
2525 CHAB3 = 1.0D0
2526
2527* common /DTFLG3/
2528 ISINGD = 0
2529 IDOUBD = 0
2530 IFLAGD = 0
2531 IDIFF = 0
2532
2533* common /DTMODL/
2534 MCGENE = 2
2535 CMODEL(1) = 'DTUNUC '
2536 CMODEL(2) = 'PHOJET '
2537 CMODEL(3) = 'LEPTO '
2538 CMODEL(4) = 'QNEUTRIN'
2539 LPHOIN = .TRUE.
2540 ELOJET = 5.0D0
2541
2542* common /DTLCUT/
2543 ECMIN = 3.5D0
2544 ECMAX = 1.0D10
2545 XBJMIN = ZERO
2546 ELMIN = ZERO
2547 EGMIN = ZERO
2548 EGMAX = 1.0D10
2549 YMIN = TINY10
2550 YMAX = 0.999D0
2551 Q2MIN = TINY10
2552 Q2MAX = 10.0D0
2553 THMIN = ZERO
2554 THMAX = TWOPI
2555 Q2LI = ZERO
2556 Q2HI = 1.0D10
2557 ECMLI = ZERO
2558 ECMHI = 1.0D10
2559
2560* common /DTVDMP/
2561 RL2 = 2.0D0
2562 INTRGE(1) = 1
2563 INTRGE(2) = 3
2564 IDPDF = 2212
2565 MODEGA = 4
2566 ISHAD(1) = 1
2567 ISHAD(2) = 1
2568 ISHAD(3) = 1
2569 EPSPOL = ZERO
2570
2571* common /DTGLGP/
2572 JSTATB = 1000
2573 JBINSB = 49
2574 CGLB = ' '
2575 IF (ITRSPT.EQ.1) THEN
2576 IOGLB = 100
2577 ELSE
2578 IOGLB = 0
2579 ENDIF
2580 LPROD = .TRUE.
2581
2582* common /DTHIS3/
2583 DO 16 I=1,50
2584 IHISPP(I) = 0
2585 IHISXS(I) = 0
2586 16 CONTINUE
2587 IXSTBL = 0
2588
2589* common /DTVARE/
2590 VARELO = ZERO
2591 VAREHI = ZERO
2592 VARCLO = ZERO
2593 VARCHI = ZERO
2594
2595* common /DTDIHA/
2596 DIBETA = -1.0D0
2597 DIALPH = ZERO
2598
2599* common /LEPTOI/
2600 RPPN = 0.0
2601 LEPIN = 0
2602 INTER = 0
2603
2604* common /QNEUTO/
2605 NEUTYP = 1
2606 NEUDEC = 0
2607
2608* common /DTEVNO/
2609 NEVENT = 1
2610 IF (ITRSPT.EQ.1) THEN
2611 ICASCA = 1
2612 ELSE
2613 ICASCA = 0
2614 ENDIF
2615
2616* default Lab.-energy
2617 EPN = 200.0D0
2618 PPN = SQRT((EPN-AAM(IJPROJ))*(EPN+AAM(IJPROJ)))
2619
2620 RETURN
2621 END
2622
2623*$ CREATE DT_AAEVT.FOR
2624*COPY DT_AAEVT
2625*
2626*===aaevt==============================================================*
2627*
2628 SUBROUTINE DT_AAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
2629 & IDP,IGLAU)
2630
2631************************************************************************
2632* This version dated 22.03.96 is written by S. Roesler. *
2633************************************************************************
2634
2635 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2636 SAVE
2637 PARAMETER ( LINP = 10 ,
2638 & LOUT = 6 ,
2639 & LDAT = 9 )
2640
2641 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2642* emulsion treatment
2643 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2644 & NCOMPO,IEMUL
2645* event flag
2646 COMMON /DTEVNO/ NEVENT,ICASCA
9aaba0d6 2647 CHARACTER*8 DATE,HHMMSS
2648 DIMENSION IDMNYR(3)
09b429a4 2649 NSD1 = 0
2650 NSD2 = 0
2651 NDD = 0
9aaba0d6 2652 KKMAT = 1
2653 NMSG = MAX(NEVTS/100,1)
2654
2655* initialization of run-statistics and histograms
2656 CALL DT_STATIS(1)
2657 CALL PHO_PHIST(1000,DUM)
2658
2659* initialization of Glauber-formalism
2660 IF (NCOMPO.LE.0) THEN
2661 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
2662 ELSE
2663 DO 1 I=1,NCOMPO
2664 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2665 1 CONTINUE
2666 ENDIF
2667 CALL DT_SIGEMU
2668
2669 CALL IDATE(IDMNYR)
2670 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2671 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
2672 CALL ITIME(IDMNYR)
2673 WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
2674 & IDMNYR(1),IDMNYR(2),IDMNYR(3)
2675 WRITE(LOUT,1001) DATE,HHMMSS
2676 1001 FORMAT(/,' DT_AAEVT: Initialisation finished. ( Date: ',A8,
2677 & ' Time: ',A8,' )')
2678
2679* generate NEVTS events
2680 DO 2 IEVT=1,NEVTS
2681
2682* print run-status message
2683 IF (MOD(IEVT,NMSG).EQ.0) THEN
2684 CALL IDATE(IDMNYR)
2685 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2686 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
2687 CALL ITIME(IDMNYR)
2688 WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
2689 & IDMNYR(1),IDMNYR(2),IDMNYR(3)
2690 WRITE(LOUT,1000) IEVT-1,NEVTS,DATE,HHMMSS
2691 1000 FORMAT(/,1X,I8,' out of ',I8,' events sampled ( Date: ',A,
2692 & ' Time: ',A,' )',/)
2693C WRITE(LOUT,1000) IEVT-1
2694C1000 FORMAT(1X,I8,' events sampled')
2695 ENDIF
2696 NEVENT = IEVT
2697* treat nuclear emulsions
2698 IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
2699* composite targets only
2700 KKMAT = -KKMAT
2701* sample this event
2702 CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,IREJ)
2703
2704 CALL PHO_PHIST(2000,DUM)
09b429a4 2705
2706 write(6,*) "Diffractive collisions", NSD1, NSD2, NDD
9aaba0d6 2707
2708 2 CONTINUE
2709
2710* print run-statistics and histograms to output-unit 6
2711 CALL PHO_PHIST(3000,DUM)
2712 CALL DT_STATIS(2)
9aaba0d6 2713 RETURN
2714 END
2715
2716*$ CREATE DT_LAEVT.FOR
2717*COPY DT_LAEVT
2718*
2719*===laevt==============================================================*
2720*
2721 SUBROUTINE DT_LAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
2722 & IDP,IGLAU)
2723
2724************************************************************************
2725* Interface to run DPMJET for lepton-nucleus interactions. *
2726* Kinematics is sampled using the equivalent photon approximation *
2727* Based on GPHERA-routine by R. Engel. *
2728* This version dated 23.03.96 is written by S. Roesler. *
2729************************************************************************
2730
2731 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2732 SAVE
2733 PARAMETER ( LINP = 10 ,
2734 & LOUT = 6 ,
2735 & LDAT = 9 )
2736 PARAMETER (TINY10=1.0D-10,TINY4=1.0D-4,
2737 & ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,THREE=3.0D0)
2738 PARAMETER (TWOPI = 6.283185307179586454D+00,
2739 & PI = TWOPI/TWO,
2740 & ALPHEM = ONE/137.0D0)
2741
2742C CHARACTER*72 HEADER
2743
2744* particle properties (BAMJET index convention)
2745 CHARACTER*8 ANAME
2746 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2747 & IICH(210),IIBAR(210),K1(210),K2(210)
2748* event history
2749 PARAMETER (NMXHKK=200000)
2750 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
2751 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
2752 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
2753* extended event history
2754 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
2755 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
2756 & IHIST(2,NMXHKK)
2757* kinematical cuts for lepton-nucleus interactions
2758 COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
2759 & Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
2760* properties of interacting particles
2761 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2762* properties of photon/lepton projectiles
2763 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
2764* kinematics at lepton-gamma vertex
2765 COMMON /DTLGVX/ PPL0(4),PPL1(4),PPG(4),PPA(4)
2766* flags for activated histograms
2767 COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
2768 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2769* emulsion treatment
2770 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2771 & NCOMPO,IEMUL
2772* Glauber formalism: cross sections
2773 COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
2774 & XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
2775 & XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
2776 & XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
2777 & XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
2778 & XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
2779 & XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
2780 & XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
2781 & XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
2782 & BSLOPE,NEBINI,NQBINI
2783* nucleon-nucleon event-generator
2784 CHARACTER*8 CMODEL
2785 LOGICAL LPHOIN
2786 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
2787* flags for input different options
2788 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2789 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2790 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2791* event flag
2792 COMMON /DTEVNO/ NEVENT,ICASCA
2793
2794 DIMENSION XDUMB(40),BGTA(4)
2795
2796* LEPTO
2797 IF (MCGENE.EQ.3) THEN
2798 STOP ' This version does not contain LEPTO !'
2799 ENDIF
2800
2801 KKMAT = 1
2802 NMSG = MAX(NEVTS/10,1)
2803
2804* mass of incident lepton
2805 AMLPT = AAM(IDP)
2806 AMLPT2 = AMLPT**2
2807 IDPPDG = IDT_IPDGHA(IDP)
2808
2809* consistency of kinematical limits
2810 Q2MIN = MAX(Q2MIN,TINY10)
2811 Q2MAX = MAX(Q2MAX,TINY10)
2812 YMIN = MIN(MAX(YMIN,TINY10),0.999D0)
2813 YMAX = MIN(MAX(YMAX,TINY10),0.999D0)
2814
2815* total energy of the lepton-nucleon system
2816 PTOTLN = SQRT( (PLEPT0(1)+PNUCL(1))**2+(PLEPT0(2)+PNUCL(2))**2
2817 & +(PLEPT0(3)+PNUCL(3))**2 )
2818 ETOTLN = PLEPT0(4)+PNUCL(4)
2819 ECMLN = SQRT((ETOTLN-PTOTLN)*(ETOTLN+PTOTLN))
2820 ECMAX = MIN(ECMAX,ECMLN)
2821 WRITE(LOUT,1003) ECMIN,ECMAX,YMIN,YMAX,Q2MIN,Q2MAX,EGMIN,
2822 & THMIN,THMAX,ELMIN
2823 1003 FORMAT(1X,'LAEVT:',16X,'kinematical cuts',/,22X,
2824 & '------------------',/,9X,'W (min) =',
2825 & F7.1,' GeV (max) =',F7.1,' GeV',/,9X,'y (min) =',
2826 & F7.3,8X,'(max) =',F7.3,/,9X,'Q^2 (min) =',F7.1,
2827 & ' GeV^2 (max) =',F7.1,' GeV^2',/,' (Lab) E_g (min) ='
2828 & ,F7.1,' GeV',/,' (Lab) theta (min) =',F7.4,8X,'(max) =',
2829 & F7.4,' for E_lpt >',F7.1,' GeV',/)
2830
2831* Lorentz-parameter for transf. into Lab
2832 BGTA(1) = PNUCL(1)/AAM(1)
2833 BGTA(2) = PNUCL(2)/AAM(1)
2834 BGTA(3) = PNUCL(3)/AAM(1)
2835 BGTA(4) = PNUCL(4)/AAM(1)
2836* LT of incident lepton into Lab and dump it in DTEVT1
2837 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
2838 & PLEPT0(1),PLEPT0(2),PLEPT0(3),PLEPT0(4),
2839 & PLTOT,PPL0(1),PPL0(2),PPL0(3),PPL0(4))
2840 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
2841 & PNUCL(1),PNUCL(2),PNUCL(3),PNUCL(4),
2842 & PLTOT,PPA(1),PPA(2),PPA(3),PPA(4))
2843* maximum energy of photon nucleon system
2844 PTOTGN = SQRT((YMAX*PPL0(1)+PPA(1))**2+(YMAX*PPL0(2)+PPA(2))**2
2845 & +(YMAX*PPL0(3)+PPA(3))**2)
2846 ETOTGN = YMAX*PPL0(4)+PPA(4)
2847 EGNMAX = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
2848 EGNMAX = MIN(EGNMAX,ECMAX)
2849* minimum energy of photon nucleon system
2850 PTOTGN = SQRT((YMIN*PPL0(1)+PPA(1))**2+(YMIN*PPL0(2)+PPA(2))**2
2851 & +(YMIN*PPL0(3)+PPA(3))**2)
2852 ETOTGN = YMIN*PPL0(4)+PPA(4)
2853 EGNMIN = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
2854 EGNMIN = MAX(EGNMIN,ECMIN)
2855
2856* limits for Glauber-initialization
2857 Q2LI = Q2MIN
2858 Q2HI = MAX(Q2LI,MIN(Q2HI,Q2MAX))
2859 ECMLI = MAX(EGNMIN,THREE)
2860 ECMHI = EGNMAX
2861 WRITE(LOUT,1004) EGNMIN,EGNMAX,ECMLI,ECMHI,Q2LI,Q2HI
2862 1004 FORMAT(1X,'resulting limits:',/,9X,'W (min) =',F7.1,
2863 & ' GeV (max) =',F7.1,' GeV',/,/,' limits for ',
2864 & 'Glauber-initialization:',/,9X,'W (min) =',F7.1,
2865 & ' GeV (max) =',F7.1,' GeV',/,9X,'Q^2 (min) =',F7.1,
2866 & ' GeV^2 (max) =',F7.1,' GeV^2',/)
2867* initialization of Glauber-formalism
2868 IF (NCOMPO.LE.0) THEN
2869 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
2870 ELSE
2871 DO 9 I=1,NCOMPO
2872 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2873 9 CONTINUE
2874 ENDIF
2875 CALL DT_SIGEMU
2876
2877* initialization of run-statistics and histograms
2878 CALL DT_STATIS(1)
2879 CALL PHO_PHIST(1000,DUM)
2880
2881* maximum photon-nucleus cross section
2882 I1 = 1
2883 I2 = 1
2884 RAT = ONE
2885 IF (EGNMAX.GE.ECMNN(NEBINI)) THEN
2886 I1 = NEBINI
2887 I2 = NEBINI
2888 RAT = ONE
2889 ELSEIF (EGNMAX.GT.ECMNN(1)) THEN
2890 DO 5 I=2,NEBINI
2891 IF (EGNMAX.LT.ECMNN(I)) THEN
2892 I1 = I-1
2893 I2 = I
2894 RAT = (EGNMAX-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
2895 GOTO 6
2896 ENDIF
2897 5 CONTINUE
2898 6 CONTINUE
2899 ENDIF
2900 SIGMAX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
2901 EGNXX = EGNMAX
2902 I1 = 1
2903 I2 = 1
2904 RAT = ONE
2905 IF (EGNMIN.GE.ECMNN(NEBINI)) THEN
2906 I1 = NEBINI
2907 I2 = NEBINI
2908 RAT = ONE
2909 ELSEIF (EGNMIN.GT.ECMNN(1)) THEN
2910 DO 7 I=2,NEBINI
2911 IF (EGNMIN.LT.ECMNN(I)) THEN
2912 I1 = I-1
2913 I2 = I
2914 RAT = (EGNMIN-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
2915 GOTO 8
2916 ENDIF
2917 7 CONTINUE
2918 8 CONTINUE
2919 ENDIF
2920 SIGXX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
2921 IF (SIGXX.GT.SIGMAX) EGNXX = EGNMIN
2922 SIGMAX = MAX(SIGMAX,SIGXX)
2923 WRITE(LOUT,'(9X,A,F8.3,A)') 'Sigma_tot (max) =',SIGMAX,' mb'
2924
2925* plot photon flux table
2926 AYMIN = LOG(YMIN)
2927 AYMAX = LOG(YMAX)
2928 AYRGE = AYMAX-AYMIN
2929 MAXTAB = 50
2930 ADY = LOG(YMAX/YMIN)/DBLE(MAXTAB-1)
2931C WRITE(LOUT,'(/,1X,A)') 'LAEVT: photon flux '
2932 DO 1 I=1,MAXTAB
2933 Y = EXP(AYMIN+ADY*DBLE(I-1))
2934 Q2LOW = MAX(Q2MIN,AMLPT2*Y**2/(ONE-Y))
2935 FF1 = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2936 & -TWO*AMLPT2*Y*(ONE/Q2LOW-ONE/Q2MAX))
2937 FF2 = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2938 & -TWO*(ONE-Y)/Y*(ONE-Q2LOW/Q2MAX))
2939C WRITE(LOUT,'(5X,3E15.4)') Y,FF1,FF2
2940 1 CONTINUE
2941
2942* maximum residual weight for flux sampling (dy/y)
2943 YY = YMIN
2944 Q2LOW = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
2945 WGHMAX = (ONE+(ONE-YY)**2)*LOG(Q2MAX/Q2LOW)
2946 & -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
2947
2948 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY0)
2949 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY1)
2950 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY2)
2951 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ0)
2952 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ1)
2953 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ2)
2954 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE0)
2955 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE1)
2956 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE2)
2957 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU0)
2958 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU1)
2959 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU2)
2960 XBLOW = 0.001D0
2961 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX0)
2962 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX1)
2963 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX2)
2964
2965 ITRY = 0
2966 ITRW = 0
2967 NC0 = 0
2968 NC1 = 0
2969
2970* generate events
2971 DO 2 IEVT=1,NEVTS
2972 IF (MOD(IEVT,NMSG).EQ.0) THEN
2973C OPEN(LDAT,FILE='/scrtch3/hr/sroesler/statusd5.out',
2974C & STATUS='UNKNOWN')
2975 WRITE(LOUT,'(1X,I8,A)') IEVT-1,' events sampled'
2976C CLOSE(LDAT)
2977 ENDIF
2978 NEVENT = IEVT
2979
2980 100 CONTINUE
2981 ITRY = ITRY+1
2982
2983* sample y
2984 101 CONTINUE
2985 ITRW = ITRW+1
2986 YY = EXP(AYRGE*DT_RNDM(RAT)+AYMIN)
2987 Q2LOW = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
2988 Q2LOG = LOG(Q2MAX/Q2LOW)
2989 WGH = (ONE+(ONE-YY)**2)*Q2LOG
2990 & -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
2991 IF (WGHMAX.LT.WGH) WRITE(LOUT,1000) YY,WGHMAX,WGH
2992 1000 FORMAT(1X,'LAEVT: weight error!',3E12.5)
2993 IF (DT_RNDM(YY)*WGHMAX.GT.WGH) GOTO 101
2994
2995* sample Q2
2996 YEFF = ONE+(ONE-YY)**2
2997 102 CONTINUE
2998 Q2 = Q2LOW*EXP(Q2LOG*DT_RNDM(YY))
2999 WGH = (YEFF-TWO*(ONE-YY)*Q2LOW/Q2)/YEFF
3000 IF (WGH.LT.DT_RNDM(Q2)) GOTO 102
3001
3002c NC0 = NC0+1
3003c CALL DT_FILHGR(YY,ONE,IHFLY0,NC0)
3004c CALL DT_FILHGR(Q2,ONE,IHFLQ0,NC0)
3005
3006* kinematics at lepton-photon vertex
3007* scattered electron
3008 YQ2 = SQRT((ONE-YY)*Q2)
3009 Q2E = Q2/(4.0D0*PLEPT0(4))
3010 E1Y = (ONE-YY)*PLEPT0(4)
3011 CALL DT_DSFECF(SIF,COF)
3012 PLEPT1(1) = YQ2*COF
3013 PLEPT1(2) = YQ2*SIF
3014 PLEPT1(3) = E1Y-Q2E
3015 PLEPT1(4) = E1Y+Q2E
3016C THETA = ACOS( (E1Y-Q2E)/(E1Y+Q2E) )
3017* radiated photon
3018 PGAMM(1) = -PLEPT1(1)
3019 PGAMM(2) = -PLEPT1(2)
3020 PGAMM(3) = PLEPT0(3)-PLEPT1(3)
3021 PGAMM(4) = PLEPT0(4)-PLEPT1(4)
3022* E_cm cut
3023 PTOTGN = SQRT( (PGAMM(1)+PNUCL(1))**2+(PGAMM(2)+PNUCL(2))**2
3024 & +(PGAMM(3)+PNUCL(3))**2 )
3025 ETOTGN = PGAMM(4)+PNUCL(4)
3026 ECMGN = (ETOTGN-PTOTGN)*(ETOTGN+PTOTGN)
3027 IF (ECMGN.LT.0.1D0) GOTO 101
3028 ECMGN = SQRT(ECMGN)
3029 IF ((ECMGN.LT.ECMIN).OR.(ECMGN.GT.ECMAX)) GOTO 101
3030
3031* Lorentz-transformation into nucleon-rest system
3032 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
3033 & PGAMM(1),PGAMM(2),PGAMM(3),PGAMM(4),
3034 & PGTOT,PPG(1),PPG(2),PPG(3),PPG(4))
3035 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
3036 & PLEPT1(1),PLEPT1(2),PLEPT1(3),PLEPT1(4),
3037 & PLTOT,PPL1(1),PPL1(2),PPL1(3),PPL1(4))
3038* temporary checks..
3039 Q2TMP = ABS(PPG(4)**2-PGTOT**2)
3040 IF (ABS(Q2-Q2TMP).GT.0.01D0) WRITE(LOUT,1001) Q2,Q2TMP
3041 1001 FORMAT(1X,'LAEVT: inconsistent kinematics (Q2,Q2TMP) ',
3042 & 2F10.4)
3043 ECMTMP = SQRT((PPG(4)+AAM(1)-PGTOT)*(PPG(4)+AAM(1)+PGTOT))
3044 IF (ABS(ECMGN-ECMTMP).GT.TINY10) WRITE(LOUT,1002) ECMGN,ECMTMP
3045 1002 FORMAT(1X,'LAEVT: inconsistent kinematics (ECMGN,ECMTMP) ',
3046 & 2F10.2)
3047 YYTMP = PPG(4)/PPL0(4)
3048 IF (ABS(YY-YYTMP).GT.0.01D0) WRITE(LOUT,1005) YY,YYTMP
3049 1005 FORMAT(1X,'LAEVT: inconsistent kinematics (YY,YYTMP) ',
3050 & 2F10.4)
3051
3052* lepton tagger (Lab)
3053 THETA = ACOS( PPL1(3)/PLTOT )
3054 IF (PPL1(4).GT.ELMIN) THEN
3055 IF ((THETA.LT.THMIN).OR.(THETA.GT.THMAX)) GOTO 101
3056 ENDIF
3057* photon energy-cut (Lab)
3058 IF (PPG(4).LT.EGMIN) GOTO 101
3059 IF (PPG(4).GT.EGMAX) GOTO 101
3060* x_Bj cut
3061 XBJ = ABS(Q2/(1.876D0*PPG(4)))
3062 IF (XBJ.LT.XBJMIN) GOTO 101
3063
3064 NC0 = NC0+1
3065 CALL DT_FILHGR( Q2,ONE,IHFLQ0,NC0)
3066 CALL DT_FILHGR( YY,ONE,IHFLY0,NC0)
3067 CALL DT_FILHGR( XBJ,ONE,IHFLX0,NC0)
3068 CALL DT_FILHGR(PPG(4),ONE,IHFLU0,NC0)
3069 CALL DT_FILHGR( ECMGN,ONE,IHFLE0,NC0)
3070
3071* rotation angles against z-axis
3072 COD = PPG(3)/PGTOT
3073C SID = SQRT((ONE-COD)*(ONE+COD))
3074 PPT = SQRT(PPG(1)**2+PPG(2)**2)
3075 SID = PPT/PGTOT
3076 COF = ONE
3077 SIF = ZERO
3078 IF (PGTOT*SID.GT.TINY10) THEN
3079 COF = PPG(1)/(SID*PGTOT)
3080 SIF = PPG(2)/(SID*PGTOT)
3081 ANORF = SQRT(COF*COF+SIF*SIF)
3082 COF = COF/ANORF
3083 SIF = SIF/ANORF
3084 ENDIF
3085
3086 IF (IXSTBL.EQ.0) THEN
3087* change to photon projectile
3088 IJPROJ = 7
3089* set virtuality
3090 VIRT = Q2
3091* re-initialize LTs with new kinematics
3092* !!PGAMM ist set in cms (ECMGN) along z
3093 EPN = ZERO
3094 PPN = ZERO
3095 CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,ECMGN,0)
3096* force Lab-system
3097 IFRAME = 1
3098* get emulsion component if requested
3099 IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
3100* convolute with cross section
3101 CALL DT_SIGGAT(Q2LOW,EGNXX,STOTX,KKMAT)
3102 CALL DT_SIGGAT(Q2,ECMGN,STOT,KKMAT)
3103 IF (STOTX.LT.STOT) WRITE(LOUT,'(1X,A,/,6E12.3)')
3104 & 'LAEVT: warning STOTX<STOT ! ',Q2LOW,EGNMAX,STOTX,
3105 & Q2,ECMGN,STOT
3106 IF (DT_RNDM(Q2)*STOTX.GT.STOT) GOTO 100
3107 NC1 = NC1+1
3108 CALL DT_FILHGR( Q2,ONE,IHFLQ1,NC1)
3109 CALL DT_FILHGR( YY,ONE,IHFLY1,NC1)
3110 CALL DT_FILHGR( XBJ,ONE,IHFLX1,NC1)
3111 CALL DT_FILHGR(PPG(4),ONE,IHFLU1,NC1)
3112 CALL DT_FILHGR( ECMGN,ONE,IHFLE1,NC1)
3113* composite targets only
3114 KKMAT = -KKMAT
3115* sample this event
3116 CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IJPROJ,EPN,KKMAT,
3117 & IREJ)
3118* rotate momenta of final state particles back in photon-nucleon syst.
3119 DO 4 I=NPOINT(4),NHKK
3120 IF ((ABS(ISTHKK(I)).EQ.1).OR.(ISTHKK(I).EQ.1000).OR.
3121 & (ISTHKK(I).EQ.1001)) THEN
3122 PX = PHKK(1,I)
3123 PY = PHKK(2,I)
3124 PZ = PHKK(3,I)
3125 CALL DT_MYTRAN(1,PX,PY,PZ,COD,SID,COF,SIF,
3126 & PHKK(1,I),PHKK(2,I),PHKK(3,I))
3127 ENDIF
3128 4 CONTINUE
3129 ENDIF
3130
3131 CALL DT_FILHGR( Q2,ONE,IHFLQ2,NC1)
3132 CALL DT_FILHGR( YY,ONE,IHFLY2,NC1)
3133 CALL DT_FILHGR( XBJ,ONE,IHFLX2,NC1)
3134 CALL DT_FILHGR(PPG(4),ONE,IHFLU2,NC1)
3135 CALL DT_FILHGR( ECMGN,ONE,IHFLE2,NC1)
3136
3137* dump this event to histograms
3138 CALL PHO_PHIST(2000,DUM)
3139
3140 2 CONTINUE
3141
3142 WGY = ALPHEM/TWOPI*WGHMAX*DBLE(ITRY)/DBLE(ITRW)
3143 WGY = WGY*LOG(YMAX/YMIN)
3144 WEIGHT = WGY*SIGMAX*DBLE(NEVTS)/DBLE(ITRY)
3145
3146C HEADER = ' LAEVT: Q^2 distribution 0'
3147C CALL DT_OUTHGR(IHFLQ0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3148C HEADER = ' LAEVT: Q^2 distribution 1'
3149C CALL DT_OUTHGR(IHFLQ1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3150C HEADER = ' LAEVT: Q^2 distribution 2'
3151C CALL DT_OUTHGR(IHFLQ2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3152C HEADER = ' LAEVT: y distribution 0'
3153C CALL DT_OUTHGR(IHFLY0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3154C HEADER = ' LAEVT: y distribution 1'
3155C CALL DT_OUTHGR(IHFLY1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3156C HEADER = ' LAEVT: y distribution 2'
3157C CALL DT_OUTHGR(IHFLY2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3158C HEADER = ' LAEVT: x distribution 0'
3159C CALL DT_OUTHGR(IHFLX0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3160C HEADER = ' LAEVT: x distribution 1'
3161C CALL DT_OUTHGR(IHFLX1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3162C HEADER = ' LAEVT: x distribution 2'
3163C CALL DT_OUTHGR(IHFLX2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3164C HEADER = ' LAEVT: E_g distribution 0'
3165C CALL DT_OUTHGR(IHFLU0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3166C HEADER = ' LAEVT: E_g distribution 1'
3167C CALL DT_OUTHGR(IHFLU1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3168C HEADER = ' LAEVT: E_g distribution 2'
3169C CALL DT_OUTHGR(IHFLU2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3170C HEADER = ' LAEVT: E_c distribution 0'
3171C CALL DT_OUTHGR(IHFLE0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3172C HEADER = ' LAEVT: E_c distribution 1'
3173C CALL DT_OUTHGR(IHFLE1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3174C HEADER = ' LAEVT: E_c distribution 2'
3175C CALL DT_OUTHGR(IHFLE2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3176
3177* print run-statistics and histograms to output-unit 6
3178 CALL PHO_PHIST(3000,DUM)
3179 IF (IXSTBL.EQ.0) CALL DT_STATIS(2)
3180
3181 RETURN
3182 END
3183
3184*$ CREATE DT_DTUINI.FOR
3185*COPY DT_DTUINI
3186*
3187*===dtuini=============================================================*
3188*
3189 SUBROUTINE DT_DTUINI(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
3190 & IDP,IEMU)
3191
3192 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3193 SAVE
3194
3195 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
3196* emulsion treatment
3197 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
3198 & NCOMPO,IEMUL
3199* Glauber formalism: flags and parameters for statistics
3200 LOGICAL LPROD
3201 CHARACTER*8 CGLB
3202 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
3203
3204 CALL DT_INIT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
3205 CALL DT_STATIS(1)
3206 CALL PHO_PHIST(1000,DUM)
3207 IF (NCOMPO.LE.0) THEN
3208 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
3209 ELSE
3210 DO 1 I=1,NCOMPO
3211 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
3212 1 CONTINUE
3213 ENDIF
3214 IF (IOGLB.NE.100) CALL DT_SIGEMU
3215 IEMU = IEMUL
3216
3217 RETURN
3218 END
3219
3220*$ CREATE DT_DTUOUT.FOR
3221*COPY DT_DTUOUT
3222*
3223*===dtuout=============================================================*
3224*
3225 SUBROUTINE DT_DTUOUT
3226
3227 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3228 SAVE
3229
3230 CALL PHO_PHIST(3000,DUM)
3231 CALL DT_STATIS(2)
3232
3233 RETURN
3234 END
3235
3236*$ CREATE DT_BEAMPR.FOR
3237*COPY DT_BEAMPR
3238*
3239*===beampr=============================================================*
3240*
3241 SUBROUTINE DT_BEAMPR(WHAT,PLAB,MODE)
3242
3243************************************************************************
3244* Initialization of event generation *
3245* This version dated 7.4.98 is written by S. Roesler. *
3246************************************************************************
3247
3248 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3249 SAVE
3250
3251 PARAMETER ( LINP = 10 ,
3252 & LOUT = 6 ,
3253 & LDAT = 9 )
3254 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
3255 PARAMETER (TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
3256
3257 LOGICAL LBEAM
3258
3259* event history
3260 PARAMETER (NMXHKK=200000)
3261 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3262 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3263 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3264* extended event history
3265 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3266 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3267 & IHIST(2,NMXHKK)
3268* properties of interacting particles
3269 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
3270* particle properties (BAMJET index convention)
3271 CHARACTER*8 ANAME
3272 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
3273 & IICH(210),IIBAR(210),K1(210),K2(210)
3274* beam momenta
3275 COMMON /DTBEAM/ P1(4),P2(4)
3276
3277C DIMENSION WHAT(6),P1(4),P2(4),P1CMS(4),P2CMS(4)
3278 DIMENSION WHAT(6),P1CMS(4),P2CMS(4)
3279
3280 DATA LBEAM /.FALSE./
3281
3282 GOTO (1,2) MODE
3283
3284 1 CONTINUE
3285
3286 E1 = WHAT(1)
3287 IF (E1.LT.ZERO) E1 = DBLE(IPZ)/DBLE(IP)*ABS(WHAT(1))
3288 E2 = WHAT(2)
3289 IF (E2.LT.ZERO) E2 = DBLE(ITZ)/DBLE(IT)*ABS(WHAT(2))
3290 PP1 = SQRT( (E1+AAM(IJPROJ))*(E1-AAM(IJPROJ)) )
3291 PP2 = SQRT( (E2+AAM(IJTARG))*(E2-AAM(IJTARG)) )
3292 TH = 1.D-6*WHAT(3)/2.D0
3293 PH = WHAT(4)*BOG
3294 P1(1) = PP1*SIN(TH)*COS(PH)
3295 P1(2) = PP1*SIN(TH)*SIN(PH)
3296 P1(3) = PP1*COS(TH)
3297 P1(4) = E1
3298 P2(1) = PP2*SIN(TH)*COS(PH)
3299 P2(2) = PP2*SIN(TH)*SIN(PH)
3300 P2(3) = -PP2*COS(TH)
3301 P2(4) = E2
3302 ECM = SQRT( (P1(4)+P2(4))**2-(P1(1)+P2(1))**2-(P1(2)+P2(2))**2
3303 & -(P1(3)+P2(3))**2 )
3304 ELAB = (ECM**2-AAM(IJPROJ)**2-AAM(IJTARG)**2)/(2.0D0*AAM(IJTARG))
3305 PLAB = SQRT( (ELAB+AAM(IJPROJ))*(ELAB-AAM(IJPROJ)) )
3306 BGX = (P1(1)+P2(1))/ECM
3307 BGY = (P1(2)+P2(2))/ECM
3308 BGZ = (P1(3)+P2(3))/ECM
3309 BGE = (P1(4)+P2(4))/ECM
3310 CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P1(1),P1(2),P1(3),P1(4),
3311 & P1TOT,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4))
3312 CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P2(1),P2(2),P2(3),P2(4),
3313 & P2TOT,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4))
3314 COD = P1CMS(3)/P1TOT
3315C SID = SQRT((ONE-COD)*(ONE+COD))
3316 PPT = SQRT(P1CMS(1)**2+P1CMS(2)**2)
3317 SID = PPT/P1TOT
3318 COF = ONE
3319 SIF = ZERO
3320 IF (P1TOT*SID.GT.TINY10) THEN
3321 COF = P1CMS(1)/(SID*P1TOT)
3322 SIF = P1CMS(2)/(SID*P1TOT)
3323 ANORF = SQRT(COF*COF+SIF*SIF)
3324 COF = COF/ANORF
3325 SIF = SIF/ANORF
3326 ENDIF
3327**check
3328C WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
3329C WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
3330C WRITE(LOUT,'(5E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),P1TOT
3331C WRITE(LOUT,'(5E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),P2TOT
3332C PAX = ZERO
3333C PAY = ZERO
3334C PAZ = P1TOT
3335C PAE = SQRT(AAM(IJPROJ)**2+PAZ**2)
3336C PBX = ZERO
3337C PBY = ZERO
3338C PBZ = -P2TOT
3339C PBE = SQRT(AAM(IJTARG)**2+PBZ**2)
3340C WRITE(LOUT,'(4E15.4)') PAX,PAY,PAZ,PAE
3341C WRITE(LOUT,'(4E15.4)') PBX,PBY,PBZ,PBE
3342C CALL DT_MYTRAN(1,PAX,PAY,PAZ,COD,SID,COF,SIF,
3343C & P1CMS(1),P1CMS(2),P1CMS(3))
3344C CALL DT_MYTRAN(1,PBX,PBY,PBZ,COD,SID,COF,SIF,
3345C & P2CMS(1),P2CMS(2),P2CMS(3))
3346C WRITE(LOUT,'(4E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4)
3347C WRITE(LOUT,'(4E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4)
3348C CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),
3349C & P1TOT,P1(1),P1(2),P1(3),P1(4))
3350C CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),
3351C & P2TOT,P2(1),P2(2),P2(3),P2(4))
3352C WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
3353C WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
3354C STOP
3355**
3356
3357 LBEAM = .TRUE.
3358
3359 RETURN
3360
3361 2 CONTINUE
3362
3363 IF (LBEAM) THEN
3364 IF ( (NPOINT(4).EQ.0).OR.(NHKK.LT.NPOINT(4)) ) RETURN
3365 DO 20 I=NPOINT(4),NHKK
430525dd 3366 IF ((ABS(ISTHKK(I)).EQ.1) .OR.
3367 & (ABS(ISTHKK(I)).EQ.2) .OR.
3368 & (ISTHKK(I).EQ.1000) .OR.
3369 & (ISTHKK(I).EQ.1001)) THEN
3370
9aaba0d6 3371 CALL DT_MYTRAN(1,PHKK(1,I),PHKK(2,I),PHKK(3,I),
3372 & COD,SID,COF,SIF,PXCMS,PYCMS,PZCMS)
3373 PECMS = PHKK(4,I)
3374 CALL DT_DALTRA(BGE,BGX,BGY,BGZ,PXCMS,PYCMS,PZCMS,PECMS,
3375 & PTOT,PHKK(1,I),PHKK(2,I),PHKK(3,I),PHKK(4,I))
3376 ENDIF
3377 20 CONTINUE
3378 ELSE
3379 MODE = -1
3380 ENDIF
3381
3382 RETURN
3383 END
3384
3385*$ CREATE DT_REJUCO.FOR
3386*COPY DT_REJUCO
3387*
3388*===rejuco=============================================================*
3389*
3390 SUBROUTINE DT_REJUCO(MODE,IREJ)
3391
3392************************************************************************
3393* REJection of Unphysical COnfigurations *
3394* MODE = 1 rejection of particles with unphysically large energy *
3395* *
3396* This version dated 27.12.2006 is written by S. Roesler. *
3397************************************************************************
3398
3399 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3400 SAVE
3401
3402 PARAMETER ( LINP = 10 ,
3403 & LOUT = 6 ,
3404 & LDAT = 9 )
3405 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
3406 PARAMETER (TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
3407
3408* maximum x_cms of final state particle
3409 PARAMETER (XCMSMX = 1.4D0)
3410
3411* event history
3412 PARAMETER (NMXHKK=200000)
3413 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3414 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3415 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3416* extended event history
3417 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3418 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3419 & IHIST(2,NMXHKK)
3420* Lorentz-parameters of the current interaction
3421 COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
3422 & UMO,PPCM,EPROJ,PPROJ
3423
3424 IREJ = 0
3425
3426 IF (MODE.EQ.1) THEN
3427 IF ( (NPOINT(4).EQ.0).OR.(NHKK.LT.NPOINT(4)) ) RETURN
3428 ECMHLF = UMO/2.0D0
3429 DO 10 I=NPOINT(4),NHKK
3430 IF ((ABS(ISTHKK(I)).EQ.1).AND.(IDHKK(I).NE.80000)) THEN
3431 XCMS = ABS(PHKK(4,I))/ECMHLF
3432 IF (XCMS.GT.XCMSMX) GOTO 9999
3433 ENDIF
3434 10 CONTINUE
3435 ENDIF
3436
3437 RETURN
3438 9999 CONTINUE
3439 IREJ = 1
3440 RETURN
3441 END
3442
3443*$ CREATE DT_EVENTB.FOR
3444*COPY DT_EVENTB
3445*
3446*===eventb=============================================================*
3447*
3448 SUBROUTINE DT_EVENTB(NCSY,IREJ)
3449
3450************************************************************************
3451* Treatment of nucleon-nucleon interactions with full two-component *
3452* Dual Parton Model. *
3453* NCSY number of nucleon-nucleon interactions *
3454* IREJ rejection flag *
3455* This version dated 14.01.2000 is written by S. Roesler *
3456************************************************************************
3457
3458 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3459 SAVE
3460 PARAMETER ( LINP = 10 ,
3461 & LOUT = 6 ,
3462 & LDAT = 9 )
3463 PARAMETER (TINY10=1.0D-10,ZERO=0.0D0,ONE=1.0D0)
3464
3465* event history
3466 PARAMETER (NMXHKK=200000)
3467 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3468 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3469 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3470* extended event history
3471 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3472 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3473 & IHIST(2,NMXHKK)
3474*! uncomment this line for internal phojet-fragmentation
3475C #include "dtu_dtevtp.inc"
3476* particle properties (BAMJET index convention)
3477 CHARACTER*8 ANAME
3478 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
3479 & IICH(210),IIBAR(210),K1(210),K2(210)
3480* flags for input different options
3481 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
3482 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
3483 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
3484* rejection counter
3485 COMMON /DTREJC/ IRPT,IRHHA,IRRES(2),LOMRES,LOBRES,
3486 & IRCHKI(2),IRFRAG,IRCRON(3),IREVT,
3487 & IREXCI(3),IRDIFF(2),IRINC
3488* properties of interacting particles
3489 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
3490* properties of photon/lepton projectiles
3491 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
3492* various options for treatment of partons (DTUNUC 1.x)
3493* (chain recombination, Cronin,..)
3494 LOGICAL LCO2CR,LINTPT
3495 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
3496 & LCO2CR,LINTPT
3497* statistics
3498 COMMON /DTSTA1/ ICREQU,ICSAMP,ICCPRO,ICDPR,ICDTA,
3499 & ICRJSS,ICVV2S,ICCHAI(2,9),ICRES(9),ICDIFF(5),
3500 & ICEVTG(8,0:30)
3501* DTUNUC-PHOJET interface, Lorentz-param. of n-n subsystem
3502 COMMON /DTLTSU/ BGX,BGY,BGZ,GAM
3503* Glauber formalism: collision properties
3504 COMMON /DTGLCP/ RPROJ,RTARG,BIMPAC,
e3f546f5 3505 & NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC,
3506 & NCP,NCT
9aaba0d6 3507* flags for diffractive interactions (DTUNUC 1.x)
3508 COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
3509* statistics: double-Pomeron exchange
3510 COMMON /DTFLG2/ INTFLG,IPOPO
3511* flags for particle decays
3512 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
3513 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
3514 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
3515* nucleon-nucleon event-generator
3516 CHARACTER*8 CMODEL
3517 LOGICAL LPHOIN
3518 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
3519C nucleon-nucleus / nucleus-nucleus interface to DPMJET
3520 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
3521 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
3522 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
3523 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
3524C model switches and parameters
3525 CHARACTER*8 MDLNA
3526 INTEGER ISWMDL,IPAMDL
3527 DOUBLE PRECISION PARMDL
3528 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
3529C initial state parton radiation (internal part)
3530 INTEGER MXISR3,MXISR4
3531 PARAMETER ( MXISR3 = 50, MXISR4 = 100 )
3532 INTEGER IFL1,IFL2,IBRA,IFANO,ISH,NACC
3533 DOUBLE PRECISION Q2SH,PT2SH,XPSH,ZPSH,THSH,SHAT
3534 COMMON /POINT6/ Q2SH(2,MXISR3),PT2SH(2,MXISR3),XPSH(2,MXISR3),
3535 & ZPSH(2,MXISR3),THSH(2,MXISR3),SHAT(MXISR3),
3536 & IFL1(2,MXISR3),IFL2(2,MXISR3),
3537 & IBRA(2,MXISR4),IFANO(2),ISH(2),NACC
3538C event debugging information
3539 INTEGER NMAXD
3540 PARAMETER (NMAXD=100)
3541 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
3542 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3543 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
3544 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3545C general process information
3546 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
3547 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
3548
3549 DIMENSION PP(4),PT(4),PTOT(4),PP1(4),PP2(4),PT1(4),PT2(4),
3550 & PPNN(4),PTNN(4),PTOTNN(4),PPSUB(4),PTSUB(4),
3551 & PPTCMS(4),PTTCMS(4),PPTMP(4),PTTMP(4),
3552 & KPRON(15),ISINGL(2000)
3553
3554* initial values for max. number of phojet scatterings and dtunuc chains
3555* to be fragmented with one pyexec call
3556 DATA MXPHFR,MXDTFR /10,100/
3557
3558 IREJ = 0
3559* pointer to first parton of the first chain in dtevt common
3560 NPOINT(3) = NHKK+1
3561* special flag for double-Pomeron statistics
3562 IPOPO = 1
3563* counter for low-mass (DTUNUC) interactions
3564 NDTUSC = 0
3565* counter for interactions treated by PHOJET
3566 NPHOSC = 0
3567
3568* scan interactions for single nucleon-nucleon interactions
3569* (this has to be checked here because Cronin modifies parton momenta)
3570 NC = NPOINT(2)
3571 IF (NCSY.GT.2000) STOP ' DT_EVENTB: NCSY > 2000 ! '
3572 DO 8 I=1,NCSY
3573 ISINGL(I) = 0
3574 MOP = JMOHKK(1,NC)
3575 MOT = JMOHKK(1,NC+1)
3576 DIFF1 = ABS(PHKK(4,MOP)-PHKK(4, NC)-PHKK(4,NC+2))
3577 DIFF2 = ABS(PHKK(4,MOT)-PHKK(4,NC+1)-PHKK(4,NC+3))
3578 IF ((DIFF1.LT.TINY10).AND.(DIFF2.LT.TINY10)) ISINGL(I) = 1
3579 NC = NC+4
3580 8 CONTINUE
3581
3582* multiple scattering of chain ends
3583 IF ((IP.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(1)
3584 IF ((IT.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(2)
3585
3586* switch to PHOJET-settings for JETSET parameter
3587 CALL DT_INITJS(1)
3588
3589* loop over nucleon-nucleon interaction
3590 NC = NPOINT(2)
3591 DO 2 I=1,NCSY
3592*
3593* pick up one nucleon-nucleon interaction from DTEVT1
3594* ppnn / ptnn - momenta of the interacting nucleons (cms)
3595* ptotnn - total momentum of the interacting nucleons (cms)
3596* pp1,2 / pt1,2 - momenta of the four partons
3597* pp / pt - total momenta of the proj / targ partons
3598* ptot - total momentum of the four partons
3599 MOP = JMOHKK(1,NC)
3600 MOT = JMOHKK(1,NC+1)
3601 DO 3 K=1,4
3602 PPNN(K) = PHKK(K,MOP)
3603 PTNN(K) = PHKK(K,MOT)
3604 PTOTNN(K) = PPNN(K)+PTNN(K)
3605 PP1(K) = PHKK(K,NC)
3606 PT1(K) = PHKK(K,NC+1)
3607 PP2(K) = PHKK(K,NC+2)
3608 PT2(K) = PHKK(K,NC+3)
3609 PP(K) = PP1(K)+PP2(K)
3610 PT(K) = PT1(K)+PT2(K)
3611 PTOT(K) = PP(K)+PT(K)
3612 3 CONTINUE
3613*
3614*-----------------------------------------------------------------------
3615* this is a complete nucleon-nucleon interaction
3616*
3617 IF (ISINGL(I).EQ.1) THEN
3618*
3619* initialize PHOJET-variables for remnant/valence-partons
3620 IHFLD(1,1) = 0
3621 IHFLD(1,2) = 0
3622 IHFLD(2,1) = 0
3623 IHFLD(2,2) = 0
3624 IHFLS(1) = 1
3625 IHFLS(2) = 1
3626* save current settings of PHOJET process and min. bias flags
3627 DO 9 K=1,11
3628 KPRON(K) = IPRON(K,1)
3629 9 CONTINUE
3630 ISWSAV = ISWMDL(2)
3631*
3632* check if forced sampling of diffractive interaction requested
3633 IF (ISINGD.LT.-1) THEN
3634 DO 90 K=1,11
3635 IPRON(K,1) = 0
3636 90 CONTINUE
3637 IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-3)) IPRON(5,1) = 1
3638 IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-4)) IPRON(6,1) = 1
3639 IF (ISINGD.EQ.-5) IPRON(4,1) = 1
3640 ENDIF
3641*
3642* for photons: a direct/anomalous interaction is not sampled
3643* in PHOJET but already in Glauber-formalism. Here we check if such
3644* an interaction is requested
3645 IF (IJPROJ.EQ.7) THEN
3646* first switch off direct interactions
3647 IPRON(8,1) = 0
3648* this is a direct interactions
3649 IF (IDIREC.EQ.1) THEN
3650 DO 12 K=1,11
3651 IPRON(K,1) = 0
3652 12 CONTINUE
3653 IPRON(8,1) = 1
3654* this is an anomalous interactions
3655* (iswmdl(2) = 0 only hard int. generated ( = 1 min. bias) )
3656 ELSEIF (IDIREC.EQ.2) THEN
3657 ISWMDL(2) = 0
3658 ENDIF
3659 ELSE
3660 IF (IDIREC.NE.0) STOP ' DT_EVENTB: IDIREC > 0 ! '
3661 ENDIF
3662*
3663* make sure that total momenta of partons, pp and pt, are on mass
3664* shell (Cronin may have srewed this up..)
3665 CALL DT_MASHEL(PP,PT,PHKK(5,MOP),PHKK(5,MOT),PPNN,PTNN,IR1)
3666 IF (IR1.NE.0) THEN
3667 IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A)')
3668 & 'EVENTB: mass shell correction rejected'
3669 GOTO 9999
3670 ENDIF
3671*
3672* initialize the incoming particles in PHOJET
3673 IF ((IP.EQ.1).AND.(IJPROJ.EQ.7)) THEN
3674 CALL PHO_SETPAR(1,22,0,VIRT)
3675 ELSE
3676 CALL PHO_SETPAR(1,IDHKK(MOP),0,ZERO)
3677 ENDIF
3678 CALL PHO_SETPAR(2,IDHKK(MOT),0,ZERO)
3679*
3680* initialize rejection loop counter for anomalous processes
3681 IRJANO = 0
3682 800 CONTINUE
3683 IRJANO = IRJANO+1
3684*
3685* temporary fix for ifano problem
3686 IFANO(1) = 0
3687 IFANO(2) = 0
3688*
3689* generate complete hadron/nucleon/photon-nucleon event with PHOJET
3690 CALL PHO_EVENT(2,PPNN,PTNN,DUM,IREJ1)
3691*
3692* for photons: special consistency check for anomalous interactions
3693 IF (IJPROJ.EQ.7) THEN
3694 IF (IRJANO.LT.30) THEN
3695 IF (IFANO(1).NE.0) THEN
3696* here, an anomalous interaction was generated. Check if it
3697* was also requested. Otherwise reject this event.
3698 IF (IDIREC.EQ.0) GOTO 800
3699 ELSE
3700* here, an anomalous interaction was not generated. Check if it
3701* was requested in which case we need to reject this event.
3702 IF (IDIREC.EQ.2) GOTO 800
3703 ENDIF
3704 ELSE
3705 WRITE(LOUT,*) ' DT_EVENTB: Warning! IRJANO > 30 ',
3706 & IRJANO,IDIREC,NEVHKK
3707 ENDIF
3708 ENDIF
3709*
3710* copy back original settings of PHOJET process and min. bias flags
3711 DO 10 K=1,11
3712 IPRON(K,1) = KPRON(K)
3713 10 CONTINUE
3714 ISWMDL(2) = ISWSAV
3715*
3716* check if PHOJET has rejected this event
3717 IF (IREJ1.NE.0) THEN
3718C IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A,I4)')
3719 WRITE(LOUT,'(1X,A,I4)')
3720 & 'EVENTB: chain system rejected',IDIREC
3721 CALL PHO_PREVNT(0)
3722 GOTO 9999
3723 ENDIF
3724*
3725* copy partons and strings from PHOJET common back into DTEVT for
3726* external fragmentation
3727 MO1 = NC
3728 MO2 = NC+3
3729*! uncomment this line for internal phojet-fragmentation
3730C CALL DT_GETFSP(MO1,MO2,PPNN,PTNN,-1)
3731 NPHOSC = NPHOSC+1
3732 CALL DT_GETPJE(MO1,MO2,PPNN,PTNN,-1,NPHOSC,IREJ1)
3733 IF (IREJ1.NE.0) THEN
3734 IF (IOULEV(1).GT.0)
3735 & WRITE(LOUT,'(1X,A,I4)') 'EVENTB: chain system rejected 1'
3736 GOTO 9999
3737 ENDIF
3738*