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