1 *CMZ : 19/11/98 21.32.03 by Federico Carminati
3 C***********************************************************************
5 SUBROUTINE ZDC_ADDANG(THETA1,PHI1,THETA2,PHI2,THETA3,PHI3)
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.
14 C 27-AUG-1992 S. Coutu
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
23 CT1=COS(THETA1*DEGRAD)
24 ST1=SIN(THETA1*DEGRAD)
27 CT2=COS(THETA2*DEGRAD)
28 ST2=SIN(THETA2*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
37 IF(THETA3.EQ.0.OR.THETA3.EQ.180.)THEN
42 IF(TEMP.LT.-1.)TEMP=-1.
43 PHI3=ACOS(TEMP)*RADDEG
44 IF(CY.LT.0.)PHI3=360.-PHI3
49 *CMZ : 21/11/98 17.33.14 by Federico Carminati
51 C==================================================================
52 SUBROUTINE ZDC_EFERMI(ID)
54 #include "zdc_common.inc"
55 #undef CERNLIB_GEANT321_GCONSP_INC
56 #include "geant321/gconsp.inc"
63 IF(RNDM(3).GE.PROBINTP(I-1).AND.RNDM(3).LT.PROBINTP(I))GO TO
66 ELSE IF(ID.EQ.13) THEN
68 IF(RNDM(3).GE.PROBINTN(I-1).AND.RNDM(3).LT.PROBINTN(I))GO TO
77 DDP(1)=PIPPO*SIN(TET)*COS(PHI)
78 DDP(2)=PIPPO*SIN(TET)*SIN(PHI)
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.
88 C *** TERMINATION OF THE ZDC SIMULATION AT THE END OF A RUN ***
91 C ORIGIN : E. SCOMPARIN
94 * CALL HROUT(0,ICYCLE,' ')
99 10000 FORMAT(1H /1H ,36('*'),' ZDC_END ',36('*')/
100 $ 1H ,20X,'ZDC simulation package terminated'/
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.
110 C *** TERMINATION OF THE ZDC SIMULATION AFTER EACH EVENT ***
113 C ORIGIN : E. SCOMPARIN
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"
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)
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)
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
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)
170 IF(ITRP(I).NE.ITCTP)THEN
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))
180 EZPH(ITCTP,J1,J2)=HITP(7,I)
181 NPPH(ITCTP,J1,J2)=HITP(6,I)
185 IF(ITRN(I).NE.ITCTN)THEN
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))
195 EZNH(ITCTN,J1,J2)=HITN(7,I)
196 NPNH(ITCTN,J1,J2)=HITN(6,I)
226 *CMZ : 21/11/98 17.33.14 by Federico Carminati
228 C==================================================================
229 SUBROUTINE ZDC_GFERMI(A,Z)
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
242 ALFA=0.18*((A/12.)**0.333)
243 XK=(4.*PI)/((1.+ALFA)*TWOPI**1.5)
248 E1=(P*P)/(2.*SIG1*SIG1)
249 E2=(P*P)/(2.*SIG2*SIG2)
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)
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
262 C-- Author : E.SCOMPARIN, 07/05/1996
265 C *** INITIALISATION FOR THE ZDC SIMULATION ***
268 C ORIGIN : E. SCOMPARIN
271 #undef CERNLIB_GEANT321_GCFLAG_INC
272 #include "geant321/gcflag.inc"
273 #undef CERNLIB_GEANT321_GCKINE_INC
274 #include "geant321/gckine.inc"
276 #include "zdc_common.inc"
278 * Initialize COMMON block ZDC_CGEOM
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/
292 * Initialize COMMON block ZDC_HITS
295 DATA ZSET,ZN1,ZP1 /'ZDC ', 'ZN1 ','ZP1 '/
296 DATA CHNMSZN/'ZNTX','ZN1 '/
297 DATA CHNMSZP/'ZPTX','ZP1 '/
300 DATA IDTYPZN,IDTYPZP /1, 2/
301 DATA NWHI,NWDI /1000, 1000/
303 * Data for Detectors ZN1,ZP1
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./
311 * Read tables of Cerenkov light produced in the fibers
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')
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)
327 10000 FORMAT(7(1X,F10.5))
335 * Set Fermi flag according to SWITCH 3
336 IF(ISWIT(3).EQ.1) IFERMI=1
337 * define the LHC energy/nucleon
340 * Open data file to read in particles from event generators
342 OPEN(27,FILE=FILEH,STATUS='OLD')
343 ELSE IF(IVNU.GE.1)THEN
344 OPEN(27,FILE=FILEV,STATUS='OLD')
346 C GET FERMI MOMENTUM DISTRIBUTIONS
347 CALL ZDC_GFERMI(207.,82.)
349 * CALL HROPEN (1,'ALIZ',FILEA,'N',1024,ISTAT)
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)')
360 10100 FORMAT(1H /1H ,36('*'),' ZDC_INIT ',36('*')/
361 + 1H ,20X,'ZDC simulation package initialised'/
366 *CMZ : 21/11/98 17.33.14 by Federico Carminati
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"
377 DIMENSION KATT(10000,2),PATT(10000,4)
379 INTEGER ITA(49),ITB(49)
380 REAL PLAB(3),UBUF(1),UB(1)
382 REAL BALP(4),DDDP(4),PLABS(4)
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/
396 CHARACTER*20 PARTNAME
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
407 VERTEX(1)= VERTEX(1)+FX+RNDM(1)*SIGX
408 VERTEX(2)= VERTEX(2)+FY+RNDM(2)*SIGY
410 ELSEIF(BEAMTY.EQ.2) THEN
412 VERTEX(1)= VERTEX(1)+FX+(2.*RNDM(1)-1.)*SIGX
413 VERTEX(2)= VERTEX(2)+FY+(2.*RNDM(2)-1.)*SIGY
417 * Fill the vertex bank
419 CALL GSVERT(VERTEX,0,0,0.,0,NVTX)
420 IF(NVTX.EQ.0)WRITE(6,*)'Error defining vertex'
422 * Read HIJING event file
425 10 READ(27,*,END=70) NMUL,EATT,JATT,NT,NP,N0,N01,N10,N11,B,NATT
427 WRITE(6,*)'HIJING EVENT HEADER= ',NMUL,EATT,JATT,NT,NP,
428 + N0,N01,N10,N11,B,NATT
430 READ(27,*) (KATT(I,1),KATT(I,2),(PATT(I,J),J=1,4),I=1,NATT)
432 IF(IVR.LT.IHIJF) GOTO 10
435 PARTI=FLOAT(NT)+FLOAT(NP)
439 * -- Remove particles with negative z-direction
440 IF(PATT(I,3).LT.0.0) GOTO 30
441 * -- Remove (if required) the spectators
443 IF(KATT(I,2).EQ.0..OR.KATT(I,2).EQ.10.) GOTO 30
445 * -- Define participant flag : UBUF(1)
446 IF(KATT(I,2).EQ.0..OR.KATT(I,2).EQ.10.) THEN
451 * -- Translate the particle ID from HIJING to GEANT
453 IF(KATT(I,1).EQ.ITA(J)) THEN
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)
468 * Fill the tracks bank
469 CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
471 WRITE(6,*)'IPART= ',IPART,'PLAB= ',PLAB
473 IF(NT.EQ.0)WRITE(6,*)'Error defining track'
478 * Read VENUS event file
480 ELSE IF(IVNU.GE.1)THEN
481 40 READ(27,*,END=80)LEVT,NREVT,NMUL,B,KOLEVT,COLEVT, PMXEVT,
484 WRITE(6,*) ' Error in reading VENUS data'
488 WRITE(6,*)'VENUS EVENT HEADER= ',LEVT,NREVT,NATT,B,KOLEVT,
489 + COLEVT,PMXEVT,EGYEVT,NP,NT
493 *-- read produced particles
494 READ(27,*,END=80) LPTL,NREVT,NRPTL,KATT(I,1), (PATT(I,J),J=
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
505 WRITE(6,*) ' ERROR IN READING VENUS DATA'
511 IF(IVR.LT.IVNUF) GOTO 40
514 PARTI=(FLOAT(NT)+FLOAT(NP))/2.
518 * -- remove particles in negative direction
519 IF(PATT(I,3).LT.0.0) GOTO 60
520 * -- if required remove the spectators
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)))
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
533 * -- Translate the particle ID from VENUS to GEANT
535 IF(KATT(I,1).EQ.ITB(J))THEN
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)
549 * Fill the tracks bank
550 CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
552 WRITE(6,*)'IPART= ',IPART,'PLAB= ',PLAB
554 IF(NT.EQ.0)WRITE(6,*)'Error defining track'
560 * Generation of test particles
562 * Definition of the incident particle
564 IPART = NINT(PKINE(1))
566 CALL GFPART(IPART,PARTNAME,ITRTYP,AMASS,CHARGE,TLIFE,UB,NB)
569 * --- In case of a ion , find out energy and momentum
570 * --- NAAP=A=Z+N - atomic number of the initial particle
572 NAAP = INT(0.5 + AMASS/0.93149432)
576 IF((IPART.GE.45).OR.(IPART.EQ.13).OR.(IPART.EQ.14)) THEN
582 ELSE IF(IPART.EQ.13)THEN
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)
593 SCANG=2*ATAN(EXP(-(PKINE(6))))
594 PLAB(1)=-PKINE(2)*SIN(SCANG)
596 PLAB(3)=PKINE(2)*COS(SCANG)
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)
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
609 GOLDHABER=SQRT(AAP*(207-AAP)/(207-1.))
611 DDP(JJ)=DDP(JJ)*GOLDHABER
617 BALP(4)=SQRT(BALP(1)**2+BALP(2)**2+BALP(3)**2+AMA**2)
621 DDDP(4)=SQRT(DDDP(1)**2+DDDP(2)**2+DDDP(3)**2+AMA**2)
622 CALL LORENF(AMA,BALP,DDDP,PLABS)
626 c WRITE(6,*)'BALP ',BALP,'DDDP ',DDDP,'PLABS ',PLABS
631 * Fill the tracks bank
632 CALL GSKINE(PLAB,IPART,NVTX,UBUF,1,NT)
633 IF(NT.EQ.0)WRITE(6,*)'Error defining track'
637 * Some debug printings
638 IF(IHIJ.EQ.2.OR.IVNU.EQ.2)THEN
647 70 WRITE(6,*) 'END OF HIJING FILE'
650 80 WRITE(6,*) 'END OF VENUS FILE'
655 *CMZ : 19/11/98 09.54.52 by Federico Carminati
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"
663 #undef CERNLIB_GEANT321_GCONSP_INC
664 #include "geant321/gconsp.inc"
665 #include "zdc_common.inc"
667 DOUBLE PRECISION TETPART,FIPART,TETDIV,FIDIV,TETSUM,FISUM
668 DOUBLE PRECISION DPLAB(3)
673 C TAKE NON-ZERO DIVERGENCY
674 c write(6,*)'PLAB (NODIV) = ',plab
678 PMOD=PMOD+PLAB(IJ)*PLAB(IJ)
682 C--- TETDIV is given in radiants and it is read by data cards
685 TETDIV=PKINE(10)*ABS(RVEC(1))
686 C--- PHIDIV is randomly chosen between 0 and 360
689 ELSE IF(ICROSS.EQ.1)THEN
690 IF(ISWIT(7).EQ.0)THEN
693 ELSE IF(ISWIT(7).EQ.1)THEN
696 ELSE IF(ISWIT(7).EQ.2)THEN
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))
707 IF (FIPART.LT.0)FIPART=FIPART+TWOPI
708 C write(6,*)' TETPART= ',tetpart,' FIPART= ',fipart
711 TETPART=TETPART*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
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
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.
730 C *** RECORDING OF ZDC HITS AFTER EACH TRACKING STEP ***
733 C ORIGIN : E. SCOMPARIN
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"
758 DIMENSION LNAM(15),LNUM(15)
761 DIMENSION VERT0(3),PVERT0(4)
762 CHARACTER*20 CHNPART0
764 DIMENSION XLOCM(3),XANGM(3),XLOCL(3),XANGL(3),FNULA(3),ZENTR(3)
769 IF(ITRA.NE.ITOLD)THEN
770 IF(INWVOL.EQ.1.AND.IDTYPE.GT.0)THEN
772 * Impact point on the ZDC (either ZPRO or ZNEU)
775 CALL GFPATH(ISET,IDET,NUMBV,NLEV,LNAM,LNUM)
776 CALL UCTOH('ZNEU',IZNEU,4,4)
777 CALL UCTOH('ZPRO',IZPRO,4,4)
779 IF(LNAM(I).EQ.IZNEU.OR.LNAM(I).EQ.IZPRO)THEN
781 CALL GLVOLU(I,LNAM,LNUM,IER)
783 CALL GMTOD(VECT,HITSBUPN(3),1)
786 ELSE IF(IDTYPE.EQ.2)THEN
787 CALL GMTOD(VECT,HITSBUPP(3),1)
797 * Particle type and energy for the primary track
799 CALL GFKINE(ITRA,VERT0,PVERT0,IPART0,NVERT0,UB,NB)
800 CALL GFPART(IPART0,CHNPART0,ITRTYP0,AMASS0,CHARGE0,TLIFE0,
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)
808 HITSBUPN(1)=FLOAT(IPART0)
812 HITSBUPP(1)=FLOAT(IPART)
815 HITSBUPN(1)=FLOAT(IPART)
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)
832 * Energy deposition in the passive material
834 IF(INWVOL.EQ.0.AND.IDTYPE.GT.0)THEN
843 HITSBUPN(KCOPY)=HITS(KCOPY)
845 ELSE IF(IDTYPE.EQ.2)THEN
847 HITSBUPP(KCOPY)=HITS(KCOPY)
850 CALL GSCHIT(ISET,IDET,ITRA,NUMBV,HITS,NHSUM,IHIT)
853 * Select charged relativistic particles
855 IF(ABS(CHARGE).EQ.0)RETURN
856 IF(GEKIN.GT.0.00001) THEN
861 IF(BETA.LT.0.67) RETURN
864 * Particle inside fiber
866 IF(NUMED.EQ.803.AND.IDTYPE.GT.0) THEN
875 * --- Calculating the position of the centre of the fiber in MARS
880 CALL GDTOM(FNULA,ZENTR,1)
882 * --- Angle between the particle trajectory and the fiber axis
884 CALL GMTOD(XANGM,XANGL,2)
885 IF(XANGL(3).GT.1.0) XANGL(3)=SIGN(1.0,XANGL(3))
889 IF (ALFAD.GT.ALFMAX) RETURN
890 IALFA = INT(1.+ALFAD/2.)
892 * --- Distance between the particle trajectory and the fiber axis
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.)
903 ELSEIF(IDTYPE.EQ.2)THEN
906 IF(BE.GT.RADIUS) THEN
907 PRINT*,'DISTANCE GREATER THAN FIBER RADIUS'
910 IBE = 1+INT(1000.*BE)
911 IF(IBE.GT.NBE)IBE=NBE
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
919 * --- Looking into the tables
920 OUT = CHARGE*CHARGE*TABLE(IBETA,IALFA,IBE)
921 CALL RNPSSN(OUT,N,IERROR)
925 CALL UCOPY(HITSBUPN,HITS,4)
926 ELSE IF(IDTYPE.EQ.2)THEN
927 CALL UCOPY(HITSBUPP,HITS,4)
929 CALL GSCHIT(ISET,IDET,ITRA,NUMBV,HITS,NHSUM,IHIT)
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)
939 C *** DEFINITION OF THE GEOMETRY OF THE PHOS ***
940 C *** NVE 24-SEP-1990 CERN GENEVA ***
943 C ORIGIN : E. SCOMPARIN
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"
963 * Fetch new track (VENUS or HIJING)
964 IF(ITRA.NE.ITOLDG)THEN
965 CALL GFKINE(ITRA,VERT,PVERT,IPART,NVERT,UBUF,NUBUF)
967 IF(IHIJ.EQ.2.OR.IVNU.EQ.2)WRITE(6,*)'ITRA,IFLPART ',ITRA,
974 *CMZ : 15/02/99 14.40.30 by Federico Carminati
975 *-- Author : Rene Brun
978 C *** DEFINITION OF THE ZDC parameters
983 subroutine zdc_setbeam(beam, fxoff, fyoff, sx, sy,
985 * beam : 1 = gaussian 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"
997 real fxoff,fyoff,sx,sy,div,angle
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
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
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
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"
1066 * Vector of random numbers distributed like a gaussian
1067 DIMENSION RNDM(2), RVEC(N)
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)