This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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