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