Update to version 3.0-5
[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)
172 COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
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
2102 IF ((IP.GT.1).AND.(MKCRON.GT.0)) THEN
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
2169 PARAMETER (NMXHKK=200000)
2170 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
2171 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
2172 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
2173* extended event history
2174 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
2175 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
2176 & IHIST(2,NMXHKK)
2177* particle properties (BAMJET index convention)
2178 CHARACTER*8 ANAME
2179 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2180 & IICH(210),IIBAR(210),K1(210),K2(210)
2181* properties of interacting particles
2182 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2183* Lorentz-parameters of the current interaction
2184 COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
2185 & UMO,PPCM,EPROJ,PPROJ
2186* flags for input different options
2187 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2188 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2189 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2190* flags for particle decays
2191 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
2192 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
2193 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
2194* cuts for variable energy runs
2195 COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
2196* Glauber formalism: flags and parameters for statistics
2197 LOGICAL LPROD
2198 CHARACTER*8 CGLB
2199 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
2200
2201 DIMENSION WHAT(6)
2202
2203 IREJ = 0
2204 ILOOP = 0
2205 100 CONTINUE
2206 IF (ILOOP.EQ.4) THEN
2207 WRITE(LOUT,1000) NEVHKK
2208 1000 FORMAT(1X,'KKINC: event ',I8,' rejected!')
2209 GOTO 9999
2210 ENDIF
2211 ILOOP = ILOOP+1
2212
2213* variable energy-runs, recalculate parameters for LT's
2214 IF ((ABS(VAREHI).GT.ZERO).OR.(IOGLB.EQ.100)) THEN
2215 PDUM = ZERO
2216 CDUM = ZERO
2217 CALL DT_LTINI(IDP,1,EPN,PDUM,CDUM,1)
2218 ENDIF
2219 IF (EPN.GT.EPROJ) THEN
2220 WRITE(LOUT,'(A,E9.3,2A,E9.3,A)')
2221 & ' Requested energy (',EPN,'GeV) exceeds',
2222 & ' initialization energy (',EPROJ,'GeV) !'
2223 STOP
2224 ENDIF
2225
2226* re-initialize /DTPRTA/
2227 IP = NPMASS
2228 IPZ = NPCHAR
2229 IT = NTMASS
2230 ITZ = NTCHAR
2231 IJPROJ = IDP
2232 IBPROJ = IIBAR(IJPROJ)
2233
2234* calculate nuclear potentials (common /DTNPOT/)
2235 CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
2236
2237* initialize treatment for residual nuclei
2238 CALL DT_RESNCL(EPN,NLOOP,1)
2239
2240* sample hadron/nucleus-nucleus interaction
2241 CALL DT_KKEVNT(KKMAT,IREJ1)
2242 IF (IREJ1.GT.0) THEN
2243 IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 1 in KKINC'
2244 GOTO 9999
2245 ENDIF
2246
2247 IF ((NPMASS.GT.1).OR.(NTMASS.GT.1)) THEN
2248
2249* intranuclear cascade of final state particles for KTAUGE generations
2250* of secondaries
2251 CALL DT_FOZOCA(LFZC,IREJ1)
2252 IF (IREJ1.GT.0) THEN
2253 IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2 in KKINC'
2254 GOTO 9999
2255 ENDIF
2256
2257* baryons unable to escape the nuclear potential are treated as
2258* excited nucleons (ISTHKK=15,16)
2259 CALL DT_SCN4BA
2260
2261* decay of resonances produced in intranuclear cascade processes
2262**sr 15-11-95 should be obsolete
2263C IF (LFZC) CALL DT_DECAY1
2264
2265 101 CONTINUE
2266* treatment of residual nuclei
2267 CALL DT_RESNCL(EPN,NLOOP,2)
2268
2269* evaporation / fission / fragmentation
2270* (if intranuclear cascade was sampled only)
2271 IF (LFZC) THEN
2272 CALL DT_FICONF(IJPROJ,IP,IPZ,IT,ITZ,NLOOP,IREJ1)
2273 IF (IREJ1.GT.1) GOTO 101
2274 IF (IREJ1.EQ.1) GOTO 100
2275 ENDIF
2276
2277 ENDIF
2278
2279* rejection of unphysical configurations
2280 CALL DT_REJUCO(1,IREJ1)
2281 IF (IREJ1.GT.0) THEN
2282 IF (IOULEV(1).GT.0)
2283 & WRITE(LOUT,*) 'rejected 3 in KKINC: too large x'
2284 GOTO 100
2285 ENDIF
2286
2287* transform finale state into Lab.
2288 IFLAG = 2
2289 CALL DT_BEAMPR(WHAT,DUM,IFLAG)
2290 IF ((IFRAME.EQ.1).AND.(IFLAG.EQ.-1)) CALL DT_LT2LAB
2291
2292 IF (IPI0.EQ.1) CALL DT_DECPI0
2293
2294C IF (NEVHKK.EQ.5) CALL DT_EVTOUT(4)
2295
2296 RETURN
2297 9999 CONTINUE
2298 IREJ = 1
2299 RETURN
2300 END
2301
2302*$ CREATE DT_DEFAUL.FOR
2303*COPY DT_DEFAUL
2304*
2305*===defaul=============================================================*
2306*
2307 SUBROUTINE DT_DEFAUL(EPN,PPN)
2308
2309************************************************************************
2310* Variables are set to default values. *
2311* This version dated 8.5.95 is written by S. Roesler. *
2312************************************************************************
2313
2314 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2315 SAVE
2316 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
2317 PARAMETER (TWOPI = 6.283185307179586454D+00)
2318
2319* particle properties (BAMJET index convention)
2320 CHARACTER*8 ANAME
2321 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2322 & IICH(210),IIBAR(210),K1(210),K2(210)
2323* nuclear potential
2324 LOGICAL LFERMI
2325 COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
2326 & EBINDP(2),EBINDN(2),EPOT(2,210),
2327 & ETACOU(2),ICOUL,LFERMI
2328* interface HADRIN-DPM
2329 COMMON /HNTHRE/ EHADTH,EHADLO,EHADHI,INTHAD,IDXTA
2330* central particle production, impact parameter biasing
2331 COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
2332* properties of interacting particles
2333 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2334* properties of photon/lepton projectiles
2335 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
2336 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2337* emulsion treatment
2338 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2339 & NCOMPO,IEMUL
2340* parameter for intranuclear cascade
2341 LOGICAL LPAULI
2342 COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
2343* various options for treatment of partons (DTUNUC 1.x)
2344* (chain recombination, Cronin,..)
2345 LOGICAL LCO2CR,LINTPT
2346 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
2347 & LCO2CR,LINTPT
2348* threshold values for x-sampling (DTUNUC 1.x)
2349 COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
2350 & SSMIMQ,VVMTHR
2351* flags for input different options
2352 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2353 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2354 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2355* n-n cross section fluctuations
2356 PARAMETER (NBINS = 1000)
2357 COMMON /DTXSFL/ FLUIXX(NBINS),IFLUCT
2358* flags for particle decays
2359 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
2360 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
2361 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
2362* diquark-breaking mechanism
2363 COMMON /DTDIQB/ DBRKR(3,8),DBRKA(3,8),CHAM1,CHAM3,CHAB1,CHAB3
2364* nucleon-nucleon event-generator
2365 CHARACTER*8 CMODEL
2366 LOGICAL LPHOIN
2367 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
2368* flags for diffractive interactions (DTUNUC 1.x)
2369 COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
2370* VDM parameter for photon-nucleus interactions
2371 COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
2372* Glauber formalism: flags and parameters for statistics
2373 LOGICAL LPROD
2374 CHARACTER*8 CGLB
2375 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
2376* kinematical cuts for lepton-nucleus interactions
2377 COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
2378 & Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
2379* flags for activated histograms
2380 COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
2381* cuts for variable energy runs
2382 COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
2383* parameters for hA-diffraction
2384 COMMON /DTDIHA/ DIBETA,DIALPH
2385* LEPTO
2386 REAL RPPN
2387 COMMON /LEPTOI/ RPPN,LEPIN,INTER
2388* steering flags for qel neutrino scattering modules
2389 COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
2390* event flag
2391 COMMON /DTEVNO/ NEVENT,ICASCA
2392
2393 DATA POTMES /0.002D0/
2394
2395* common /DTNPOT/
2396 DO 10 I=1,2
2397 PFERMP(I) = ZERO
2398 PFERMN(I) = ZERO
2399 EBINDP(I) = ZERO
2400 EBINDN(I) = ZERO
2401 DO 11 J=1,210
2402 EPOT(I,J) = ZERO
2403 11 CONTINUE
2404* nucleus independent meson potential
2405 EPOT(I,13) = POTMES
2406 EPOT(I,14) = POTMES
2407 EPOT(I,15) = POTMES
2408 EPOT(I,16) = POTMES
2409 EPOT(I,23) = POTMES
2410 EPOT(I,24) = POTMES
2411 EPOT(I,25) = POTMES
2412 10 CONTINUE
2413 FERMOD = 0.55D0
2414 ETACOU(1) = ZERO
2415 ETACOU(2) = ZERO
2416 ICOUL = 1
2417 LFERMI = .TRUE.
2418
2419* common /HNTHRE/
2420 EHADTH = -99.0D0
2421 EHADLO = 4.06D0
2422 EHADHI = 6.0D0
2423 INTHAD = 1
2424 IDXTA = 2
2425
2426* common /DTIMPA/
2427 ICENTR = 0
2428 BIMIN = ZERO
2429 BIMAX = 1.0D10
2430 XSFRAC = 1.0D0
2431
2432* common /DTPRTA/
2433 IP = 1
2434 IPZ = 1
2435 IT = 1
2436 ITZ = 1
2437 IJPROJ = 1
2438 IBPROJ = 1
2439 IJTARG = 1
2440 IBTARG = 1
2441* common /DTGPRO/
2442 VIRT = ZERO
2443 DO 14 I=1,4
2444 PGAMM(I) = ZERO
2445 PLEPT0(I) = ZERO
2446 PLEPT1(I) = ZERO
2447 PNUCL(I) = ZERO
2448 14 CONTINUE
2449 IDIREC = 0
2450
2451* common /DTFOTI/
2452**sr 7.4.98: changed after corrected B-sampling
2453C TAUFOR = 4.4D0
2454 TAUFOR = 3.5D0
2455 KTAUGE = 25
2456 ITAUVE = 1
2457 INCMOD = 1
2458 LPAULI = .TRUE.
2459
2460* common /DTCHAI/
2461 SEASQ = ONE
2462 MKCRON = 1
2463 CRONCO = 0.64D0
2464 ISICHA = 0
2465 CUTOF = 100.0D0
2466 LCO2CR = .FALSE.
2467 IRECOM = 1
2468 LINTPT = .TRUE.
2469
2470* common /DTXCUT/
2471* definition of soft quark distributions
2472 XSEACU = 0.05D0
2473 UNON = 2.0D0
2474 UNOM = 1.5D0
2475 UNOSEA = 5.0D0
2476* cutoff parameters for x-sampling
2477 CVQ = 1.0D0
2478 CDQ = 2.0D0
2479C CSEA = 0.3D0
2480 CSEA = 0.1D0
2481 SSMIMA = 1.2D0
2482 SSMIMQ = SSMIMA**2
2483 VVMTHR = 2.0D0
2484
2485* common /DTXSFL/
2486 IFLUCT = 0
2487
2488* common /DTFRPA/
2489 PDB = 0.15D0
2490 PDBSEA(1) = 0.0D0
2491 PDBSEA(2) = 0.0D0
2492 PDBSEA(3) = 0.0D0
2493 ISIG0 = 0
2494 IPI0 = 0
2495 NMSTU = 0
2496 NPARU = 0
2497 NMSTJ = 0
2498 NPARJ = 0
2499
2500* common /DTDIQB/
2501 DO 15 I=1,8
2502 DBRKR(1,I) = 5.0D0
2503 DBRKR(2,I) = 5.0D0
2504 DBRKR(3,I) = 10.0D0
2505 DBRKA(1,I) = ZERO
2506 DBRKA(2,I) = ZERO
2507 DBRKA(3,I) = ZERO
2508 15 CONTINUE
2509 CHAM1 = 0.2D0
2510 CHAM3 = 0.5D0
2511 CHAB1 = 0.7D0
2512 CHAB3 = 1.0D0
2513
2514* common /DTFLG3/
2515 ISINGD = 0
2516 IDOUBD = 0
2517 IFLAGD = 0
2518 IDIFF = 0
2519
2520* common /DTMODL/
2521 MCGENE = 2
2522 CMODEL(1) = 'DTUNUC '
2523 CMODEL(2) = 'PHOJET '
2524 CMODEL(3) = 'LEPTO '
2525 CMODEL(4) = 'QNEUTRIN'
2526 LPHOIN = .TRUE.
2527 ELOJET = 5.0D0
2528
2529* common /DTLCUT/
2530 ECMIN = 3.5D0
2531 ECMAX = 1.0D10
2532 XBJMIN = ZERO
2533 ELMIN = ZERO
2534 EGMIN = ZERO
2535 EGMAX = 1.0D10
2536 YMIN = TINY10
2537 YMAX = 0.999D0
2538 Q2MIN = TINY10
2539 Q2MAX = 10.0D0
2540 THMIN = ZERO
2541 THMAX = TWOPI
2542 Q2LI = ZERO
2543 Q2HI = 1.0D10
2544 ECMLI = ZERO
2545 ECMHI = 1.0D10
2546
2547* common /DTVDMP/
2548 RL2 = 2.0D0
2549 INTRGE(1) = 1
2550 INTRGE(2) = 3
2551 IDPDF = 2212
2552 MODEGA = 4
2553 ISHAD(1) = 1
2554 ISHAD(2) = 1
2555 ISHAD(3) = 1
2556 EPSPOL = ZERO
2557
2558* common /DTGLGP/
2559 JSTATB = 1000
2560 JBINSB = 49
2561 CGLB = ' '
2562 IF (ITRSPT.EQ.1) THEN
2563 IOGLB = 100
2564 ELSE
2565 IOGLB = 0
2566 ENDIF
2567 LPROD = .TRUE.
2568
2569* common /DTHIS3/
2570 DO 16 I=1,50
2571 IHISPP(I) = 0
2572 IHISXS(I) = 0
2573 16 CONTINUE
2574 IXSTBL = 0
2575
2576* common /DTVARE/
2577 VARELO = ZERO
2578 VAREHI = ZERO
2579 VARCLO = ZERO
2580 VARCHI = ZERO
2581
2582* common /DTDIHA/
2583 DIBETA = -1.0D0
2584 DIALPH = ZERO
2585
2586* common /LEPTOI/
2587 RPPN = 0.0
2588 LEPIN = 0
2589 INTER = 0
2590
2591* common /QNEUTO/
2592 NEUTYP = 1
2593 NEUDEC = 0
2594
2595* common /DTEVNO/
2596 NEVENT = 1
2597 IF (ITRSPT.EQ.1) THEN
2598 ICASCA = 1
2599 ELSE
2600 ICASCA = 0
2601 ENDIF
2602
2603* default Lab.-energy
2604 EPN = 200.0D0
2605 PPN = SQRT((EPN-AAM(IJPROJ))*(EPN+AAM(IJPROJ)))
2606
2607 RETURN
2608 END
2609
2610*$ CREATE DT_AAEVT.FOR
2611*COPY DT_AAEVT
2612*
2613*===aaevt==============================================================*
2614*
2615 SUBROUTINE DT_AAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
2616 & IDP,IGLAU)
2617
2618************************************************************************
2619* This version dated 22.03.96 is written by S. Roesler. *
2620************************************************************************
2621
2622 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2623 SAVE
2624 PARAMETER ( LINP = 10 ,
2625 & LOUT = 6 ,
2626 & LDAT = 9 )
2627
2628 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2629* emulsion treatment
2630 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2631 & NCOMPO,IEMUL
2632* event flag
2633 COMMON /DTEVNO/ NEVENT,ICASCA
2634
2635 CHARACTER*8 DATE,HHMMSS
2636 DIMENSION IDMNYR(3)
2637
2638 KKMAT = 1
2639 NMSG = MAX(NEVTS/100,1)
2640
2641* initialization of run-statistics and histograms
2642 CALL DT_STATIS(1)
2643 CALL PHO_PHIST(1000,DUM)
2644
2645* initialization of Glauber-formalism
2646 IF (NCOMPO.LE.0) THEN
2647 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
2648 ELSE
2649 DO 1 I=1,NCOMPO
2650 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2651 1 CONTINUE
2652 ENDIF
2653 CALL DT_SIGEMU
2654
2655 CALL IDATE(IDMNYR)
2656 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2657 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
2658 CALL ITIME(IDMNYR)
2659 WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
2660 & IDMNYR(1),IDMNYR(2),IDMNYR(3)
2661 WRITE(LOUT,1001) DATE,HHMMSS
2662 1001 FORMAT(/,' DT_AAEVT: Initialisation finished. ( Date: ',A8,
2663 & ' Time: ',A8,' )')
2664
2665* generate NEVTS events
2666 DO 2 IEVT=1,NEVTS
2667
2668* print run-status message
2669 IF (MOD(IEVT,NMSG).EQ.0) THEN
2670 CALL IDATE(IDMNYR)
2671 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2672 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
2673 CALL ITIME(IDMNYR)
2674 WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
2675 & IDMNYR(1),IDMNYR(2),IDMNYR(3)
2676 WRITE(LOUT,1000) IEVT-1,NEVTS,DATE,HHMMSS
2677 1000 FORMAT(/,1X,I8,' out of ',I8,' events sampled ( Date: ',A,
2678 & ' Time: ',A,' )',/)
2679C WRITE(LOUT,1000) IEVT-1
2680C1000 FORMAT(1X,I8,' events sampled')
2681 ENDIF
2682 NEVENT = IEVT
2683* treat nuclear emulsions
2684 IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
2685* composite targets only
2686 KKMAT = -KKMAT
2687* sample this event
2688 CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,IREJ)
2689
2690 CALL PHO_PHIST(2000,DUM)
2691
2692 2 CONTINUE
2693
2694* print run-statistics and histograms to output-unit 6
2695 CALL PHO_PHIST(3000,DUM)
2696 CALL DT_STATIS(2)
2697
2698 RETURN
2699 END
2700
2701*$ CREATE DT_LAEVT.FOR
2702*COPY DT_LAEVT
2703*
2704*===laevt==============================================================*
2705*
2706 SUBROUTINE DT_LAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
2707 & IDP,IGLAU)
2708
2709************************************************************************
2710* Interface to run DPMJET for lepton-nucleus interactions. *
2711* Kinematics is sampled using the equivalent photon approximation *
2712* Based on GPHERA-routine by R. Engel. *
2713* This version dated 23.03.96 is written by S. Roesler. *
2714************************************************************************
2715
2716 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2717 SAVE
2718 PARAMETER ( LINP = 10 ,
2719 & LOUT = 6 ,
2720 & LDAT = 9 )
2721 PARAMETER (TINY10=1.0D-10,TINY4=1.0D-4,
2722 & ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,THREE=3.0D0)
2723 PARAMETER (TWOPI = 6.283185307179586454D+00,
2724 & PI = TWOPI/TWO,
2725 & ALPHEM = ONE/137.0D0)
2726
2727C CHARACTER*72 HEADER
2728
2729* particle properties (BAMJET index convention)
2730 CHARACTER*8 ANAME
2731 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2732 & IICH(210),IIBAR(210),K1(210),K2(210)
2733* event history
2734 PARAMETER (NMXHKK=200000)
2735 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
2736 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
2737 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
2738* extended event history
2739 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
2740 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
2741 & IHIST(2,NMXHKK)
2742* kinematical cuts for lepton-nucleus interactions
2743 COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
2744 & Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
2745* properties of interacting particles
2746 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
2747* properties of photon/lepton projectiles
2748 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
2749* kinematics at lepton-gamma vertex
2750 COMMON /DTLGVX/ PPL0(4),PPL1(4),PPG(4),PPA(4)
2751* flags for activated histograms
2752 COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
2753 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2754* emulsion treatment
2755 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2756 & NCOMPO,IEMUL
2757* Glauber formalism: cross sections
2758 COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
2759 & XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
2760 & XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
2761 & XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
2762 & XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
2763 & XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
2764 & XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
2765 & XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
2766 & XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
2767 & BSLOPE,NEBINI,NQBINI
2768* nucleon-nucleon event-generator
2769 CHARACTER*8 CMODEL
2770 LOGICAL LPHOIN
2771 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
2772* flags for input different options
2773 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
2774 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
2775 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
2776* event flag
2777 COMMON /DTEVNO/ NEVENT,ICASCA
2778
2779 DIMENSION XDUMB(40),BGTA(4)
2780
2781* LEPTO
2782 IF (MCGENE.EQ.3) THEN
2783 STOP ' This version does not contain LEPTO !'
2784 ENDIF
2785
2786 KKMAT = 1
2787 NMSG = MAX(NEVTS/10,1)
2788
2789* mass of incident lepton
2790 AMLPT = AAM(IDP)
2791 AMLPT2 = AMLPT**2
2792 IDPPDG = IDT_IPDGHA(IDP)
2793
2794* consistency of kinematical limits
2795 Q2MIN = MAX(Q2MIN,TINY10)
2796 Q2MAX = MAX(Q2MAX,TINY10)
2797 YMIN = MIN(MAX(YMIN,TINY10),0.999D0)
2798 YMAX = MIN(MAX(YMAX,TINY10),0.999D0)
2799
2800* total energy of the lepton-nucleon system
2801 PTOTLN = SQRT( (PLEPT0(1)+PNUCL(1))**2+(PLEPT0(2)+PNUCL(2))**2
2802 & +(PLEPT0(3)+PNUCL(3))**2 )
2803 ETOTLN = PLEPT0(4)+PNUCL(4)
2804 ECMLN = SQRT((ETOTLN-PTOTLN)*(ETOTLN+PTOTLN))
2805 ECMAX = MIN(ECMAX,ECMLN)
2806 WRITE(LOUT,1003) ECMIN,ECMAX,YMIN,YMAX,Q2MIN,Q2MAX,EGMIN,
2807 & THMIN,THMAX,ELMIN
2808 1003 FORMAT(1X,'LAEVT:',16X,'kinematical cuts',/,22X,
2809 & '------------------',/,9X,'W (min) =',
2810 & F7.1,' GeV (max) =',F7.1,' GeV',/,9X,'y (min) =',
2811 & F7.3,8X,'(max) =',F7.3,/,9X,'Q^2 (min) =',F7.1,
2812 & ' GeV^2 (max) =',F7.1,' GeV^2',/,' (Lab) E_g (min) ='
2813 & ,F7.1,' GeV',/,' (Lab) theta (min) =',F7.4,8X,'(max) =',
2814 & F7.4,' for E_lpt >',F7.1,' GeV',/)
2815
2816* Lorentz-parameter for transf. into Lab
2817 BGTA(1) = PNUCL(1)/AAM(1)
2818 BGTA(2) = PNUCL(2)/AAM(1)
2819 BGTA(3) = PNUCL(3)/AAM(1)
2820 BGTA(4) = PNUCL(4)/AAM(1)
2821* LT of incident lepton into Lab and dump it in DTEVT1
2822 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
2823 & PLEPT0(1),PLEPT0(2),PLEPT0(3),PLEPT0(4),
2824 & PLTOT,PPL0(1),PPL0(2),PPL0(3),PPL0(4))
2825 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
2826 & PNUCL(1),PNUCL(2),PNUCL(3),PNUCL(4),
2827 & PLTOT,PPA(1),PPA(2),PPA(3),PPA(4))
2828* maximum energy of photon nucleon system
2829 PTOTGN = SQRT((YMAX*PPL0(1)+PPA(1))**2+(YMAX*PPL0(2)+PPA(2))**2
2830 & +(YMAX*PPL0(3)+PPA(3))**2)
2831 ETOTGN = YMAX*PPL0(4)+PPA(4)
2832 EGNMAX = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
2833 EGNMAX = MIN(EGNMAX,ECMAX)
2834* minimum energy of photon nucleon system
2835 PTOTGN = SQRT((YMIN*PPL0(1)+PPA(1))**2+(YMIN*PPL0(2)+PPA(2))**2
2836 & +(YMIN*PPL0(3)+PPA(3))**2)
2837 ETOTGN = YMIN*PPL0(4)+PPA(4)
2838 EGNMIN = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
2839 EGNMIN = MAX(EGNMIN,ECMIN)
2840
2841* limits for Glauber-initialization
2842 Q2LI = Q2MIN
2843 Q2HI = MAX(Q2LI,MIN(Q2HI,Q2MAX))
2844 ECMLI = MAX(EGNMIN,THREE)
2845 ECMHI = EGNMAX
2846 WRITE(LOUT,1004) EGNMIN,EGNMAX,ECMLI,ECMHI,Q2LI,Q2HI
2847 1004 FORMAT(1X,'resulting limits:',/,9X,'W (min) =',F7.1,
2848 & ' GeV (max) =',F7.1,' GeV',/,/,' limits for ',
2849 & 'Glauber-initialization:',/,9X,'W (min) =',F7.1,
2850 & ' GeV (max) =',F7.1,' GeV',/,9X,'Q^2 (min) =',F7.1,
2851 & ' GeV^2 (max) =',F7.1,' GeV^2',/)
2852* initialization of Glauber-formalism
2853 IF (NCOMPO.LE.0) THEN
2854 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
2855 ELSE
2856 DO 9 I=1,NCOMPO
2857 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2858 9 CONTINUE
2859 ENDIF
2860 CALL DT_SIGEMU
2861
2862* initialization of run-statistics and histograms
2863 CALL DT_STATIS(1)
2864 CALL PHO_PHIST(1000,DUM)
2865
2866* maximum photon-nucleus cross section
2867 I1 = 1
2868 I2 = 1
2869 RAT = ONE
2870 IF (EGNMAX.GE.ECMNN(NEBINI)) THEN
2871 I1 = NEBINI
2872 I2 = NEBINI
2873 RAT = ONE
2874 ELSEIF (EGNMAX.GT.ECMNN(1)) THEN
2875 DO 5 I=2,NEBINI
2876 IF (EGNMAX.LT.ECMNN(I)) THEN
2877 I1 = I-1
2878 I2 = I
2879 RAT = (EGNMAX-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
2880 GOTO 6
2881 ENDIF
2882 5 CONTINUE
2883 6 CONTINUE
2884 ENDIF
2885 SIGMAX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
2886 EGNXX = EGNMAX
2887 I1 = 1
2888 I2 = 1
2889 RAT = ONE
2890 IF (EGNMIN.GE.ECMNN(NEBINI)) THEN
2891 I1 = NEBINI
2892 I2 = NEBINI
2893 RAT = ONE
2894 ELSEIF (EGNMIN.GT.ECMNN(1)) THEN
2895 DO 7 I=2,NEBINI
2896 IF (EGNMIN.LT.ECMNN(I)) THEN
2897 I1 = I-1
2898 I2 = I
2899 RAT = (EGNMIN-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
2900 GOTO 8
2901 ENDIF
2902 7 CONTINUE
2903 8 CONTINUE
2904 ENDIF
2905 SIGXX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
2906 IF (SIGXX.GT.SIGMAX) EGNXX = EGNMIN
2907 SIGMAX = MAX(SIGMAX,SIGXX)
2908 WRITE(LOUT,'(9X,A,F8.3,A)') 'Sigma_tot (max) =',SIGMAX,' mb'
2909
2910* plot photon flux table
2911 AYMIN = LOG(YMIN)
2912 AYMAX = LOG(YMAX)
2913 AYRGE = AYMAX-AYMIN
2914 MAXTAB = 50
2915 ADY = LOG(YMAX/YMIN)/DBLE(MAXTAB-1)
2916C WRITE(LOUT,'(/,1X,A)') 'LAEVT: photon flux '
2917 DO 1 I=1,MAXTAB
2918 Y = EXP(AYMIN+ADY*DBLE(I-1))
2919 Q2LOW = MAX(Q2MIN,AMLPT2*Y**2/(ONE-Y))
2920 FF1 = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2921 & -TWO*AMLPT2*Y*(ONE/Q2LOW-ONE/Q2MAX))
2922 FF2 = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
2923 & -TWO*(ONE-Y)/Y*(ONE-Q2LOW/Q2MAX))
2924C WRITE(LOUT,'(5X,3E15.4)') Y,FF1,FF2
2925 1 CONTINUE
2926
2927* maximum residual weight for flux sampling (dy/y)
2928 YY = YMIN
2929 Q2LOW = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
2930 WGHMAX = (ONE+(ONE-YY)**2)*LOG(Q2MAX/Q2LOW)
2931 & -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
2932
2933 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY0)
2934 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY1)
2935 CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY2)
2936 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ0)
2937 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ1)
2938 CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ2)
2939 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE0)
2940 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE1)
2941 CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE2)
2942 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU0)
2943 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU1)
2944 CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU2)
2945 XBLOW = 0.001D0
2946 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX0)
2947 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX1)
2948 CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX2)
2949
2950 ITRY = 0
2951 ITRW = 0
2952 NC0 = 0
2953 NC1 = 0
2954
2955* generate events
2956 DO 2 IEVT=1,NEVTS
2957 IF (MOD(IEVT,NMSG).EQ.0) THEN
2958C OPEN(LDAT,FILE='/scrtch3/hr/sroesler/statusd5.out',
2959C & STATUS='UNKNOWN')
2960 WRITE(LOUT,'(1X,I8,A)') IEVT-1,' events sampled'
2961C CLOSE(LDAT)
2962 ENDIF
2963 NEVENT = IEVT
2964
2965 100 CONTINUE
2966 ITRY = ITRY+1
2967
2968* sample y
2969 101 CONTINUE
2970 ITRW = ITRW+1
2971 YY = EXP(AYRGE*DT_RNDM(RAT)+AYMIN)
2972 Q2LOW = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
2973 Q2LOG = LOG(Q2MAX/Q2LOW)
2974 WGH = (ONE+(ONE-YY)**2)*Q2LOG
2975 & -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
2976 IF (WGHMAX.LT.WGH) WRITE(LOUT,1000) YY,WGHMAX,WGH
2977 1000 FORMAT(1X,'LAEVT: weight error!',3E12.5)
2978 IF (DT_RNDM(YY)*WGHMAX.GT.WGH) GOTO 101
2979
2980* sample Q2
2981 YEFF = ONE+(ONE-YY)**2
2982 102 CONTINUE
2983 Q2 = Q2LOW*EXP(Q2LOG*DT_RNDM(YY))
2984 WGH = (YEFF-TWO*(ONE-YY)*Q2LOW/Q2)/YEFF
2985 IF (WGH.LT.DT_RNDM(Q2)) GOTO 102
2986
2987c NC0 = NC0+1
2988c CALL DT_FILHGR(YY,ONE,IHFLY0,NC0)
2989c CALL DT_FILHGR(Q2,ONE,IHFLQ0,NC0)
2990
2991* kinematics at lepton-photon vertex
2992* scattered electron
2993 YQ2 = SQRT((ONE-YY)*Q2)
2994 Q2E = Q2/(4.0D0*PLEPT0(4))
2995 E1Y = (ONE-YY)*PLEPT0(4)
2996 CALL DT_DSFECF(SIF,COF)
2997 PLEPT1(1) = YQ2*COF
2998 PLEPT1(2) = YQ2*SIF
2999 PLEPT1(3) = E1Y-Q2E
3000 PLEPT1(4) = E1Y+Q2E
3001C THETA = ACOS( (E1Y-Q2E)/(E1Y+Q2E) )
3002* radiated photon
3003 PGAMM(1) = -PLEPT1(1)
3004 PGAMM(2) = -PLEPT1(2)
3005 PGAMM(3) = PLEPT0(3)-PLEPT1(3)
3006 PGAMM(4) = PLEPT0(4)-PLEPT1(4)
3007* E_cm cut
3008 PTOTGN = SQRT( (PGAMM(1)+PNUCL(1))**2+(PGAMM(2)+PNUCL(2))**2
3009 & +(PGAMM(3)+PNUCL(3))**2 )
3010 ETOTGN = PGAMM(4)+PNUCL(4)
3011 ECMGN = (ETOTGN-PTOTGN)*(ETOTGN+PTOTGN)
3012 IF (ECMGN.LT.0.1D0) GOTO 101
3013 ECMGN = SQRT(ECMGN)
3014 IF ((ECMGN.LT.ECMIN).OR.(ECMGN.GT.ECMAX)) GOTO 101
3015
3016* Lorentz-transformation into nucleon-rest system
3017 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
3018 & PGAMM(1),PGAMM(2),PGAMM(3),PGAMM(4),
3019 & PGTOT,PPG(1),PPG(2),PPG(3),PPG(4))
3020 CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
3021 & PLEPT1(1),PLEPT1(2),PLEPT1(3),PLEPT1(4),
3022 & PLTOT,PPL1(1),PPL1(2),PPL1(3),PPL1(4))
3023* temporary checks..
3024 Q2TMP = ABS(PPG(4)**2-PGTOT**2)
3025 IF (ABS(Q2-Q2TMP).GT.0.01D0) WRITE(LOUT,1001) Q2,Q2TMP
3026 1001 FORMAT(1X,'LAEVT: inconsistent kinematics (Q2,Q2TMP) ',
3027 & 2F10.4)
3028 ECMTMP = SQRT((PPG(4)+AAM(1)-PGTOT)*(PPG(4)+AAM(1)+PGTOT))
3029 IF (ABS(ECMGN-ECMTMP).GT.TINY10) WRITE(LOUT,1002) ECMGN,ECMTMP
3030 1002 FORMAT(1X,'LAEVT: inconsistent kinematics (ECMGN,ECMTMP) ',
3031 & 2F10.2)
3032 YYTMP = PPG(4)/PPL0(4)
3033 IF (ABS(YY-YYTMP).GT.0.01D0) WRITE(LOUT,1005) YY,YYTMP
3034 1005 FORMAT(1X,'LAEVT: inconsistent kinematics (YY,YYTMP) ',
3035 & 2F10.4)
3036
3037* lepton tagger (Lab)
3038 THETA = ACOS( PPL1(3)/PLTOT )
3039 IF (PPL1(4).GT.ELMIN) THEN
3040 IF ((THETA.LT.THMIN).OR.(THETA.GT.THMAX)) GOTO 101
3041 ENDIF
3042* photon energy-cut (Lab)
3043 IF (PPG(4).LT.EGMIN) GOTO 101
3044 IF (PPG(4).GT.EGMAX) GOTO 101
3045* x_Bj cut
3046 XBJ = ABS(Q2/(1.876D0*PPG(4)))
3047 IF (XBJ.LT.XBJMIN) GOTO 101
3048
3049 NC0 = NC0+1
3050 CALL DT_FILHGR( Q2,ONE,IHFLQ0,NC0)
3051 CALL DT_FILHGR( YY,ONE,IHFLY0,NC0)
3052 CALL DT_FILHGR( XBJ,ONE,IHFLX0,NC0)
3053 CALL DT_FILHGR(PPG(4),ONE,IHFLU0,NC0)
3054 CALL DT_FILHGR( ECMGN,ONE,IHFLE0,NC0)
3055
3056* rotation angles against z-axis
3057 COD = PPG(3)/PGTOT
3058C SID = SQRT((ONE-COD)*(ONE+COD))
3059 PPT = SQRT(PPG(1)**2+PPG(2)**2)
3060 SID = PPT/PGTOT
3061 COF = ONE
3062 SIF = ZERO
3063 IF (PGTOT*SID.GT.TINY10) THEN
3064 COF = PPG(1)/(SID*PGTOT)
3065 SIF = PPG(2)/(SID*PGTOT)
3066 ANORF = SQRT(COF*COF+SIF*SIF)
3067 COF = COF/ANORF
3068 SIF = SIF/ANORF
3069 ENDIF
3070
3071 IF (IXSTBL.EQ.0) THEN
3072* change to photon projectile
3073 IJPROJ = 7
3074* set virtuality
3075 VIRT = Q2
3076* re-initialize LTs with new kinematics
3077* !!PGAMM ist set in cms (ECMGN) along z
3078 EPN = ZERO
3079 PPN = ZERO
3080 CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,ECMGN,0)
3081* force Lab-system
3082 IFRAME = 1
3083* get emulsion component if requested
3084 IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
3085* convolute with cross section
3086 CALL DT_SIGGAT(Q2LOW,EGNXX,STOTX,KKMAT)
3087 CALL DT_SIGGAT(Q2,ECMGN,STOT,KKMAT)
3088 IF (STOTX.LT.STOT) WRITE(LOUT,'(1X,A,/,6E12.3)')
3089 & 'LAEVT: warning STOTX<STOT ! ',Q2LOW,EGNMAX,STOTX,
3090 & Q2,ECMGN,STOT
3091 IF (DT_RNDM(Q2)*STOTX.GT.STOT) GOTO 100
3092 NC1 = NC1+1
3093 CALL DT_FILHGR( Q2,ONE,IHFLQ1,NC1)
3094 CALL DT_FILHGR( YY,ONE,IHFLY1,NC1)
3095 CALL DT_FILHGR( XBJ,ONE,IHFLX1,NC1)
3096 CALL DT_FILHGR(PPG(4),ONE,IHFLU1,NC1)
3097 CALL DT_FILHGR( ECMGN,ONE,IHFLE1,NC1)
3098* composite targets only
3099 KKMAT = -KKMAT
3100* sample this event
3101 CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IJPROJ,EPN,KKMAT,
3102 & IREJ)
3103* rotate momenta of final state particles back in photon-nucleon syst.
3104 DO 4 I=NPOINT(4),NHKK
3105 IF ((ABS(ISTHKK(I)).EQ.1).OR.(ISTHKK(I).EQ.1000).OR.
3106 & (ISTHKK(I).EQ.1001)) THEN
3107 PX = PHKK(1,I)
3108 PY = PHKK(2,I)
3109 PZ = PHKK(3,I)
3110 CALL DT_MYTRAN(1,PX,PY,PZ,COD,SID,COF,SIF,
3111 & PHKK(1,I),PHKK(2,I),PHKK(3,I))
3112 ENDIF
3113 4 CONTINUE
3114 ENDIF
3115
3116 CALL DT_FILHGR( Q2,ONE,IHFLQ2,NC1)
3117 CALL DT_FILHGR( YY,ONE,IHFLY2,NC1)
3118 CALL DT_FILHGR( XBJ,ONE,IHFLX2,NC1)
3119 CALL DT_FILHGR(PPG(4),ONE,IHFLU2,NC1)
3120 CALL DT_FILHGR( ECMGN,ONE,IHFLE2,NC1)
3121
3122* dump this event to histograms
3123 CALL PHO_PHIST(2000,DUM)
3124
3125 2 CONTINUE
3126
3127 WGY = ALPHEM/TWOPI*WGHMAX*DBLE(ITRY)/DBLE(ITRW)
3128 WGY = WGY*LOG(YMAX/YMIN)
3129 WEIGHT = WGY*SIGMAX*DBLE(NEVTS)/DBLE(ITRY)
3130
3131C HEADER = ' LAEVT: Q^2 distribution 0'
3132C CALL DT_OUTHGR(IHFLQ0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3133C HEADER = ' LAEVT: Q^2 distribution 1'
3134C CALL DT_OUTHGR(IHFLQ1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3135C HEADER = ' LAEVT: Q^2 distribution 2'
3136C CALL DT_OUTHGR(IHFLQ2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3137C HEADER = ' LAEVT: y distribution 0'
3138C CALL DT_OUTHGR(IHFLY0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3139C HEADER = ' LAEVT: y distribution 1'
3140C CALL DT_OUTHGR(IHFLY1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3141C HEADER = ' LAEVT: y distribution 2'
3142C CALL DT_OUTHGR(IHFLY2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3143C HEADER = ' LAEVT: x distribution 0'
3144C CALL DT_OUTHGR(IHFLX0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3145C HEADER = ' LAEVT: x distribution 1'
3146C CALL DT_OUTHGR(IHFLX1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3147C HEADER = ' LAEVT: x distribution 2'
3148C CALL DT_OUTHGR(IHFLX2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3149C HEADER = ' LAEVT: E_g distribution 0'
3150C CALL DT_OUTHGR(IHFLU0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3151C HEADER = ' LAEVT: E_g distribution 1'
3152C CALL DT_OUTHGR(IHFLU1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3153C HEADER = ' LAEVT: E_g distribution 2'
3154C CALL DT_OUTHGR(IHFLU2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3155C HEADER = ' LAEVT: E_c distribution 0'
3156C CALL DT_OUTHGR(IHFLE0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3157C HEADER = ' LAEVT: E_c distribution 1'
3158C CALL DT_OUTHGR(IHFLE1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3159C HEADER = ' LAEVT: E_c distribution 2'
3160C CALL DT_OUTHGR(IHFLE2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
3161
3162* print run-statistics and histograms to output-unit 6
3163 CALL PHO_PHIST(3000,DUM)
3164 IF (IXSTBL.EQ.0) CALL DT_STATIS(2)
3165
3166 RETURN
3167 END
3168
3169*$ CREATE DT_DTUINI.FOR
3170*COPY DT_DTUINI
3171*
3172*===dtuini=============================================================*
3173*
3174 SUBROUTINE DT_DTUINI(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
3175 & IDP,IEMU)
3176
3177 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3178 SAVE
3179
3180 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
3181* emulsion treatment
3182 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
3183 & NCOMPO,IEMUL
3184* Glauber formalism: flags and parameters for statistics
3185 LOGICAL LPROD
3186 CHARACTER*8 CGLB
3187 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
3188
3189 CALL DT_INIT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
3190 CALL DT_STATIS(1)
3191 CALL PHO_PHIST(1000,DUM)
3192 IF (NCOMPO.LE.0) THEN
3193 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
3194 ELSE
3195 DO 1 I=1,NCOMPO
3196 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
3197 1 CONTINUE
3198 ENDIF
3199 IF (IOGLB.NE.100) CALL DT_SIGEMU
3200 IEMU = IEMUL
3201
3202 RETURN
3203 END
3204
3205*$ CREATE DT_DTUOUT.FOR
3206*COPY DT_DTUOUT
3207*
3208*===dtuout=============================================================*
3209*
3210 SUBROUTINE DT_DTUOUT
3211
3212 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3213 SAVE
3214
3215 CALL PHO_PHIST(3000,DUM)
3216 CALL DT_STATIS(2)
3217
3218 RETURN
3219 END
3220
3221*$ CREATE DT_BEAMPR.FOR
3222*COPY DT_BEAMPR
3223*
3224*===beampr=============================================================*
3225*
3226 SUBROUTINE DT_BEAMPR(WHAT,PLAB,MODE)
3227
3228************************************************************************
3229* Initialization of event generation *
3230* This version dated 7.4.98 is written by S. Roesler. *
3231************************************************************************
3232
3233 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3234 SAVE
3235
3236 PARAMETER ( LINP = 10 ,
3237 & LOUT = 6 ,
3238 & LDAT = 9 )
3239 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
3240 PARAMETER (TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
3241
3242 LOGICAL LBEAM
3243
3244* event history
3245 PARAMETER (NMXHKK=200000)
3246 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3247 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3248 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3249* extended event history
3250 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3251 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3252 & IHIST(2,NMXHKK)
3253* properties of interacting particles
3254 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
3255* particle properties (BAMJET index convention)
3256 CHARACTER*8 ANAME
3257 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
3258 & IICH(210),IIBAR(210),K1(210),K2(210)
3259* beam momenta
3260 COMMON /DTBEAM/ P1(4),P2(4)
3261
3262C DIMENSION WHAT(6),P1(4),P2(4),P1CMS(4),P2CMS(4)
3263 DIMENSION WHAT(6),P1CMS(4),P2CMS(4)
3264
3265 DATA LBEAM /.FALSE./
3266
3267 GOTO (1,2) MODE
3268
3269 1 CONTINUE
3270
3271 E1 = WHAT(1)
3272 IF (E1.LT.ZERO) E1 = DBLE(IPZ)/DBLE(IP)*ABS(WHAT(1))
3273 E2 = WHAT(2)
3274 IF (E2.LT.ZERO) E2 = DBLE(ITZ)/DBLE(IT)*ABS(WHAT(2))
3275 PP1 = SQRT( (E1+AAM(IJPROJ))*(E1-AAM(IJPROJ)) )
3276 PP2 = SQRT( (E2+AAM(IJTARG))*(E2-AAM(IJTARG)) )
3277 TH = 1.D-6*WHAT(3)/2.D0
3278 PH = WHAT(4)*BOG
3279 P1(1) = PP1*SIN(TH)*COS(PH)
3280 P1(2) = PP1*SIN(TH)*SIN(PH)
3281 P1(3) = PP1*COS(TH)
3282 P1(4) = E1
3283 P2(1) = PP2*SIN(TH)*COS(PH)
3284 P2(2) = PP2*SIN(TH)*SIN(PH)
3285 P2(3) = -PP2*COS(TH)
3286 P2(4) = E2
3287 ECM = SQRT( (P1(4)+P2(4))**2-(P1(1)+P2(1))**2-(P1(2)+P2(2))**2
3288 & -(P1(3)+P2(3))**2 )
3289 ELAB = (ECM**2-AAM(IJPROJ)**2-AAM(IJTARG)**2)/(2.0D0*AAM(IJTARG))
3290 PLAB = SQRT( (ELAB+AAM(IJPROJ))*(ELAB-AAM(IJPROJ)) )
3291 BGX = (P1(1)+P2(1))/ECM
3292 BGY = (P1(2)+P2(2))/ECM
3293 BGZ = (P1(3)+P2(3))/ECM
3294 BGE = (P1(4)+P2(4))/ECM
3295 CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P1(1),P1(2),P1(3),P1(4),
3296 & P1TOT,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4))
3297 CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P2(1),P2(2),P2(3),P2(4),
3298 & P2TOT,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4))
3299 COD = P1CMS(3)/P1TOT
3300C SID = SQRT((ONE-COD)*(ONE+COD))
3301 PPT = SQRT(P1CMS(1)**2+P1CMS(2)**2)
3302 SID = PPT/P1TOT
3303 COF = ONE
3304 SIF = ZERO
3305 IF (P1TOT*SID.GT.TINY10) THEN
3306 COF = P1CMS(1)/(SID*P1TOT)
3307 SIF = P1CMS(2)/(SID*P1TOT)
3308 ANORF = SQRT(COF*COF+SIF*SIF)
3309 COF = COF/ANORF
3310 SIF = SIF/ANORF
3311 ENDIF
3312**check
3313C WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
3314C WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
3315C WRITE(LOUT,'(5E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),P1TOT
3316C WRITE(LOUT,'(5E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),P2TOT
3317C PAX = ZERO
3318C PAY = ZERO
3319C PAZ = P1TOT
3320C PAE = SQRT(AAM(IJPROJ)**2+PAZ**2)
3321C PBX = ZERO
3322C PBY = ZERO
3323C PBZ = -P2TOT
3324C PBE = SQRT(AAM(IJTARG)**2+PBZ**2)
3325C WRITE(LOUT,'(4E15.4)') PAX,PAY,PAZ,PAE
3326C WRITE(LOUT,'(4E15.4)') PBX,PBY,PBZ,PBE
3327C CALL DT_MYTRAN(1,PAX,PAY,PAZ,COD,SID,COF,SIF,
3328C & P1CMS(1),P1CMS(2),P1CMS(3))
3329C CALL DT_MYTRAN(1,PBX,PBY,PBZ,COD,SID,COF,SIF,
3330C & P2CMS(1),P2CMS(2),P2CMS(3))
3331C WRITE(LOUT,'(4E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4)
3332C WRITE(LOUT,'(4E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4)
3333C CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),
3334C & P1TOT,P1(1),P1(2),P1(3),P1(4))
3335C CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),
3336C & P2TOT,P2(1),P2(2),P2(3),P2(4))
3337C WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
3338C WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
3339C STOP
3340**
3341
3342 LBEAM = .TRUE.
3343
3344 RETURN
3345
3346 2 CONTINUE
3347
3348 IF (LBEAM) THEN
3349 IF ( (NPOINT(4).EQ.0).OR.(NHKK.LT.NPOINT(4)) ) RETURN
3350 DO 20 I=NPOINT(4),NHKK
3351 IF ((ABS(ISTHKK(I)).EQ.1).OR.(ISTHKK(I).EQ.1000).OR.
3352 & (ISTHKK(I).EQ.1001)) THEN
3353 CALL DT_MYTRAN(1,PHKK(1,I),PHKK(2,I),PHKK(3,I),
3354 & COD,SID,COF,SIF,PXCMS,PYCMS,PZCMS)
3355 PECMS = PHKK(4,I)
3356 CALL DT_DALTRA(BGE,BGX,BGY,BGZ,PXCMS,PYCMS,PZCMS,PECMS,
3357 & PTOT,PHKK(1,I),PHKK(2,I),PHKK(3,I),PHKK(4,I))
3358 ENDIF
3359 20 CONTINUE
3360 ELSE
3361 MODE = -1
3362 ENDIF
3363
3364 RETURN
3365 END
3366
3367*$ CREATE DT_REJUCO.FOR
3368*COPY DT_REJUCO
3369*
3370*===rejuco=============================================================*
3371*
3372 SUBROUTINE DT_REJUCO(MODE,IREJ)
3373
3374************************************************************************
3375* REJection of Unphysical COnfigurations *
3376* MODE = 1 rejection of particles with unphysically large energy *
3377* *
3378* This version dated 27.12.2006 is written by S. Roesler. *
3379************************************************************************
3380
3381 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3382 SAVE
3383
3384 PARAMETER ( LINP = 10 ,
3385 & LOUT = 6 ,
3386 & LDAT = 9 )
3387 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
3388 PARAMETER (TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
3389
3390* maximum x_cms of final state particle
3391 PARAMETER (XCMSMX = 1.4D0)
3392
3393* event history
3394 PARAMETER (NMXHKK=200000)
3395 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3396 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3397 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3398* extended event history
3399 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3400 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3401 & IHIST(2,NMXHKK)
3402* Lorentz-parameters of the current interaction
3403 COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
3404 & UMO,PPCM,EPROJ,PPROJ
3405
3406 IREJ = 0
3407
3408 IF (MODE.EQ.1) THEN
3409 IF ( (NPOINT(4).EQ.0).OR.(NHKK.LT.NPOINT(4)) ) RETURN
3410 ECMHLF = UMO/2.0D0
3411 DO 10 I=NPOINT(4),NHKK
3412 IF ((ABS(ISTHKK(I)).EQ.1).AND.(IDHKK(I).NE.80000)) THEN
3413 XCMS = ABS(PHKK(4,I))/ECMHLF
3414 IF (XCMS.GT.XCMSMX) GOTO 9999
3415 ENDIF
3416 10 CONTINUE
3417 ENDIF
3418
3419 RETURN
3420 9999 CONTINUE
3421 IREJ = 1
3422 RETURN
3423 END
3424
3425*$ CREATE DT_EVENTB.FOR
3426*COPY DT_EVENTB
3427*
3428*===eventb=============================================================*
3429*
3430 SUBROUTINE DT_EVENTB(NCSY,IREJ)
3431
3432************************************************************************
3433* Treatment of nucleon-nucleon interactions with full two-component *
3434* Dual Parton Model. *
3435* NCSY number of nucleon-nucleon interactions *
3436* IREJ rejection flag *
3437* This version dated 14.01.2000 is written by S. Roesler *
3438************************************************************************
3439
3440 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3441 SAVE
3442 PARAMETER ( LINP = 10 ,
3443 & LOUT = 6 ,
3444 & LDAT = 9 )
3445 PARAMETER (TINY10=1.0D-10,ZERO=0.0D0,ONE=1.0D0)
3446
3447* event history
3448 PARAMETER (NMXHKK=200000)
3449 COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
3450 & JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
3451 & PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
3452* extended event history
3453 COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
3454 & IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
3455 & IHIST(2,NMXHKK)
3456*! uncomment this line for internal phojet-fragmentation
3457C #include "dtu_dtevtp.inc"
3458* particle properties (BAMJET index convention)
3459 CHARACTER*8 ANAME
3460 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
3461 & IICH(210),IIBAR(210),K1(210),K2(210)
3462* flags for input different options
3463 LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
3464 COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
3465 & LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
3466* rejection counter
3467 COMMON /DTREJC/ IRPT,IRHHA,IRRES(2),LOMRES,LOBRES,
3468 & IRCHKI(2),IRFRAG,IRCRON(3),IREVT,
3469 & IREXCI(3),IRDIFF(2),IRINC
3470* properties of interacting particles
3471 COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
3472* properties of photon/lepton projectiles
3473 COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
3474* various options for treatment of partons (DTUNUC 1.x)
3475* (chain recombination, Cronin,..)
3476 LOGICAL LCO2CR,LINTPT
3477 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
3478 & LCO2CR,LINTPT
3479* statistics
3480 COMMON /DTSTA1/ ICREQU,ICSAMP,ICCPRO,ICDPR,ICDTA,
3481 & ICRJSS,ICVV2S,ICCHAI(2,9),ICRES(9),ICDIFF(5),
3482 & ICEVTG(8,0:30)
3483* DTUNUC-PHOJET interface, Lorentz-param. of n-n subsystem
3484 COMMON /DTLTSU/ BGX,BGY,BGZ,GAM
3485* Glauber formalism: collision properties
3486 COMMON /DTGLCP/ RPROJ,RTARG,BIMPAC,
3487 & NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC
3488* flags for diffractive interactions (DTUNUC 1.x)
3489 COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
3490* statistics: double-Pomeron exchange
3491 COMMON /DTFLG2/ INTFLG,IPOPO
3492* flags for particle decays
3493 COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
3494 & IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
3495 & NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
3496* nucleon-nucleon event-generator
3497 CHARACTER*8 CMODEL
3498 LOGICAL LPHOIN
3499 COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
3500C nucleon-nucleus / nucleus-nucleus interface to DPMJET
3501 INTEGER IDEQP,IDEQB,IHFLD,IHFLS
3502 DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
3503 COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
3504 & IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
3505C model switches and parameters
3506 CHARACTER*8 MDLNA
3507 INTEGER ISWMDL,IPAMDL
3508 DOUBLE PRECISION PARMDL
3509 COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
3510C initial state parton radiation (internal part)
3511 INTEGER MXISR3,MXISR4
3512 PARAMETER ( MXISR3 = 50, MXISR4 = 100 )
3513 INTEGER IFL1,IFL2,IBRA,IFANO,ISH,NACC
3514 DOUBLE PRECISION Q2SH,PT2SH,XPSH,ZPSH,THSH,SHAT
3515 COMMON /POINT6/ Q2SH(2,MXISR3),PT2SH(2,MXISR3),XPSH(2,MXISR3),
3516 & ZPSH(2,MXISR3),THSH(2,MXISR3),SHAT(MXISR3),
3517 & IFL1(2,MXISR3),IFL2(2,MXISR3),
3518 & IBRA(2,MXISR4),IFANO(2),ISH(2),NACC
3519C event debugging information
3520 INTEGER NMAXD
3521 PARAMETER (NMAXD=100)
3522 INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
3523 & KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3524 COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
3525 & KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
3526C general process information
3527 INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
3528 COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
3529
3530 DIMENSION PP(4),PT(4),PTOT(4),PP1(4),PP2(4),PT1(4),PT2(4),
3531 & PPNN(4),PTNN(4),PTOTNN(4),PPSUB(4),PTSUB(4),
3532 & PPTCMS(4),PTTCMS(4),PPTMP(4),PTTMP(4),
3533 & KPRON(15),ISINGL(2000)
3534
3535* initial values for max. number of phojet scatterings and dtunuc chains
3536* to be fragmented with one pyexec call
3537 DATA MXPHFR,MXDTFR /10,100/
3538
3539 IREJ = 0
3540* pointer to first parton of the first chain in dtevt common
3541 NPOINT(3) = NHKK+1
3542* special flag for double-Pomeron statistics
3543 IPOPO = 1
3544* counter for low-mass (DTUNUC) interactions
3545 NDTUSC = 0
3546* counter for interactions treated by PHOJET
3547 NPHOSC = 0
3548
3549* scan interactions for single nucleon-nucleon interactions
3550* (this has to be checked here because Cronin modifies parton momenta)
3551 NC = NPOINT(2)
3552 IF (NCSY.GT.2000) STOP ' DT_EVENTB: NCSY > 2000 ! '
3553 DO 8 I=1,NCSY
3554 ISINGL(I) = 0
3555 MOP = JMOHKK(1,NC)
3556 MOT = JMOHKK(1,NC+1)
3557 DIFF1 = ABS(PHKK(4,MOP)-PHKK(4, NC)-PHKK(4,NC+2))
3558 DIFF2 = ABS(PHKK(4,MOT)-PHKK(4,NC+1)-PHKK(4,NC+3))
3559 IF ((DIFF1.LT.TINY10).AND.(DIFF2.LT.TINY10)) ISINGL(I) = 1
3560 NC = NC+4
3561 8 CONTINUE
3562
3563* multiple scattering of chain ends
3564 IF ((IP.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(1)
3565 IF ((IT.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(2)
3566
3567* switch to PHOJET-settings for JETSET parameter
3568 CALL DT_INITJS(1)
3569
3570* loop over nucleon-nucleon interaction
3571 NC = NPOINT(2)
3572 DO 2 I=1,NCSY
3573*
3574* pick up one nucleon-nucleon interaction from DTEVT1
3575* ppnn / ptnn - momenta of the interacting nucleons (cms)
3576* ptotnn - total momentum of the interacting nucleons (cms)
3577* pp1,2 / pt1,2 - momenta of the four partons
3578* pp / pt - total momenta of the proj / targ partons
3579* ptot - total momentum of the four partons
3580 MOP = JMOHKK(1,NC)
3581 MOT = JMOHKK(1,NC+1)
3582 DO 3 K=1,4
3583 PPNN(K) = PHKK(K,MOP)
3584 PTNN(K) = PHKK(K,MOT)
3585 PTOTNN(K) = PPNN(K)+PTNN(K)
3586 PP1(K) = PHKK(K,NC)
3587 PT1(K) = PHKK(K,NC+1)
3588 PP2(K) = PHKK(K,NC+2)
3589 PT2(K) = PHKK(K,NC+3)
3590 PP(K) = PP1(K)+PP2(K)
3591 PT(K) = PT1(K)+PT2(K)
3592 PTOT(K) = PP(K)+PT(K)
3593 3 CONTINUE
3594*
3595*-----------------------------------------------------------------------
3596* this is a complete nucleon-nucleon interaction
3597*
3598 IF (ISINGL(I).EQ.1) THEN
3599*
3600* initialize PHOJET-variables for remnant/valence-partons
3601 IHFLD(1,1) = 0
3602 IHFLD(1,2) = 0
3603 IHFLD(2,1) = 0
3604 IHFLD(2,2) = 0
3605 IHFLS(1) = 1
3606 IHFLS(2) = 1
3607* save current settings of PHOJET process and min. bias flags
3608 DO 9 K=1,11
3609 KPRON(K) = IPRON(K,1)
3610 9 CONTINUE
3611 ISWSAV = ISWMDL(2)
3612*
3613* check if forced sampling of diffractive interaction requested
3614 IF (ISINGD.LT.-1) THEN
3615 DO 90 K=1,11
3616 IPRON(K,1) = 0
3617 90 CONTINUE
3618 IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-3)) IPRON(5,1) = 1
3619 IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-4)) IPRON(6,1) = 1
3620 IF (ISINGD.EQ.-5) IPRON(4,1) = 1
3621 ENDIF
3622*
3623* for photons: a direct/anomalous interaction is not sampled
3624* in PHOJET but already in Glauber-formalism. Here we check if such
3625* an interaction is requested
3626 IF (IJPROJ.EQ.7) THEN
3627* first switch off direct interactions
3628 IPRON(8,1) = 0
3629* this is a direct interactions
3630 IF (IDIREC.EQ.1) THEN
3631 DO 12 K=1,11
3632 IPRON(K,1) = 0
3633 12 CONTINUE
3634 IPRON(8,1) = 1
3635* this is an anomalous interactions
3636* (iswmdl(2) = 0 only hard int. generated ( = 1 min. bias) )
3637 ELSEIF (IDIREC.EQ.2) THEN
3638 ISWMDL(2) = 0
3639 ENDIF
3640 ELSE
3641 IF (IDIREC.NE.0) STOP ' DT_EVENTB: IDIREC > 0 ! '
3642 ENDIF
3643*
3644* make sure that total momenta of partons, pp and pt, are on mass
3645* shell (Cronin may have srewed this up..)
3646 CALL DT_MASHEL(PP,PT,PHKK(5,MOP),PHKK(5,MOT),PPNN,PTNN,IR1)
3647 IF (IR1.NE.0) THEN
3648 IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A)')
3649 & 'EVENTB: mass shell correction rejected'
3650 GOTO 9999
3651 ENDIF
3652*
3653* initialize the incoming particles in PHOJET
3654 IF ((IP.EQ.1).AND.(IJPROJ.EQ.7)) THEN
3655 CALL PHO_SETPAR(1,22,0,VIRT)
3656 ELSE
3657 CALL PHO_SETPAR(1,IDHKK(MOP),0,ZERO)
3658 ENDIF
3659 CALL PHO_SETPAR(2,IDHKK(MOT),0,ZERO)
3660*
3661* initialize rejection loop counter for anomalous processes
3662 IRJANO = 0
3663 800 CONTINUE
3664 IRJANO = IRJANO+1
3665*
3666* temporary fix for ifano problem
3667 IFANO(1) = 0
3668 IFANO(2) = 0
3669*
3670* generate complete hadron/nucleon/photon-nucleon event with PHOJET
3671 CALL PHO_EVENT(2,PPNN,PTNN,DUM,IREJ1)
3672*
3673* for photons: special consistency check for anomalous interactions
3674 IF (IJPROJ.EQ.7) THEN
3675 IF (IRJANO.LT.30) THEN
3676 IF (IFANO(1).NE.0) THEN
3677* here, an anomalous interaction was generated. Check if it
3678* was also requested. Otherwise reject this event.
3679 IF (IDIREC.EQ.0) GOTO 800
3680 ELSE
3681* here, an anomalous interaction was not generated. Check if it
3682* was requested in which case we need to reject this event.
3683 IF (IDIREC.EQ.2) GOTO 800
3684 ENDIF
3685 ELSE
3686 WRITE(LOUT,*) ' DT_EVENTB: Warning! IRJANO > 30 ',
3687 & IRJANO,IDIREC,NEVHKK
3688 ENDIF
3689 ENDIF
3690*
3691* copy back original settings of PHOJET process and min. bias flags
3692 DO 10 K=1,11
3693 IPRON(K,1) = KPRON(K)
3694 10 CONTINUE
3695 ISWMDL(2) = ISWSAV
3696*
3697* check if PHOJET has rejected this event
3698 IF (IREJ1.NE.0) THEN
3699C IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A,I4)')
3700 WRITE(LOUT,'(1X,A,I4)')
3701 & 'EVENTB: chain system rejected',IDIREC
3702 CALL PHO_PREVNT(0)
3703 GOTO 9999
3704 ENDIF
3705*
3706* copy partons and strings from PHOJET common back into DTEVT for
3707* external fragmentation
3708 MO1 = NC
3709 MO2 = NC+3
3710*! uncomment this line for internal phojet-fragmentation
3711C CALL DT_GETFSP(MO1,MO2,PPNN,PTNN,-1)
3712 NPHOSC = NPHOSC+1
3713 CALL DT_GETPJE(MO1,MO2,PPNN,PTNN,-1,NPHOSC,IREJ1)
3714 IF (IREJ1.NE.0) THEN
3715 IF (IOULEV(1).GT.0)
3716 & WRITE(LOUT,'(1X,A,I4)') 'EVENTB: chain system rejected 1'
3717 GOTO 9999
3718 ENDIF
3719*
3720* update statistics counter
3721 ICEVTG(IDCH(NC),29) = ICEVTG(IDCH(NC),29)+1
3722*
3723*-----------------------------------------------------------------------
3724* this interaction involves "remnants"
3725*
3726 ELSE
3727*
3728* total mass of this system
3729 PPTOT = SQRT(PTOT(1)**2+PTOT(2)**2+PTOT(3)**2)
3730 AMTOT2 = (PTOT(4)-PPTOT)*(PTOT(4)+PPTOT)
3731 IF (AMTOT2.LT.ZERO) THEN
3732 AMTOT = ZERO
3733 ELSE
3734 AMTOT = SQRT(AMTOT2)
3735 ENDIF
3736*
3737* systems with masses larger than elojet are treated with PHOJET
3738 IF (AMTOT.GT.ELOJET) THEN
3739*
3740* initialize PHOJET-variables for remnant/valence-partons
3741* projectile parton flavors and valence flag
3742 IHFLD(1,1) = IDHKK(NC)
3743 IHFLD(1,2) = IDHKK(NC+2)
3744 IHFLS(1) = 0
3745 IF ((IDCH(NC).EQ.6).OR.(IDCH(NC).EQ.7)
3746 & .OR.(IDCH(NC).EQ.8)) IHFLS(1) = 1
3747* target parton flavors and valence flag
3748 IHFLD(2,1) = IDHKK(NC+1)
3749 IHFLD(2,2) = IDHKK(NC+3)
3750 IHFLS(2) = 0
3751 IF ((IDCH(NC).EQ.4).OR.(IDCH(NC).EQ.5)
3752 & .OR.(IDCH(NC).EQ.8)) IHFLS(2) = 1
3753* flag signalizing PHOJET how to treat the remnant:
3754* iremn = -1 sea-quark remnant: PHOJET takes flavors from ihfld
3755* iremn > -1 valence remnant: PHOJET assumes flavors according
3756* to mother particle
3757 IREMN1 = IHFLS(1)-1
3758 IREMN2 = IHFLS(2)-1
3759*
3760* initialize the incoming particles in PHOJET
3761 IF ((