2 * +-------------------------------------------------------------+
8 * | S. Roesler+), R. Engel#), J. Ranft*) |
11 * | CH-1211 Geneva 23, Switzerland |
12 * | Email: Stefan.Roesler@cern.ch |
14 * | #) University of Delaware, BRI |
15 * | Newark, DE 19716, USA |
17 * | *) University of Siegen, Dept. of Physics |
18 * | D-57068 Siegen, Germany |
21 * | http://home.cern.ch/sroesler/dpmjet3.html |
24 * | Monte Carlo models used for event generation: |
25 * | PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1 |
27 * +-------------------------------------------------------------+
30 *===init===============================================================*
33 SUBROUTINE DT_INIT(NCASES,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
36 ************************************************************************
37 * Initialization of event generation *
38 * This version dated 7.4.98 is written by S. Roesler. *
39 ************************************************************************
41 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
44 PARAMETER ( LINP = 5 ,
48 PARAMETER (ZERO=0.0D0,ONE=1.0D0)
50 * particle properties (BAMJET index convention)
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
56 COMMON /DTPAIN/ BTYPE(30)
63 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
66 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
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),
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
90 COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
91 * various options for treatment of partons (DTUNUC 1.x)
92 * (chain recombination, Cronin,..)
94 COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
96 * threshold values for x-sampling (DTUNUC 1.x)
97 COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
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
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
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
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
146 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
148 COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
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
157 COMMON /LEPTOI/ RPPN,LEPIN,INTER
158 * steering flags for qel neutrino scattering modules
159 COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
161 COMMON /DTEVNO/ NEVENT,ICASCA
166 DIMENSION XDUMB(40),IPRANG(5)
168 PARAMETER (MXCARD=58)
169 CHARACTER*78 CLINE,CTITLE
171 CHARACTER*8 BLANK,SDUM
172 CHARACTER*10 CODE,CODEWD
174 LOGICAL LSTART,LEINP,LXSTAB
175 DIMENSION WHAT(6),CODE(MXCARD)
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',
194 DATA LSTART,LXSTAB,IFIRST /.TRUE.,.FALSE.,1/
197 * --- Added by Chiara
199 CHARACTER*100 ALIROOT
205 *---------------------------------------------------------------------
206 * at the first call of INIT: initialize event generation
210 * initialization and test of the random number generator
211 IF (ITRSPT.NE.1) THEN
213 CALL FL48UT (ISRM48,ISEED1,ISEED2)
214 CALL FL48IN (54217137,ISEED1,ISEED2)
217 * initialization of BAMJET, DECAY and HADRIN
222 * set default values for input variables
223 CALL DT_DEFAUL(EPN,PPN)
226 * flag for collision energy input
231 *---------------------------------------------------------------------
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
250 * read control card from input-unit LINP
251 C READ(LINP,'(A78)',END=9999) CLINE
252 * ### Read control card from specified file
253 * ### Changed by Chiara (original version LINP=5)
255 * + FILE='/home/oppedisa/AliRoot/new/DPMJET/inp/PbPbLHC.inp',
258 CALL GETENVF('ALICE_ROOT',ALIROOT)
259 LNROOT = LNBLNK(ALIROOT)
261 FILNAM=ALIROOT(1:LNROOT)//'/DPMJET/inp/ppLHC.inp'
262 OPEN(UNIT=7,FILE=FILNAM,STATUS='OLD')
263 OPEN(UNIT=14,FILE="nuclear.bin",STATUS='OLD')
264 * OPEN(UNIT=6,FILE="dpm.out",STATUS='UNKNOWN')
266 READ(7,'(A78)',END=9999) CLINE
268 IF (CLINE(1:1).EQ.'*') THEN
270 C WRITE(LOUT,'(A78)') CLINE
273 C READ(CLINE,1000,END=9999) CODEWD,(WHAT(I),I=1,6),SDUM
274 C1000 FORMAT(A10,6E10.0,A8)
278 READ(CLINE,1006,END=9999) CODEWD,CWHAT,SDUM
279 1006 FORMAT(A10,A60,A8)
280 READ(CWHAT,*,END=1007) (WHAT(I),I=1,6)
282 WRITE(LOUT,1001) CODEWD,(WHAT(I),I=1,6),SDUM
283 1001 FORMAT(A10,6G10.3,A8)
287 * check for valid control card and get card index
290 IF (CODEWD.EQ.CODE(I)) ICW = I
293 WRITE(LOUT,1002) CODEWD
294 1002 FORMAT(/,1X,'---> ',A10,': invalid control-card !',/)
298 *------------------------------------------------------------
299 * TITLE , PROJPAR , TARPAR , ENERGY , MOMENTUM,
300 & 100 , 110 , 120 , 130 , 140 ,
302 *------------------------------------------------------------
303 * CMENERGY, EMULSION, FERMI , TAUFOR , PAULI ,
304 & 150 , 160 , 170 , 180 , 190 ,
306 *------------------------------------------------------------
307 * COULOMB , HADRIN , EVAP , EMCCHECK, MODEL ,
308 & 200 , 210 , 220 , 230 , 240 ,
310 *------------------------------------------------------------
311 * PHOINPUT, GLAUBERI, FLUCTUAT, CENTRAL , RECOMBIN,
312 & 250 , 260 , 270 , 280 , 290 ,
314 *------------------------------------------------------------
315 * COMBIJET, XCUTS , INTPT , CRONINPT, SEADISTR,
316 & 300 , 310 , 320 , 330 , 340 ,
318 *------------------------------------------------------------
319 * SEASU3 , DIQUARKS, RESONANC, DIFFRACT, SINGLECH,
320 & 350 , 360 , 370 , 380 , 390 ,
322 *------------------------------------------------------------
323 * NOFRAGME, HADRONIZE, POPCORN , PARDECAY, BEAM ,
324 & 400 , 410 , 420 , 430 , 440 ,
326 *------------------------------------------------------------
327 * LUND-MSTU, LUND-MSTJ, LUND-MDCY, LUND-PARJ, LUND-PARU,
328 & 450 , 451 , 452 , 460 , 470 ,
330 *------------------------------------------------------------
331 * OUTLEVEL, FRAME , L-TAG , L-ETAG , ECMS-CUT,
332 & 480 , 490 , 500 , 510 , 520 ,
334 *------------------------------------------------------------
335 * VDM-PAR1, HISTOGRAM, XS-TABLE , GLAUB-PAR, GLAUB-INI,
336 & 530 , 540 , 550 , 560 , 565 ,
338 *------------------------------------------------------------
339 * , , VDM-PAR2, XS-QELPRO, RNDMINIT ,
342 *------------------------------------------------------------
343 * LEPTO-CUT, LEPTO-LST,LEPTO-PARL, START , STOP )
344 & 600 , 610 , 620 , 630 , 640 ) , ICW
346 *------------------------------------------------------------
350 *********************************************************************
352 * control card: codewd = TITLE *
354 * what (1..6), sdum no meaning *
356 * Note: The control-card following this must consist of *
357 * a string of characters usually giving the title of *
360 *********************************************************************
363 C READ(LINP,'(A78)') CTITLE
364 * ### Read control card from specified file
365 * ### Changed by Chiara (original version LINP=5)
366 READ(7,'(A78)') CTITLE
368 WRITE(LOUT,'(//,5X,A78,//)') CTITLE
371 *********************************************************************
373 * control card: codewd = PROJPAR *
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 *
380 * Note: If sdum is defined what (1..2) have no meaning. *
382 *********************************************************************
385 IF (SDUM.EQ.BLANK) THEN
393 IF (SDUM.EQ.BTYPE(II)) THEN
398 ELSEIF (II.EQ.27) THEN
400 ELSEIF (II.EQ.28) THEN
402 ELSEIF (II.EQ.29) THEN
407 IBPROJ = IIBAR(IJPROJ)
409 IF ((IJPROJ.EQ.7).AND.(WHAT(1).GT.ZERO)) VIRT = WHAT(1)
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)
416 IF (IJPROJ.EQ.0) THEN
418 1110 FORMAT(/,1X,'invalid PROJPAR card !',/)
424 *********************************************************************
426 * control card: codewd = TARPAR *
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 *
433 * Note: If sdum is defined what (1..2) have no meaning. *
435 *********************************************************************
438 IF (SDUM.EQ.BLANK) THEN
446 IF (SDUM.EQ.BTYPE(II)) THEN
450 IBTARG = IIBAR(IJTARG)
453 IF (IJTARG.EQ.0) THEN
455 1120 FORMAT(/,1X,'invalid TARPAR card !',/)
461 *********************************************************************
463 * control card: codewd = ENERGY *
465 * what (1) = energy (GeV) of projectile in Lab. *
466 * if what(1) < 0: |what(1)| = kinetic energy *
468 * if |what(2)| > 0: min. energy for variable *
470 * what (2) = max. energy for variable energy runs *
471 * if what(2) < 0: |what(2)| = kinetic energy *
473 *********************************************************************
479 IF ((ABS(WHAT(2)).GT.ZERO).AND.
480 & (ABS(WHAT(2)).GT.ABS(WHAT(1)))) THEN
488 *********************************************************************
490 * control card: codewd = MOMENTUM *
492 * what (1) = momentum (GeV/c) of projectile in Lab. *
493 * default: 200 GeV/c *
494 * what (2..6), sdum no meaning *
496 *********************************************************************
505 *********************************************************************
507 * control card: codewd = CMENERGY *
509 * what (1) = energy in nucleon-nucleon cms. *
511 * what (2..6), sdum no meaning *
513 *********************************************************************
522 *********************************************************************
524 * control card: codewd = EMULSION *
526 * definition of nuclear emulsions *
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 *
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. *
543 *********************************************************************
546 IF ((WHAT(1).GT.ZERO).AND.(WHAT(2).GT.ZERO)
547 & .AND.(ABS(WHAT(3)).GT.ZERO)) THEN
549 IF (NCOMPO.GT.NCOMPX) THEN
553 IEMUMA(NCOMPO) = INT(WHAT(1))
554 IEMUCH(NCOMPO) = INT(WHAT(2))
555 EMUFRA(NCOMPO) = WHAT(3)
557 C CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
559 IF ((WHAT(4).GT.ZERO).AND.(WHAT(5).GT.ZERO)
560 & .AND.(ABS(WHAT(6)).GT.ZERO)) THEN
562 IF (NCOMPO.GT.NCOMPX) THEN
566 IEMUMA(NCOMPO) = INT(WHAT(4))
567 IEMUCH(NCOMPO) = INT(WHAT(5))
568 EMUFRA(NCOMPO) = WHAT(6)
569 C CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
571 1600 FORMAT(1X,'too many emulsion components - program stopped')
574 *********************************************************************
576 * control card: codewd = FERMI *
578 * what (1) = -1 Fermi-motion of nucleons not treated *
580 * what (2) = scale factor for Fermi-momentum *
582 * what (3..6), sdum no meaning *
584 *********************************************************************
587 IF (WHAT(1).EQ.-1.0D0) THEN
593 IF (XMOD.GE.ZERO) FERMOD = XMOD
596 *********************************************************************
598 * control card: codewd = TAUFOR *
600 * formation time supressed intranuclear cascade *
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 *
606 * what (2) number of generations followed *
608 * what (3) = 1. p_t-dependent formation zone *
609 * = 2. constant formation zone *
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 *
618 * what (5..6), sdum no meaning *
620 *********************************************************************
624 KTAUGE = INT(WHAT(2))
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))
632 *********************************************************************
634 * control card: codewd = PAULI *
636 * what (1) = -1 Pauli's principle for secondary *
637 * interactions not treated *
639 * what (2..6), sdum no meaning *
641 *********************************************************************
644 IF (WHAT(1).EQ.-1.0D0) THEN
651 *********************************************************************
653 * control card: codewd = COULOMB *
655 * what (1) = -1. Coulomb-energy treatment switched off *
657 * what (2..6), sdum no meaning *
659 *********************************************************************
663 IF (WHAT(1).EQ.-1.0D0) THEN
670 *********************************************************************
672 * control card: codewd = HADRIN *
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 *
681 * what (2) upper threshold in total energy (GeV) below *
682 * which interactions are sampled by HADRIN *
684 * what (3..6), sdum no meaning *
686 *********************************************************************
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)
694 *********************************************************************
696 * control card: codewd = EVAP *
698 * evaporation module *
700 * what (1) =< -1 ==> evaporation is switched off *
701 * >= 1 ==> evaporation is performed *
703 * what (1) = i1 + i2*10 + i3*100 + i4*10000 *
704 * (i1, i2, i3, i4 >= 0 ) *
706 * i1 is the flag for selecting the T=0 level density option used *
707 * = 1: standard EVAP level densities with Cook pairing *
709 * = 2: Z,N-dependent Gilbert & Cameron level densities *
711 * = 3: Julich A-dependent level densities *
712 * = 4: Z,N-dependent Brancazio & Cameron level densities *
714 * i2 >= 1: high energy fission activated *
715 * (default high energy fission activated) *
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 *
729 * = 6: Second Ignyatuk (1975, 2nd) energy dependence *
730 * for level densities with fit 2 Iljinov & Mebel set of *
732 * = 7: Second Ignyatuk (1975, 2nd) energy dependence *
733 * for level densities with fit 3 Iljinov & Mebel set of *
735 * = 8: Second Ignyatuk (1975, 2nd) energy dependence *
736 * for level densities with fit 4 Iljinov & Mebel set of *
739 * i4 >= 1: Original Gilbert and Cameron pairing energies used *
740 * (default Cook's modified pairing energies) *
742 * what (2) = ig + 10 * if (ig and if must have the same sign) *
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 *
755 *********************************************************************
759 IF (WHAT(1).LE.-1.0D0) THEN
766 IF ( NINT (WHAT (1)) .GE. 10000 ) THEN
768 JLVHLP = NINT (WHAT (1)) / 10000
769 WHAT (1) = WHAT (1) - 10000.D+00 * JLVHLP
771 IF ( NINT (WHAT (1)) .GE. 100 ) THEN
772 JLVMOD = NINT (WHAT (1)) / 100
773 WHAT (1) = WHAT (1) - 100.D+00 * JLVMOD
775 IF ( NINT (WHAT (1)) .GE. 10 ) THEN
777 JLVHLP = NINT (WHAT (1)) / 10
778 WHAT (1) = WHAT (1) - 10.D+00 * JLVHLP
779 ELSE IF ( NINT (WHTSAV) .NE. 0 ) THEN
782 IF ( NINT (WHAT (1)) .GE. 0 ) THEN
784 ILVMOD = NINT (WHAT(1))
785 IF ( ABS (NINT (WHAT (2))) .GE. 10 ) THEN
787 JLVHLP = NINT (WHAT (2)) / 10
788 WHAT (2) = WHAT (2) - 10.D+00 * JLVHLP
789 ELSE IF ( NINT (WHAT (2)) .NE. 0 ) THEN
792 IF ( NINT (WHAT (2)) .GE. 0 ) THEN
797 **sr heavies are always put to /FKFHVY/
798 C IF ( NINT (WHAT(3)) .GE. 1 ) THEN
814 *********************************************************************
816 * control card: codewd = EMCCHECK *
818 * extended energy-momentum / quantum-number conservation check *
820 * what (1) = -1 extended check not performed *
822 * what (2..6), sdum no meaning *
824 *********************************************************************
827 IF (WHAT(1).EQ.-1) THEN
834 *********************************************************************
836 * control card: codewd = MODEL *
838 * Model to be used to treat nucleon-nucleon interactions *
840 * sdum = DTUNUC two-chain model *
841 * = PHOJET multiple chains including minijets *
843 * = QNEUTRIN quasi-elastic neutrino scattering *
847 * what (1) (variable INTER) *
848 * = 1 gamma exchange *
851 * = 4 gamma/Z0 exchange *
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 *
861 * what (2..6) no meaning *
863 *********************************************************************
866 IF (SDUM.EQ.CMODEL(1)) THEN
868 ELSEIF (SDUM.EQ.CMODEL(2)) THEN
870 ELSEIF (SDUM.EQ.CMODEL(3)) THEN
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
877 IF ((IWHAT.EQ.1 ).OR.(IWHAT.EQ.2 ).OR.
878 & (IWHAT.EQ.10).OR.(IWHAT.EQ.11))
881 STOP ' Unknown model !'
885 *********************************************************************
887 * control card: codewd = PHOINPUT *
889 * Start of input-section for PHOJET-specific input-cards *
890 * Note: This section will not be finished before giving *
892 * what (1..6), sdum no meaning *
894 *********************************************************************
899 C 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)
905 WRITE(LOUT,'(1X,A)')'INIT: reading PHOJET-input failed'
912 *********************************************************************
914 * control card: codewd = GLAUBERI *
916 * Pre-initialization of impact parameter selection *
918 * what (1..6), sdum no meaning *
920 *********************************************************************
923 IF (IFIRST.NE.99) THEN
924 CALL DT_RNDMST(12,34,56,78)
926 OPEN(40,FILE='shm.out',STATUS='UNKNOWN')
927 C OPEN(11,FILE='shm.dbg',STATUS='UNKNOWN')
938 ADP = (APHI-APLOW)/DBLE(IPPN)
959 IT = ITLOW+(NCIT-1)*IDIT
962 C IIP = (IPHI-IPLOW)/IDIP
963 C IF (IIP.EQ.0) IIP = 1
964 C IF (IT.EQ.IPLOW) IIP = 0
968 CC IF (NCIP.LE.IIP) THEN
969 C IP = IPLOW+(NCIP-1)*IDIP
973 IF (IP.GT.IT) GOTO 472
976 APPN = APLOW+DBLE(NCP-1)*ADP
979 OPEN(12,FILE='shm.sta',STATUS='UNKNOWN')
980 WRITE(12,'(1X,2I5,E15.3)') IP,IT,PPN
987 CALL DT_NEWHGR(XDUM,XDUM,XDUM,XDUMB,-1,IHDUM)
988 CALL DT_NEWHGR(XLIM1,XLIM2,XLIM3,XDUMB,IBIN,IHSHMA)
991 C IF ((IP.GT.10).OR.(IT.GT.10)) THEN
999 CALL DT_SHMAKI(IP,IDUM1,IT,IDUM1,IJPROJ,PPN,99)
1000 SIGAV = SIGAV+XSPRO(1,1,1)
1003 CALL DT_FILHGR(XC,BSITE(1,1,1,J),IHSHMA,I)
1007 CALL DT_EVTHIS(IDUM)
1009 C CALL OUTGEN(IHSHMA,0,0,0,0,0,HEADER,0,NEVFIT,ONE,0,1,-1)
1011 C CALL GENFIT(XPARA)
1012 C WRITE(40,'(2I4,E11.3,F6.0,5E11.3)')
1013 C & IP,IT,PPN,SIGAV/DBLE(NEVFIT),XPARA
1023 *********************************************************************
1025 * control card: codewd = FLUCTUAT *
1027 * Treatment of cross section fluctuations *
1029 * what (1) = 1 treat cross section fluctuations *
1031 * what (1..6), sdum no meaning *
1033 *********************************************************************
1037 IF (WHAT(1).EQ.ONE) THEN
1043 *********************************************************************
1045 * control card: codewd = CENTRAL *
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 *
1056 * what (4..6), sdum no meaning *
1058 *********************************************************************
1061 ICENTR = INT(WHAT(1))
1062 IF (ICENTR.LT.0) THEN
1063 IF (ICENTR.GT.-100) THEN
1072 *********************************************************************
1074 * control card: codewd = RECOMBIN *
1076 * Chain recombination *
1077 * (recombine S-S and V-V chains to V-S chains) *
1079 * what (1) = -1. recombination switched off default: 1 *
1080 * what (2..6), sdum no meaning *
1082 *********************************************************************
1086 IF (WHAT(1).EQ.-1.0D0) IRECOM = 0
1089 *********************************************************************
1091 * control card: codewd = COMBIJET *
1093 * chain fusion (2 q-aq --> qq-aqaq) *
1095 * what (1) = 1 fusion treated *
1097 * what (2) minimum number of uncombined chains from *
1098 * single projectile or target nucleons *
1100 * what (3..6), sdum no meaning *
1102 *********************************************************************
1106 IF (INT(WHAT(1)).EQ.1) LCO2CR = .TRUE.
1107 IF (WHAT(2).GE.ZERO) CUTOF = WHAT(2)
1110 *********************************************************************
1112 * control card: codewd = XCUTS *
1114 * thresholds for x-sampling *
1116 * what (1) defines lower threshold for val.-q x-value (CVQ) *
1118 * what (2) defines lower threshold for val.-qq x-value (CDQ) *
1120 * what (3) defines lower threshold for sea-q x-value (CSEA) *
1122 * what (4) sea-q x-values in S-S chains (SSMIMA) *
1124 * what (5) not used *
1126 * what (6), sdum no meaning *
1128 * Note: Lower thresholds (what(1..3)) are def. as x_thr=CXXX/ECM *
1130 *********************************************************************
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
1140 IF (WHAT(5).GT.2.0D0) VVMTHR = WHAT(5)
1143 *********************************************************************
1145 * control card: codewd = INTPT *
1147 * what (1) = -1 intrinsic transverse momenta of partons *
1148 * not treated default: 1 *
1149 * what (2..6), sdum no meaning *
1151 *********************************************************************
1154 IF (WHAT(1).EQ.-1.0D0) THEN
1161 *********************************************************************
1163 * control card: codewd = CRONINPT *
1165 * Cronin effect (multiple scattering of partons at chain ends) *
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 *
1171 *********************************************************************
1174 IF (WHAT(1).EQ.-1.0D0) THEN
1182 *********************************************************************
1184 * control card: codewd = SEADISTR *
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 *
1193 *********************************************************************
1197 XSEACU = 1.05D0-XSEACO
1199 IF (UNON.LT.0.1D0) UNON = 2.0D0
1201 IF (UNOM.LT.0.1D0) UNOM = 1.5D0
1203 IF (UNOSEA.LT.0.1D0) UNOSEA = 5.0D0
1206 *********************************************************************
1208 * control card: codewd = SEASU3 *
1210 * Treatment of strange-quarks at chain ends *
1212 * what (1) (SEASQ) strange-quark supression factor *
1213 * iflav = 1.+rndm*(2.+SEASQ) *
1215 * what (2..6), sdum no meaning *
1217 *********************************************************************
1223 *********************************************************************
1225 * control card: codewd = DIQUARKS *
1227 * what (1) = -1. sea-diquark/antidiquark-pairs not treated *
1229 * what (2..6), sdum no meaning *
1231 *********************************************************************
1234 IF (WHAT(1).EQ.-1.0D0) THEN
1241 *********************************************************************
1243 * control card: codewd = RESONANC *
1245 * treatment of low mass chains *
1247 * what (1) = -1 low chain masses are not corrected for resonance *
1248 * masses (obsolete for BAMJET-fragmentation) *
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 *
1257 *********************************************************************
1263 IF (WHAT(1).EQ.-ONE) IRESCO = 0
1264 IF (WHAT(2).EQ.-ONE) IMSHL = 0
1265 IF (WHAT(3).EQ.-ONE) IRESRJ = 1
1268 *********************************************************************
1270 * control card: codewd = DIFFRACT *
1272 * Treatment of diffractive events *
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) *
1283 * what (2) = (IDOUBD) 0 no double diffraction *
1284 * 1 double diffraction included *
1285 * 2 double diffractive events only *
1287 * what (3) = 1 projectile diffraction treated (2-channel form.) *
1289 * what (4) = alpha-parameter in projectile diffraction *
1291 * what (5..6), sdum no meaning *
1293 *********************************************************************
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
1300 1380 FORMAT(1X,'INIT: inconsistent DIFFRACT - input !',/,
1301 & 11X,'IDOUBD is reset to zero')
1304 IF (WHAT(3).GT.ZERO) DIBETA = WHAT(3)
1305 IF (WHAT(4).GT.ZERO) DIALPH = WHAT(4)
1308 *********************************************************************
1310 * control card: codewd = SINGLECH *
1312 * what (1) = 1. Regge contribution (one chain) included *
1314 * what (2..6), sdum no meaning *
1316 *********************************************************************
1320 IF (WHAT(1).EQ.ONE) ISICHA = 1
1323 *********************************************************************
1325 * control card: codewd = NOFRAGME *
1327 * biased chain hadronization *
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 *
1341 *********************************************************************
1345 ICHAIN = INT(WHAT(I))
1346 IF ((ICHAIN.LE.-1).AND.(ICHAIN.GE.-9))
1347 & LHADRO(ABS(ICHAIN)) = .FALSE.
1351 *********************************************************************
1353 * control card: codewd = HADRONIZE *
1355 * hadronization model and parameter switch *
1357 * what (1) = 1 hadronization via BAMJET *
1358 * = 2 hadronization via JETSET *
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 *
1365 * what (3..6), sdum no meaning *
1367 *********************************************************************
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))
1377 *********************************************************************
1379 * control card: codewd = POPCORN *
1381 * "Popcorn-effect" in fragmentation and diquark breaking diagrams *
1383 * what (1) = (PDB) frac. of diquark fragmenting directly into *
1384 * baryons (PYTHIA/JETSET fragmentation) *
1385 * (JETSET: = 0. Popcorn mechanism switched off) *
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 *
1392 * what (4..6), sdum no meaning *
1394 *********************************************************************
1397 IF (WHAT(1).GE.0.0D0) PDB = WHAT(1)
1398 IF (WHAT(2).GE.0.0D0) THEN
1402 IF (WHAT(3).GE.0.0D0) PDBSEA(3) = WHAT(3)
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))
1410 *********************************************************************
1412 * control card: codewd = PARDECAY *
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 *
1419 *********************************************************************
1422 IF (WHAT(1).EQ.ONE) ISIG0 = 1
1423 IF (WHAT(1).EQ.2.0D0) IPI0 = 1
1426 *********************************************************************
1428 * control card: codewd = BEAM *
1430 * definition of beam parameters *
1432 * what (1/2) > 0 : energy of beam 1/2 (GeV) *
1433 * < 0 : abs(what(1/2)) energy per charge of *
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 *
1441 * Note: this card requires previously defined projectile and *
1442 * target identities (PROJPAR, TARPAR) *
1444 *********************************************************************
1447 CALL DT_BEAMPR(WHAT,PPN,1)
1453 *********************************************************************
1455 * control card: codewd = LUND-MSTU *
1457 * set parameter MSTU in JETSET-common /LUDAT1/ *
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 *
1464 * default: default-Lund or corresponding to *
1465 * the set given in HADRONIZE *
1467 *********************************************************************
1470 IF (WHAT(1).GT.ZERO) THEN
1472 IMSTU(NMSTU) = INT(WHAT(1))
1473 MSTUX(NMSTU) = INT(WHAT(2))
1475 IF (WHAT(3).GT.ZERO) THEN
1477 IMSTU(NMSTU) = INT(WHAT(3))
1478 MSTUX(NMSTU) = INT(WHAT(4))
1480 IF (WHAT(5).GT.ZERO) THEN
1482 IMSTU(NMSTU) = INT(WHAT(5))
1483 MSTUX(NMSTU) = INT(WHAT(6))
1487 *********************************************************************
1489 * control card: codewd = LUND-MSTJ *
1491 * set parameter MSTJ in JETSET-common /LUDAT1/ *
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 *
1498 * default: default-Lund or corresponding to *
1499 * the set given in HADRONIZE *
1501 *********************************************************************
1504 IF (WHAT(1).GT.ZERO) THEN
1506 IMSTJ(NMSTJ) = INT(WHAT(1))
1507 MSTJX(NMSTJ) = INT(WHAT(2))
1509 IF (WHAT(3).GT.ZERO) THEN
1511 IMSTJ(NMSTJ) = INT(WHAT(3))
1512 MSTJX(NMSTJ) = INT(WHAT(4))
1514 IF (WHAT(5).GT.ZERO) THEN
1516 IMSTJ(NMSTJ) = INT(WHAT(5))
1517 MSTJX(NMSTJ) = INT(WHAT(6))
1521 *********************************************************************
1523 * control card: codewd = LUND-MDCY *
1525 * set parameter MDCY(I,1) for particle decays in JETSET-common *
1528 * what (1-6) = PDG particle index of particle which should *
1530 * default: default-Lund or forced in *
1533 *********************************************************************
1537 IF (WHAT(I).NE.ZERO) THEN
1539 KC = PYCOMP(INT(WHAT(I)))
1546 *********************************************************************
1548 * control card: codewd = LUND-PARJ *
1550 * set parameter PARJ in JETSET-common /LUDAT1/ *
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 *
1557 * default: default-Lund or corresponding to *
1558 * the set given in HADRONIZE *
1560 *********************************************************************
1563 IF (WHAT(1).NE.ZERO) THEN
1565 IPARJ(NPARJ) = INT(WHAT(1))
1566 PARJX(NPARJ) = WHAT(2)
1568 IF (WHAT(3).NE.ZERO) THEN
1570 IPARJ(NPARJ) = INT(WHAT(3))
1571 PARJX(NPARJ) = WHAT(4)
1573 IF (WHAT(5).NE.ZERO) THEN
1575 IPARJ(NPARJ) = INT(WHAT(5))
1576 PARJX(NPARJ) = WHAT(6)
1580 *********************************************************************
1582 * control card: codewd = LUND-PARU *
1584 * set parameter PARJ in JETSET-common /LUDAT1/ *
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 *
1591 * default: default-Lund or corresponding to *
1592 * the set given in HADRONIZE *
1594 *********************************************************************
1597 IF (WHAT(1).GT.ZERO) THEN
1599 IPARU(NPARU) = INT(WHAT(1))
1600 PARUX(NPARU) = WHAT(2)
1602 IF (WHAT(3).GT.ZERO) THEN
1604 IPARU(NPARU) = INT(WHAT(3))
1605 PARUX(NPARU) = WHAT(4)
1607 IF (WHAT(5).GT.ZERO) THEN
1609 IPARU(NPARU) = INT(WHAT(5))
1610 PARUX(NPARU) = WHAT(6)
1614 *********************************************************************
1616 * control card: codewd = OUTLEVEL *
1618 * output control switches *
1620 * what (1) = internal rejection informations default: 0 *
1621 * what (2) = energy-momentum conservation check output *
1623 * what (3) = internal warning messages default: 0 *
1624 * what (4..6), sdum not yet used *
1626 *********************************************************************
1630 IOULEV(K) = INT(WHAT(K))
1634 *********************************************************************
1636 * control card: codewd = FRAME *
1638 * frame in which final state is given in DTEVT1 *
1640 * what (1) = 1 target rest frame (laboratory) *
1641 * = 2 nucleon-nucleon cms *
1644 *********************************************************************
1647 KFRAME = INT(WHAT(1))
1648 IF ((KFRAME.GE.1).AND.(KFRAME.LE.2)) IFRAME = KFRAME
1651 *********************************************************************
1653 * control card: codewd = L-TAG *
1656 * definition of kinematical cuts for radiated photon and *
1657 * outgoing lepton detection in lepton-nucleus interactions *
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 *
1668 *********************************************************************
1679 *********************************************************************
1681 * control card: codewd = L-ETAG *
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 *
1690 *********************************************************************
1693 ELMIN = MAX(WHAT(1),ZERO)
1694 EGMIN = MAX(WHAT(2),ZERO)
1695 EGMAX = MAX(WHAT(3),ZERO)
1698 *********************************************************************
1700 * control card: codewd = ECMS-CUT *
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 *
1708 *********************************************************************
1713 IF (ECMIN.GT.ECMAX) ECMIN = ECMAX
1714 XBJMIN = MAX(WHAT(3),ZERO)
1717 *********************************************************************
1719 * control card: codewd = VDM-PAR1 *
1721 * parameters in gamma-nucleus cross section calculation *
1723 * what (1) = Lambda^2 default: 2. *
1724 * what (2) lower limit in M^2 integration *
1727 * = 3 (m_phi)^2 default: 1 *
1728 * what (3) upper limit in M^2 integration *
1731 * = 3 s default: 3 *
1732 * what (4) CKMT F_2 structure function *
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 *
1739 * = 4 PHOJET cross sections default: 4 *
1741 * what (6), sdum no meaning *
1743 *********************************************************************
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))
1753 *********************************************************************
1755 * control card: codewd = HISTOGRAM *
1757 * activate different classes of histograms *
1759 * default: no histograms *
1761 *********************************************************************
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
1774 *********************************************************************
1776 * control card: codewd = XS-TABLE *
1778 * output of cross section table for requested interaction *
1779 * - particle production deactivated ! - *
1781 * what (1) lower energy limit for tabulation *
1783 * < 0 nucleon-nucleon cms *
1784 * what (2) upper energy limit for tabulation *
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 *
1794 *********************************************************************
1797 IF (WHAT(1).EQ.99999.0D0) THEN
1798 IRATIO = INT(WHAT(2))
1801 CMENER = ABS(WHAT(2))
1802 IF (.NOT.LXSTAB) THEN
1808 IF ((.NOT.LXSTAB).OR.(CMENER.NE.CMEOLD)) THEN
1810 IF (WHAT(2).GT.ZERO)
1811 & CMENER = SQRT(2.0D0*AAM(1)**2+2.0D0*WHAT(2)*AAM(1))
1814 C WRITE(LOUT,*) 'CMENER = ',CMENER
1815 CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,CMENER,1)
1818 CALL DT_XSTABL(WHAT,IXSQEL,IRATIO)
1823 *********************************************************************
1825 * control card: codewd = GLAUB-PAR *
1827 * parameters in Glauber-formalism *
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 *
1835 * what (4) = 1 read pre-calculated impact-parameter distrib. *
1837 * =-1 dump pre-calculated impact-parameter distrib. *
1839 * = 100 read pre-calculated impact-parameter distrib. *
1840 * for variable projectile/target/energy runs *
1843 * what (5..6) no meaning *
1844 * sdum if |what (4)| = 1 name of in/output-file (sdum.glb) *
1846 *********************************************************************
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))
1858 *********************************************************************
1860 * control card: codewd = GLAUB-INI *
1862 * pre-initialization of profile function *
1864 * what (1) lower energy limit for initialization *
1866 * < 0 nucleon-nucleon cms *
1867 * what (2) upper energy limit for initialization *
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) *
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 *
1881 * what (6) no meaning *
1884 *********************************************************************
1888 CALL DT_GLBINI(WHAT)
1891 *********************************************************************
1893 * control card: codewd = VDM-PAR2 *
1895 * parameters in gamma-nucleus cross section calculation *
1897 * what (1) = 0 no suppression of shadowing by direct photon *
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 *
1908 * eps*R*Q^2/M^2 = what(4)*Q^2/M^2 default: 0 *
1909 * what (5..6), sdum no meaning *
1911 *********************************************************************
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))
1920 *********************************************************************
1922 * control card: XS-QELPRO *
1924 * what (1..6), sdum no meaning *
1926 *********************************************************************
1929 IXSQEL = ABS(WHAT(1))
1932 *********************************************************************
1934 * control card: RNDMINIT *
1936 * initialization of random number generator *
1938 * what (1..4) values for initialization (= 1..168) *
1939 * what (5..6), sdum no meaning *
1941 *********************************************************************
1944 IF ((WHAT(1).LT.1.0D0).OR.(WHAT(1).GT.168.0D0)) THEN
1949 IF ((WHAT(2).LT.1.0D0).OR.(WHAT(2).GT.168.0D0)) THEN
1954 IF ((WHAT(3).LT.1.0D0).OR.(WHAT(3).GT.168.0D0)) THEN
1959 IF ((WHAT(4).LT.1.0D0).OR.(WHAT(4).GT.168.0D0)) THEN
1964 CALL DT_RNDMST(NA1,NA2,NA3,NA4)
1967 *********************************************************************
1969 * control card: codewd = LEPTO-CUT *
1971 * set parameter CUT in LEPTO-common /LEPTOU/ *
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 *
1978 * default: default-LEPTO parameters *
1980 *********************************************************************
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)
1988 *********************************************************************
1990 * control card: codewd = LEPTO-LST *
1992 * set parameter LST in LEPTO-common /LEPTOU/ *
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 *
1999 * default: default-LEPTO parameters *
2001 *********************************************************************
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))
2009 *********************************************************************
2011 * control card: codewd = LEPTO-PARL *
2013 * set parameter PARL in LEPTO-common /LEPTOU/ *
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 *
2020 * default: default-LEPTO parameters *
2022 *********************************************************************
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)
2030 *********************************************************************
2032 * control card: codewd = START *
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 *
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 *
2047 *********************************************************************
2051 * check for cross-section table output only
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))
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..)
2070 C IF (IDP.EQ.26) IDP = 5
2071 C IF (IDP.EQ.27) IDP = 6
2073 * redefine collision energy
2075 IF (ABS(VAREHI).GT.ZERO) THEN
2077 IF (VARELO.LT.EHADLO) VARELO = EHADLO
2078 CALL DT_LTINI(IDP,IDT,VARELO,PDUM,VARCLO,1)
2080 CALL DT_LTINI(IDP,IDT,VAREHI,PDUM,VARCHI,1)
2082 CALL DT_LTINI(IDP,IDT,EPN,PPN,CMENER,1)
2085 1003 FORMAT(1X,'INIT: collision energy not defined!',/,
2086 & 1X,' -program stopped- ')
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
2094 1004 FORMAT(1X,/,'Warning! Evaporation request rejected since',
2095 & ' central collisions forced.')
2102 * initialization of evaporation-module
2104 * initialize evaporation if the code is not used as Fluka event generator
2105 IF (ITRSPT.NE.1) THEN
2109 IF (LEVPRT) LHEAVY = .TRUE.
2112 * save the default JETSET-parameter
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
2122 STOP ' This version does not contain LEPTO !'
2126 * initialization of quasi-elastic neutrino scattering
2127 IF (MCGENE.EQ.4) THEN
2128 IF (IJPROJ.EQ.5) THEN
2130 ELSEIF (IJPROJ.EQ.6) THEN
2132 ELSEIF (IJPROJ.EQ.135) THEN
2134 ELSEIF (IJPROJ.EQ.136) THEN
2136 ELSEIF (IJPROJ.EQ.133) THEN
2138 ELSEIF (IJPROJ.EQ.134) THEN
2143 * normalize fractions of emulsion components
2144 IF (NCOMPO.GT.0) THEN
2147 SUMFRA = SUMFRA+EMUFRA(I)
2149 IF (SUMFRA.GT.ZERO) THEN
2151 EMUFRA(I) = EMUFRA(I)/SUMFRA
2156 * disallow Cronin's multiple scattering for nucleus-nucleus interactions
2157 IF ((IP.GT.1).AND.(MKCRON.GT.0)) THEN
2159 1005 FORMAT(/,1X,'INIT: multiple scattering disallowed',/)
2163 * initialization of Glauber-formalism (moved to xAEVT, sr 26.3.96)
2164 C IF (NCOMPO.LE.0) THEN
2165 C CALL DT_SHMAKI(IP,IPZ,IT,ITZ,IDP,PPN,IGLAU)
2168 C CALL DT_SHMAKI(IP,IPZ,IEMUMA(I),IEMUCH(I),IDP,PPN,0)
2172 * pre-tabulation of elastic cross-sections
2173 CALL DT_SIGTBL(JDUM,JDUM,DUM,DUM,-1)
2179 *********************************************************************
2181 * control card: codewd = STOP *
2183 * stop of the event generation *
2185 * what (1..6) no meaning *
2187 *********************************************************************
2191 9000 FORMAT(1X,'---> unexpected end of input !')
2198 *===kkinc==============================================================*
2201 SUBROUTINE DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,
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 ************************************************************************
2211 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2214 PARAMETER ( LINP = 5 ,
2218 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY5=1.0D-5,
2219 & TINY2=1.0D-2,TINY3=1.0D-3)
2225 PARAMETER (NMXHKK=200000)
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),
2234 * particle properties (BAMJET index convention)
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
2256 COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
2263 IF (ILOOP.EQ.4) THEN
2264 WRITE(LOUT,1000) NEVHKK
2265 1000 FORMAT(1X,'KKINC: event ',I8,' rejected!')
2270 * variable energy-runs, recalculate parameters for LT's
2271 IF ((ABS(VAREHI).GT.ZERO).OR.(IOGLB.EQ.100)) THEN
2274 CALL DT_LTINI(IDP,1,EPN,PDUM,CDUM,1)
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) !'
2283 * re-initialize /DTPRTA/
2289 IBPROJ = IIBAR(IJPROJ)
2291 * calculate nuclear potentials (common /DTNPOT/)
2292 CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
2294 * initialize treatment for residual nuclei
2295 CALL DT_RESNCL(EPN,NLOOP,1)
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'
2304 IF ((NPMASS.GT.1).OR.(NTMASS.GT.1)) THEN
2306 * intranuclear cascade of final state particles for KTAUGE generations
2308 CALL DT_FOZOCA(LFZC,IREJ1)
2309 IF (IREJ1.GT.0) THEN
2310 IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2 in KKINC'
2314 * baryons unable to escape the nuclear potential are treated as
2315 * excited nucleons (ISTHKK=15,16)
2318 * decay of resonances produced in intranuclear cascade processes
2319 **sr 15-11-95 should be obsolete
2320 C IF (LFZC) CALL DT_DECAY1
2323 * treatment of residual nuclei
2324 CALL DT_RESNCL(EPN,NLOOP,2)
2326 * evaporation / fission / fragmentation
2327 * (if intranuclear cascade was sampled only)
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
2336 * transform finale state into Lab.
2338 CALL DT_BEAMPR(WHAT,DUM,IFLAG)
2339 IF ((IFRAME.EQ.1).AND.(IFLAG.EQ.-1)) CALL DT_LT2LAB
2341 IF (IPI0.EQ.1) CALL DT_DECPI0
2343 C IF (NEVHKK.EQ.5) CALL DT_EVTOUT(4)
2351 *===defaul=============================================================*
2353 CDECK ID>, DT_DEFAUL
2354 SUBROUTINE DT_DEFAUL(EPN,PPN)
2356 ************************************************************************
2357 * Variables are set to default values. *
2358 * This version dated 8.5.95 is written by S. Roesler. *
2359 ************************************************************************
2361 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2363 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
2364 PARAMETER (TWOPI = 6.283185307179586454D+00)
2366 * particle properties (BAMJET index convention)
2368 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2369 & IICH(210),IIBAR(210),K1(210),K2(210)
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
2384 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2386 * emulsion treatment
2387 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2389 * parameter for intranuclear cascade
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,
2397 * threshold values for x-sampling (DTUNUC 1.x)
2398 COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
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
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
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
2436 COMMON /LEPTOI/ RPPN,LEPIN,INTER
2437 * steering flags for qel neutrino scattering modules
2438 COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
2440 COMMON /DTEVNO/ NEVENT,ICASCA
2442 DATA POTMES /0.002D0/
2453 * nucleus independent meson potential
2462 **sr 7.4.98: changed after corrected B-sampling
2503 **sr 7.4.98: changed after corrected B-sampling
2522 * definition of soft quark distributions
2527 * cutoff parameters for x-sampling
2573 CMODEL(1) = 'DTUNUC '
2574 CMODEL(2) = 'PHOJET '
2575 CMODEL(3) = 'LEPTO '
2576 CMODEL(4) = 'QNEUTRIN'
2613 IF (ITRSPT.EQ.1) THEN
2648 IF (ITRSPT.EQ.1) THEN
2654 * default Lab.-energy
2656 PPN = SQRT((EPN-AAM(IJPROJ))*(EPN+AAM(IJPROJ)))
2661 *===aaevt==============================================================*
2664 SUBROUTINE DT_AAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
2667 ************************************************************************
2668 * This version dated 22.03.96 is written by S. Roesler. *
2669 ************************************************************************
2671 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2674 PARAMETER ( LINP = 5 ,
2678 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2680 * emulsion treatment
2681 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
2684 COMMON /DTEVNO/ NEVENT,ICASCA
2686 CHARACTER*8 DATE,HHMMSS
2690 NMSG = MAX(NEVTS/100,1)
2692 * initialization of run-statistics and histograms
2695 CALL PHO_PHIST(1000,DUM)
2697 * initialization of Glauber-formalism
2698 IF (NCOMPO.LE.0) THEN
2699 CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
2702 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2708 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2709 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
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,' )')
2717 * generate NEVTS events
2720 * print run-status message
2721 IF (MOD(IEVT,NMSG).EQ.0) THEN
2723 WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
2724 & IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
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,' )',/)
2731 C WRITE(LOUT,1000) IEVT-1
2732 C1000 FORMAT(1X,I8,' events sampled')
2735 * treat nuclear emulsions
2736 IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
2737 * composite targets only
2740 CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,IREJ)
2742 CALL PHO_PHIST(2000,DUM)
2746 * print run-statistics and histograms to output-unit 6
2748 CALL PHO_PHIST(3000,DUM)
2755 *===laevt==============================================================*
2758 SUBROUTINE DT_LAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
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 ************************************************************************
2768 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2771 PARAMETER ( LINP = 5 ,
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,
2779 & ALPHEM = ONE/137.0D0)
2781 C CHARACTER*72 HEADER
2783 * particle properties (BAMJET index convention)
2785 COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
2786 & IICH(210),IIBAR(210),K1(210),K2(210)
2789 PARAMETER (NMXHKK=200000)
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),
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
2810 PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
2812 * emulsion treatment
2813 COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
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
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
2835 COMMON /DTEVNO/ NEVENT,ICASCA
2837 DIMENSION XDUMB(40),BGTA(4)
2840 IF (MCGENE.EQ.3) THEN
2842 STOP ' This version does not contain LEPTO !'
2847 NMSG = MAX(NEVTS/10,1)
2849 * mass of incident lepton
2852 IDPPDG = IDT_IPDGHA(IDP)
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)
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,
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',/)
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)
2901 * limits for Glauber-initialization
2903 Q2HI = MAX(Q2LI,MIN(Q2HI,Q2MAX))
2904 ECMLI = MAX(EGNMIN,THREE)
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)
2917 CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
2922 * initialization of run-statistics and histograms
2925 CALL PHO_PHIST(1000,DUM)
2927 * maximum photon-nucleus cross section