| 0b5dd071 |
1 | CDECK ID>, GGCDES. |
| 2 | *======================================================================= |
| 3 | CDECK ID>, GGINIT. |
| 4 | SUBROUTINE gginit |
| 5 | *----------------------------------------------------------------------- |
| 6 | * GGHIC initialization |
| 7 | * Author: Yu.Kharlov |
| 8 | *----------------------------------------------------------------------- |
| 9 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 10 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 11 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 12 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 13 | CHARACTER lumfil*80 |
| 14 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 15 | & gvconst(4,10) |
| 16 | REAL nc |
| 17 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 18 | & xstote, ssbr(10) |
| 19 | COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5), |
| 20 | & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5) |
| 21 | COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr, |
| 22 | & xmll, xmlr, xmnl, xtanb, xmha, xmu, |
| 23 | & xmt, xat, xmbr, xab, u11, v11 |
| 24 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 25 | INTEGER lz,la |
| 26 | REAL Rion,Gamma |
| 27 | COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) |
| 28 | DOUBLE PRECISION CKIN |
| 29 | SAVE /PYSUBS/ |
| 30 | COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) |
| 31 | DOUBLE PRECISION PARP,PARI |
| 32 | SAVE /PYPARS/ |
| 33 | COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) |
| 34 | DOUBLE PRECISION P,V |
| 35 | SAVE /PYJETS/ |
| 36 | COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) |
| 37 | DOUBLE PRECISION PARU,PARJ |
| 38 | SAVE /PYDAT1/ |
| 39 | COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) |
| 40 | DOUBLE PRECISION PMAS,PARF,VCKM |
| 41 | SAVE /PYDAT2/ |
| 42 | COMMON/PYDAT4/CHAF(500,2) |
| 43 | CHARACTER CHAF*16 |
| 44 | SAVE /PYDAT4/ |
| 45 | C MXSS = maximum number of modes |
| 46 | C NSSMOD = number of modes |
| 47 | C ISSMOD = initial particle |
| 48 | C JSSMOD = final particles |
| 49 | C GSSMOD = width |
| 50 | C BSSMOD = branching ratio |
| 51 | INTEGER MXSS |
| 52 | PARAMETER (MXSS=1000) |
| 53 | COMMON/SSMODE/NSSMOD,ISSMOD(MXSS),JSSMOD(5,MXSS),GSSMOD(MXSS) |
| 54 | $,BSSMOD(MXSS) |
| 55 | INTEGER NSSMOD,ISSMOD,JSSMOD |
| 56 | REAL GSSMOD,BSSMOD |
| 57 | SAVE /SSMODE/ |
| 58 | C SUSY parameters |
| 59 | C AMGLSS = gluino mass |
| 60 | C AMQKSS = squark mass |
| 61 | C AMLLSS = left-slepton mass |
| 62 | C AMLRSS = right-slepton mass |
| 63 | C AMNLSS = sneutrino mass |
| 64 | C TWOM1 = Higgsino mass = - mu |
| 65 | C RV2V1 = ratio v2/v1 of vev's |
| 66 | C AMTLSS,AMTRSS = left,right stop masses |
| 67 | C AMT1SS,AMT2SS = light,heavy stop masses |
| 68 | C AMBLSS,AMBRSS = left,right sbottom masses |
| 69 | C AMB1SS,AMB2SS = light,heavy sbottom masses |
| 70 | C AMZiSS = signed mass of Zi |
| 71 | C ZMIXSS = Zi mixing matrix |
| 72 | C AMWiSS = signed Wi mass |
| 73 | C GAMMAL,GAMMAR = Wi left, right mixing angles |
| 74 | C AMHL,AMHH,AMHA = neutral Higgs h0, H0, A0 masses |
| 75 | C AMHC = charged Higgs H+ mass |
| 76 | C ALFAH = Higgs mixing angle |
| 77 | C AAT = stop trilinear term |
| 78 | C THETAT = stop mixing angle |
| 79 | C AAB = sbottom trilinear term |
| 80 | C THETAB = sbottom mixing angle |
| 81 | COMMON/SSPAR/AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS |
| 82 | $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS |
| 83 | $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS |
| 84 | $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS |
| 85 | $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS(4,4) |
| 86 | $,AMW1SS,AMW2SS |
| 87 | $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT |
| 88 | $,AAB,THETAB,AAL,THETAL,AMGVSS |
| 89 | REAL AMGLSS,AMULSS,AMURSS,AMDLSS,AMDRSS,AMSLSS |
| 90 | $,AMSRSS,AMCLSS,AMCRSS,AMBLSS,AMBRSS,AMB1SS,AMB2SS |
| 91 | $,AMTLSS,AMTRSS,AMT1SS,AMT2SS,AMELSS,AMERSS,AMMLSS,AMMRSS |
| 92 | $,AMLLSS,AMLRSS,AML1SS,AML2SS,AMN1SS,AMN2SS,AMN3SS |
| 93 | $,TWOM1,RV2V1,AMZ1SS,AMZ2SS,AMZ3SS,AMZ4SS,ZMIXSS |
| 94 | $,AMW1SS,AMW2SS |
| 95 | $,GAMMAL,GAMMAR,AMHL,AMHH,AMHA,AMHC,ALFAH,AAT,THETAT |
| 96 | $,AAB,THETAB,AAL,THETAL,AMGVSS |
| 97 | REAL AMZISS(4) |
| 98 | EQUIVALENCE (AMZISS(1),AMZ1SS) |
| 99 | SAVE /SSPAR/ |
| 82764a9c |
100 | REAL xpar(10),ypar(6) |
| 0b5dd071 |
101 | INTEGER ipar(2) |
| 102 | INTEGER pycomp |
| 103 | REAL*8 ggrnd |
| 104 | EXTERNAL pydata |
| 105 | la = ia |
| 106 | lz = iz |
| 107 | * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 108 | * |
| 109 | * Check bounds on 2 photon mass |
| 110 | * |
| 111 | IF (iproc .EQ. 1) THEN |
| 112 | IF (amin .LE. 2.31) THEN |
| 113 | WRITE (6,*) |
| 114 | & ' %GGINIT: minimal 2-photon mass is too low, ', |
| 115 | & 'set to 3m(pho)' |
| 116 | amin = 2.31 |
| 117 | END IF |
| 118 | * |
| 119 | ELSE IF (iproc .EQ. 2) THEN |
| 120 | kc = pycomp (kf_onium) |
| 121 | xmres = pmas(kc,1) |
| 122 | IF (kf_onium .EQ. 441) THEN |
| 123 | xgtres = 0.013 |
| 124 | ELSE |
| 125 | xgtres = pmas(kc,2) |
| 126 | END IF |
| 127 | c. xggres = ggwid(kf_onium) |
| 128 | IF (xgtres .EQ. 0.) xgtres = 0.1 |
| 129 | IF (xmres.LE.amin .OR. xmres.GE.amax) THEN |
| 130 | WRITE(6,*) |
| 131 | & ' %GGINIT: Xmres is out of the Luminosity mass range' |
| 132 | STOP |
| 133 | ENDIF |
| 134 | * |
| 135 | ELSE IF (iproc .EQ. 3) THEN |
| 136 | kc = pycomp (kferm) |
| 137 | amprt(1) = pmas(kc,1) |
| 138 | IF (amin .LT. 2*amprt(1)) THEN |
| 139 | WRITE (6,*) |
| 140 | & ' %GGINIT: minimal 2-photon mass is too low, ', |
| 141 | & 'set to 2m(fermion)' |
| 142 | amin = 2*amprt(1) |
| 143 | END IF |
| 144 | * |
| 145 | ELSE IF (iproc .EQ. 4) THEN |
| 146 | IF (amprt(4) .LE. 0.) amprt(4) = pmas(pycomp(24),1) |
| 147 | IF (amin .LT. 2*amprt(4)) then |
| 148 | WRITE (6,*) |
| 149 | & ' %GGINIT: minimal 2-photon mass is too low, ', |
| 150 | & 'set to 2m(W+-)' |
| 151 | amin = 2*amprt(4) |
| 152 | END IF |
| 153 | * |
| 154 | ELSE IF (iproc .EQ. 5) THEN |
| 155 | WRITE (6,*) 'Process 5 is not implemented yet' |
| 156 | STOP |
| 157 | C CALL ssmssm (xm1, xm2, xmg, xms, xmtl, xmtr, xmll, xmlr, xmnl, |
| 158 | C & xtanb, xmha, xmu, xmt, xat, xmbr, xab, iallow) |
| 159 | C IF (iallow .NE. 0) THEN |
| 160 | C WRITE (6,*) 'Lightest neutralino is not LSP, ', |
| 161 | C & 'input other MSSM parameters' |
| 162 | C STOP |
| 163 | C END IF |
| 164 | C amprt(1) = abs (amw1ss) |
| 165 | C amprt(2) = abs (amz1ss) |
| 166 | C WRITE (6,'('' %GGINIT'','// |
| 167 | C & ''', M(chi+) = '',f5.1,'// |
| 168 | C & ''', M(chi0) = '',f5.1)') amprt(1), amprt(2) |
| 169 | C * |
| 170 | C IF (amin .LT. 2*amprt(1)) THEN |
| 171 | C WRITE (6,*) |
| 172 | C & ' %GGINIT: minimal 2-photon mass is too low, ', |
| 173 | C & 'set to 2M(chi+)' |
| 174 | C amin = 2*amprt(1) |
| 175 | C END IF |
| 176 | * |
| 177 | ELSE IF (iproc .EQ. 6) THEN |
| 178 | kc1 = pycomp (kv1) |
| 179 | amprt(1) = pmas(kc1,1) + pmas(kc1,3) |
| 180 | kc2 = pycomp (kv2) |
| 181 | amprt(2) = pmas(kc2,1) + pmas(kc2,3) |
| 182 | IF (amin .LT. amprt(1)+amprt(2)) THEN |
| 183 | WRITE (6,*) |
| 184 | & ' %GGINIT: minimal 2-photon mass is too low, ', |
| 185 | & 'set to m(kv1)+m(kv2)' |
| 186 | amin = amprt(1)+amprt(2) |
| 187 | END IF |
| 188 | * |
| 189 | END IF |
| 190 | * |
| 191 | IF (amin .GT. amax) THEN |
| 192 | WRITE (6,*) ' %GGINIT: AMIN > AMAX, execution stopped' |
| 193 | STOP |
| 194 | END IF |
| 195 | * |
| 196 | * ----- TWOGAM Initialization ----- |
| 197 | * |
| 198 | gamma = ebmn*ia/amas |
| 199 | xpar(1)= ebmn*ia ! Beam energy |
| 200 | xpar(2)= amin ! Minimum gg mass |
| 201 | xpar(3)= amax ! Maximum gg mass |
| 202 | xpar(8)= 100. ! Weight, but not too long |
| 203 | * |
| 204 | ipar(1)= 1 ! Unweighted events |
| 205 | ipar(2)= 0 ! Continuum generation |
| 206 | * |
| 207 | IF (iproc .EQ. 2) THEN |
| 208 | xpar(9) = xmres ! Resonance mass |
| 209 | xpar(10)= xgtres ! Resonance total width |
| 210 | ipar(2) = 1 ! Resonance generation |
| 211 | END IF |
| 212 | * |
| 213 | ypar(1)= iz ! Ion charge |
| 214 | ypar(2)= ia ! Ion atomic number |
| 215 | ypar(3)= amas ! Ion mass |
| 216 | ypar(4)= gamma ! Gamma factor of ion |
| 217 | ypar(5)= ymin ! Minimum gg-rapidity |
| 218 | ypar(6)= ymax ! Maximum gg-rapidity |
| 219 | eb = xpar(1) |
| 220 | * |
| 221 | CALL twoini (6,ipar,xpar,ypar) |
| 222 | * |
| 223 | * The Luminosity function initialization |
| 224 | * |
| 225 | * ilumf = +1 ! Reading the Luminosity from file |
| 226 | * ilumf = -1 ! New calculation of the Luminosity function |
| 227 | * |
| 228 | CALL lumfun (ilumf,lumfil,ymin,ymax,ny,amin,amax,nmas) |
| 229 | * |
| 230 | * Integral of lum. function over the (M,y) region |
| 231 | xlumint = 0.0 |
| 232 | hmas = (amax-amin)/nmas |
| 233 | f1 = dldmx(amin ,ny) |
| 234 | DO i=1,nmas |
| 235 | xm = amin + i*hmas |
| 236 | f2 = dldmx(xm-hmas/2.,ny) |
| 237 | f3 = dldmx(xm ,ny) |
| 238 | xlumint = xlumint +(f1 + 4.*f2 + f3)/6. |
| 239 | f1 = f3 |
| 240 | END DO |
| 241 | xlumint = xlumint * hmas |
| 242 | * ----- end of TWOGAM initialization ----- |
| 243 | * |
| 244 | IF (iproc .EQ. 1) THEN |
| 245 | * ----- PYTHIA initialization ----- |
| 246 | parp(2) = dble(amin) |
| 247 | mstp(14) = 20 |
| 248 | mstp(82) = 1 |
| 249 | mstp(171) = 1 |
| 250 | mstp(172) = 1 |
| 251 | DO ip = 1,2 |
| 252 | p(1,ip) = 0.D0 |
| 253 | p(2,ip) = 0.D0 |
| 254 | END DO |
| 255 | p(1,3) = dble( amax/2) |
| 256 | p(2,3) = dble(-amax/2) |
| 257 | p(1,4) = dble( amax/2) |
| 258 | p(2,4) = dble( amax/2) |
| 259 | p(1,5) = 0.D0 |
| 260 | p(2,5) = 0.D0 |
| 261 | CALL pyinit ('5MOM','gamma','gamma',0.0D0) |
| 262 | * ----- end of PYTHIA initialization ----- |
| 263 | * |
| 264 | ELSE IF (iproc .EQ. 2) THEN |
| 265 | * --- Cross section for narrow Resonance production --- |
| 266 | totxsc = dldmx(xmres,ny) ! dL/dMx in GeV^-1 |
| 267 | njres = kf_onium / 10 |
| 268 | njres = kf_onium - 10 * njres |
| 269 | xnorm = 4. * pi**2 * njres * xggres * gev2nb / xmres**2 |
| 270 | xstot = xnorm * totxsc |
| 271 | xstote = 0. |
| 272 | * |
| 273 | ELSE IF (iproc .EQ. 3) THEN |
| 274 | qf = abs(kchg(kc,1))/3. |
| 275 | nc = abs(kchg(kc,2))*2. + 1. |
| 276 | ELSE IF (iproc .EQ. 5) THEN |
| 277 | * Chargino storing in JETSET |
| 278 | kc = pycomp (41) |
| 279 | chaf(kc,1) = 'chi_1' |
| 280 | kchg(kc,1) = -3 |
| 281 | kchg(kc,2) = 0 |
| 282 | kchg(kc,3) = 1 |
| 283 | pmas(kc,1) = amprt(1) |
| 284 | qf = 1. |
| 285 | nc = 1. |
| 286 | * Neutralino storing in JETSET |
| 287 | kc = pycomp (42) |
| 288 | chaf(kc,1) = 'chi_1' |
| 289 | kchg(kc,1) = 0 |
| 290 | kchg(kc,2) = 0 |
| 291 | kchg(kc,3) = 1 |
| 292 | pmas(kc,1) = amprt(2) |
| 293 | * |
| 294 | * Gaugino mixing parameters |
| 295 | u11 = -zmixss(4,1)*sin(gammar)/sqrt(2.) + |
| 296 | & zmixss(2,1)*cos(gammar) |
| 297 | v11 = zmixss(3,1)*sin(gammal)/sqrt(2.) + |
| 298 | & zmixss(2,1)*cos(gammal) |
| 299 | DO idc = 1,nssmod |
| 300 | IF (issmod(idc) .EQ. 39) THEN |
| 301 | IF (moddcy .EQ. 1) THEN |
| 302 | DO jmod = 1,5 |
| 303 | IF (jssmod(jmod,idc) .EQ. -12) THEN |
| 304 | xsbra = bssmod(idc) |
| 305 | xsbra = xsbra*xsbra*4 |
| 306 | GOTO 10 |
| 307 | END IF |
| 308 | END DO |
| 309 | ELSE |
| 310 | DO jmod = 1,5 |
| 311 | IF (jssmod(jmod,idc) .EQ. -12) THEN ! e nu_e |
| 312 | ssbr(1) = bssmod(idc) |
| 313 | GOTO 20 |
| 314 | ELSE IF (jssmod(jmod,idc) .EQ. -14) THEN ! mu nu_mu |
| 315 | ssbr(2) = bssmod(idc) |
| 316 | GOTO 20 |
| 317 | ELSE IF (jssmod(jmod,idc) .EQ. -16) THEN ! tau nu_tau |
| 318 | ssbr(3) = bssmod(idc) |
| 319 | GOTO 20 |
| 320 | ELSE IF (jssmod(jmod,idc) .EQ. -2) THEN ! u dbar |
| 321 | ssbr(4) = bssmod(idc) |
| 322 | GOTO 20 |
| 323 | ELSE IF (jssmod(jmod,idc) .EQ. 4) THEN ! c sbar |
| 324 | ssbr(5) = bssmod(idc) |
| 325 | GOTO 20 |
| 326 | END IF |
| 327 | END DO |
| 328 | 20 CONTINUE |
| 329 | END IF |
| 330 | END IF |
| 331 | END DO |
| 332 | * |
| 333 | 10 CONTINUE |
| 334 | ELSE IF (iproc .EQ. 6) THEN |
| 335 | * Fix parameters of process gg -> V1V2 |
| 336 | IF (kv1 .GT. kv2) THEN |
| 337 | kvtmp = kv1 |
| 338 | kv1 = kv2 |
| 339 | kv2 = kvtmp |
| 340 | END IF |
| 341 | IF (kv1.EQ.113 .AND. kv2.EQ.113) THEN |
| 342 | idx_vv = 1 |
| 343 | ELSE IF (kv1.EQ.113 .AND. kv2.EQ.223) THEN |
| 344 | idx_vv = 2 |
| 345 | ELSE IF (kv1.EQ.113 .AND. kv2.EQ.333) THEN |
| 346 | idx_vv = 3 |
| 347 | ELSE IF (kv1.EQ.113 .AND. kv2.EQ.443) THEN |
| 348 | idx_vv = 4 |
| 349 | ELSE IF (kv1.EQ.223 .AND. kv2.EQ.223) THEN |
| 350 | idx_vv = 5 |
| 351 | ELSE IF (kv1.EQ.223 .AND. kv2.EQ.333) THEN |
| 352 | idx_vv = 6 |
| 353 | ELSE IF (kv1.EQ.223 .AND. kv2.EQ.443) THEN |
| 354 | idx_vv = 7 |
| 355 | ELSE IF (kv1.EQ.333 .AND. kv2.EQ.333) THEN |
| 356 | idx_vv = 8 |
| 357 | ELSE IF (kv1.EQ.333 .AND. kv2.EQ.443) THEN |
| 358 | idx_vv = 9 |
| 359 | ELSE IF (kv1.EQ.443 .AND. kv2.EQ.443) THEN |
| 360 | idx_vv = 10 |
| 361 | ELSE |
| 362 | WRITE (6,*) ' %GGINIT: unknown (V1 V2) combination' |
| 363 | WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2 |
| 364 | STOP |
| 365 | END IF |
| 366 | DO 101 igvpar=1,4 |
| 367 | 101 gvpar(igvpar) = gvconst(igvpar,idx_vv) |
| 368 | IF (gvpar(1) .EQ. 0.0) THEN |
| 369 | WRITE (6,*) ' %GGINIT: (V1 V2) combination not implemented' |
| 370 | WRITE (6,'(3x,''KV1 = '',i3,'', KV2 = '',i3)') kv1, kv2 |
| 371 | STOP |
| 372 | END IF |
| 373 | END IF |
| 374 | * |
| 375 | IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR. |
| 376 | & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN |
| 377 | * Max cross section estimation |
| 378 | xsmax0 = 0. |
| 379 | DO nest = 1,100 |
| 380 | xm = amin + ggrnd(0) * (amax-amin) |
| 381 | wsq = xm*xm |
| 382 | CALL ggxsec |
| 383 | END DO |
| 384 | xssum = 0. |
| 385 | ntry = 0 |
| 386 | END IF |
| 387 | * |
| 388 | RETURN |
| 389 | END |
| 390 | * |
| 391 | *======================================================================= |
| 392 | CDECK ID>, GGWID. |
| 393 | FUNCTION ggwid (kf_res) |
| 394 | *----------------------------------------------------------------------- |
| 395 | * Routine to define 2-gamma width of a resonance |
| 396 | * with JETSET code KF_RES |
| 397 | * Author: Yu.Kharlov |
| 398 | *----------------------------------------------------------------------- |
| 399 | REAL wid(5) |
| 400 | DATA wid /9.E-09, 1.E-06, 5.E-06, 6.3E-06, 0.41E-06/ |
| 401 | IF (kf_res .EQ. 111) THEN |
| 402 | ggwid = wid(1) |
| 403 | ELSE IF (kf_res .EQ. 221) THEN |
| 404 | ggwid = wid(2) |
| 405 | ELSE IF (kf_res .EQ. 331) THEN |
| 406 | ggwid = wid(3) |
| 407 | ELSE IF (kf_res .EQ. 441 |
| 408 | & .OR. kf_res .EQ. 10441 |
| 409 | & .OR. kf_res .EQ. 445) THEN |
| 410 | ggwid = wid(4) |
| 411 | ELSE IF (kf_res .EQ. 551 |
| 412 | & .OR. kf_res .EQ. 10551 |
| 413 | & .OR. kf_res .EQ. 555) THEN |
| 414 | ggwid = wid(5) |
| 415 | END IF |
| 416 | RETURN |
| 417 | END |
| 418 | * |
| 419 | *======================================================================= |
| 420 | CDECK ID>, GGRUN. |
| 421 | SUBROUTINE ggrun |
| 422 | *----------------------------------------------------------------------- |
| 423 | * Single event generation |
| 424 | * Author: Yu.Kharlov |
| 425 | *----------------------------------------------------------------------- |
| 426 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 427 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 428 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 429 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 430 | CHARACTER lumfil*80 |
| 431 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 432 | & gvconst(4,10) |
| 433 | REAL nc |
| 434 | COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5), |
| 435 | & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5) |
| 436 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 437 | & xstote, ssbr(10) |
| 438 | COMMON /ggkin/ xkin(10) |
| 439 | COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) |
| 440 | DOUBLE PRECISION P,V |
| 441 | SAVE /PYJETS/ |
| 442 | COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) |
| 443 | DOUBLE PRECISION PARU,PARJ |
| 444 | SAVE /PYDAT1/ |
| 445 | COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) |
| 446 | DOUBLE PRECISION PMAS,PARF,VCKM |
| 447 | SAVE /PYDAT2/ |
| 82764a9c |
448 | REAL pgam1(4),pgam2(4) |
| 0b5dd071 |
449 | REAL p3(4), p4(4) |
| 450 | REAL*8 dbetaz |
| 451 | REAL*8 ggrnd |
| 452 | * |
| 453 | 1 CALL gg2gam (amas,wsq,ygg,q1sq,q2sq, |
| 454 | & pgam1,pgam2,ptag1,ptag2,weight) |
| 455 | * |
| 456 | * Cross section calculation |
| 457 | IF (iproc.EQ.1 .OR. iproc.EQ.3 .OR. |
| 458 | & iproc.EQ.4 .OR. iproc.EQ.5 .OR. iproc.EQ.6) THEN |
| 459 | CALL ggxsec |
| 460 | xssum = xssum + xscur |
| 461 | ntry = ntry + 1 |
| 462 | * Select event on cross section of process gg -> XX |
| 463 | IF (xscur0/xsmax0 .LT. sngl (ggrnd(0))) GOTO 1 |
| 464 | * Total cross section calculation due to current statistics |
| 465 | xstot = xssum/ntry |
| 466 | xstote = xstot / sqrt(float(ntry)) |
| 467 | END IF |
| 468 | * |
| 469 | xmgg = sqrt (wsq) |
| 470 | xmg1 = sqrt (q1sq) |
| 471 | xmg2 = sqrt (q2sq) |
| 472 | dbetaz = dtanh(1.d0*ygg) |
| 473 | * |
| 474 | * calculation of the 4-vector of two-photon system |
| 475 | * and filling arrays pgg(1:2,1:5), kgg(1:2) |
| 476 | DO ip=1,4 |
| 477 | p2g(ip) = pgam1(ip) + pgam2(ip) |
| 478 | pgg(1,ip) = pgam1(ip) |
| 479 | pgg(2,ip) = pgam2(ip) |
| 480 | END DO |
| 481 | p2g(5) = xmgg |
| 482 | pgg(1,5) = -xmg1 |
| 483 | pgg(2,5) = -xmg2 |
| 484 | kgg(1) = 22 |
| 485 | kgg(2) = 22 |
| 486 | ngg = 2 |
| 487 | * |
| 488 | * =======> PROCESS # 1 |
| 489 | IF (iproc .EQ. 1) THEN |
| 490 | DO ip = 1,4 |
| 491 | p(1,ip) = dble (pgam1(ip)) |
| 492 | p(2,ip) = dble (pgam2(ip)) |
| 493 | END DO |
| 494 | p(1,5) = dble (-xmg1) |
| 495 | p(2,5) = dble (-xmg2) |
| 496 | CALL pyevnt |
| 497 | * |
| 498 | * =======> PROCESS # 2 |
| 499 | ELSE IF (iproc .EQ. 2) THEN |
| 500 | ngg = ngg + 1 |
| 501 | DO ip = 1,5 |
| 502 | pgg(ngg,ip) = p2g(ip) |
| 503 | END DO |
| 504 | kgg(ngg) = kf_onium |
| 505 | * Fill the JETSET common block |
| 506 | CALL gg2py (ngg) |
| 507 | CALL pyexec |
| 508 | * |
| 509 | * =======> PROCESS # 3, # 4, # 5, # 6 |
| 510 | ELSE IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. |
| 511 | & iproc.EQ.5 .OR. iproc.EQ.6) THEN |
| 512 | phi = 2*pi*ggrnd(0) |
| 513 | IF (iproc.EQ.3 .OR. iproc.EQ.4 .OR. iproc.EQ.5) THEN |
| 514 | e3 = 0.5*sqrt(wsq) |
| 515 | e4 = e3 |
| 516 | beta = xkin(1) |
| 517 | costhe = xkin(2) |
| 518 | am1 = xkin(3) |
| 519 | am2 = xkin(3) |
| 520 | p3abs = e3 * beta |
| 521 | ELSE IF (iproc.EQ.6) THEN |
| 522 | am1 = xkin(1) |
| 523 | am2 = xkin(2) |
| 524 | costhe = xkin(3) |
| 525 | p3abs = xkin(4) |
| 526 | e3 = xkin(5) |
| 527 | e4 = xmgg - e3 |
| 528 | END IF |
| 529 | sinthe = sqrt (1 - costhe*costhe) |
| 530 | p3(1) = p3abs * cos(phi) * sinthe |
| 531 | p3(2) = p3abs * sin(phi) * sinthe |
| 532 | p3(3) = p3abs * costhe |
| 533 | p3(4) = e3 |
| 534 | p4(4) = e4 |
| 535 | DO ip = 1,3 |
| 536 | p4(ip) = -p3(ip) |
| 537 | END DO |
| 538 | igg1 = ngg + 1 |
| 539 | igg2 = ngg + 2 |
| 540 | DO ip = 1,4 |
| 541 | pgg(igg1,ip) = p3(ip) |
| 542 | pgg(igg2,ip) = p4(ip) |
| 543 | END DO |
| 544 | pgg(igg1,5) = am1 |
| 545 | pgg(igg2,5) = am2 |
| 546 | ngg = ngg + 2 |
| 547 | * |
| 548 | IF (iproc .EQ. 3) THEN |
| 549 | kgg(igg1) = kferm |
| 550 | kgg(igg2) = -kferm |
| 551 | CALL gg2py (ngg) |
| 552 | ELSE IF (iproc .EQ. 4) THEN |
| 553 | kgg(igg1) = 24 |
| 554 | kgg(igg2) = -24 |
| 555 | CALL gg2py (ngg) |
| 556 | ELSE IF (iproc .EQ. 5) THEN |
| 557 | kgg(igg1) = 41 |
| 558 | kgg(igg2) = -41 |
| 559 | * Decay of charginos into W+neutralino |
| 560 | CALL ggdecy (igg1) |
| 561 | CALL ggdecy (igg2) |
| 562 | * |
| 563 | CALL gg2py (ngg) |
| 564 | k(igg1,1) = 11 |
| 565 | k(igg2,1) = 11 |
| 566 | ELSE IF (iproc.EQ.6) THEN |
| 567 | kgg(igg1) = kv1 |
| 568 | kgg(igg2) = kv2 |
| 569 | CALL gg2py (ngg) |
| 570 | END IF |
| 571 | * |
| 572 | * Boost event into initial frame |
| 573 | mstu(33) = 1 |
| 574 | CALL pyrobo(3,ngg, 0.d0,0.d0, 0.d0,0.d0, dbetaz) |
| 575 | CALL pyexec |
| 576 | * |
| 577 | END IF |
| 578 | * |
| 579 | RETURN |
| 580 | END |
| 581 | * |
| 582 | *======================================================================= |
| 583 | CDECK ID>, GGEXIT. |
| 584 | SUBROUTINE ggexit |
| 585 | *----------------------------------------------------------------------- |
| 586 | * Cross section calculation and end of run and |
| 587 | * print out cross section and related variables |
| 588 | * Author: Yu.Kharlov |
| 589 | *----------------------------------------------------------------------- |
| 590 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 591 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 592 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 593 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 594 | CHARACTER lumfil*80 |
| 595 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 596 | & xstote, ssbr(10) |
| 597 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 598 | & gvconst(4,10) |
| 599 | REAL nc |
| 600 | * |
| d5692e9e |
601 | WRITE (6,'(/1x,79(''=''))') |
| 602 | WRITE (6,'(1x,''I'',36x,''TPHIC'',36x, ''I'')') |
| 603 | WRITE (6,'(1x,''I'',34x,''Process '',i1,34x,''I'')') iproc |
| 604 | WRITE (6,'(1x,''I'',34x,''---------'', 34x,''I'')') |
| 605 | WRITE (6,'(1x,''I'',77x, ''I'')') |
| 606 | WRITE (6,'(1x,''I'',20x,''Events thrown : '',i8,31x,''I'')') |
| 0b5dd071 |
607 | & ntry |
| d5692e9e |
608 | WRITE (6,'(1x,''I'',20x,''Events accepted : '',i8,31x,''I'')') |
| 0b5dd071 |
609 | & nevent |
| d5692e9e |
610 | WRITE (6,'(1x,''I'',20x,''Cross section : '',e10.3,'' +- '','// |
| 0b5dd071 |
611 | & 'e10.3,'' nb'',12x ''I'')') xstot, xstote |
| d5692e9e |
612 | WRITE (6,'(1x,''I'',77x, ''I'')') |
| 613 | WRITE (6,'(1x,79(''=''))') |
| 0b5dd071 |
614 | * |
| 615 | IF (iproc .EQ. 1) THEN |
| 616 | * --- Alternative cross section for minimum bias events --- |
| 617 | totxsc = 0.0 |
| 618 | hmas = (amax-amin)/nmas |
| 619 | f1 = dldmx(amin ,ny) * dchad(amin) |
| 620 | DO i=1,nmas |
| 621 | xm = amin + i*hmas |
| 622 | f2 = dldmx(xm-hmas/2.,ny) * dchad(xm-hmas/2.) |
| 623 | f3 = dldmx(xm ,ny) * dchad(xm) |
| 624 | totxsc = totxsc +(f1 + 4.*f2 + f3)/6. |
| 625 | f1 = f3 |
| 626 | END DO |
| 627 | sigtt = totxsc * hmas |
| 628 | sigmx = sigtt/sqrt(float(nevent)) |
| 629 | WRITE (6,*) ' CrosSec =',sigtt,' +/-',sigmx,' mb' |
| 630 | END IF |
| 631 | * |
| 632 | RETURN |
| 633 | END |
| 634 | * |
| 635 | REAL FUNCTION dchad(xm) |
| 636 | *----------------------------------------------------------------------- |
| 637 | * Total cross section for gamma-gamma collisions vs. Xm in mb |
| 638 | * Parametrization is given by G.A.Schuler, T.Sjostrand, |
| 639 | * CERN-TH.7193/94, formula (9). |
| 640 | * Coded by Yu.Kharlov 20.03.1995. |
| 641 | *----------------------------------------------------------------------- |
| 642 | DATA epsilon, eta /0.0808, 0.4525/ |
| 643 | dchad = 0.211e-3 * (xm**2)**epsilon + 0.297e-3/(xm**2)**eta |
| 644 | RETURN |
| 645 | END |
| 646 | * |
| 647 | *======================================================================= |
| 648 | CDECK ID>, GGXSEC. |
| 649 | SUBROUTINE ggxsec |
| 650 | *----------------------------------------------------------------------- |
| 651 | * Cross section calculation in current event for iproc=1,3,4,5,6 |
| 652 | * Author: Yu.Kharlov |
| 653 | *----------------------------------------------------------------------- |
| 654 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 655 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 656 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 657 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 658 | CHARACTER lumfil*80 |
| 659 | COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5), |
| 660 | & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5) |
| 661 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 662 | & gvconst(4,10) |
| 663 | REAL nc |
| 664 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 665 | & xstote, ssbr(10) |
| 666 | COMMON /ggkin/ xkin(10) |
| 667 | COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) |
| 668 | DOUBLE PRECISION PMAS,PARF,VCKM |
| 669 | SAVE /PYDAT2/ |
| 670 | PARAMETER (epsilon = 0.0808, eta = -0.4525) |
| 671 | REAL*8 ggrnd |
| 672 | REAL*8 pymass |
| 673 | * |
| 674 | IF (iproc .EQ. 1) THEN |
| 675 | * Parametrization by G.A.Schuler, T.Sjostrand, |
| 676 | * CERN-TH.7193/94, formula (9). |
| 677 | xscur0 = 211. * wsq**epsilon + 297. * wsq**eta |
| 678 | xscur = xscur0 * xlumint |
| 679 | ELSE IF (iproc .EQ. 3 .OR. iproc .EQ. 5) THEN |
| 680 | beta2 = 1.0 - 4.0*amprt(1)**2/wsq |
| 681 | beta2 = min (beta2, 1.0) |
| 682 | beta = sqrt (beta2) |
| 683 | IF (beta2 .LT. 0.9) THEN |
| 684 | 1 costhe = 2.*ggrnd(0) - 1. |
| 685 | costh2 = costhe*costhe |
| 686 | xnum = 1. + 2.*beta2 * (1.-beta2) * (1.-costh2) - |
| 687 | & (beta2*costh2)**2 |
| 688 | xden = (1 - beta2*costh2)**2 |
| 689 | ffcur = xnum / xden |
| 690 | ffmax = (1. + beta2) / (1. - beta2) |
| 691 | rff = ggrnd(0) |
| 692 | IF (ffcur .LE. rff*ffmax) GOTO 1 |
| 693 | ELSE |
| 694 | small = 1.0 - cos(thetamin) |
| 695 | 2 zz = exp (ggrnd(0) * alog(2./small)) |
| 696 | C sgn = sign (1., 2.*ggrnd(0)-1.) |
| 697 | sgn = 2.*ggrnd(0)-1. |
| 698 | sgn = sign (1., sgn) |
| 699 | costhe = sgn * (zz-1.) / (zz+1.) |
| 700 | ffcur = (1.+costhe*costhe) / (1.-beta2*costhe*costhe) |
| 701 | ffmax = 2. / (1.- costhe*costhe) |
| 702 | rff = ggrnd(0) |
| 703 | IF (ffcur .LE. rff*ffmax) GOTO 2 |
| 704 | END IF |
| 705 | * |
| 706 | xkin(1) = beta |
| 707 | xkin(2) = costhe |
| 708 | xkin(3) = amprt(1) |
| 709 | * |
| 710 | * Yu.Kharlov, 3.06.1995 |
| 711 | totfac = 4. * pi * alpha**2 * gev2nb * qf**4 * nc / wsq |
| 712 | IF (beta2 .LT. 0.99) THEN |
| 713 | xscur0 =((3.-beta2**2)/(2.*beta) * alog ((1.+beta)/(1.-beta))- |
| 714 | & 2. + beta2) * beta |
| 715 | ELSE |
| 716 | xscur0 = 2. * alog (sqrt(wsq)/amprt(1)) - 1.0 |
| 717 | END IF |
| 718 | xscur0 = totfac * xscur0 * xsbra |
| 719 | xscur = xscur0 * xlumint |
| 720 | ELSE IF (iproc .EQ. 4) THEN |
| 721 | xmw = amprt(4) |
| 722 | beta2 = 1. - 4.*xmw**2/wsq |
| 723 | beta = sqrt (beta2) |
| 724 | 3 costhe = 2.*ggrnd(0) - 1. |
| 725 | costh2 = costhe*costhe |
| 726 | xnum = 19. - 6.*beta2 * (1.-beta2) + |
| 727 | & 2.*beta2 * (8. - 3.*beta2) * costh2 + |
| 728 | & 3.*(beta2*costh2)**2 |
| 729 | xden = (1. - beta2*costh2)**2 |
| 730 | ffcur = xnum / xden |
| 731 | ffmax = (19. + 10.*beta2 + beta2**2) / (1. - beta2)**2 |
| 732 | rff = ggrnd(0) |
| 733 | IF (ffcur .LE. rff*ffmax) GOTO 3 |
| 734 | * |
| 735 | xkin(1) = beta |
| 736 | xkin(2) = costhe |
| 737 | xkin(3) = xmw |
| 738 | * |
| 739 | * Yu.Kharlov, 3.06.1995 |
| 740 | totfac = pi * alpha**2 * beta * gev2nb/ wsq |
| 741 | xscur0 = 2. * (22. - 9.*beta2 + 3.*beta2**2) / (1. - beta2) - |
| 742 | & 3. * (1.-beta2**2)/beta * alog ((1.+beta)/(1.-beta)) |
| 743 | xscur0 = totfac * xscur0 * xsbra |
| 744 | xscur = xscur0 * xlumint |
| 745 | ELSE IF (iproc .EQ. 6) THEN |
| 746 | xmgg = sqrt (wsq) |
| 747 | bb = gvpar(3) + gvpar(4)*alog(xmgg/50.0) |
| 748 | am1 = sngl (pymass(kv1)) |
| 749 | am2 = sngl (pymass(kv2)) |
| 750 | p1 = pfork (xmgg,am1,am2) |
| 751 | e1 = sqrt (p1*p1 + am1*am1) |
| 752 | tmin = am1*am1 - xmgg*(e1+p1*costhv1) |
| 753 | tmax = am1*am1 - xmgg*(e1-p1*costhv1) |
| 754 | fmin = exp (bb*tmin) |
| 755 | fmax = exp (bb*tmax) |
| 756 | IF (fmax .GT. 1.e-10) THEN |
| 757 | tvar = 1./bb * dlog (fmax - ggrnd(0)*(fmax-fmin)) |
| 758 | ELSE |
| 759 | tvar = tmax + 1./bb * dlog (1. - ggrnd(0)) |
| 760 | END IF |
| 761 | costhe = (e1*xmgg - am1*am1 + tvar) / xmgg / p1 |
| 762 | IF (costhe .GT. 1.0) costhe = 1.0 |
| 763 | IF (costhe .LT.-1.0) costhe = -1.0 |
| 764 | IF (ggrnd(0) .LT. 0.5) costhe = -costhe |
| 765 | * |
| 766 | xkin(1) = am1 |
| 767 | xkin(2) = am2 |
| 768 | xkin(3) = costhe |
| 769 | xkin(4) = p1 |
| 770 | xkin(5) = e1 |
| 771 | * |
| 772 | IF (kv1 .EQ. 443 .AND. kv2 .EQ. 443) THEN |
| 773 | factor = 1.0 / alog(xmgg/3.097) |
| 774 | ELSE |
| 775 | factor = 1.0 |
| 776 | END IF |
| 777 | xscur0 = factor * gvpar(1) * xmgg**gvpar(2) / bb * exp(bb*tmax) |
| 778 | xscur = xscur0 * xlumint |
| 779 | * |
| 780 | END IF |
| 781 | * |
| 782 | * Find max cross section of gg-process |
| 783 | * |
| 784 | IF (xscur0 .GT. xsmax0) xsmax0 = xscur0 |
| 785 | * |
| 786 | RETURN |
| 787 | END |
| 788 | * |
| 789 | *======================================================================= |
| 790 | CDECK ID>, GGDECY. |
| 791 | SUBROUTINE ggdecy (igg) |
| 792 | *----------------------------------------------------------------------- |
| 793 | * Routine to decay a particle number I in GGEVNT common block |
| 794 | * Input : IGG - particle number in common GGEVNT arrays |
| 795 | * Author : Yu.Kharlov |
| 796 | *----------------------------------------------------------------------- |
| 797 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 798 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 799 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 800 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 801 | CHARACTER lumfil*80 |
| 802 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 803 | & gvconst(4,10) |
| 804 | REAL nc |
| 805 | COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5), |
| 806 | & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5) |
| 807 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 808 | & xstote, ssbr(10) |
| 809 | COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) |
| 810 | DOUBLE PRECISION PMAS,PARF,VCKM |
| 811 | SAVE /PYDAT2/ |
| 812 | COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr, |
| 813 | & xmll, xmlr, xmnl, xtanb, xmha, xmu, |
| 814 | & xmt, xat, xmbr, xab, u11, v11 |
| 815 | REAL p1(4), pcm(4), p2(4), p3(4), p4(4), ptmp(3) |
| 816 | REAL o1(3,3), o2(3,3) |
| 817 | INTEGER pycomp |
| 818 | DATA ifl /0/ |
| 819 | REAL*8 ggrnd |
| 820 | * |
| 821 | DO ip = 1,4 |
| 822 | p1(ip) = pgg(igg,ip) |
| 823 | END DO |
| 824 | * |
| 825 | IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 3) THEN |
| 826 | * Decay chi_1^+- --> chi_1^0 W^+- in CMS, uniformly |
| 827 | eb = (amprt(1)**2 + amprt(2)**2 - amprt(4)**2) / (2 * amprt(1)) |
| 828 | ec = (amprt(1)**2 - amprt(2)**2 + amprt(4)**2) / (2 * amprt(1)) |
| 829 | pb = sqrt (eb*eb - amprt(2)*amprt(2)) |
| 830 | pc = pb |
| 831 | costhe = 2*ggrnd(0) - 1. |
| 832 | sinthe = sqrt (1 - costhe*costhe) |
| 833 | phi = 2*pi*ggrnd(0) |
| 834 | pcm(1) = pb * sinthe * cos(phi) |
| 835 | pcm(2) = pb * sinthe * sin(phi) |
| 836 | pcm(3) = pb * costhe |
| 837 | pcm(4) = eb |
| 838 | CALL lorenb (amprt(1),p1,pcm,p2) |
| 839 | DO ip = 1,3 |
| 840 | pcm(ip) = -pcm(ip) |
| 841 | END DO |
| 842 | pcm(ip) = ec |
| 843 | CALL lorenb (amprt(1),p1,pcm,p3) |
| 844 | DO ip = 1,4 |
| 845 | pgg(ngg+1,ip) = p2(ip) |
| 846 | pgg(ngg+2,ip) = p3(ip) |
| 847 | END DO |
| 848 | pgg(ngg+1,5) = amprt(2) |
| 849 | pgg(ngg+2,5) = amprt(4) |
| 850 | kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg)) |
| 851 | kgg(ngg+2) = 24 * sign(1.,1.*kgg(igg)) |
| 852 | ngg = ngg + 2 |
| 853 | * |
| 854 | ELSE IF (abs(kgg(igg)) .EQ. 41 .AND. iproc .EQ. 5) THEN |
| 855 | * Decay chi_1^+- --> e^+- nu_e chi_1^0 in CMS, |
| 856 | * according to the exact matrix element |
| 857 | * |
| 858 | xm23max = amprt(1) - amprt(2) |
| 859 | xmw = sngl (pmas(pycomp(24),1)) |
| 860 | gwt = sngl (pmas(pycomp(24),2)) |
| 861 | atmin = - atan (xmw/gwt) |
| 862 | atmax = atan ((xm23max**2 - xmw**2) / (xmw*gwt)) |
| 863 | am1 = amprt(1) |
| 864 | am4 = amprt(2) |
| 865 | * |
| 866 | 1 CONTINUE |
| 867 | * |
| 868 | * Random mass of (2+3) according to Breit-Wigner |
| 869 | * |
| 870 | xm23 = xmw**2 + xmw*gwt * |
| 871 | * tan (atmin + ggrnd(0)*(atmax - atmin)) |
| 872 | xm23 = max (0., xm23) |
| 873 | xm23 = sqrt (xm23) |
| 874 | IF (xm23 .GT. xm23max) GOTO 1 ! to avoid boundary effect |
| 875 | e23 = (am1**2 + xm23**2 - am4**2) / 2. / am1 |
| 876 | * |
| 877 | * Rest frame of (2+3) |
| 878 | * |
| 879 | phi = 2.*pi * ggrnd(0) |
| 880 | costhe = 2.*ggrnd(0) - 1. |
| 881 | sinthe = sqrt (1. - costhe*costhe) |
| 882 | * |
| 883 | pp = xm23 / 2. |
| 884 | p2(1) = pp * sinthe * cos(phi) |
| 885 | p2(2) = pp * sinthe * sin(phi) |
| 886 | p2(3) = pp * costhe |
| 887 | p2(4) = pp |
| 888 | * |
| 889 | DO ip = 1,3 |
| 890 | p3(ip) = -p2(ip) |
| 891 | END DO |
| 892 | p3(4) = pp |
| 893 | * |
| 894 | * Transformation from rest frame of (2+3) to rest frame of (1), |
| 895 | * 3-momenta p2, p3 are along z-axis |
| 896 | * |
| 897 | gamma = e23 / xm23 |
| 898 | gambe = sqrt (gamma*gamma - 1.) |
| 899 | * |
| 900 | p2z = gamma*p2(3) + gambe*p2(4) |
| 901 | e2 = gamma*p2(4) + gambe*p2(3) |
| 902 | p3z = gamma*p3(3) + gambe*p3(4) |
| 903 | e3 = gamma*p3(4) + gambe*p3(3) |
| 904 | * |
| 905 | p2(3) = p2z |
| 906 | p2(4) = e2 |
| 907 | p3(3) = p3z |
| 908 | p3(4) = e3 |
| 909 | * |
| 910 | * Arbitrary rotation of rest frame of (1) |
| 911 | * |
| 912 | phi = 2.*pi * ggrnd(0) |
| 913 | costhe = 2.*ggrnd(0) - 1. |
| 914 | sinthe = sqrt (1. - costhe*costhe) |
| 915 | cosphi = cos(phi) |
| 916 | sinphi = sin(phi) |
| 917 | * |
| 918 | o1(1,1) = costhe |
| 919 | o1(1,2) = 0. |
| 920 | o1(1,3) = -sinthe |
| 921 | o1(2,1) = 0. |
| 922 | o1(2,2) = 1. |
| 923 | o1(2,3) = 0. |
| 924 | o1(3,1) = sinthe |
| 925 | o1(3,2) = 0. |
| 926 | o1(3,3) = costhe |
| 927 | * |
| 928 | o2(1,1) = cosphi |
| 929 | o2(1,2) = -sinphi |
| 930 | o2(1,3) = 0. |
| 931 | o2(2,1) = sinphi |
| 932 | o2(2,2) = cosphi |
| 933 | o2(2,3) = 0. |
| 934 | o2(3,1) = 0. |
| 935 | o2(3,2) = 0. |
| 936 | o2(3,3) = 1. |
| 937 | * |
| 938 | DO ix = 1,3 |
| 939 | ptmp(ix) = 0. |
| 940 | DO jx = 1,3 |
| 941 | ptmp(ix) = ptmp(ix) + o1(ix,jx) * p2(jx) |
| 942 | END DO |
| 943 | END DO |
| 944 | DO ix = 1,3 |
| 945 | p2(ix) = ptmp(ix) |
| 946 | END DO |
| 947 | DO ix = 1,3 |
| 948 | ptmp(ix) = 0. |
| 949 | DO jx = 1,3 |
| 950 | ptmp(ix) = ptmp(ix) + o2(ix,jx) * p2(jx) |
| 951 | END DO |
| 952 | END DO |
| 953 | DO ix = 1,3 |
| 954 | p2(ix) = ptmp(ix) |
| 955 | END DO |
| 956 | * |
| 957 | DO ix = 1,3 |
| 958 | ptmp(ix) = 0. |
| 959 | DO jx = 1,3 |
| 960 | ptmp(ix) = ptmp(ix) + o1(ix,jx) * p3(jx) |
| 961 | END DO |
| 962 | END DO |
| 963 | DO ix = 1,3 |
| 964 | p3(ix) = ptmp(ix) |
| 965 | END DO |
| 966 | DO ix = 1,3 |
| 967 | ptmp(ix) = 0. |
| 968 | DO jx = 1,3 |
| 969 | ptmp(ix) = ptmp(ix) + o2(ix,jx) * p3(jx) |
| 970 | END DO |
| 971 | END DO |
| 972 | DO ix = 1,3 |
| 973 | p3(ix) = ptmp(ix) |
| 974 | END DO |
| 975 | * |
| 976 | p12 = am1 * p2(4) |
| 977 | * |
| 978 | e4 = am1 - e23 |
| 979 | pp =sqrt (e4*e4 - am4*am4) |
| 980 | p4(1) = pp * sinthe * cos(phi) |
| 981 | p4(2) = pp * sinthe * sin(phi) |
| 982 | p4(3) = -pp * costhe |
| 983 | p4(4) = e4 |
| 984 | * |
| 985 | width = (- 4.*p12**2*u11**2 - 4.*p12**2*v11**2 + |
| 986 | & 4.*p12*xm23**2*u11**2 - 2.*p12*am4**2*u11**2 - |
| 987 | & 2.*p12*am4**2*v11**2 + 2.*p12*am1**2*u11**2 + |
| 988 | & 2.*p12*am1**2*v11**2 - xm23**4*u11**2 + |
| 989 | & xm23**2*am4**2*u11**2 - 2.*xm23**2*am4*am1*u11*v11 - |
| 990 | & xm23**2*am1**2*u11**2) / |
| 991 | & ((xm23**2-xmw**2)**2 + (xmw*gwt)**2) |
| 992 | width = width * pp * xm23 |
| 993 | * |
| 994 | widmax = (2.*am1*xm23max * |
| 995 | & (2.*(xm23max*u11)**2 + am1**2*(u11**2 + v11**2)) + |
| 996 | & (xm23max*am4*u11)**2) * (am1**2 - am4**2) / 2./am1 * |
| 997 | & xm23 / ((xm23**2-xmw**2)**2 + (xmw*gwt)**2) |
| 998 | * |
| 999 | * Select current phase space configuration |
| 1000 | * |
| 1001 | IF (width .GT. widmax) WRITE (6,*) |
| 1002 | & '%GGDECY: Maximum chi+ width violated' |
| 1003 | rwid = widmax * ggrnd(0) |
| 1004 | IF (width .LE. rwid) GOTO 1 |
| 1005 | * |
| 1006 | * Pick up final lepton flavour |
| 1007 | rflav = ggrnd(0) |
| 1008 | IF (moddcy .EQ. 1) THEN |
| 1009 | IF (rflav .LE. 0.5) THEN |
| 1010 | lflav1 = 11 |
| 1011 | lflav2 = -12 |
| 1012 | ELSE |
| 1013 | lflav1 = 13 |
| 1014 | lflav2 = -14 |
| 1015 | END IF |
| 1016 | ELSE |
| 1017 | rr1 = ssbr(1) |
| 1018 | rr2 = ssbr(1)+ssbr(2) |
| 1019 | rr3 = ssbr(1)+ssbr(2)+ssbr(3) |
| 1020 | rr4 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4) |
| 1021 | rr5 = ssbr(1)+ssbr(2)+ssbr(3)+ssbr(4)+ssbr(5) |
| 1022 | IF (rflav .LE. rr1) THEN |
| 1023 | lflav1 = 11 |
| 1024 | lflav2 = -12 |
| 1025 | ELSE IF (rflav .LE. rr2) THEN |
| 1026 | lflav1 = 13 |
| 1027 | lflav2 = -14 |
| 1028 | ELSE IF (rflav .LE. rr3) THEN |
| 1029 | lflav1 = 15 |
| 1030 | lflav2 = -16 |
| 1031 | ELSE IF (rflav .LE. rr4) THEN |
| 1032 | lflav1 = -2 |
| 1033 | lflav2 = 1 |
| 1034 | ELSE IF (rflav .LE. rr5) THEN |
| 1035 | lflav1 = -4 |
| 1036 | lflav2 = 3 |
| 1037 | ELSE |
| 1038 | GOTO 1 |
| 1039 | END IF |
| 1040 | END IF |
| 1041 | * |
| 1042 | * Transformation to CMS of 2-gamma |
| 1043 | CALL lorenb (amprt(1),p1,p2,pcm) |
| 1044 | DO ip = 1,4 |
| 1045 | pgg(ngg+1,ip) = pcm(ip) |
| 1046 | END DO |
| 1047 | pgg(ngg+1,5) = 0. |
| 1048 | kgg(ngg+1) = lflav1 * sign(1.,1.*kgg(igg)) |
| 1049 | ngg = ngg + 1 |
| 1050 | * |
| 1051 | CALL lorenb (amprt(1),p1,p3,pcm) |
| 1052 | DO ip = 1,4 |
| 1053 | pgg(ngg+1,ip) = pcm(ip) |
| 1054 | END DO |
| 1055 | pgg(ngg+1,5) = 0. |
| 1056 | kgg(ngg+1) = lflav2 * sign(1.,1.*kgg(igg)) |
| 1057 | ngg = ngg + 1 |
| 1058 | * |
| 1059 | CALL lorenb (amprt(1),p1,p4,pcm) |
| 1060 | DO ip = 1,4 |
| 1061 | pgg(ngg+1,ip) = pcm(ip) |
| 1062 | END DO |
| 1063 | pgg(ngg+1,5) = amprt(2) |
| 1064 | kgg(ngg+1) = 42 * sign(1.,1.*kgg(igg)) |
| 1065 | ngg = ngg + 1 |
| 1066 | END IF |
| 1067 | * |
| 1068 | RETURN |
| 1069 | END |
| 1070 | * |
| 1071 | *======================================================================= |
| 1072 | CDECK ID>, GG2PY. |
| 1073 | SUBROUTINE gg2py (nngg) |
| 1074 | *----------------------------------------------------------------------- |
| 1075 | * Routine to fill PYTHIA event record by particles from GG |
| 1076 | * Input: NNGG : number of particles |
| 1077 | * Author: Yu.Kharlov |
| 1078 | *----------------------------------------------------------------------- |
| 1079 | COMMON /ggevnt/ nrun,ievent,wsq,ygg,xmg1,xmg2, p2g(5), |
| 1080 | & ptag1(4),ptag2(4), ngg, kgg(10),pgg(20,5) |
| 1081 | COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) |
| 1082 | DOUBLE PRECISION P,V |
| 1083 | SAVE /PYJETS/ |
| 1084 | DO igg = 1,nngg |
| 1085 | IF (igg .LE. 2) THEN |
| 1086 | k(igg,1) = 21 |
| 1087 | ELSE |
| 1088 | k(igg,1)=1 |
| 1089 | END IF |
| 1090 | k(igg,2) = kgg(igg) |
| 1091 | DO ip = 1,5 |
| 1092 | p(igg,ip) = pgg(igg,ip) |
| 1093 | END DO |
| 1094 | END DO |
| 1095 | n = nngg |
| 1096 | RETURN |
| 1097 | END |
| 1098 | * |
| 1099 | *======================================================================= |
| 1100 | CDECK ID>, PFORK. |
| 1101 | FUNCTION pfork (am0,am1,am2) |
| 1102 | *----------------------------------------------------------------------- |
| 1103 | * Routine to find a momentum of a secondary particle in the decay |
| 1104 | * of a particle with mass AM0 in the rest to 2 particles with |
| 1105 | * masses AM1 and AM2. |
| 1106 | *----------------------------------------------------------------------- |
| 1107 | denom = am0**4 + am1**4 + am2**4 - |
| 1108 | & 2.*(am0*am1)**2 - 2.*(am0*am2)**2 - 2.*(am1*am2)**2 |
| 1109 | IF (denom .GT. 0.) THEN |
| 1110 | pfork = sqrt (denom) / 2./am0 |
| 1111 | ELSE |
| 1112 | pfork = -1. |
| 1113 | END IF |
| 1114 | RETURN |
| 1115 | END |
| 1116 | *======================================================================= |
| 1117 | CDECK ID>, GG2GAM. |
| 1118 | SUBROUTINE gg2gam (amas,wsq,y3,qs1,qs2,pgam1,pgam2,ptag1,ptag2,wt) |
| 1119 | C*********************************************************************** |
| 1120 | C |
| 1121 | C Generator of the Direct Two Photon Processes in S.A.Sadovsky |
| 1122 | C coherent heavy ion collisions: A+B --> 1+2+3 17.01.1995 |
| 1123 | C with double differencial gg-luminosity function Yu.Kharlov |
| 1124 | C 20.03.1995 |
| 1125 | C |
| 1126 | C LUMFUN must be called first to initialise the generator |
| 1127 | C |
| 1128 | C*********************************************************************** |
| 1129 | C |
| 1130 | C DETAILED DESCRIPTION |
| 1131 | C ======== =========== |
| 1132 | C |
| 1133 | C INPUTS: only through common /TWOGENC/ |
| 1134 | C ======= |
| 1135 | C |
| 1136 | C OUTPUTS: AMAS The mass of ion. |
| 1137 | C ======== WSQ Square mass of the gamma-gamma system |
| 1138 | C Y3 The rapidity of the gamma-gamma system |
| 1139 | C QS1 Virtual mass squared of photon 1. |
| 1140 | C QS2 Virtual mass squared of photon 2. |
| 1141 | C PGAM1 Four-vector of photon 1 |
| 1142 | C PGAM2 Four-vector of photon 2 |
| 1143 | C PTAG1 Four-vector of tag 1 |
| 1144 | C PTAG2 Four-vector of tag 2 |
| 1145 | C WT The weight of the event (dummy so far) |
| 1146 | C |
| 1147 | C The remaining variables are internal to the program. |
| 1148 | C |
| 1149 | C*********************************************************************** |
| 1150 | INTEGER icnter,ncntrs,nevts,lout,debug |
| 1151 | PARAMETER (ncntrs =7) |
| 1152 | COMMON /twgenc/s,eb,wmin,wmax,th1pmn,th1pmx,th2pmn,th2pmx, |
| 1153 | + st12mn,st12mx,st12rt,st22mn,st22mx,st22rt,ommin,omrat,eps, |
| 1154 | + wghtsm,wghtgt,wghtmx,wtmxfd,xmres,xgres, |
| 1155 | + icnter(ncntrs),nevts,lout,lres,lwate,debug |
| 1156 | SAVE /twgenc/ |
| 1157 | REAL st12mn,st12mx,st12rt,st22mn,st22mx,st22rt |
| 1158 | REAL wmin,wmax,ommin,omrat,eps,th1pmn,th1pmx |
| 1159 | REAL th2pmn,th2pmx,s,eb,wghtsm,wghtmx,wghtgt,wtmxfd |
| 1160 | REAL xmres,xgres |
| 1161 | LOGICAL lres,lwate |
| 1162 | C |
| 82764a9c |
1163 | REAL alpha,pi |
| 0b5dd071 |
1164 | PARAMETER (alpha=1./137.036,pi=3.14159265) |
| 1165 | C |
| 1166 | REAL ptag1(10),ptag2(10),pgam1(10),pgam2(10),wt |
| 82764a9c |
1167 | REAL amas,wsq |
| 0b5dd071 |
1168 | C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 82764a9c |
1169 | REAL*4 y3,xm,q0,qs,qs1,qs2 |
| 1170 | REAL*8 pa(5),pb(5) |
| 0b5dd071 |
1171 | C |
| 82764a9c |
1172 | REAL*8 pi2 |
| 0b5dd071 |
1173 | C |
| 1174 | LOGICAL lstr |
| 1175 | DATA q0,pi2,lstr / 0.060, 6.28318530718 d0, .false. / |
| 1176 | REAL*8 ggrnd |
| 1177 | C |
| 1178 | 999 CONTINUE |
| 1179 | icnter(1) = icnter(1) + 1 |
| 1180 | C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 1181 | IF (lstr) go to 10 |
| 1182 | lstr =.true. |
| 1183 | qs = q0/sqrt(2.) |
| 1184 | C |
| 1185 | pa(1)= 0D0 |
| 1186 | pa(2)= 0D0 |
| 1187 | pa(3)= dsqrt(1D0*(eb*eb-amas*amas)) |
| 1188 | pa(4)= eb*1D0 |
| 1189 | pa(5)= amas*1D0 |
| 1190 | C |
| 1191 | pb(1)= 0D0 |
| 1192 | pb(2)= 0D0 |
| 1193 | pb(3)=-pa(3) |
| 1194 | pb(4)= pa(4) |
| 1195 | pb(5)= amas*1D0 |
| 1196 | C |
| 1197 | 10 xm = xmres |
| 1198 | * Generate 2gamma rapidity y3 and mass xm |
| 1199 | CALL ranlum(y3,xm,wt) |
| 1200 | C |
| 1201 | wsq = xm*xm |
| 1202 | C |
| 1203 | * Energy and z-momentum of 2gamma system |
| 1204 | egg = xm*cosh(y3) |
| 1205 | pgg = xm*sinh(y3) |
| 1206 | * |
| 1207 | * Squared masses of the photons |
| 1208 | qs1 = - qs**2 * dlog(ggrnd(0)) ! is it true that squared mass |
| 1209 | qs2 = - qs**2 * dlog(ggrnd(0)) ! is distributed in exponent? |
| 1210 | xm1 = -qs1 |
| 1211 | xm2 = -qs2 |
| 1212 | * |
| 1213 | * Energies and z-momenta of the photons |
| 1214 | xfun12 = wsq + xm1 - xm2 |
| 1215 | * |
| 1216 | eg1 = (egg + pgg*sqrt(1 - 4*wsq*xm1/(xfun12*xfun12))) * |
| 1217 | & xfun12/(2*wsq) |
| 1218 | eg2 = egg - eg1 |
| 1219 | pg1 = sqrt (eg1**2 - xm1) |
| 1220 | pg2 = -sqrt (eg2**2 - xm2) |
| 1221 | * |
| 1222 | * 4-momenta of the photons |
| 1223 | DO ip=1,2 |
| 1224 | pgam1(ip) = 0. |
| 1225 | pgam2(ip) = 0. |
| 1226 | END DO |
| 1227 | pgam1(3) = pg1 |
| 1228 | pgam1(4) = eg1 |
| 1229 | pgam2(3) = pg2 |
| 1230 | pgam2(4) = eg2 |
| 1231 | * |
| 1232 | * 4-momenta of the recoil ions |
| 1233 | DO ip=1,4 |
| 1234 | ptag1(ip) = pa(ip) - pgam1(ip) |
| 1235 | ptag2(ip) = pb(ip) - pgam2(ip) |
| 1236 | END DO |
| 1237 | * |
| 1238 | RETURN |
| 1239 | END |
| 1240 | * |
| 1241 | *======================================================================= |
| 1242 | CDECK ID>, RANLUM. |
| 1243 | SUBROUTINE ranlum (ry,rm,wlum) |
| 1244 | C G.V.Khaoustov |
| 1245 | C 12.12.1994 |
| 1246 | C Corrected: Yu.Kharlov |
| 1247 | C 23.05.1995 |
| 1248 | C 24-JUN-1997 |
| 1249 | C |
| 1250 | C generate two random numbers (ry and rm) with probability proportional |
| 1251 | C to luminosity function in ymin-ymax rapidity and xmin-xmax mass range |
| 1252 | C |
| 1253 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 1254 | & gvconst(4,10) |
| 1255 | REAL nc |
| 1256 | PARAMETER (isize=100000) |
| 1257 | COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny, |
| 1258 | + xmin,xmax,amn,amx,nnms,alfun(isize) |
| 1259 | PARAMETER (NCNTRS =7) |
| 1260 | COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX, |
| 1261 | + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS, |
| 1262 | + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES, |
| 1263 | + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG |
| 1264 | LOGICAL lres,lwate |
| 1265 | REAL*4 wlum |
| 1266 | DATA ifl/0/, stpy/0./, stpm/0./ |
| 1267 | REAL*8 ggrnd |
| 1268 | C |
| 1269 | IF (ifl.EQ.0) THEN |
| 1270 | IF (ymin.LT.ymn .OR. ymin.GT.ymx .OR. |
| 1271 | * ymax.LT.ymn .OR. ymax.GT.ymx .OR. |
| 1272 | * xmin.LT.amn .OR. xmin.GT.amx .OR. |
| 1273 | * xmax.LT.amn .OR. xmax.GT.amx) THEN |
| 1274 | WRITE (6,*) 'Luminosity function is not defined in this range' |
| 1275 | STOP |
| 1276 | ENDIF |
| 1277 | aky = ymax-ymin |
| 1278 | akm = xmax-xmin |
| 1279 | stpy = (ymx-ymn)/nny |
| 1280 | stpm = (amx-amn)/nnms |
| 1281 | rmin = exp(cexp*(xmax)) |
| 1282 | rmax = exp(cexp*(xmin)) |
| 1283 | dr = rmax-rmin |
| 1284 | ifl = 1 |
| 1285 | ENDIF |
| 1286 | C |
| 1287 | 1 CONTINUE |
| 1288 | IF (.NOT.lres) THEN |
| 1289 | rm = ggrnd(0) ! random mass |
| 1290 | rm = rmin + rm*dr |
| 1291 | rm = alog(rm)/cexp ! exponetial law |
| 1292 | ELSE |
| 1293 | wsq = xmres**2 + xmres*xgres * tan(pi*(ggrnd(0)-0.5)) |
| 1294 | IF (wsq .LT. xmin*xmin) GOTO 1 |
| 1295 | IF (wsq .GT. xmax*xmax) GOTO 1 |
| 1296 | rm = sqrt (wsq) |
| 1297 | END IF |
| 1298 | C |
| 1299 | ry = ggrnd(0) ! random Y |
| 1300 | * |
| 1301 | * Find y-limits for generated gg-mass: y_max = log(E_max/M_gg) |
| 1302 | * |
| 1303 | ym_max = alog (egg_max/rm) |
| 1304 | ygg_max = min (ymax, ym_max) |
| 1305 | ygg_min = max (ymin,-ym_max) |
| 1306 | ry = ygg_min + ry*(ygg_max-ygg_min) |
| 1307 | iy = (ry-ymn)/stpy+1. |
| 1308 | C |
| 1309 | im = (rm-amn)/stpm |
| 1310 | r = anorm*exp(cexp*rm)*ggrnd(0) ! normalization |
| 1311 | p = (ry-(iy-1)*stpy-ymn)/stpy ! four point interpolation |
| 1312 | q = (rm-im*stpm-amn)/stpm |
| 1313 | k = iy+im*(nny+1) |
| 1314 | f = (1.-p)*(1.-q)*alfun(k) + p*(1.-q)*alfun(k+1) + |
| 1315 | * q*(1.-p)*alfun(k+nny+1) + q*p*alfun(k+2+nny) |
| 1316 | IF (r .GT. f) GOTO 1 |
| 1317 | wlum = 1000.*f ! differential luminosity in GeV^-1 |
| 1318 | * |
| 1319 | RETURN |
| 1320 | END |
| 1321 | * |
| 1322 | *======================================================================= |
| 1323 | CDECK ID>, GGDATA. |
| 1324 | BLOCK DATA ggdata |
| 1325 | *----------------------------------------------------------------------- |
| 1326 | * Some useful constants and masses of MSSM particles |
| 1327 | * Author: Yu.Kharlov |
| 1328 | *----------------------------------------------------------------------- |
| 1329 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 1330 | & gvconst(4,10) |
| 1331 | REAL nc |
| 1332 | COMMON /ggini/ iproc,nevent,ilumf,lumfil,ebmn,eb,iz,ia,amas, |
| 1333 | & amin,amax,ymin,ymax,nmas,ny, kferm, |
| 1334 | & kf_onium,xmres,xgtres,xggres, xlumint, moddcy, |
| 1335 | & thetamin, costhv1, kv1,kv2,gvpar(4) |
| 1336 | CHARACTER lumfil*80 |
| 1337 | COMMON /ggxs/ xsmax0, xscur0, xscur, xsbra, xssum, ntry, xstot, |
| 1338 | & xstote, ssbr(10) |
| 1339 | COMMON /ggmssm/ xm1, xm2, xmg, xms, xmtl, xmtr, |
| 1340 | & xmll, xmlr, xmnl, xtanb, xmha, xmu, |
| 1341 | & xmt, xat, xmbr, xab, u11, v11 |
| 1342 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 1343 | INTEGER lz,la |
| 1344 | REAL Rion,Gamma |
| 1345 | PARAMETER (al = 1./128.) |
| 1346 | DATA pi /3.14159276/, alpha /al/, gev2nb /389379./, |
| 1347 | & hbarc /0.197327/ |
| 1348 | DATA iproc, nevent, ilumf /1, 10000, -1/ |
| 1349 | DATA lumfil /'gglum.dat'/ |
| 1350 | DATA Rion /0.0/ |
| 1351 | DATA amprt /110.,10.,150.,-1.,0.106/ |
| 1352 | DATA thetamin /1.e-03/ |
| 1353 | DATA costhv1 /1.0/ |
| 1354 | DATA xsbra /1./ |
| 1355 | DATA ny, nmas /50, 50/ |
| 1356 | DATA xm1, xm2, xmg, xms, xmtl, xmtr, |
| 1357 | & xmll, xmlr, xmnl, xtanb, xmha, xmu, |
| 1358 | & xmt, xat, xmbr, xab |
| 1359 | & /30., 50., 250., 700., 600., 600., |
| 1360 | & 400., 400., 400., 4., 300., -300., |
| 1361 | & 174., 300., 300., 300./ |
| 1362 | DATA gvconst / 63.0, 0.22, 6.0, 1.0, ! rho rho |
| 1363 | & 14.0, 0.22, 6.0, 1.0, ! rho omega |
| 1364 | & 7.0, 0.22, 4.0, 1.0, ! rho phi |
| 1365 | & 0.05, 0.80, 2.5, 0.0, ! rho J/psi |
| 1366 | & 0.0, 0.00, 0.0, 0.0, ! omega omega |
| 1367 | & 0.8, 0.22, 4.0, 1.0, ! omega phi |
| 1368 | & 6.E-3, 0.80, 2.5, 0.0, ! omega J/psi |
| 1369 | & 0.2, 0.22, 2.0, 1.0, ! phi phi |
| 1370 | & 3.E-3, 0.80, 1.5, 0.0, ! phi J/psi |
| 1371 | & 6.E-5, 2.00, 1.5, 0.0/ ! J/psi J/psi |
| 1372 | END |
| 1373 | * |
| 1374 | *======================================================================= |
| 1375 | CDECK ID>, EVPRINT. |
| 1376 | SUBROUTINE evprint(Lun,Iev,P3,Ygg,ifl,Eb,ia,iz,amas,ptag1,ptag2) |
| 1377 | *----------------------------------------------------------------------- |
| 1378 | * Event output for internal use by GGHIC authors |
| 1379 | *----------------------------------------------------------------------- |
| 1380 | COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) |
| 1381 | DOUBLE PRECISION P,V |
| 1382 | SAVE /PYJETS/ |
| 1383 | REAL ptag1(4),ptag2(4),p3(5) |
| 1384 | C |
| 1385 | WRITE(lun,100) iev,n,Eb,p3(5),ygg,ifl,-ifl |
| 1386 | WRITE(lun,101) 90000+ia,iz,ptag1,amas |
| 1387 | WRITE(lun,101) 90000+ia,iz,ptag2,amas |
| 1388 | WRITE(lun,101) 90000+ 3, 0,(p3(j),j=1,5) |
| 1389 | IF (n.NE.0) THEN |
| 1390 | DO i=1,n |
| 1391 | ichg=pychge(k(i,2))/3 |
| 1392 | WRITE(lun,101) k(i,2),ichg,(p(i,j),j=1,5) |
| 1393 | ENDDO |
| 1394 | ENDIF |
| 1395 | RETURN |
| 1396 | 100 FORMAT(i6,i3,3e16.8,2x,2i8) |
| 1397 | 101 FORMAT(i6,i3,3(1x,e14.8),e13.8,e13.8) |
| 1398 | END |
| 1399 | * |
| 1400 | *======================================================================= |
| 1401 | C |
| 1402 | CDECK ID>, TWOINI. |
| 1403 | SUBROUTINE TWOINI(LOUTX,IPAR,XPAR,YPAR) TWOG0965 |
| 1404 | C***********************************************************************TWOG0966 |
| 1405 | C TWOG0967 |
| 1406 | C Initialises the gamma-gamma generator of Ion collisions TWOGAM. TWOG0968 |
| 1407 | C TWOG0969 |
| 1408 | C INPUTS: LOUTX Logical unit for print-out TWOG0970 |
| 1409 | C ======= XPAR(1) EBM Beam energy GeV TWOG0971 |
| 1410 | C XPAR(2) WMN Minimum two-gamma mass GeV TWOG0972 |
| 1411 | C XPAR(3) WMX Maximum two-gamma mass GeV TWOG0973 |
| 1412 | C TWOG0974 |
| 1413 | C XPAR(8) WGHTMX Maximum weight TWOG0978 |
| 1414 | C XPAR(9) XMRES Mass of resonance GeV TWOG0979 |
| 1415 | C XPAR(10) XGRES Total width of resonance GeV TWOG0980 |
| 1416 | C TWOG0981 |
| 1417 | C IPAR(1) 0 Unweighted events TWOG0982 |
| 1418 | C 1 Weighted events TWOG0983 |
| 1419 | C IPAR(2) 0 Continuum production TWOG0984 |
| 1420 | C 1 Resonance formation TWOG0985 |
| 1421 | C TWOG0986 |
| 1422 | C EXTERNALS: none TWOG0987 |
| 1423 | C ========== TWOG0988 |
| 1424 | C TWOG0989 |
| 1425 | C AUTHOR/DATE: S.A. Sadovsky, 17 January 1995 TWOG0990 |
| 1426 | C =========== TWOG0991 |
| 1427 | C TWOG0992 |
| 1428 | C***********************************************************************TWOG0993 |
| 1429 | C IMPLICIT NONE TWOG0994 |
| 1430 | INTEGER ICNTER,NCNTRS,NEVTS,LOUT,DEBUG TWOG0995 |
| 1431 | PARAMETER (NCNTRS =7) TWOG0996 |
| 1432 | COMMON /TWGENC/S,EB,WMIN,WMAX,TH1PMN,TH1PMX,TH2PMN,TH2PMX, TWOG0997 |
| 1433 | + ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT,OMMIN,OMRAT,EPS, TWOG0998 |
| 1434 | + WGHTSM,WGHTGT,WGHTMX,WTMXFD,XMRES,XGRES, TWOG0999 |
| 1435 | + ICNTER(NCNTRS),NEVTS,LOUT,LRES,LWATE,DEBUG TWOG1000 |
| 1436 | SAVE /TWGENC/ TWOG1001 |
| 1437 | REAL ST12MN,ST12MX,ST12RT,ST22MN,ST22MX,ST22RT TWOG1002 |
| 1438 | REAL WMIN,WMAX,OMMIN,OMRAT,EPS,TH1PMN,TH1PMX TWOG1003 |
| 1439 | REAL TH2PMN,TH2PMX,S,EB,WGHTSM,WGHTMX,WGHTGT,WTMXFD TWOG1004 |
| 1440 | REAL XMRES,XGRES TWOG1005 |
| 1441 | LOGICAL LRES,LWATE TWOG1006 |
| 1442 | C TWOG1007 |
| 1443 | REAL ALPHA,PI,TWOPI,XME,XMU TWOG1008 |
| 1444 | PARAMETER (ALPHA=1./137.036,PI=3.14159265) TWOG1009 |
| 1445 | PARAMETER (TWOPI=2.*PI,XME=0.000511,XMU=0.105658) TWOG1010 |
| 1446 | C TWOG1011 |
| 1447 | REAL XPAR(11),YPAR(6) TWOG1012 |
| 1448 | INTEGER LOUTX,IPAR(10),I TWOG1013 |
| 1449 | C TWOG1014 |
| 1450 | EB = XPAR(1) TWOG1015 |
| 1451 | WMIN = XPAR(2) TWOG1016 |
| 1452 | WMAX = XPAR(3) TWOG1017 |
| 1453 | WGHTMX= XPAR(8) TWOG1022 |
| 1454 | WTMXFD= 0 TWOG1023 |
| 1455 | XMRES = XPAR(9) TWOG1024 |
| 1456 | XGRES = XPAR(10) TWOG1025 |
| 1457 | LOUT = LOUTX TWOG1026 |
| 1458 | LWATE = IPAR(1) .EQ. 1 TWOG1027 |
| 1459 | LRES = IPAR(2) .EQ. 1 TWOG1028 |
| 1460 | EPS = 1./64.*(XME*WMIN**2/EB**3)**2/100. TWOG1029 |
| 1461 | C+ TWOG1030 |
| 1462 | C EPS is a small constant. We throw 1/2 COS(th/2)/(eps+SIN(th/2)) TWOG1031 |
| 1463 | C rather than 1/2 COT(th/2) to be able to start theta from zero. TWOG1032 |
| 1464 | C- TWOG1033 |
| 1465 | S = (2.*EB)**2 TWOG1054 |
| 1466 | OMMIN = WMIN**2/(4.*EB) TWOG1055 |
| 1467 | OMRAT = EB/OMMIN TWOG1056 |
| 1468 | C+ TWOG1057 |
| 1469 | C Print out initial conditions. TWOG1058 |
| 1470 | C- TWOG1059 |
| 1471 | WRITE(LOUT,1000) TWOG1060 |
| 1472 | 1000 FORMAT(1H ,10X,27('-'),' TWINIT ',27('-')) TWOG1061 |
| 1473 | WRITE(LOUT,1001) TWOG1062 |
| 1474 | 1001 FORMAT(1H ,10X,'I',60X,'I') TWOG1063 |
| 1475 | C |
| 1476 | WRITE(LOUT,1002)EB,YPAR(3),YPAR(4),YPAR(1),YPAR(2),WMIN,WMAX, |
| 1477 | + YPAR(5),YPAR(6) |
| 1478 | C TWOG1065 |
| 1479 | 1002 FORMAT TWOG1066 |
| 1480 | + (1H ,10X,'I Beam energy ',E15.5, TWOG1067 |
| 1481 | + ' GeV',11X,'I',/, TWOG1068 |
| 1482 | + 1H ,10X,'I Mass of ion ',F12.5, TWOG1067 |
| 1483 | + ' GeV',11X,'I',/, TWOG1068 |
| 1484 | + 1H ,10X,'I Gamma of ion ',F12.5, TWOG1067 |
| 1485 | + ' ',11X,'I',/, TWOG1068 |
| 1486 | + 1H ,10X,'I Charge of ion ',F12.5, TWOG1067 |
| 1487 | + ' ',11X,'I',/, TWOG1068 |
| 1488 | + 1H ,10X,'I Atomic number of ion ',F12.5, TWOG1067 |
| 1489 | + ' ',11X,'I',/, TWOG1068 |
| 1490 | + 1H ,10X,'I Minimum gamma-gamma mass ',F12.5, TWOG1069 |
| 1491 | + ' GeV',11X,'I',/, TWOG1070 |
| 1492 | + 1H ,10X,'I Maximum gamma-gamma mass ',F12.5, TWOG1071 |
| 1493 | + ' GeV',11X,'I',/, TWOG1080 |
| 1494 | + 1H ,10X,'I Minimum rapidity ',F12.5, TWOG1069 |
| 1495 | + ' ',11X,'I',/, TWOG1070 |
| 1496 | + 1H ,10X,'I Maximum rapidity ',F12.5, TWOG1071 |
| 1497 | + ' ',11X,'I' ) TWOG1080 |
| 1498 | IF (LRES) THEN TWOG1081 |
| 1499 | WRITE(LOUT,1001) TWOG1082 |
| 1500 | WRITE(LOUT,1003)XMRES,XGRES TWOG1083 |
| 1501 | ENDIF TWOG1084 |
| 1502 | 1003 FORMAT TWOG1085 |
| 1503 | + (1H ,10X,'I Resonance mass ',F12.5, TWOG1086 |
| 1504 | + ' GeV',11X,'I',/, TWOG1087 |
| 1505 | + 1H ,10X,'I Resonance total width ',F12.5, TWOG1088 |
| 1506 | + ' GeV',11X,'I') TWOG1089 |
| 1507 | C |
| 1508 | IF (LWATE) THEN TWOG1090 |
| 1509 | WRITE(LOUT,1001) TWOG1091 |
| 1510 | WRITE(LOUT,1004) TWOG1092 |
| 1511 | ELSE TWOG1095 |
| 1512 | WRITE(LOUT,1001) TWOG1096 |
| 1513 | WRITE(LOUT,1005)WGHTMX TWOG1097 |
| 1514 | 1004 FORMAT(1H ,10X,'I Produce weighted events ',27X,'I') TWOG1094 |
| 1515 | 1005 FORMAT(1H ,10X,'I Maximum weight ',F12.5, TWOG1099 |
| 1516 | + ' ',11X,'I') TWOG1100 |
| 1517 | ENDIF TWOG1101 |
| 1518 | C TWOG1102 |
| 1519 | WRITE(LOUT,1001) TWOG1103 |
| 1520 | WRITE(LOUT,1000) TWOG1104 |
| 1521 | C TWOG1105 |
| 1522 | DO 10 I=1,NCNTRS TWOG1106 |
| 1523 | 10 ICNTER(I)=0 TWOG1107 |
| 1524 | WGHTGT = 0. TWOG1108 |
| 1525 | WGHTSM = 0. TWOG1109 |
| 1526 | RETURN TWOG1110 |
| 1527 | END TWOG1111 |
| 1528 | C========================================================================= |
| 1529 | CDECK ID>, DLDMX. |
| 1530 | FUNCTION DLDMX(Rm,Ny) |
| 1531 | PARAMETER (isize=100000) |
| 1532 | C S.A.Sadovsky |
| 1533 | C Calculate the Luminosity integral over dY 5.02.1995 |
| 1534 | C in the rapidity range ymin-ymax |
| 1535 | C using luminosity interpolation table |
| 1536 | C Ny/2 - number of subintervals for Simpson integration |
| 1537 | C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 1538 | COMMON/lumfunct/anorm,cexp,ymin,ymax,ymn,ymx,nny, |
| 1539 | + xmin,xmax,amn,amx,nnms,alfun(isize) |
| 1540 | stpy=(ymx-ymn)/nny |
| 1541 | stpm=(amx-amn)/nnms ! Interpolation steps |
| 1542 | C |
| 1543 | im=(rm-amn)/stpm |
| 1544 | q =(rm-im*stpm-amn)/stpm |
| 1545 | C |
| 1546 | My = Ny/2 |
| 1547 | Hy =(ymax-ymin)/My ! Integration step over Y |
| 1548 | ry = ymin |
| 1549 | iy =(ry-ymn)/stpy+1. |
| 1550 | p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 1 |
| 1551 | k = iy+im*(nny+1) |
| 1552 | F1 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+ |
| 1553 | * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny) |
| 1554 | C |
| 1555 | SIM = 0.0 |
| 1556 | DO 10 i=1,My |
| 1557 | C |
| 1558 | ry = ymin + Hy*i |
| 1559 | iy =(ry-ymn)/stpy+1. |
| 1560 | p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 3 |
| 1561 | k = iy+im*(nny+1) |
| 1562 | F3 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+ |
| 1563 | * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny) |
| 1564 | C |
| 1565 | ry = ry - Hy*0.5 |
| 1566 | iy =(ry-ymn)/stpy+1. |
| 1567 | p =(ry-(iy-1)*stpy-ymn)/stpy ! Four point interpolation: 2 |
| 1568 | k = iy+im*(nny+1) |
| 1569 | F2 =(1.-p)*(1.-q)*alfun(k)+p*(1.-q)*alfun(k+1)+ |
| 1570 | * q*(1.-p)*alfun(k+nny+1)+q*p*alfun(k+2+nny) |
| 1571 | C |
| 1572 | SIM = SIM + (F1 + 4.*F2 + F3)/6. |
| 1573 | 10 F1 = F3 |
| 1574 | C |
| 1575 | DLDMX = 1000.*SIM*HY ! dL/dMx in GeV^-1 |
| 1576 | RETURN |
| 1577 | END |
| 1578 | *======================================================================= |
| 1579 | CDECK ID>, LUMFUN. |
| 1580 | SUBROUTINE LUMFUN(ifl,file,ymin,ymax,ny,xmn,xmx,nms) |
| 1581 | C----------------------------------------------------------------------- |
| 1582 | C preparation of luminosity function for generator |
| 1583 | C |
| 1584 | C ifl=-1 - calculate function + save in file 'file' |
| 1585 | C ifl= 0 - calculate function |
| 1586 | C ifl= 1 - read function from file |
| 1587 | C ymin,ymax,ny - rapidity range & number of points |
| 1588 | C xmn,xmx,nms - mass range & number of points |
| 1589 | C Attention: dimensions of alfun are alfun(ny+1,nms+1) |
| 1590 | C |
| 1591 | C Called by GGINIT |
| 1592 | C |
| 1593 | C Author: G.V.Khaoustov 12.12.1994 |
| 1594 | C Modified: Yu.Kharlov 18-JUN-1997 |
| 1595 | C----------------------------------------------------------------------- |
| 1596 | COMMON /ggpar/ pi,hbarc,gev2nb,alpha, amprt(5), qf,nc, egg_max, |
| 1597 | & gvconst(4,10) |
| 1598 | REAL nc |
| 1599 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 1600 | INTEGER lz,la |
| 1601 | REAL Rion,Gamma |
| 1602 | PARAMETER (isize=100000) |
| 1603 | CHARACTER*(*) file |
| 1604 | REAL work(1000) |
| 1605 | COMMON/lumfunct/anorm,cexp,umn,umx,ymn,ymx,nny, |
| 1606 | + wmn,wmx,amn,amx,nnms,alfun(isize) |
| 1607 | C |
| 1608 | umn = ymin ! parameters for luminosity generation |
| 1609 | umx = ymax |
| 1610 | wmn = xmn |
| 1611 | wmx = xmx |
| 1612 | C |
| 1613 | C Define R : radii of the ions (in fm) (if set non-zero, use it) |
| 1614 | IF (Rion .LE. 0) THEN |
| 1615 | rion = 1.2*float(la)**(1.0/3.0) |
| 1616 | ENDIF |
| 1617 | * |
| 1618 | * Fing max energy of gg-system (with safity factor 2) |
| 1619 | * |
| 1620 | egg_max = 2.*hbarc * gamma / rion |
| 1621 | * |
| 1622 | IF(ifl.EQ.1) THEN ! read lum. function from file |
| 1623 | OPEN(unit=10,err=1,file=file,status='old', |
| 1624 | * form='formatted') |
| 1625 | READ(10, *,err=3) anorm,cexp |
| 1626 | READ(10,99,err=3) ymn,ymx,nny,amn,amx,nnms |
| 1627 | READ(10,100,err=3) (alfun(I),i=1,(nny+1)*(nnms+1)) |
| 1628 | CLOSE(unit=10) |
| 1629 | PRINT *,'The Luminosity has been read from File ' |
| 1630 | RETURN |
| 1631 | ENDIF |
| 1632 | c |
| 1633 | IF (ny*nms .GT. isize) GOTO 2 |
| 1634 | PRINT *,'Start the Luminosity calculation' |
| 1635 | nny = ny |
| 1636 | ymn = ymin |
| 1637 | ymx = ymax |
| 1638 | nnms= nms |
| 1639 | amn = xmn |
| 1640 | amx = xmx |
| 1641 | c |
| 1642 | stepy=(ymx-ymn)/nny |
| 1643 | stepm=(amx-amn)/nnms |
| 1644 | * |
| 1645 | * Lum. function calculation |
| 1646 | * |
| 1647 | DO IM = 0,nnms |
| 1648 | am=amn+im*stepm |
| 1649 | * Find Y-limits for mass AM (Yu.Kharlov) |
| 1650 | ym_max = alog (egg_max/am) |
| 1651 | ym_max = max(ym_max,5.*stepy) |
| 1652 | ygg_max = min (ymax, ym_max) |
| 1653 | ygg_min = max (ymin,-ym_max) |
| 1654 | DO IY = 1,nny+1 |
| 1655 | Y = ymn+(iy-1)*stepy |
| 1656 | IF (Y.LT.ygg_max .AND. Y.GT.ygg_min) THEN |
| 1657 | dl = D2LDMDY(1000.*am,Y) ! Differential Luminosity in MeV^-1 |
| 1658 | ELSE |
| 1659 | dl = 0.0 |
| 1660 | END IF |
| 1661 | alfun(iy+(nny+1)*im) = DL |
| 1662 | ENDDO |
| 1663 | ENDDO |
| 1664 | c |
| 1665 | s=0. |
| 1666 | DO i=1,nnms+1 ! find max values for exp. parameters calcul. |
| 1667 | am=0. |
| 1668 | DO j=1,nny+1 |
| 1669 | IF(am.LT.alfun(j+(nny+1)*(i-1))) am=alfun(j+(nny+1)*(i-1)) |
| 1670 | ENDDO |
| 1671 | work(i)=am |
| 1672 | s=s+am |
| 1673 | ENDDO |
| 1674 | s=(s-work(nnms+1))*stepm |
| 1675 | c |
| 1676 | smax=0. ! calculation optimized parameters for exp. generator |
| 1677 | iflag=-1 |
| 1678 | DO 101 i=1,nnms |
| 1679 | DO 102 j=i+1,nnms+1 |
| 1680 | IF (work(j) .LE. 0.0) GOTO 102 |
| 1681 | c=alog(work(i)/work(j))/((i-j)*stepm) |
| 1682 | am=work(i)/exp(c*(i-1)*stepm) |
| 1683 | se=am/c*(exp(c*(amx-amn))-1.) |
| 1684 | IF(s/se.lt.smax) GOTO 102 |
| 1685 | DO k=1,nnms+1 |
| 1686 | amp=am*exp(c*(k-1)*stepm) |
| 1687 | IF((work(k)-amp)/amp.GT.0.00001) THEN ! Accuracy problem |
| 1688 | GOTO 102 |
| 1689 | ENDIF |
| 1690 | ENDDO |
| 1691 | smax =s/se |
| 1692 | anorm=am |
| 1693 | cexp =c |
| 1694 | iflag=1 |
| 1695 | 102 CONTINUE |
| 1696 | 101 CONTINUE |
| 1697 | IF(iflag.LT.0) THEN |
| 1698 | PRINT *,'Fault to find best exponential approximation' |
| 1699 | STOP |
| 1700 | ENDIF |
| 1701 | c |
| 1702 | IF(ifl.EQ.0) RETURN ! save lum. function |
| 1703 | OPEN(unit=10,err=1,file=file,status='unknown', |
| 1704 | * form='formatted') |
| 1705 | WRITE(10, *,err=3) anorm,cexp |
| 1706 | WRITE(10,99,err=3) ymn,ymx,nny,amn,amx,nnms |
| 1707 | WRITE(10,100,err=3) (alfun(i),i=1,(nny+1)*(nnms+1)) |
| 1708 | CLOSE(unit=10) |
| 1709 | PRINT *,'The Luminosity calculation is finished' |
| 1710 | RETURN |
| 1711 | 1 PRINT *,' file open error, file:',file |
| 1712 | STOP |
| 1713 | 2 PRINT *,' array "ALFUN" to small to contain luminosity function' |
| 1714 | STOP |
| 1715 | 3 PRINT *,' file read/writr error, file:',file |
| 1716 | STOP |
| 1717 | 100 FORMAT(5e13.5) |
| 1718 | 99 FORMAT(2e13.5,i5,2e13.5,i5) |
| 1719 | END |
| 1720 | C====================================================================== |
| 1721 | CDECK ID>, D2LDMDY. |
| 1722 | REAL FUNCTION D2LDMDY(M,Y) |
| 1723 | C----------------------------------------------------------------------- |
| 1724 | C Calculate the Luminosity-function d^2L/dMdY |
| 1725 | C where M is the invariant mass and Y the rapidity |
| 1726 | C of the photon-photon-system |
| 1727 | C |
| 1728 | C Calculations are based on the formulae of |
| 1729 | C N. Baron, G. Baur, Phys. Rev. C |
| 1730 | C C. Cahn, J. D. Jackson, Nucl. Phys. |
| 1731 | C |
| 1732 | C Kai Hencken, 29. 9. 1994 |
| 1733 | C |
| 1734 | C M : invariant mass (MeV) |
| 1735 | C Y : rapidity |
| 1736 | C |
| 1737 | C D2LDMDY : diff. Lum. (MeV**-1) |
| 1738 | C |
| 1739 | C Called by LUMFUN |
| 1740 | C----------------------------------------------------------------------- |
| 1741 | IMPLICIT NONE |
| 1742 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 1743 | INTEGER lz,la |
| 1744 | REAL Rion,Gamma |
| 1745 | REAL M,Y |
| 1746 | REAL RM,W1,W2 |
| 1747 | COMMON/D2LIParam/RM,W1,W2 |
| 1748 | |
| 1749 | REAL NClass,DeltaL |
| 1750 | EXTERNAL NClass,DeltaL |
| 1751 | |
| 1752 | C |
| 1753 | C physical constants hbar*c and 1/alpha |
| 1754 | C |
| 1755 | REAL hbarc,falphai |
| 1756 | PARAMETER (hbarc=197.327053,falphai=137.0359895) |
| 1757 | |
| 1758 | REAL LClass,DL |
| 1759 | C |
| 1760 | C calculate photon frequencies from M,Y |
| 1761 | C |
| 1762 | W1 = M/2.0 * EXP ( Y) |
| 1763 | W2 = M/2.0 * EXP (-Y) |
| 1764 | C |
| 1765 | C calculate Rho in MeV^-1 instead of fm |
| 1766 | C |
| 1767 | RM = Rion / hbarc |
| 1768 | C |
| 1769 | C calculate the "classical" value first. |
| 1770 | C |
| 1771 | |
| 1772 | LClass = NClass (W1,RM,Gamma) * NClass (W2,RM,Gamma) |
| 1773 | |
| 1774 | C |
| 1775 | C substract from this the Delta-value |
| 1776 | C |
| 1777 | DL = DeltaL () |
| 1778 | |
| 1779 | D2LDMDY = 2.0 / M * float(LZ)**4 / falphai**2 * (LClass - DL) |
| 1780 | |
| 1781 | RETURN |
| 1782 | END |
| 1783 | *======================================================================= |
| 1784 | CDECK ID>, NCLASS. |
| 1785 | REAL FUNCTION NClass (W, R, Gamma) |
| 1786 | C |
| 1787 | C the classical photon-number without impact parameter cutoff: |
| 1788 | C |
| 1789 | IMPLICIT NONE |
| 1790 | REAL W,R, Gamma |
| 1791 | |
| 1792 | REAL Pi |
| 1793 | PARAMETER (Pi = 3.141592) |
| 1794 | |
| 1795 | REAL BesselK0,BesselK1 |
| 1796 | EXTERNAL BesselK0,BesselK1 |
| 1797 | |
| 1798 | REAL Xi |
| 1799 | |
| 1800 | Xi = W*R/Gamma |
| 1801 | |
| 1802 | NClass = 2.0/PI * (Xi*BesselK0 (Xi)*BesselK1 (Xi) - |
| 1803 | & Xi**2/2.0 *(BesselK1(Xi)**2 - BesselK0(Xi)**2)) |
| 1804 | |
| 1805 | RETURN |
| 1806 | END |
| 1807 | *======================================================================= |
| 1808 | CDECK ID>, DELTAL. |
| 1809 | REAL FUNCTION DeltaL () |
| 1810 | C |
| 1811 | C the difference to the classical value |
| 1812 | C Called by d2LdMdY |
| 1813 | C |
| 1814 | IMPLICIT NONE |
| 1815 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 1816 | INTEGER lz,la |
| 1817 | REAL Rion,Gamma |
| 1818 | REAL RM,W1,W2 |
| 1819 | COMMON/D2LIParam/RM,W1,W2 |
| 1820 | |
| 1821 | DOUBLE PRECISION XGAUSS(126),WGAUSS(126) |
| 1822 | COMMON /GAUSSD/XGAUSS,WGAUSS |
| 1823 | |
| 1824 | INTEGER N,I |
| 1825 | REAL b1 |
| 1826 | REAL*8 t,tmin,tmax |
| 1827 | REAL*8 Sum,Int,Int2 |
| 1828 | |
| 1829 | REAL Itgb1 |
| 1830 | EXTERNAL Itgb1 |
| 1831 | |
| 1832 | REAL Pi |
| 1833 | REAL*8 Accur |
| 1834 | PARAMETER (Pi = 3.141592,Accur = 1D -2) |
| 1835 | |
| 1836 | C |
| 1837 | C integrate first over b1 |
| 1838 | C |
| 1839 | |
| 1840 | C |
| 1841 | C Loop incrementing the boundary |
| 1842 | C |
| 1843 | tmin = 0.00 |
| 1844 | tmax = 0.25 |
| 1845 | Sum = 0.00 |
| 1846 | |
| 1847 | 100 CONTINUE |
| 1848 | C |
| 1849 | C Loop for the Gauss integration |
| 1850 | C |
| 1851 | Int=0.0 |
| 1852 | DO N=1,6 |
| 1853 | Int2 = Int |
| 1854 | Int=0.0 |
| 1855 | DO I=2**N-1,2**(N+1)-2 |
| 1856 | t = (tmax-tmin)/2.0*XGAUSS(I) + (tmax+tmin)/2.0 |
| 1857 | b1 = RM * DEXP (t) |
| 1858 | Int = Int + WGAUSS(I) * Itgb1(b1) * b1**2 |
| 1859 | ENDDO |
| 1860 | Int = (tmax-tmin)/2.0*Int |
| 1861 | IF (Int .NE. 0.) THEN |
| 1862 | IF (DABS ((Int2-Int)/Int) .LT. Accur) GOTO 200 |
| 1863 | END IF |
| 1864 | ENDDO |
| 1865 | IF (int .LT. 1.D-20) THEN |
| 1866 | int = 0.0D0 |
| 1867 | ELSE |
| 1868 | WRITE(*,*) ' (b1) GAUSS MAY BE INACCURATE' |
| 1869 | END IF |
| 1870 | |
| 1871 | 200 CONTINUE |
| 1872 | Sum = Sum + Int |
| 1873 | IF (sum .NE. 0.0) THEN |
| 1874 | IF (ABS (Int2/Sum) .LT. Accur) GOTO 300 |
| 1875 | ELSE |
| 1876 | IF (b1 .GT. 1.E+06) GOTO 300 |
| 1877 | END IF |
| 1878 | tmin = tmax |
| 1879 | tmax = tmax + 0.5 |
| 1880 | GOTO 100 |
| 1881 | |
| 1882 | 300 CONTINUE |
| 1883 | |
| 1884 | DeltaL = 4.0*Pi * Sum |
| 1885 | RETURN |
| 1886 | END |
| 1887 | |
| 1888 | CDECK ID>, ITGB1. |
| 1889 | REAL FUNCTION Itgb1(b) |
| 1890 | C |
| 1891 | C integration then over b2 |
| 1892 | C Called by DeltaL |
| 1893 | C |
| 1894 | C Corrected by S.A.Sadovsky |
| 1895 | C 25.10.1994 |
| 1896 | IMPLICIT NONE |
| 1897 | c +seq,d2lpar. Not used Yu.Kh. 18-JUN-1997 |
| 1898 | REAL b |
| 1899 | REAL*8 b1,bmin,bmax |
| 1900 | REAL RM,W1,W2 |
| 1901 | COMMON/D2LIParam/RM,W1,W2 |
| 1902 | DOUBLE PRECISION XGAUSS(126),WGAUSS(126) |
| 1903 | COMMON /GAUSSD/XGAUSS,WGAUSS |
| 1904 | INTEGER N,I |
| 1905 | REAL*8 b2,Int,Int2 |
| 1906 | REAL*8 Itgb2 |
| 1907 | EXTERNAL Itgb2 |
| 1908 | REAL*8 Accur |
| 1909 | PARAMETER (Accur = 1.D-2) |
| 1910 | |
| 1911 | b1 = b |
| 1912 | bmin = b1 - 2.0*RM |
| 1913 | IF (RM .GT. bmin) THEN |
| 1914 | bmin = RM |
| 1915 | ENDIF |
| 1916 | bmax = b1 + 2.0 * RM |
| 1917 | Int = 0.0 |
| 1918 | DO N=1,6 |
| 1919 | Int2 = Int |
| 1920 | Int = 0.0 |
| 1921 | DO I=2**N-1,2**(N+1)-2 |
| 1922 | b2 = (bmax-bmin)/2.0*XGAUSS(I) + (bmax+bmin)/2.0 |
| 1923 | Int = Int + WGAUSS(I) * Itgb2(b1,b2) * b2 |
| 1924 | END DO |
| 1925 | Int = (bmax-bmin)/2.0*Int |
| 1926 | IF (Int.NE.0.) THEN |
| 1927 | IF (DABS((Int2 - Int)/Int) .LT. Accur) GOTO 100 |
| 1928 | ENDIF |
| 1929 | ENDDO |
| 1930 | IF (int .LT. 1.D-20) THEN |
| 1931 | int = 0.0D0 |
| 1932 | ELSE |
| 1933 | WRITE(*,*) ' (b2) GAUSS MAY BE INACCURATE, Itgb1=',Int |
| 1934 | END IF |
| 1935 | 100 CONTINUE |
| 1936 | Itgb1 = Int |
| 1937 | RETURN |
| 1938 | END |
| 1939 | |
| 1940 | CDECK ID>, ITGB2. |
| 1941 | REAL*8 FUNCTION Itgb2 (b1,b2) |
| 1942 | C |
| 1943 | C The function to be integrated over |
| 1944 | C Called by Itgb1 |
| 1945 | C |
| 1946 | C Corrected by S.A.Sadovsky |
| 1947 | C 25.10.1994 |
| 1948 | IMPLICIT NONE |
| 1949 | COMMON/D2LParam/lz,la,Rion,Gamma |
| 1950 | INTEGER lz,la |
| 1951 | REAL Rion,Gamma |
| 1952 | REAL*8 b1,b2,arg |
| 1953 | REAL RM,W1,W2 |
| 1954 | COMMON/D2LIParam/RM,W1,W2 |
| 1955 | |
| 1956 | REAL Nphoton |
| 1957 | EXTERNAL Nphoton |
| 1958 | |
| 1959 | arg = (b1**2 + b2**2 - 4.0 * RM**2) / (2.0*b1*b2) |
| 1960 | IF (arg .GT. 1.D0) arg = 1.D0 |
| 1961 | Itgb2 = Nphoton(W1,sngl(b1),Gamma) * |
| 1962 | & Nphoton(W2,sngl(b2),Gamma) * |
| 1963 | & DACOS (arg) |
| 1964 | |
| 1965 | RETURN |
| 1966 | END |
| 1967 | |
| 1968 | CDECK ID>, NPHOTON. |
| 1969 | REAL FUNCTION Nphoton (W,Rho,Gamma) |
| 1970 | C |
| 1971 | C The differential photonnumber for a nucleus |
| 1972 | C (without form factor) |
| 1973 | C |
| 1974 | IMPLICIT NONE |
| 1975 | REAL W,Rho,Gamma |
| 1976 | |
| 1977 | REAL BESSELK1 |
| 1978 | EXTERNAL BESSELK1 |
| 1979 | |
| 1980 | REAL Wphib,WGamma |
| 1981 | |
| 1982 | REAL PI |
| 1983 | PARAMETER (PI=3.141592) |
| 1984 | |
| 1985 | WGamma = W/Gamma |
| 1986 | C |
| 1987 | C without form factor |
| 1988 | C |
| 1989 | Wphib = WGamma*BESSELK1 (WGamma*Rho) |
| 1990 | |
| 1991 | Nphoton = 1.0/PI**2 * Wphib**2 |
| 1992 | |
| 1993 | RETURN |
| 1994 | END |
| 1995 | |
| 1996 | CDECK ID>, BESSELK0. |
| 1997 | REAL FUNCTION BESSELK0(X) |
| 1998 | C |
| 1999 | C Special functions I0,K0,I1,K1 |
| 2000 | C see Abramowitz&Stegun |
| 2001 | C |
| 2002 | IMPLICIT NONE |
| 2003 | REAL X,Y |
| 2004 | |
| 2005 | REAL BESSELI0 |
| 2006 | EXTERNAL BESSELI0 |
| 2007 | |
| 2008 | IF (X .LT. 2.0) THEN |
| 2009 | |
| 2010 | Y = X**2/4.0 |
| 2011 | |
| 2012 | BESSELK0 = |
| 2013 | & (-LOG(X/2.0)*BESSELI0(X))+(-.57721566+Y*(0.42278420 |
| 2014 | & +Y*(0.23069756+Y*(0.3488590E-1+Y*(0.262698E-2 |
| 2015 | & +Y*(0.10750E-3+Y*0.740E-5)))))) |
| 2016 | |
| 2017 | ELSE |
| 2018 | |
| 2019 | Y = 2.0/X |
| 2020 | |
| 2021 | BESSELK0 = |
| 2022 | & (EXP(-X)/SQRT(X))*(1.25331414+Y*(-0.7832358E-1 |
| 2023 | & +Y*(0.2189568E-1+Y*(-0.1062446E-1+Y*(0.587872E-2 |
| 2024 | & +Y*(-0.251540E-2+Y*0.53208E-3)))))) |
| 2025 | |
| 2026 | ENDIF |
| 2027 | |
| 2028 | RETURN |
| 2029 | END |
| 2030 | |
| 2031 | |
| 2032 | CDECK ID>, BESSELI0. |
| 2033 | REAL FUNCTION BESSELI0(X) |
| 2034 | IMPLICIT NONE |
| 2035 | |
| 2036 | REAL X,Y,AX |
| 2037 | |
| 2038 | AX = ABS(X) |
| 2039 | |
| 2040 | IF (AX .LT. 3.75) THEN |
| 2041 | |
| 2042 | Y = (X/3.75)**2 |
| 2043 | |
| 2044 | BESSELI0 = |
| 2045 | & 1.0+Y*(3.5156229+Y*(3.0899424+Y*(1.2067492 |
| 2046 | & +Y*(0.2659732+Y*(0.360768E-1+Y*0.45813E-2))))) |
| 2047 | |
| 2048 | ELSE |
| 2049 | |
| 2050 | Y = 3.75/AX |
| 2051 | |
| 2052 | BESSELI0 = |
| 2053 | & (EXP(AX)/SQRT(AX))*(0.39894228+Y*(0.1328592E-1 |
| 2054 | & +Y*(0.225319E-2+Y*(-0.157565E-2+Y*(0.916281E-2 |
| 2055 | & +Y*(-0.2057706E-1+Y*(0.2635537E-1+Y*(-0.1647633E-1 |
| 2056 | & +Y*0.392377E-2)))))))) |
| 2057 | |
| 2058 | ENDIF |
| 2059 | |
| 2060 | RETURN |
| 2061 | END |
| 2062 | |
| 2063 | |
| 2064 | CDECK ID>, BESSELK1. |
| 2065 | REAL FUNCTION BESSELK1(X) |
| 2066 | IMPLICIT NONE |
| 2067 | |
| 2068 | REAL X,Y |
| 2069 | |
| 2070 | REAL BESSELI1 |
| 2071 | EXTERNAL BESSELI1 |
| 2072 | |
| 2073 | IF (X .LT. 2.0) THEN |
| 2074 | |
| 2075 | Y = X**2/4.0 |
| 2076 | BESSELK1 = |
| 2077 | & (LOG(X/2.0)*BESSELI1(X))+(1.0/X)*(1.0+Y*(0.15443144 |
| 2078 | & +Y*(-0.67278579+Y*(-0.18156897+Y*(-0.1919402E-1 |
| 2079 | & +Y*(-0.110404E-2+Y*(-0.4686E-4))))))) |
| 2080 | |
| 2081 | ELSE |
| 2082 | |
| 2083 | Y=2.0/X |
| 2084 | BESSELK1 = |
| 2085 | & (EXP(-X)/SQRT(X))*(1.25331414+Y*(0.23498619 |
| 2086 | & +Y*(-0.3655620E-1+Y*(0.1504268E-1+Y*(-0.780353E-2 |
| 2087 | & +Y*(0.325614E-2+Y*(-0.68245E-3))))))) |
| 2088 | |
| 2089 | ENDIF |
| 2090 | |
| 2091 | RETURN |
| 2092 | END |
| 2093 | |
| 2094 | CDECK ID>, BESSELI1. |
| 2095 | REAL FUNCTION BESSELI1(X) |
| 2096 | IMPLICIT NONE |
| 2097 | |
| 2098 | REAL X,Y,AX |
| 2099 | |
| 2100 | AX = ABS(X) |
| 2101 | |
| 2102 | IF (AX .LT. 3.75) THEN |
| 2103 | |
| 2104 | Y = (X/3.75)**2 |
| 2105 | |
| 2106 | BESSELI1 = |
| 2107 | & AX*(0.5+Y*(0.87890594+Y*(0.51498869+Y*(0.15084934 |
| 2108 | & +Y*(0.2658733E-1+Y*(0.301532E-2+Y*0.32411E-3)))))) |
| 2109 | |
| 2110 | |
| 2111 | ELSE |
| 2112 | |
| 2113 | Y = 3.75/AX |
| 2114 | |
| 2115 | BESSELI1 = |
| 2116 | & 0.2282967E-1+Y*(-0.2895312E-1+Y*(0.1787654E-1 |
| 2117 | & -Y*0.420059E-2)) |
| 2118 | |
| 2119 | BESSELI1 = |
| 2120 | & 0.39894228+Y*(-0.3988024E-1+Y*(-0.362018E-2 |
| 2121 | & +Y*(0.163801E-2+Y*(-0.1031555E-1+Y*BESSELI1)))) |
| 2122 | |
| 2123 | BESSELI1 = BESSELI1 * |
| 2124 | & EXP(AX)/SQRT(AX) |
| 2125 | |
| 2126 | ENDIF |
| 2127 | |
| 2128 | IF (X .LT. 0) THEN |
| 2129 | BESSELI1 = -BESSELI1 |
| 2130 | ENDIF |
| 2131 | |
| 2132 | RETURN |
| 2133 | END |
| 2134 | |
| 2135 | CDECK ID>, GAUSSDAT. |
| 2136 | C DATA for the Gauss integration, x-values and weight for a |
| 2137 | C N=2,4,8,16,32,64 point Gauss integration. |
| 2138 | C |
| 2139 | C based on program GAUSS64 of D. Trautmann |
| 2140 | C data from Abramowitz & Stegun |
| 2141 | C |
| 2142 | C KAI HENCKEN, FEBRUAR 1993 |
| 2143 | C |
| 2144 | BLOCKDATA GAUSSDATA |
| 2145 | IMPLICIT NONE |
| 2146 | DOUBLE PRECISION XGAUSS(126),WGAUSS(126) |
| 2147 | COMMON /GAUSSD/XGAUSS,WGAUSS |
| 2148 | |
| 2149 | DATA XGAUSS(1)/ .57735026918962576D0/ |
| 2150 | DATA XGAUSS(2)/-.57735026918962576D0/ |
| 2151 | DATA WGAUSS(1)/ 1.00000000000000000D0/ |
| 2152 | DATA WGAUSS(2)/ 1.00000000000000000D0/ |
| 2153 | |
| 2154 | DATA XGAUSS(3)/ .33998104358485627D0/ |
| 2155 | DATA XGAUSS(4)/ .86113631159405258D0/ |
| 2156 | DATA XGAUSS(5)/-.33998104358485627D0/ |
| 2157 | DATA XGAUSS(6)/-.86113631159405258D0/ |
| 2158 | DATA WGAUSS(3)/ .65214515486254613D0/ |
| 2159 | DATA WGAUSS(4)/ .34785484513745385D0/ |
| 2160 | DATA WGAUSS(5)/ .65214515486254613D0/ |
| 2161 | DATA WGAUSS(6)/ .34785484513745385D0/ |
| 2162 | |
| 2163 | DATA XGAUSS(7)/ .18343464249564981D0/ |
| 2164 | DATA XGAUSS(8)/ .52553240991632899D0/ |
| 2165 | DATA XGAUSS(9)/ .79666647741362674D0/ |
| 2166 | DATA XGAUSS(10)/ .96028985649753623D0/ |
| 2167 | DATA XGAUSS(11)/-.18343464249564981D0/ |
| 2168 | DATA XGAUSS(12)/-.52553240991632899D0/ |
| 2169 | DATA XGAUSS(13)/-.79666647741362674D0/ |
| 2170 | DATA XGAUSS(14)/-.96028985649753623D0/ |
| 2171 | DATA WGAUSS(7)/ .36268378337836198D0/ |
| 2172 | DATA WGAUSS(8)/ .31370664587788727D0/ |
| 2173 | DATA WGAUSS(9)/ .22238103445337448D0/ |
| 2174 | DATA WGAUSS(10)/ .10122853629037627D0/ |
| 2175 | DATA WGAUSS(11)/ .36268378337836198D0/ |
| 2176 | DATA WGAUSS(12)/ .31370664587788727D0/ |
| 2177 | DATA WGAUSS(13)/ .22238103445337448D0/ |
| 2178 | DATA WGAUSS(14)/ .10122853629037627D0/ |
| 2179 | |
| 2180 | DATA XGAUSS(15)/ .0950125098376374402D0/ |
| 2181 | DATA XGAUSS(16)/ .281603550779258913D0/ |
| 2182 | DATA XGAUSS(17)/ .458016777657227386D0/ |
| 2183 | DATA XGAUSS(18)/ .617876244402643748D0/ |
| 2184 | DATA XGAUSS(19)/ .755404408355003034D0/ |
| 2185 | DATA XGAUSS(20)/ .865631202387831744D0/ |
| 2186 | DATA XGAUSS(21)/ .944575023073232576D0/ |
| 2187 | DATA XGAUSS(22)/ .989400934991649933D0/ |
| 2188 | DATA XGAUSS(23)/-.0950125098376374402D0/ |
| 2189 | DATA XGAUSS(24)/-.281603550779258913D0/ |
| 2190 | DATA XGAUSS(25)/-.458016777657227386D0/ |
| 2191 | DATA XGAUSS(26)/-.617876244402643748D0/ |
| 2192 | DATA XGAUSS(27)/-.755404408355003034D0/ |
| 2193 | DATA XGAUSS(28)/-.865631202387831744D0/ |
| 2194 | DATA XGAUSS(29)/-.944575023073232576D0/ |
| 2195 | DATA XGAUSS(30)/-.989400934991649933D0/ |
| 2196 | DATA WGAUSS(15)/ .189450610455068496D0/ |
| 2197 | DATA WGAUSS(16)/ .182603415044923589D0/ |
| 2198 | DATA WGAUSS(17)/ .169156519395002538D0/ |
| 2199 | DATA WGAUSS(18)/ .149595988816576732D0/ |
| 2200 | DATA WGAUSS(19)/ .124628971255533872D0/ |
| 2201 | DATA WGAUSS(20)/ .0951585116824927848D0/ |
| 2202 | DATA WGAUSS(21)/ .0622535239386478929D0/ |
| 2203 | DATA WGAUSS(22)/ .0271524594117540949D0/ |
| 2204 | DATA WGAUSS(23)/ .189450610455068496D0/ |
| 2205 | DATA WGAUSS(24)/ .182603415044923589D0/ |
| 2206 | DATA WGAUSS(25)/ .169156519395002538D0/ |
| 2207 | DATA WGAUSS(26)/ .149595988816576732D0/ |
| 2208 | DATA WGAUSS(27)/ .124628971255533872D0/ |
| 2209 | DATA WGAUSS(28)/ .0951585116824927848D0/ |
| 2210 | DATA WGAUSS(29)/ .0622535239386478929D0/ |
| 2211 | DATA WGAUSS(30)/ .0271524594117540949D0/ |
| 2212 | |
| 2213 | DATA XGAUSS(31)/ .0483076656877383162D0/ |
| 2214 | DATA XGAUSS(32)/ .144471961582796493D0/ |
| 2215 | DATA XGAUSS(33)/ .239287362252137075D0/ |
| 2216 | DATA XGAUSS(34)/ .331868602282127650D0/ |
| 2217 | DATA XGAUSS(35)/ .421351276130635345D0/ |
| 2218 | DATA XGAUSS(36)/ .506899908932229390D0/ |
| 2219 | DATA XGAUSS(37)/ .587715757240762329D0/ |
| 2220 | DATA XGAUSS(38)/ .663044266930215201D0/ |
| 2221 | DATA XGAUSS(39)/ .732182118740289680D0/ |
| 2222 | DATA XGAUSS(40)/ .794483795967942407D0/ |
| 2223 | DATA XGAUSS(41)/ .849367613732569970D0/ |
| 2224 | DATA XGAUSS(42)/ .896321155766052124D0/ |
| 2225 | DATA XGAUSS(43)/ .934906075937739689D0/ |
| 2226 | DATA XGAUSS(44)/ .964762255587506430D0/ |
| 2227 | DATA XGAUSS(45)/ .985611511545268335D0/ |
| 2228 | DATA XGAUSS(46)/ .997263861849481564D0/ |
| 2229 | DATA XGAUSS(47)/-.0483076656877383162D0/ |
| 2230 | DATA XGAUSS(48)/-.144471961582796493D0/ |
| 2231 | DATA XGAUSS(49)/-.239287362252137075D0/ |
| 2232 | DATA XGAUSS(50)/-.331868602282127650D0/ |
| 2233 | DATA XGAUSS(51)/-.421351276130635345D0/ |
| 2234 | DATA XGAUSS(52)/-.506899908932229390D0/ |
| 2235 | DATA XGAUSS(53)/-.587715757240762329D0/ |
| 2236 | DATA XGAUSS(54)/-.663044266930215201D0/ |
| 2237 | DATA XGAUSS(55)/-.732182118740289680D0/ |
| 2238 | DATA XGAUSS(56)/-.794483795967942407D0/ |
| 2239 | DATA XGAUSS(57)/-.849367613732569970D0/ |
| 2240 | DATA XGAUSS(58)/-.896321155766052124D0/ |
| 2241 | DATA XGAUSS(59)/-.934906075937739689D0/ |
| 2242 | DATA XGAUSS(60)/-.964762255587506430D0/ |
| 2243 | DATA XGAUSS(61)/-.985611511545268335D0/ |
| 2244 | DATA XGAUSS(62)/-.997263861849481564D0/ |
| 2245 | DATA WGAUSS(31)/ .0965400885147278006D0/ |
| 2246 | DATA WGAUSS(32)/ .0956387200792748594D0/ |
| 2247 | DATA WGAUSS(33)/ .0938443990808045654D0/ |
| 2248 | DATA WGAUSS(34)/ .0911738786957638847D0/ |
| 2249 | DATA WGAUSS(35)/ .0876520930044038111D0/ |
| 2250 | DATA WGAUSS(36)/ .0833119242269467552D0/ |
| 2251 | DATA WGAUSS(37)/ .0781938957870703065D0/ |
| 2252 | DATA WGAUSS(38)/ .0723457941088485062D0/ |
| 2253 | DATA WGAUSS(39)/ .0658222227763618468D0/ |
| 2254 | DATA WGAUSS(40)/ .0586840934785355471D0/ |
| 2255 | DATA WGAUSS(41)/ .0509980592623761762D0/ |
| 2256 | DATA WGAUSS(42)/ .0428358980222266807D0/ |
| 2257 | DATA WGAUSS(43)/ .0342738629130214331D0/ |
| 2258 | DATA WGAUSS(44)/ .0253920653092620595D0/ |
| 2259 | DATA WGAUSS(45)/ .0162743947309056706D0/ |
| 2260 | DATA WGAUSS(46)/ .00701861000947009660D0/ |
| 2261 | DATA WGAUSS(47)/ .0965400885147278006D0/ |
| 2262 | DATA WGAUSS(48)/ .0956387200792748594D0/ |
| 2263 | DATA WGAUSS(49)/ .0938443990808045654D0/ |
| 2264 | DATA WGAUSS(50)/ .0911738786957638847D0/ |
| 2265 | DATA WGAUSS(51)/ .0876520930044038111D0/ |
| 2266 | DATA WGAUSS(52)/ .0833119242269467552D0/ |
| 2267 | DATA WGAUSS(53)/ .0781938957870703065D0/ |
| 2268 | DATA WGAUSS(54)/ .0723457941088485062D0/ |
| 2269 | DATA WGAUSS(55)/ .0658222227763618468D0/ |
| 2270 | DATA WGAUSS(56)/ .0586840934785355471D0/ |
| 2271 | DATA WGAUSS(57)/ .0509980592623761762D0/ |
| 2272 | DATA WGAUSS(58)/ .0428358980222266807D0/ |
| 2273 | DATA WGAUSS(59)/ .0342738629130214331D0/ |
| 2274 | DATA WGAUSS(60)/ .0253920653092620595D0/ |
| 2275 | DATA WGAUSS(61)/ .0162743947309056706D0/ |
| 2276 | DATA WGAUSS(62)/ .00701861000947009660D0/ |
| 2277 | |
| 2278 | DATA XGAUSS(63)/ .02435029266342443250D0/ |
| 2279 | DATA XGAUSS(64)/ .0729931217877990394D0/ |
| 2280 | DATA XGAUSS(65)/ .121462819296120554D0/ |
| 2281 | DATA XGAUSS(66)/ .169644420423992818D0/ |
| 2282 | DATA XGAUSS(67)/ .217423643740007084D0/ |
| 2283 | DATA XGAUSS(68)/ .264687162208767416D0/ |
| 2284 | DATA XGAUSS(69)/ .311322871990210956D0/ |
| 2285 | DATA XGAUSS(70)/ .357220158337668116D0/ |
| 2286 | DATA XGAUSS(71)/ .402270157963991604D0/ |
| 2287 | DATA XGAUSS(72)/ .446366017253464088D0/ |
| 2288 | DATA XGAUSS(73)/ .489403145707052957D0/ |
| 2289 | DATA XGAUSS(74)/ .531279464019894546D0/ |
| 2290 | DATA XGAUSS(75)/ .571895646202634034D0/ |
| 2291 | DATA XGAUSS(76)/ .611155355172393250D0/ |
| 2292 | DATA XGAUSS(77)/ .648965471254657340D0/ |
| 2293 | DATA XGAUSS(78)/ .685236313054233243D0/ |
| 2294 | DATA XGAUSS(79)/ .719881850171610827D0/ |
| 2295 | DATA XGAUSS(80)/ .752819907260531897D0/ |
| 2296 | DATA XGAUSS(81)/ .783972358943341408D0/ |
| 2297 | DATA XGAUSS(82)/ .813265315122797560D0/ |
| 2298 | DATA XGAUSS(83)/ .840629296252580363D0/ |
| 2299 | DATA XGAUSS(84)/ .865999398154092820D0/ |
| 2300 | DATA XGAUSS(85)/ .889315445995114106D0/ |
| 2301 | DATA XGAUSS(86)/ .910522137078502806D0/ |
| 2302 | DATA XGAUSS(87)/ .929569172131939576D0/ |
| 2303 | DATA XGAUSS(88)/ .946411374858402816D0/ |
| 2304 | DATA XGAUSS(89)/ .961008799652053719D0/ |
| 2305 | DATA XGAUSS(90)/ .973326827789910964D0/ |
| 2306 | DATA XGAUSS(91)/ .983336253884625957D0/ |
| 2307 | DATA XGAUSS(92)/ .991013371476744321D0/ |
| 2308 | DATA XGAUSS(93)/ .996340116771955279D0/ |
| 2309 | DATA XGAUSS(94)/ .999305041735772139D0/ |
| 2310 | DATA XGAUSS(95)/-.02435029266342443250D0/ |
| 2311 | DATA XGAUSS(96)/-.0729931217877990394D0/ |
| 2312 | DATA XGAUSS(97)/-.121462819296120554D0/ |
| 2313 | DATA XGAUSS(98)/-.169644420423992818D0/ |
| 2314 | DATA XGAUSS(99)/-.217423643740007084D0/ |
| 2315 | DATA XGAUSS(100)/-.264687162208767416D0/ |
| 2316 | DATA XGAUSS(101)/-.311322871990210956D0/ |
| 2317 | DATA XGAUSS(102)/-.357220158337668116D0/ |
| 2318 | DATA XGAUSS(103)/-.402270157963991604D0/ |
| 2319 | DATA XGAUSS(104)/-.446366017253464088D0/ |
| 2320 | DATA XGAUSS(105)/-.489403145707052957D0/ |
| 2321 | DATA XGAUSS(106)/-.531279464019894546D0/ |
| 2322 | DATA XGAUSS(107)/-.571895646202634034D0/ |
| 2323 | DATA XGAUSS(108)/-.611155355172393250D0/ |
| 2324 | DATA XGAUSS(109)/-.648965471254657340D0/ |
| 2325 | DATA XGAUSS(110)/-.685236313054233243D0/ |
| 2326 | DATA XGAUSS(111)/-.719881850171610827D0/ |
| 2327 | DATA XGAUSS(112)/-.752819907260531897D0/ |
| 2328 | DATA XGAUSS(113)/-.783972358943341408D0/ |
| 2329 | DATA XGAUSS(114)/-.813265315122797560D0/ |
| 2330 | DATA XGAUSS(115)/-.840629296252580363D0/ |
| 2331 | DATA XGAUSS(116)/-.865999398154092820D0/ |
| 2332 | DATA XGAUSS(117)/-.889315445995114106D0/ |
| 2333 | DATA XGAUSS(118)/-.910522137078502806D0/ |
| 2334 | DATA XGAUSS(119)/-.929569172131939576D0/ |
| 2335 | DATA XGAUSS(120)/-.946411374858402816D0/ |
| 2336 | DATA XGAUSS(121)/-.961008799652053719D0/ |
| 2337 | DATA XGAUSS(122)/-.973326827789910964D0/ |
| 2338 | DATA XGAUSS(123)/-.983336253884625957D0/ |
| 2339 | DATA XGAUSS(124)/-.991013371476744321D0/ |
| 2340 | DATA XGAUSS(125)/-.996340116771955279D0/ |
| 2341 | DATA XGAUSS(126)/-.999305041735772139D0/ |
| 2342 | DATA WGAUSS(63)/ .0486909570091397204D0/ |
| 2343 | DATA WGAUSS(64)/ .0485754674415034269D0/ |
| 2344 | DATA WGAUSS(65)/ .0483447622348029572D0/ |
| 2345 | DATA WGAUSS(66)/ .0479993885964583077D0/ |
| 2346 | DATA WGAUSS(67)/ .0475401657148303087D0/ |
| 2347 | DATA WGAUSS(68)/ .0469681828162100173D0/ |
| 2348 | DATA WGAUSS(69)/ .0462847965813144172D0/ |
| 2349 | DATA WGAUSS(70)/ .0454916279274181445D0/ |
| 2350 | DATA WGAUSS(71)/ .0445905581637565631D0/ |
| 2351 | DATA WGAUSS(72)/ .0435837245293234534D0/ |
| 2352 | DATA WGAUSS(73)/ .0424735151236535890D0/ |
| 2353 | DATA WGAUSS(74)/ .0412625632426235286D0/ |
| 2354 | DATA WGAUSS(75)/ .0399537411327203414D0/ |
| 2355 | DATA WGAUSS(76)/ .0385501531786156291D0/ |
| 2356 | DATA WGAUSS(77)/ .0370551285402400460D0/ |
| 2357 | DATA WGAUSS(78)/ .0354722132568823838D0/ |
| 2358 | DATA WGAUSS(79)/ .0338051618371416094D0/ |
| 2359 | DATA WGAUSS(80)/ .0320579283548515535D0/ |
| 2360 | DATA WGAUSS(81)/ .0302346570724024789D0/ |
| 2361 | DATA WGAUSS(82)/ .0283396726142594832D0/ |
| 2362 | DATA WGAUSS(83)/ .0263774697150546587D0/ |
| 2363 | DATA WGAUSS(84)/ .0243527025687108733D0/ |
| 2364 | DATA WGAUSS(85)/ .0222701738083832542D0/ |
| 2365 | DATA WGAUSS(86)/ .0201348231535302094D0/ |
| 2366 | DATA WGAUSS(87)/ .0179517157756973431D0/ |
| 2367 | DATA WGAUSS(88)/ .0157260304760247193D0/ |
| 2368 | DATA WGAUSS(89)/ .0134630478967186426D0/ |
| 2369 | DATA WGAUSS(90)/ .0111681394601311288D0/ |
| 2370 | DATA WGAUSS(91)/ .00884675982636394772D0/ |
| 2371 | DATA WGAUSS(92)/ .00650445796897836286D0/ |
| 2372 | DATA WGAUSS(93)/ .00414703326056246764D0/ |
| 2373 | DATA WGAUSS(94)/ .00178328072169643295D0/ |
| 2374 | DATA WGAUSS(95)/ .0486909570091397204D0/ |
| 2375 | DATA WGAUSS(96)/ .0485754674415034269D0/ |
| 2376 | DATA WGAUSS(97)/ .0483447622348029572D0/ |
| 2377 | DATA WGAUSS(98)/ .0479993885964583077D0/ |
| 2378 | DATA WGAUSS(99)/ .0475401657148303087D0/ |
| 2379 | DATA WGAUSS(100)/ .0469681828162100173D0/ |
| 2380 | DATA WGAUSS(101)/ .0462847965813144172D0/ |
| 2381 | DATA WGAUSS(102)/ .0454916279274181445D0/ |
| 2382 | DATA WGAUSS(103)/ .0445905581637565631D0/ |
| 2383 | DATA WGAUSS(104)/ .0435837245293234534D0/ |
| 2384 | DATA WGAUSS(105)/ .0424735151236535890D0/ |
| 2385 | DATA WGAUSS(106)/ .0412625632426235286D0/ |
| 2386 | DATA WGAUSS(107)/ .0399537411327203414D0/ |
| 2387 | DATA WGAUSS(108)/ .0385501531786156291D0/ |
| 2388 | DATA WGAUSS(109)/ .0370551285402400460D0/ |
| 2389 | DATA WGAUSS(110)/ .0354722132568823838D0/ |
| 2390 | DATA WGAUSS(111)/ .0338051618371416094D0/ |
| 2391 | DATA WGAUSS(112)/ .0320579283548515535D0/ |
| 2392 | DATA WGAUSS(113)/ .0302346570724024789D0/ |
| 2393 | DATA WGAUSS(114)/ .0283396726142594832D0/ |
| 2394 | DATA WGAUSS(115)/ .0263774697150546587D0/ |
| 2395 | DATA WGAUSS(116)/ .0243527025687108733D0/ |
| 2396 | DATA WGAUSS(117)/ .0222701738083832542D0/ |
| 2397 | DATA WGAUSS(118)/ .0201348231535302094D0/ |
| 2398 | DATA WGAUSS(119)/ .0179517157756973431D0/ |
| 2399 | DATA WGAUSS(120)/ .0157260304760247193D0/ |
| 2400 | DATA WGAUSS(121)/ .0134630478967186426D0/ |
| 2401 | DATA WGAUSS(122)/ .0111681394601311288D0/ |
| 2402 | DATA WGAUSS(123)/ .00884675982636394772D0/ |
| 2403 | DATA WGAUSS(124)/ .00650445796897836286D0/ |
| 2404 | DATA WGAUSS(125)/ .00414703326056246764D0/ |
| 2405 | DATA WGAUSS(126)/ .00178328072169643295D0/ |
| 2406 | |
| 2407 | END |
| 2408 | *======================================================================= |
| 2409 | SUBROUTINE LORENB (U,PS,PI,PF) |
| 2410 | C |
| 2411 | C CERN PROGLIB# U102 LORENB .VERSION KERNFOR 4.04 821124 |
| 2412 | C ORIG. 20/08/75 L.PAPE |
| 2413 | C |
| 2414 | DOUBLE PRECISION PF4, FN |
| 2415 | DIMENSION PS(4),PI(4),PF(4) |
| 2416 | |
| 2417 | IF (PS(4).EQ.U) GO TO 17 |
| 2418 | PF4 = (PI(4)*PS(4)+PI(3)*PS(3)+PI(2)*PS(2)+PI(1)*PS(1)) / U |
| 2419 | FN = (PF4+PI(4)) / (PS(4)+U) |
| 2420 | PF(1)= PI(1) + FN*PS(1) |
| 2421 | PF(2)= PI(2) + FN*PS(2) |
| 2422 | PF(3)= PI(3) + FN*PS(3) |
| 2423 | PF(4)= PF4 |
| 2424 | GO TO 18 |
| 2425 | C |
| 2426 | 17 PF(1)= PI(1) |
| 2427 | PF(2)= PI(2) |
| 2428 | PF(3)= PI(3) |
| 2429 | PF(4)= PI(4) |
| 2430 | C |
| 2431 | 18 CONTINUE |
| 2432 | C |
| 2433 | RETURN |
| 2434 | C |
| 2435 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame