ZDC/AliZDCf.F

   1 *CMZ :          19/11/98  21.32.03  by  Federico Carminati
   2 *-- Author :
   3 C***********************************************************************
   4 C
   5       SUBROUTINE ZDC_ADDANG(THETA1,PHI1,THETA2,PHI2,THETA3,PHI3)
   6 C
   7 C  This subroutine takes spherical polar angles THETA1 and PHI1
   8 C  and adds to them the direction PHI2 on a cone of half-opening angle
   9 C  THETA2 (= angle from cone axis to cone surface) to produce
  10 C  the new space angles THETA3 and PHI3. All angles are in
  11 C  degrees, with 0<=THETA<=180 and 0<=PHI<360. There is no
  12 C  other restriction on THETA1, PHI1, THETA2 or PHI2.
  13 C
  14 C             27-AUG-1992          S. Coutu
  15 C
  16 C***********************************************************************
  17 #undef CERNLIB_GEANT321_GCONSP_INC
  18 #include "geant321/gconsp.inc"
  19       DOUBLE PRECISION THETA1,PHI1,THETA2,PHI2,THETA3,PHI3
  20       DOUBLE PRECISION TEMP,CX,CY,CZ,CT1,ST1,CT2,ST2,
  21      +     CP1,SP1,CP2,SP2,RTHETA3
  22 C
  23       CT1=COS(THETA1*DEGRAD)
  24       ST1=SIN(THETA1*DEGRAD)
  25       CP1=COS(PHI1*DEGRAD)
  26       SP1=SIN(PHI1*DEGRAD)
  27       CT2=COS(THETA2*DEGRAD)
  28       ST2=SIN(THETA2*DEGRAD)
  29       CP2=COS(PHI2*DEGRAD)
  30       SP2=SIN(PHI2*DEGRAD)
  31       CX=CT1*CP1*ST2*CP2+ST1*CP1*CT2-SP1*ST2*SP2
  32       CY=CT1*SP1*ST2*CP2+ST1*SP1*CT2+CP1*ST2*SP2
  33       CZ=CT1*CT2-ST1*ST2*CP2
  34 C
  35       RTHETA3=ACOS(CZ)
  36       THETA3=RTHETA3*RADDEG
  37       IF(THETA3.EQ.0.OR.THETA3.EQ.180.)THEN
  38          PHI3=0.
  39       ELSE
  40          TEMP=CX/SIN(RTHETA3)
  41          IF(TEMP.GT.1.)TEMP=1.
  42          IF(TEMP.LT.-1.)TEMP=-1.
  43          PHI3=ACOS(TEMP)*RADDEG
  44          IF(CY.LT.0.)PHI3=360.-PHI3
  45       ENDIF
  46 *
  47       END
  48
  49 *CMZ :          21/11/98  17.33.14  by  Federico Carminati
  50 *-- Author :
  51 C==================================================================
  52       SUBROUTINE ZDC_EFERMI(ID)
  53 *
  54 #include "zdc_common.inc"
  55 #undef CERNLIB_GEANT321_GCONSP_INC
  56 #include "geant321/gconsp.inc"
  57       DIMENSION RNDM(3)
  58 *
  59       CALL GRNDM(RNDM,3)
  60 *
  61       IF(ID.EQ.14)THEN
  62          DO I=1,200
  63             IF(RNDM(3).GE.PROBINTP(I-1).AND.RNDM(3).LT.PROBINTP(I))GO TO
  64      $           10
  65          END DO
  66       ELSE IF(ID.EQ.13) THEN
  67          DO I=1,200
  68             IF(RNDM(3).GE.PROBINTN(I-1).AND.RNDM(3).LT.PROBINTN(I))GO TO
  69      $           10
  70          END DO
  71       ENDIF
  72 *
  73    10 PIPPO=PP(I)+0.001
  74       PHI=TWOPI*RNDM(1)
  75       COST=1.-2*RNDM(2)
  76       TET=ACOS(COST)
  77       DDP(1)=PIPPO*SIN(TET)*COS(PHI)
  78       DDP(2)=PIPPO*SIN(TET)*SIN(PHI)
  79       DDP(3)=PIPPO*COST
  80 *
  81       END
  82
  83 *CMZ :          23/02/99  07.53.48  by  Federico Carminati
  84 *CMZ :  2.03/01 06/08/98  10.13.59  by  Federico Carminati
  85 *-- Author :   E.SCOMPARIN, 13/05/1996.
  86       SUBROUTINE ZDC_END
  87 C
  88 C *** TERMINATION OF THE ZDC SIMULATION AT THE END OF A RUN ***
  89 C
  90 C CALLED BY : SXEND
  91 C ORIGIN    : E. SCOMPARIN
  92 C
  93 C
  94 *      CALL HROUT(0,ICYCLE,' ')
  95 *      CALL HREND('ALIZ')
  96 *      CLOSE(1)
  97
  98       PRINT 10000
  99 10000 FORMAT(1H /1H ,36('*'),' ZDC_END  ',36('*')/
 100      $ 1H ,20X,'ZDC simulation package terminated'/
 101      $ 1H ,80('*'))
 102 *
 103       END
 104
 105 *CMZ :          23/02/99  07.53.48  by  Federico Carminati
 106 *CMZ :  2.03/01 07/08/98  09.48.28  by  Federico Carminati
 107 *-- Author :   E.SCOMPARIN, 13/05/1996.
 108       SUBROUTINE ZDC_EVE
 109 C
 110 C *** TERMINATION OF THE ZDC SIMULATION AFTER EACH EVENT ***
 111 C
 112 C CALLED BY : SXEVE
 113 C ORIGIN    : E. SCOMPARIN
 114 C
 115 C --  RETRIEVE HITS INFORMATION FOR THE CURRENT EVENT
 116 #undef CERNLIB_GEANT321_GCFLAG_INC
 117 #include "geant321/gcflag.inc"
 118 #include "zdc_common.inc"
 119 *
 120       INTEGER NUMVSP(NVZP),NUMVSN(NVZN),
 121      +        NUMBP(NVZP,1000),NUMBN(NVZN,1000),
 122      +        ITRP(1000),ITRN(1000)
 123       REAL HITP(NHZ,1000),HITN(NHZ,1000)
 124 *
 125       REAL ENEH(1000),XZPH(1000),YZPH(1000),EZPH(1000,NZPTX,NZPTY),
 126      +                XZNH(1000),YZNH(1000),EZNH(1000,NZNTX,NZNTY)
 127       INTEGER NPAH(1000),SFLH(1000),NPPH(1000,NZPTX,NZPTY),
 128      +                              NPNH(1000,NZNTX,NZNTY)
 129 *
 130       DATA NUMVSP/0,0/
 131       DATA NUMVSN/0,0/
 132 *
 133       DO JZERO=1,1000
 134          DO KZERO=1,NHZ
 135             HITP(KZERO,JZERO)=0
 136             HITN(KZERO,JZERO)=0
 137          ENDDO
 138       ENDDO
 139 *
 140       DO I=1,1000
 141          NPAH(I)=0
 142          XZPH(I)=99.
 143          YZPH(I)=99.
 144          XZNH(I)=99.
 145          YZNH(I)=99.
 146       ENDDO
 147 *
 148       ITOLD=0
 149 *      CALL GTREVE
 150 *
 151 C --  RETRIEVE HITS INFORMATION FOR THE CURRENT EVENT
 152       IF(ISECF.EQ.1.OR.(ISECF.EQ.0.AND.MOD(IEVENT,1).EQ.0))THEN
 153          WRITE(6,*)'EVENT NO. ',IEVENT
 154          CALL GPHITS('*','*')
 155          WRITE(6,*)
 156       ENDIF
 157 *
 158       NHMAXP=1000
 159       NHMAXN=1000
 160 *
 161       CALL GFHITS(ZSET,ZP1,NVZP,NHZ,NHMAXP,0,
 162      +            NUMVSP,ITRP,NUMBP,HITP,NHITP)
 163       CALL GFHITS(ZSET,ZN1,NVZN,NHZ,NHMAXN,0,
 164      +            NUMVSN,ITRN,NUMBN,HITN,NHITN)
 165 *
 166 * Fill n-tuple
 167 *
 168       ITCTP=0
 169       DO I=1,NHITP
 170          IF(ITRP(I).NE.ITCTP)THEN
 171             ITCTP=ITRP(I)
 172             NPAH(ITCTP)=HITP(1,I)
 173             ENEH(ITCTP)=HITP(2,I)
 174             XZPH(ITCTP)=HITP(3,I)
 175             YZPH(ITCTP)=HITP(4,I)
 176             SFLH(ITCTP)=INT(HITP(5,I))
 177          ENDIF
 178          J1=NUMBP(1,I)
 179          J2=NUMBP(2,I)
 180          EZPH(ITCTP,J1,J2)=HITP(7,I)
 181          NPPH(ITCTP,J1,J2)=HITP(6,I)
 182       ENDDO
 183       ITCTN=0
 184       DO I=1,NHITN
 185          IF(ITRN(I).NE.ITCTN)THEN
 186             ITCTN=ITRN(I)
 187             NPAH(ITCTN)=HITN(1,I)
 188             ENEH(ITCTN)=HITN(2,I)
 189             XZNH(ITCTN)=HITN(3,I)
 190             YZNH(ITCTN)=HITN(4,I)
 191             SFLH(ITCTN)=INT(HITN(5,I))
 192          ENDIF
 193          J1=NUMBN(1,I)
 194          J2=NUMBN(2,I)
 195          EZNH(ITCTN,J1,J2)=HITN(7,I)
 196          NPNH(ITCTN,J1,J2)=HITN(6,I)
 197       ENDDO
 198 *
 199       DO I=1,1000
 200          IF(NPAH(I).NE.0)THEN
 201             NPA=NPAH(I)
 202             ENE=ENEH(I)
 203             XZP=XZPH(I)
 204             YZP=YZPH(I)
 205             SFL=SFLH(I)
 206             DO J=1,NZPTX
 207                DO K=1,NZPTY
 208                   EZP(J,K)=EZPH(I,J,K)
 209                   NPP(J,K)=NPPH(I,J,K)
 210                ENDDO
 211             ENDDO
 212             XZN=XZNH(I)
 213             YZN=YZNH(I)
 214             DO J=1,NZNTX
 215                DO K=1,NZNTY
 216                   EZN(J,K)=EZNH(I,J,K)
 217                   NPN(J,K)=NPNH(I,J,K)
 218                ENDDO
 219             ENDDO
 220 *            CALL HFNT(10)
 221          ENDIF
 222       ENDDO
 223 *
 224       END
 225
 226 *CMZ :          21/11/98  17.33.14  by  Federico Carminati
 227 *-- Author :
 228 C==================================================================
 229       SUBROUTINE ZDC_GFERMI(A,Z)
 230 *
 231 #include "zdc_common.inc"
 232 #undef CERNLIB_GEANT321_GCONSP_INC
 233 #include "geant321/gconsp.inc"
 234 *     CALCOLO DELLA DISTRIBUZIONE DEI MOMENTI SECONDO
 235 *     LA DISTRIBUZIONE DOPPIA GAUSSIANA (ILINOV), UGUALE
 236 *     PER NEUTRONI E PROTONI
 237 *
 238       PROBINTP(0)=0.0
 239       PROBINTN(0)=0.0
 240       SIG1=0.113
 241       SIG2=0.250
 242       ALFA=0.18*((A/12.)**0.333)
 243       XK=(4.*PI)/((1.+ALFA)*TWOPI**1.5)
 244 *
 245       DO I=1,200
 246          P=FLOAT(I)*0.005
 247          PP(I)=P
 248          E1=(P*P)/(2.*SIG1*SIG1)
 249          E2=(P*P)/(2.*SIG2*SIG2)
 250          F1=EXP(-(E1))
 251          F2=EXP(-(E2))
 252          PROBP=XK*P*P*(F1/(SIG1**3)+F2/(SIG2**3))*0.005
 253          PROBINTP(I)=PROBINTP(I-1)+PROBP
 254          PROBINTN(I)=PROBINTP(I)
 255       END DO
 256 *
 257       END
 258
 259 *CMZ :          23/02/99  07.53.48  by  Federico Carminati
 260 *CMZ :  2.03/01 25/08/98  23.39.54  by  Federico Carminati
 261 *-- Author :
 262 C-- Author :    E.SCOMPARIN, 07/05/1996
 263       SUBROUTINE ZDC_INIT
 264 C
 265 C *** INITIALISATION FOR THE ZDC  SIMULATION ***
 266 C
 267 C CALLED BY : SXINIT
 268 C ORIGIN    : E. SCOMPARIN
 269 C
 270 C
 271 #undef CERNLIB_GEANT321_GCFLAG_INC
 272 #include "geant321/gcflag.inc"
 273 #undef CERNLIB_GEANT321_GCKINE_INC
 274 #include "geant321/gckine.inc"
 275 C
 276 #include "zdc_common.inc"
 277 *
 278 * Initialize COMMON block ZDC_CGEOM
 279 *
 280       DATA HDZN/4.0,4.0,50.0/
 281       DATA HDZP/8.0,8.0,75.0/
 282 C Coordinates of the center of the ZDC front face in the MRS
 283       DATA ZNPOS/-0.5,0.,11613./
 284       DATA ZPPOS/-19.0,0.,11563./
 285       DATA FIZN/0.,0.01825,50.0/
 286       DATA FIZP/0.,0.01825,75.0/
 287       DATA GRZN/0.025,0.025,50.0/
 288       DATA GRZP/0.040,0.040,75.0/
 289       DATA NCEN/11,11,0/
 290       DATA NCEP/10,10,0/
 291 *
 292 * Initialize COMMON block ZDC_HITS
 293 *
 294 * Data for set ZDC
 295       DATA ZSET,ZN1,ZP1 /'ZDC ', 'ZN1 ','ZP1 '/
 296       DATA CHNMSZN/'ZNTX','ZN1 '/
 297       DATA CHNMSZP/'ZPTX','ZP1 '/
 298       DATA NBITSZN /2,2/
 299       DATA NBITSZP /2,2/
 300       DATA IDTYPZN,IDTYPZP /1, 2/
 301       DATA NWHI,NWDI /1000, 1000/
 302 *
 303 * Data for Detectors ZN1,ZP1
 304 *
 305       DATA CHNAMH /'IPAR','EPAR','XPAR','YPAR','SFLG','FELE','ENER'/
 306       DATA NBITSH /8,32,32,32,1,32,32/
 307       DATA ORIG /0.,0.,10.,10.,0.,0.,0./
 308       DATA FACT /1.,1.,10000.,10000.,1.,10000.,10000./
 309       DATA NHSUM /2/
 310 *
 311 * Read tables of Cerenkov light produced in the fibers
 312 *
 313       OPEN(67,FILE='ZDC/light22620362207s',
 314      +    FORM='FORMATTED', STATUS='OLD')
 315       OPEN(68,FILE='ZDC/light22620362208s',
 316      +    FORM='FORMATTED', STATUS='OLD')
 317       OPEN(69,FILE='ZDC/light22620362209s',
 318      +    FORM='FORMATTED', STATUS='OLD')
 319       OPEN(70,FILE='ZDC/light22620362210s',
 320      +    FORM='FORMATTED', STATUS='OLD')
 321       DO 10 K=1,NALFA
 322          READ (67,10000)(TABLE(1,K,J),J=1,NBE)
 323          READ (68,10000)(TABLE(2,K,J),J=1,NBE)
 324          READ (69,10000)(TABLE(3,K,J),J=1,NBE)
 325          READ (70,10000)(TABLE(4,K,J),J=1,NBE)
 326    10 CONTINUE
 327 10000 FORMAT(7(1X,F10.5))
 328       CLOSE (67)
 329       CLOSE (68)
 330       CLOSE (69)
 331       CLOSE (70)
 332
 333       IFERMI=1
 334 *
 335 *  Set Fermi flag according to SWITCH 3
 336       IF(ISWIT(3).EQ.1) IFERMI=1
 337 *  define the LHC energy/nucleon
 338       EPERN=2760.
 339 *
 340 *  Open data file to read in particles from event generators
 341       IF(IHIJ.GE.1)THEN
 342          OPEN(27,FILE=FILEH,STATUS='OLD')
 343       ELSE IF(IVNU.GE.1)THEN
 344          OPEN(27,FILE=FILEV,STATUS='OLD')
 345       ENDIF
 346 C  GET FERMI MOMENTUM DISTRIBUTIONS
 347       CALL ZDC_GFERMI(207.,82.)
 348 C
 349 *      CALL HROPEN (1,'ALIZ',FILEA,'N',1024,ISTAT)
 350 *
 351 *      CALL HBNT(10,'ALIZDC','D')
 352 *      CALL HBNAME(10,'ALI',NPA,
 353 *     +'NPA[0,112]:U,ENE,SFL[0,1]:U,'//
 354 *     +'XZP:R*4:16:[-99,99],YZP:R*4:16:[-99,99],'//
 355 *     +'EZP(2,2),NPP(2,2),'//
 356 *     +'XZN:R*4:16:[-99,99],YZN:R*4:16:[-99,99],'//
 357 *     +'EZN(2,2),NPN(2,2)')
 358 *
 359       PRINT 10100
 360 10100 FORMAT(1H /1H ,36('*'),' ZDC_INIT ',36('*')/
 361      + 1H ,20X,'ZDC simulation package initialised'/
 362      + 1H ,80('*'))
 363
 364       END
 365
 366 *CMZ :          21/11/98  17.33.14  by  Federico Carminati
 367 *-- Author :
 368       SUBROUTINE ZDC_KINE(NT)
 369 *     Generation of event kinematics                              *
 370 * ----------------------------------------------------------------*
 371 #undef CERNLIB_GEANT321_GCKINE_INC
 372 #include "geant321/gckine.inc"
 373 #undef CERNLIB_GEANT321_GCFLAG_INC
 374 #include "geant321/gcflag.inc"
 375 #include "zdc_common.inc"
 376 *
 377       DIMENSION KATT(10000,2),PATT(10000,4)
 378 *
 379       INTEGER ITA(49),ITB(49)
 380       REAL PLAB(3),UBUF(1),UB(1)
 381       DIMENSION VERTEX(3)
 382       REAL BALP(4),DDDP(4),PLABS(4)
 383 *
 384       DATA UBUF/0./
 385       DATA ITA/22,-11,11,0,-13,13,111,211,-211,130,321,-321,
 386      +2112,2212,-2212,310,221,3122,3222,3212,3112,3322,3312,
 387      +3334,-2112,-3122,-3112,-3212,-3222,-3322,-3312,-3334,
 388      +-15,15,411,-411,421,-421,431,-431,4122,0,0,0,0,0,0,0,0/
 389       DATA ITB/10,-12,12,0,-14,14,110,120,-120,230,130,-130,
 390      +1220,1120,-1120,0,220,2130,1130,1230,2230,1330,2330,3331,
 391      +-1220,-2130,-1130,-1230,-2230,-1330,-2330,-3331,-16,16,
 392      +-240,240,-140,140,0,0,2140,0,0,0,0,0,0,0,0/
 393       DATA IVR/0/
 394 *
 395       DIMENSION RNDM(2)
 396       CHARACTER*20 PARTNAME
 397 *
 398 * --- BEAMTY=1 Gaussian beam (width sigx sigy, mean values fx fy)
 399 * --- BEAMTY=2 Uniform beam  (width sigx sigy, mean values fx fy)
 400 *     Data cards BMTY, INBM
 401 *
 402       DO KZ=1,3
 403          VERTEX(KZ)=0
 404       ENDDO
 405       IF(BEAMTY.EQ.1) THEN
 406          CALL SXGAUS(RNDM,2)
 407          VERTEX(1)= VERTEX(1)+FX+RNDM(1)*SIGX
 408          VERTEX(2)= VERTEX(2)+FY+RNDM(2)*SIGY
 409          VERTEX(3)= 0.05
 410       ELSEIF(BEAMTY.EQ.2) THEN
 411          CALL GRNDM(RNDM,2)
 412          VERTEX(1)= VERTEX(1)+FX+(2.*RNDM(1)-1.)*SIGX
 413          VERTEX(2)= VERTEX(2)+FY+(2.*RNDM(2)-1.)*SIGY
 414          VERTEX(3)= 0.05
 415       ENDIF
 416 *
 417 *  Fill the vertex bank
 418 *
 419       CALL GSVERT(VERTEX,0,0,0.,0,NVTX)
 420       IF(NVTX.EQ.0)WRITE(6,*)'Error defining vertex'
 421 *
 422 * Read HIJING event file
 423 *
 424       IF(IHIJ.GE.1)THEN
 425    10    READ(27,*,END=70) NMUL,EATT,JATT,NT,NP,N0,N01,N10,N11,B,NATT
 426          IF(IHIJ.EQ.2)THEN
 427             WRITE(6,*)'HIJING EVENT HEADER= ',NMUL,EATT,JATT,NT,NP,
 428      +      N0,N01,N10,N11,B,NATT
 429          ENDIF
 430          READ(27,*) (KATT(I,1),KATT(I,2),(PATT(I,J),J=1,4),I=1,NATT)
 431          IVR=IVR+1
 432          IF(IVR.LT.IHIJF) GOTO 10
 433 *
 434          BIMP=B
 435          PARTI=FLOAT(NT)+FLOAT(NP)
 436          RMULTI=FLOAT(NATT)
 437 *
 438          DO I=1,NATT
 439 * --  Remove particles with negative z-direction
 440             IF(PATT(I,3).LT.0.0) GOTO 30
 441 * --  Remove (if required) the spectators
 442             IF(IHIJSP.EQ.1)THEN
 443                IF(KATT(I,2).EQ.0..OR.KATT(I,2).EQ.10.) GOTO 30
 444             ENDIF
 445 * --  Define participant flag : UBUF(1)
 446             IF(KATT(I,2).EQ.0..OR.KATT(I,2).EQ.10.) THEN
 447                UBUF(1)=0.
 448             ELSE
 449                UBUF(1)=1.
 450             ENDIF
 451 * --  Translate the particle ID from HIJING to GEANT
 452             DO J=1,49
 453                IF(KATT(I,1).EQ.ITA(J)) THEN
 454                   IPART=J
 455                   GOTO 20
 456                ENDIF
 457             ENDDO
 458    20       CONTINUE
 459 *
 460             DO JJ=1,3
 461                PLAB(JJ)=PATT(I,JJ)
 462             ENDDO
 463 *
 464 * Rotations to account for beam divergence and crossing angle
 465             IF(BMDIV.NE.0.) CALL ZDC_ROTP(PLAB,0)
 466             IF(BMCRA.NE.0.) CALL ZDC_ROTP(PLAB,1)
 467 *
 468 * Fill the tracks bank
 469             CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
 470             IF(IHIJ.EQ.2)THEN
 471                WRITE(6,*)'IPART= ',IPART,'PLAB= ',PLAB
 472             ENDIF
 473             IF(NT.EQ.0)WRITE(6,*)'Error defining track'
 474 *
 475    30       CONTINUE
 476          ENDDO
 477 *
 478 * Read VENUS event file
 479 *
 480       ELSE IF(IVNU.GE.1)THEN
 481    40    READ(27,*,END=80)LEVT,NREVT,NMUL,B,KOLEVT,COLEVT, PMXEVT,
 482      +   EGYEVT,NP,NT,NATT
 483          IF(LEVT.NE.1) THEN
 484             WRITE(6,*) ' Error in reading VENUS data'
 485             STOP
 486          ENDIF
 487          IF(IVNU.EQ.2)THEN
 488             WRITE(6,*)'VENUS EVENT HEADER= ',LEVT,NREVT,NATT,B,KOLEVT,
 489      +      COLEVT,PMXEVT,EGYEVT,NP,NT
 490          ENDIF
 491 *
 492          DO I=1,NATT
 493 *--  read produced particles
 494             READ(27,*,END=80) LPTL,NREVT,NRPTL,KATT(I,1), (PATT(I,J),J=
 495      +      1,4),AMASSS
 496 *
 497 *         lptl ......... record label (lptl=3)
 498 *         nrevt ........ event number
 499 *         nrptl ........ particle number
 500 *         katt(1,1)..... particle id
 501 *         patt ......... 4-momentum (px,py,pz,en) in lab
 502 *         amasss ........ mass of particle
 503 *
 504             IF(LPTL.NE.3) THEN
 505                WRITE(6,*) ' ERROR IN READING VENUS DATA'
 506                STOP
 507             ENDIF
 508          ENDDO
 509 *
 510          IVR=IVR+1
 511          IF(IVR.LT.IVNUF) GOTO 40
 512 *
 513          BIMP=B
 514          PARTI=(FLOAT(NT)+FLOAT(NP))/2.
 515          RMULTI=FLOAT(NATT)
 516 *
 517          DO I=1,NATT
 518 * --  remove particles in negative direction
 519             IF(PATT(I,3).LT.0.0) GOTO 60
 520 * --  if required remove the spectators
 521             IF(IVNSP.EQ.1)THEN
 522                IF((KATT(I,1).EQ.1120.OR.KATT(I,1).EQ.1220).AND. (PATT(I
 523      +         ,4).LT.(EPERN*1.005).AND. PATT(I,4).GT.(EPERN*0.995)))
 524      +         GOTO 60
 525             ENDIF
 526 * --  flag the spectators (UBUF(1)=0)
 527             IF((KATT(I,1).EQ.1120.OR.KATT(I,1).EQ.1220).AND. (PATT(I,4)
 528      +      .GT.(EPERN*0.995).AND. PATT(I,4).LT.(EPERN*1.005))) THEN
 529                UBUF(1)=0.
 530             ELSE
 531                UBUF(1)=1.
 532             ENDIF
 533 * --  Translate the particle ID from VENUS to GEANT
 534             DO J=1,49
 535                IF(KATT(I,1).EQ.ITB(J))THEN
 536                   IPART=J
 537                   GOTO 50
 538                ENDIF
 539             ENDDO
 540    50       CONTINUE
 541             DO JJ=1,3
 542                PLAB(JJ)=PATT(I,JJ)
 543             ENDDO
 544 *
 545 * Rotations to account for beam divergence and crossing angle
 546             IF(BMCRA.NE.0.)CALL ZDC_ROTP(PLAB,1)
 547             IF(BMDIV.NE.0.)CALL ZDC_ROTP(PLAB,0)
 548 *
 549 * Fill the tracks bank
 550             CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
 551             IF(IVNU.EQ.2)THEN
 552                WRITE(6,*)'IPART= ',IPART,'PLAB= ',PLAB
 553             ENDIF
 554             IF(NT.EQ.0)WRITE(6,*)'Error defining track'
 555    60       CONTINUE
 556          ENDDO
 557 *
 558       ELSE
 559 *
 560 * Generation of test particles
 561 *
 562 * Definition of the incident particle
 563 *
 564          IPART = NINT(PKINE(1))
 565          UBUF(1)=1
 566          CALL GFPART(IPART,PARTNAME,ITRTYP,AMASS,CHARGE,TLIFE,UB,NB)
 567          AAP=1.0
 568 *
 569 * --- In case of a ion , find out energy and momentum
 570 * --- NAAP=A=Z+N - atomic number of the initial particle
 571          IF(IPART.GE.45) THEN
 572             NAAP = INT(0.5 + AMASS/0.93149432)
 573             AAP = NAAP
 574             UBUF(1)=0.
 575          ENDIF
 576          IF((IPART.GE.45).OR.(IPART.EQ.13).OR.(IPART.EQ.14)) THEN
 577             UBUF(1)=0.
 578          ENDIF
 579 *
 580          IF(IPART.EQ.14)THEN
 581             AMA=0.938
 582          ELSE IF(IPART.EQ.13)THEN
 583             AMA=0.939
 584          ELSE
 585             AMA=AMASS
 586          ENDIF
 587 *
 588          IF(PKINE(6).EQ.0.)THEN
 589             PLAB(1)=PKINE(2)*PKINE(3)
 590             PLAB(2)=PKINE(2)*PKINE(4)
 591             PLAB(3)=PKINE(2)*PKINE(5)
 592          ELSE
 593             SCANG=2*ATAN(EXP(-(PKINE(6))))
 594             PLAB(1)=-PKINE(2)*SIN(SCANG)
 595             PLAB(2)=0.
 596             PLAB(3)=PKINE(2)*COS(SCANG)
 597          ENDIF
 598 *
 599 * Rotations to account for beam divergence and crossing angle
 600          IF(BMCRA.NE.0.)CALL ZDC_ROTP(PLAB,1)
 601          IF(BMDIV.NE.0.)CALL ZDC_ROTP(PLAB,0)
 602 *
 603 * --  in case of a nucleon, if required, apply the Fermi momentum
 604          IF(UBUF(1).EQ.0.AND.IFERMI.EQ.1) THEN
 605             IF(IPART.EQ.13.OR.IPART.EQ.14)THEN
 606                CALL ZDC_EFERMI(IPART)
 607             ELSE IF(IPART.GE.45)THEN
 608                CALL ZDC_EFERMI(13)
 609                GOLDHABER=SQRT(AAP*(207-AAP)/(207-1.))
 610                DO JJ=1,3
 611                   DDP(JJ)=DDP(JJ)*GOLDHABER
 612                ENDDO
 613             ENDIF
 614             DO II=1,3
 615                BALP(II)=-PLAB(II)
 616             ENDDO
 617             BALP(4)=SQRT(BALP(1)**2+BALP(2)**2+BALP(3)**2+AMA**2)
 618             DO II=1,3
 619                DDDP(II)=DDP(II)
 620             ENDDO
 621             DDDP(4)=SQRT(DDDP(1)**2+DDDP(2)**2+DDDP(3)**2+AMA**2)
 622             CALL LORENF(AMA,BALP,DDDP,PLABS)
 623             DO II=1,3
 624                PLAB(II)=PLABS(II)
 625             ENDDO
 626 c           WRITE(6,*)'BALP ',BALP,'DDDP ',DDDP,'PLABS ',PLABS
 627          ENDIF
 628 *
 629          BIMP=999.
 630 *
 631 *  Fill the tracks bank
 632          CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
 633          IF(NT.EQ.0)WRITE(6,*)'Error defining track'
 634
 635       ENDIF
 636 *
 637 * Some debug printings
 638       IF(IHIJ.EQ.2.OR.IVNU.EQ.2)THEN
 639          CALL GPVERT(NVTX)
 640          DO KK=1,NT
 641             CALL GPKINE(KK)
 642          ENDDO
 643       ENDIF
 644 *
 645 C
 646       RETURN
 647    70 WRITE(6,*) 'END OF HIJING FILE'
 648       IEORUN=1
 649       GOTO 999
 650    80 WRITE(6,*) 'END OF VENUS FILE'
 651       IEORUN=1
 652 *
 653   999 END
 654
 655 *CMZ :          19/11/98  09.54.52  by  Federico Carminati
 656 *-- Author :
 657       SUBROUTINE  ZDC_ROTP(PLAB,ICROSS)
 658 #undef CERNLIB_GEANT321_GCKINE_INC
 659 #include "geant321/gckine.inc"
 660 #undef CERNLIB_GEANT321_GCFLAG_INC
 661 #include "geant321/gcflag.inc"
 662 C
 663 #undef CERNLIB_GEANT321_GCONSP_INC
 664 #include "geant321/gconsp.inc"
 665 #include "zdc_common.inc"
 666 C
 667       DOUBLE PRECISION TETPART,FIPART,TETDIV,FIDIV,TETSUM,FISUM
 668       DOUBLE PRECISION DPLAB(3)
 669       DIMENSION PLAB(3)
 670       DIMENSION RVEC(1)
 671       DIMENSION RNDM(1)
 672 C
 673 C TAKE NON-ZERO DIVERGENCY
 674 c        write(6,*)'PLAB (NODIV) = ',plab
 675       PMOD=0.0
 676       DO IJ=1,3
 677          DPLAB(IJ)=PLAB(IJ)
 678          PMOD=PMOD+PLAB(IJ)*PLAB(IJ)
 679       ENDDO
 680       PMOD=SQRT(PMOD)
 681       IF(ICROSS.EQ.0)THEN
 682 C---  TETDIV is given in radiants and it is read by data cards
 683          CALL RNORML(RVEC,1)
 684 c          WRITE(6,*)RVEC
 685          TETDIV=PKINE(10)*ABS(RVEC(1))
 686 C---  PHIDIV is randomly chosen between 0 and 360
 687          CALL GRNDM(RNDM,1)
 688          FIDIV=RNDM(1)*TWOPI
 689       ELSE IF(ICROSS.EQ.1)THEN
 690          IF(ISWIT(7).EQ.0)THEN
 691             TETDIV=0.
 692             FIDIV=0.
 693          ELSE IF(ISWIT(7).EQ.1)THEN
 694             TETDIV=PKINE(4)
 695             FIDIV=0.
 696          ELSE IF(ISWIT(7).EQ.2)THEN
 697             TETDIV=PKINE(4)
 698             FIDIV=PIBY2
 699          ENDIF
 700       ENDIF
 701       TETPART=DATAN(DSQRT(DPLAB(1)**2+DPLAB(2)**2)/DPLAB(3))
 702       IF(DPLAB(2).NE.0..OR.DPLAB(1).NE.0.) THEN
 703          FIPART=DATAN2(DPLAB(2),DPLAB(1))
 704       ELSE
 705          FIPART=0.
 706       ENDIF
 707       IF (FIPART.LT.0)FIPART=FIPART+TWOPI
 708 C        write(6,*)' TETPART= ',tetpart,' FIPART= ',fipart
 709       TETDIV=TETDIV*RADDEG
 710       FIDIV=FIDIV*RADDEG
 711       TETPART=TETPART*RADDEG
 712       FIPART=FIPART*RADDEG
 713 c        CALL ZDC_ADDANG(TETDIV,FIDIV,TETPART,FIPART,TETSUM,FISUM)
 714       CALL ZDC_ADDANG(TETPART,FIPART,TETDIV,FIDIV,TETSUM,FISUM)
 715 C        write(6,*)' TETSUM= ',tetsum,'FISUM= ',fisum
 716       TETSUM=TETSUM*DEGRAD
 717       FISUM=FISUM*DEGRAD
 718       PLAB(1)=PMOD*DSIN(TETSUM)*DCOS(FISUM)
 719       PLAB(2)=PMOD*DSIN(TETSUM)*DSIN(FISUM)
 720       PLAB(3)=PMOD*DCOS(TETSUM)
 721 C        write(6,*)'PLAB (DIV) = ',plab
 722       RETURN
 723       END
 724
 725 *CMZ :          19/11/98  09.54.52  by  Federico Carminati
 726 *CMZ :  2.03/01 25/08/98  23.40.04  by  Federico Carminati
 727 *-- Author :   E.SCOMPARIN, 13/05/1996.
 728       SUBROUTINE ZDC_STEP
 729 C
 730 C *** RECORDING OF ZDC HITS AFTER EACH TRACKING STEP ***
 731 C
 732 C CALLED BY : SXSTEP
 733 C ORIGIN    : E. SCOMPARIN
 734 C
 735 #undef CERNLIB_GEANT321_GCKINE_INC
 736 #undef CERNLIB_GEANT321_GCONSP_INC
 737 #undef CERNLIB_GEANT321_GCVOLU_INC
 738 #undef CERNLIB_GEANT321_GCKMAX_INC
 739 #undef CERNLIB_GEANT321_GCKING_INC
 740 #undef CERNLIB_GEANT321_GCNUM_INC
 741 #undef CERNLIB_GEANT321_GCTRAK_INC
 742 #undef CERNLIB_GEANT321_GCFLAG_INC
 743 #undef CERNLIB_GEANT321_GCTMED_INC
 744 #undef CERNLIB_GEANT321_GCMATE_INC
 745 #undef CERNLIB_GEANT321_GCSETS_INC
 746 #include "geant321/gckine.inc"
 747 #include "geant321/gconsp.inc"
 748 #include "geant321/gcvolu.inc"
 749 #include "geant321/gcking.inc"
 750 #include "geant321/gcnum.inc"
 751 #include "geant321/gctrak.inc"
 752 #include "geant321/gcflag.inc"
 753 #include "geant321/gctmed.inc"
 754 #include "geant321/gcmate.inc"
 755 #include "geant321/gcsets.inc"
 756 #include "zdc_common.inc"
 757 *
 758       DIMENSION LNAM(15),LNUM(15)
 759       DIMENSION UB(1)
 760 *
 761       DIMENSION VERT0(3),PVERT0(4)
 762       CHARACTER*20 CHNPART0
 763 *
 764       DIMENSION XLOCM(3),XANGM(3),XLOCL(3),XANGL(3),FNULA(3),ZENTR(3)
 765       DATA ALFMAX /110.0/
 766 *
 767 * New track
 768 *
 769       IF(ITRA.NE.ITOLD)THEN
 770          IF(INWVOL.EQ.1.AND.IDTYPE.GT.0)THEN
 771 *
 772 * Impact point on the ZDC (either ZPRO or ZNEU)
 773 *
 774             NLEV0=NLEVEL
 775             CALL GFPATH(ISET,IDET,NUMBV,NLEV,LNAM,LNUM)
 776             CALL UCTOH('ZNEU',IZNEU,4,4)
 777             CALL UCTOH('ZPRO',IZPRO,4,4)
 778             DO I=NLEV,1,-1
 779                IF(LNAM(I).EQ.IZNEU.OR.LNAM(I).EQ.IZPRO)THEN
 780                   NLEVEL=0
 781                   CALL GLVOLU(I,LNAM,LNUM,IER)
 782                   IF(IDTYPE.EQ.1)THEN
 783                      CALL GMTOD(VECT,HITSBUPN(3),1)
 784                      HITSBUPP(3)=99.
 785                      HITSBUPP(4)=99.
 786                   ELSE IF(IDTYPE.EQ.2)THEN
 787                      CALL GMTOD(VECT,HITSBUPP(3),1)
 788                      HITSBUPN(3)=99.
 789                      HITSBUPN(4)=99.
 790                   ENDIF
 791                   NLEVEL=NLEV0
 792                   GOTO 10
 793                ENDIF
 794             ENDDO
 795    10       CONTINUE
 796 *
 797 * Particle type and energy for the primary track
 798 *
 799             CALL GFKINE(ITRA,VERT0,PVERT0,IPART0,NVERT0,UB,NB)
 800             CALL GFPART(IPART0,CHNPART0,ITRTYP0,AMASS0,CHARGE0,TLIFE0,
 801      +      UB,NB)
 802             PVERTM=VMOD(PVERT0,3)
 803             GEKIN0=SQRT(PVERTM**2+AMASS0**2)
 804             IF(IPART.NE.IPART0.OR.GEKIN.LT.GEKIN0*0.95)THEN
 805                HITSBUPP(1)=FLOAT(IPART0)
 806                HITSBUPP(2)=GEKIN0
 807                HITSBUPP(5)=1.
 808                HITSBUPN(1)=FLOAT(IPART0)
 809                HITSBUPN(2)=GEKIN0
 810                HITSBUPN(5)=1.
 811             ELSE
 812                HITSBUPP(1)=FLOAT(IPART)
 813                HITSBUPP(2)=GEKIN
 814                HITSBUPP(5)=0.
 815                HITSBUPN(1)=FLOAT(IPART)
 816                HITSBUPN(2)=GEKIN
 817                HITSBUPN(5)=0.
 818             ENDIF
 819             HITSBUPP(6)=0.
 820             HITSBUPP(7)=0.
 821             HITSBUPN(6)=0.
 822             HITSBUPN(7)=0.
 823             IF(IDTYPE.EQ.1)THEN
 824                CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HITSBUPN,IHIT)
 825             ELSE IF(IDTYPE.EQ.2)THEN
 826                CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HITSBUPP,IHIT)
 827             ENDIF
 828             ITOLD=ITRA
 829          ENDIF
 830       ENDIF
 831 *
 832 * Energy deposition in the passive material
 833 *
 834       IF(INWVOL.EQ.0.AND.IDTYPE.GT.0)THEN
 835          HITS(6)=0.
 836          IF(ISTOP.EQ.2)THEN
 837             HITS(7)=GEKIN
 838          ELSE
 839             HITS(7)=DESTEP
 840          ENDIF
 841          IF(IDTYPE.EQ.1)THEN
 842             DO KCOPY=1,5
 843                HITSBUPN(KCOPY)=HITS(KCOPY)
 844             ENDDO
 845          ELSE IF(IDTYPE.EQ.2)THEN
 846             DO KCOPY=1,5
 847                HITSBUPP(KCOPY)=HITS(KCOPY)
 848             ENDDO
 849          ENDIF
 850          CALL GSCHIT(ISET,IDET,ITRA,NUMBV,HITS,NHSUM,IHIT)
 851       ENDIF
 852 *
 853 * Select charged relativistic particles
 854 *
 855       IF(ABS(CHARGE).EQ.0)RETURN
 856       IF(GEKIN.GT.0.00001) THEN
 857          BETA  = VECT(7)/GETOT
 858       ELSE
 859          RETURN
 860       ENDIF
 861       IF(BETA.LT.0.67) RETURN
 862
 863 *
 864 * Particle inside fiber
 865 *
 866       IF(NUMED.EQ.803.AND.IDTYPE.GT.0) THEN
 867 *
 868          XLOCM(1)= VECT(1)
 869          XLOCM(2)= VECT(2)
 870          XLOCM(3)= VECT(3)
 871          XANGM(1)= VECT(4)
 872          XANGM(2)= VECT(5)
 873          XANGM(3)= VECT(6)
 874 *
 875 * --- Calculating the position of the centre of the fiber in MARS
 876          DO KZERO=1,3
 877             FNULA(KZERO)=0
 878             ZENTR(KZERO)=0
 879          ENDDO
 880          CALL GDTOM(FNULA,ZENTR,1)
 881 *
 882 * --- Angle between the particle trajectory and the fiber axis
 883 *
 884          CALL GMTOD(XANGM,XANGL,2)
 885          IF(XANGL(3).GT.1.0) XANGL(3)=SIGN(1.0,XANGL(3))
 886          ALFAR=ACOS(XANGL(3))
 887          ALFAD=ALFAR*RADDEG
 888 *
 889          IF (ALFAD.GT.ALFMAX) RETURN
 890          IALFA = INT(1.+ALFAD/2.)
 891 *
 892 * --- Distance between the particle trajectory and the fiber axis
 893 *
 894          CALL GMTOD(XLOCM,XLOCL,1)
 895          IF(ABS(XANGL(1)).GT.0.0001) THEN
 896             DCOEFF = XANGL(2)/XANGL(1)
 897             BE=ABS(XLOCL(2)-DCOEFF*XLOCL(1))/SQRT(DCOEFF*DCOEFF+1.)
 898          ELSE
 899             BE = ABS(XLOCL(1))
 900          ENDIF
 901          IF(IDTYPE.EQ.1)THEN
 902             RADIUS=FIZN(2)
 903          ELSEIF(IDTYPE.EQ.2)THEN
 904             RADIUS=FIZP(2)
 905          ENDIF
 906          IF(BE.GT.RADIUS) THEN
 907             PRINT*,'DISTANCE GREATER THAN FIBER RADIUS'
 908             RETURN
 909          ENDIF
 910          IBE = 1+INT(1000.*BE)
 911          IF(IBE.GT.NBE)IBE=NBE
 912 *
 913 * --- Finding out proper index for beta
 914          IF(BETA.GE.0.67 .AND. BETA.LE.0.75) IBETA=1
 915          IF(BETA.GT.0.75 .AND. BETA.LE.0.85) IBETA=2
 916          IF(BETA.GT.0.85 .AND. BETA.LE.0.95) IBETA=3
 917          IF(BETA.GT.0.95 .AND. BETA.LE.1.00) IBETA=4
 918 *
 919 * --- Looking into the tables
 920          OUT = CHARGE*CHARGE*TABLE(IBETA,IALFA,IBE)
 921          CALL RNPSSN(OUT,N,IERROR)
 922          HITS(6)=FLOAT(N)
 923          HITS(7)=0.
 924          IF(IDTYPE.EQ.1)THEN
 925             CALL UCOPY(HITSBUPN,HITS,4)
 926          ELSE IF(IDTYPE.EQ.2)THEN
 927             CALL UCOPY(HITSBUPP,HITS,4)
 928          ENDIF
 929          CALL GSCHIT(ISET,IDET,ITRA,NUMBV,HITS,NHSUM,IHIT)
 930       ENDIF
 931 *
 932   999 END
 933
 934 *CMZ :          18/11/98  16.55.53  by  Federico Carminati
 935 *CMZ :  2.03/01 06/08/98  10.13.59  by  Federico Carminati
 936 *-- Author :   E. SCOMPARIN
 937       SUBROUTINE ZDC_TRKI(IFLAG)
 938 C
 939 C *** DEFINITION OF THE GEOMETRY OF THE PHOS ***
 940 C *** NVE 24-SEP-1990 CERN GENEVA ***
 941 C
 942 C CALLED BY : SXTRKI
 943 C ORIGIN    : E. SCOMPARIN
 944 C
 945 #undef CERNLIB_GEANT321_GCKINE_INC
 946 #undef CERNLIB_GEANT321_GCKMAX_INC
 947 #undef CERNLIB_GEANT321_GCKING_INC
 948 #undef CERNLIB_GEANT321_GCTRAK_INC
 949 #undef CERNLIB_GEANT321_GCTMED_INC
 950 #undef CERNLIB_GEANT321_GCFLAG_INC
 951 #undef CERNLIB_GEANT321_GCMATE_INC
 952 #include "geant321/gckine.inc"
 953 #include "geant321/gcking.inc"
 954 #include "geant321/gctrak.inc"
 955 #include "geant321/gctmed.inc"
 956 #include "geant321/gcflag.inc"
 957 #include "geant321/gcmate.inc"
 958 #include "zdc_common.inc"
 959
 960       REAL UBUF(1)
 961       DATA ITOLDG/0/
 962
 963 * Fetch new track (VENUS or HIJING)
 964       IF(ITRA.NE.ITOLDG)THEN
 965          CALL GFKINE(ITRA,VERT,PVERT,IPART,NVERT,UBUF,NUBUF)
 966          IFLPART=INT(UBUF(1))
 967          IF(IHIJ.EQ.2.OR.IVNU.EQ.2)WRITE(6,*)'ITRA,IFLPART ',ITRA,
 968      +   IFLPART
 969          ITOLDG=ITRA
 970       ENDIF
 971 *
 972       END
 973
 974 *CMZ :          15/02/99  14.40.30  by  Federico Carminati
 975 *-- Author :  Rene Brun
 976       SUBROUTINE ZDC_DEFS
 977 C
 978 C *** DEFINITION OF THE ZDC parameters
 979 C     CALLED BY : SXDEFS
 980 C
 981       end
 982
 983       subroutine zdc_setbeam(beam, fxoff, fyoff, sx, sy,
 984      +               div, angle, cross)
 985 * beam  : 1 = gaussian beam
 986 *       : 2 = uniform beam
 987 * fxoff : x-coordinate of beam offset
 988 * fyoff : y-coordinate of beam offset
 989 * sx    : sigma-x of the beam (gaussian or uniform)
 990 * sy    : sigma-y of the beam (gaussian or uniform)
 991 * div   : divergency of the beam (32*10**-6 rad for LHC)
 992 * angle : beam crossing angle (100*10**-6 rad for LHC)
 993 * cross : 1 = horizontal beam crossing
 994 *       : 2 = vertical beam crossing
 995 #include "zdc_common.inc"
 996       integer beam, cross
 997       real fxoff,fyoff,sx,sy,div,angle
 998       beamty = type
 999       fx     = fxoff
1000       fy     = fyoff
1001       sigx   = sx
1002       sigy   = sy
1003       bmdiv  = div
1004       bmcra  = angle
1005       iflcr  = cross
1006       end
1007
1008       subroutine zdc_sethijing(hij, hijf, hijsp, file)
1009 * IHIJ  : 1 = read HIJING event file
1010 *       : 2 =  "     "      "    "    + debug
1011 * IHIJF : event number of the first event to be read from file
1012 * IHIJSP: 0 = read all particles
1013 *       : 1 = remove spectator nucleons
1014 #include "zdc_common.inc"
1015       integer hij,hijf,hijsp
1016       character *(*) file
1017       ihij   = hij
1018       ihijf  = hijf
1019       ihijsp = hijsp
1020       fileh  = file
1021       end
1022
1023       subroutine zdc_setvenus(hiv, hivf, hivsp, file)
1024 * IHIV  : 1 = read VENUS event file
1025 *       : 2 =  "     "      "    "    + debug
1026 * IHIVF : event number of the first event to be read from file
1027 * IHIVSP: 0 = read all particles
1028 *       : 1 = remove spectator nucleons
1029 #include "zdc_common.inc"
1030       integer hiv,hivf,hivsp
1031       character *(*) file
1032       ihiv   = hiv
1033       ihivf  = hivf
1034       ihivsp = hivsp
1035       filev  = file
1036       end
1037
1038       subroutine zdc_setkine(code, pmom, cx, cy, cz, type, fermi)
1039 * code     : GEANT code of the test particle
1040 * pmom     : absolute value of particle momentum
1041 * cx,cy,cz : director cosines of the track (if type)
1042 * type     :  0 = take director cosines from cx,cy,cz)
1043 *          : <>0 = pseudorapidity of the test particle
1044 * fermi    : 0 = no Fermi motion for the spectator nucleons
1045 *          : 1 = Fermi motion for the spectator nucleons
1046 #undef CERNLIB_GEANT321_GCKINE_INC
1047 #include "geant321/gckine.inc"
1048 #include "zdc_common.inc"
1049       integer code,type,fermi
1050       real pmom,cx,cy,cz
1051       pkine(1)  = code
1052       pkine(2)  = pmom
1053       pkine(3)  = cx
1054       pkine(4)  = cy
1055       pkine(5)  = cz
1056       pkine(6)  = type
1057       ifermi    = fermi
1058       end
1059 *CMZ :          19/11/98  09.31.36  by  Federico Carminati
1060 *-- Author :    Federico Carminati   19/11/98
1061       SUBROUTINE SXGAUS(RVEC,N)
1062 #undef CERNLIB_GEANT321_GCONSP_INC
1063 #include "geant321/gconsp.inc"
1064 *KEND.
1065 *
1066 * Vector of random numbers distributed like a gaussian
1067       DIMENSION RNDM(2), RVEC(N)
1068       DO K=1,(N-1)/2+1
1069          CALL GRNDM(RNDM,2)
1070          PHI = TWOPI*RNDM(1)
1071          RHO = SQRT(-2.*LOG(RNDM(2)))
1072          RVEC(2*K-1) = RHO*SIN(PHI)
1073          IF(2*K.LE.N) RVEC(2*K  ) = RHO*COS(PHI)
1074       ENDDO
1075 *
1076       END