source: trunk/MagicSoft/Simulation/Corsika/Mmcs/cerenkov.f@ 694

Last change on this file since 694 was 687, checked in by harald, 24 years ago
Changes of Ciro and Dennis to write all events of one run in only one file. Be careful, from now on you need also a new reflector.
File size: 46.0 KB
Line 
1C=======================================================================
2
3 SUBROUTINE CERENE( STEPCR )
4
5C-----------------------------------------------------------------------
6C CEREN(KOV RADIATION FROM) E(LECTRONS)
7C
8C CREATION OF CERENKOV PHOTONS ALONG A TRACK OF ELECTRONS
9C CERENKOV RADIATION IS ONLY CALCULATED FOR LOWEST OBSERVATION LEVEL
10C THE COORDINATES ON EGS-STACK ARE AT THE END OF STEP EXCEPT E(NP),
11C WHICH IS AT THE BEGINNING OF STEP
12C THIS SUBROUTINE IS CALLED FROM ELECTR
13C ARGUMENT:
14C STEPCR = STEP LENGTH FOR ELECTRON OR POSITRON (REAL*4)
15C
16C AUTHOR : M. ROZANSKA UNIVERSITY OF KRAKOW
17C S. MARTINEZ UNIVERSITY OF MADRID
18C F. ARQUEROS UNIVERSITY OF MADRID
19C CHANGES : D. HECK IK3 FZK KARLSRUHE
20C R. ATTALLAH UNIVERSITY OF PERPIGNAN
21C-----------------------------------------------------------------------
22
23c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
24 parameter (xct=1)
25 parameter (yct=2)
26 parameter (zct=3)
27 parameter (ctthet=4)
28 parameter (ctphi=5)
29 parameter (ctdiam=6)
30 parameter (ctfoc=7)
31c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
32
33c IMPLICIT NONE
34*KEEP,BUFFS.
35 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
36 INTEGER MAXBUF,MAXLEN
37 PARAMETER (MAXBUF=39*7)
38 PARAMETER (MAXLEN=12)
39 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
40 * RUNE(MAXBUF),DATAB(MAXBUF)
41 INTEGER LH
42 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
43 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
44 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
45*KEEP,CONST.
46 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
47 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
48*KEEP,EPCONT.
49 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
50 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
51 * TSCAT,IAUSFL
52 DOUBLE PRECISION EDEP,RATIO
53 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
54 * EKE,ELKE,BETA2,GLE,TSCAT
55 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
56*KEEP,LONGI.
57 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
58 * NSTEP,LLONGI,FLGFIT
59 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
60 * SPLONG(0:1040,9),THSTEP,THSTPI
61 INTEGER NSTEP
62 LOGICAL LLONGI,FLGFIT
63*KEEP,MAGANG.
64 COMMON /MAGANG/ ARRANG,ARRANR,COSANG,SINANG
65 DOUBLE PRECISION ARRANG,ARRANR,COSANG,SINANG
66*KEEP,OBSPAR.
67 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
68 * THETPR,PHIPR,NOBSLV
69 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
70 * THETAP,THETPR(2),PHIP,PHIPR(2)
71 INTEGER NOBSLV
72*KEEP,PAM.
73 COMMON /PAM/ PAMA,SIGNUM
74 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
75*KEEP,PARPAR.
76 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
77 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
78 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
79 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
80 INTEGER ITYPE,LEVL
81*KEEP,RANDPA.
82 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
83 DOUBLE PRECISION FAC,U1,U2
84 REAL RD(3000)
85 INTEGER ISEED(103,10),NSEQ
86 LOGICAL KNOR
87*KEEP,RUNPAR.
88 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
89 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
90 * MONIOU,MDEBUG,NUCNUC,
91 * CETAPE,
92 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
93 * N1STTR,MDBASE,
94 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
95 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
96 * ,GHEISH,GHESIG
97 COMMON /RUNPAC/ DSN,HOST,USER
98 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
99 REAL STEPFC
100 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
101 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
102 * N1STTR,MDBASE
103 INTEGER CETAPE
104 CHARACTER*79 DSN
105 CHARACTER*20 HOST,USER
106
107 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
108 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
109 * ,GHEISH,GHESIG
110*KEEP,STACKE.
111 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
112 DOUBLE PRECISION E(60),TIME(60)
113 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
114 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
115*KEEP,CEREN1.
116 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
117 * CERSIZ,LCERFI
118 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
119 REAL CERSIZ
120 LOGICAL LCERFI
121*KEEP,CEREN2.
122 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
123 * DCERX,DCERY,ACERX,ACERY,
124 * XCMAX,YCMAX,EPSX,EPSY,
125 * DCERXI,DCERYI,FCERX,FCERY,
126 * XSCATT,YSCATT,CERXOS,CERYOS,
127 * NCERX,NCERY,ICERML
128 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
129 * DCERX,DCERY,ACERX,ACERY,
130 * XCMAX,YCMAX,EPSX,EPSY,
131 * DCERXI,DCERYI,FCERX,FCERY,
132 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
133 INTEGER NCERX,NCERY,ICERML
134c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
135*keep,certel.
136 common /certel/ cormxd,cord,coralp,ctpars,omega,
137 + photn,photnp,phpt,pht,vphot,
138 + vchi,veta,vzeta,vchim,vetam,vzetam,
139 + lambda,mu,nu,nctels,ncph
140 double precision cormxd,cord,coralp,ctpars(20,7),omega(20,3,3),
141 + photn(3),photnp(3),phpt(3),pht,vphot(3),
142 + vchi(3),veta(3),vzeta(3),vchim,vetam,vzetam,
143 + lambda,mu,nu
144 integer nctels,ncph(5)
145 double precision xg,yg,zg,xgp,ygp,zgp,up,vp,wp,xpcut,ypcut,zpcut
146 equivalence (photn(1) ,xg) ,(photn(2) ,yg) ,(photn(3) ,zg) ,
147 + (photnp(1),xgp) ,(photnp(2),ygp) ,(photnp(3),zgp),
148 + (phpt(1) ,xpcut),(phpt(2) ,ypcut),(phpt(3) ,zpcut),
149 + (vphot(1) ,up) ,(vphot(2) ,vp) ,(vphot(3) ,wp)
150 character *72 ctfile
151*keep,graal1.
152 common /graal1/ wavelength ! (nm)
153 real wavelength
154c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
155*KEND.
156
157 COMMON /ACLOCK/ NCLOCK,JCLOCK
158 DOUBLE PRECISION BETAE,BETAF,BETAI,CTHETA,DBETA,ECR,
159 * ETA1,E1,STHETA,STHETF,STHETI
160 DOUBLE PRECISION RHOF,THICK
161 REAL A,B,CC,COSCR,COSDEL,DVCOR,DXXX,DYYY,FSTEPI,
162 * HTOP,H2,PATHCR,PHICER,PHOTCT,RADINV,
163 * SINCR,SINDEL,SINPSI,SINPS2,STEPCR,UEMIS2,US,VCOR,
164 * VEMIS2,VS,WEMIS,XCER1,XCER2,XEMIS,XXX,
165 * YCER1,YCER2,YEMIS,YYY
166 INTEGER I,ISTC,I1,JCLOCK,LPCT1,NCLOCK,NSTEPC
167 EXTERNAL RHOF,THICK
168C-----------------------------------------------------------------------
169C_____IF (NCLOCK.GT.JCLOCK) THEN
170C______WRITE(MDEBUG,*)'CERENE: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
171C______CALL AUSGB2
172C_____ENDIF
173
174C-----------------------------------------------------------------------
175
176C SKIP PARTICLES OUT OF ZENITH ANGULAR CUT
177 IF ( W(NP) .LT. C(29) ) RETURN
178C E(NP) IS ENERGY AT BEGIN OF STEP
179 E1 = E(NP)
180
181C LOOK WETHER CERENKOV CONDITION IS FULFILLED AT BEGIN OF STEP
182 BETAI = SQRT( 1.D0 - (PAMA(2)*1.D3/E1)**2 )
183C REFRACTIVE INDEX PARAMETRISATION: N=1+ETA = ETA1
184 H2 = -Z(NP)
185 HTOP = H2 + VSTEP * W(NP)
186 ETA1 = 1.D0 + ETADSN * RHOF(DBLE(HTOP))
187 CTHETA = 1.D0 /( ETA1 * BETAI)
188 STHETI = 1.D0 - CTHETA**2
189 IF ( STHETI .GT. 0.D0 ) THEN
190C PARTICLE IS ABOVE ENERGY THRESHOLD IF EMISSION ANGLE IS >0
191 PHOTCT = CYIELD * STEPCR * STHETI
192 NSTEPC = PHOTCT / CERSIZ + 1
193 IF ( NSTEPC .LT. 1 ) RETURN
194 FSTEPI = 1. / REAL(NSTEPC)
195C CALCULATE INCREMENTS AND START VALUES FOR POSITION AND VELOCITY
196 DVCOR = -VSTEP * FSTEPI
197 VCOR = VSTEP - 0.5 * DVCOR
198 DBETA = -2.D0*FSTEPI*EDEP*(PAMA(2)*1.D3)**2 / (E1**3*BETAI)
199 BETAE = BETAI - 0.5D0 * DBETA
200 ELSE
201
202C LOOK WETHER CERENKOV CONDITION IS FULFILLED AT END OF STEP, BUT NOT
203C AT THE BEGINNING. THIS MAY HAPPEN ONLY ABOVE ABOUT 22 KM
204 IF ( HTOP .LT. 22.E5 ) RETURN
205C ENERGY AT END OF STEP IS ENEW (FROM COMMON EPCONT)
206 BETAF = SQRT( 1.D0 - (PAMA(2)*1.D3/ENEW)**2 )
207C REFRACTIVE INDEX PARAMETRISATION: N=1+ETA = ETA1
208 ETA1 = 1.D0 + ETADSN * RHOF( DBLE(H2) )
209 CTHETA = 1.D0 /( ETA1 * BETAF)
210 STHETF = 1.D0 - CTHETA**2
211C PARTICLE IS BELOW ENERGY THRESHOLD IF EMISSION ANGLE IS 0
212 IF ( STHETF .LE. 0.D0 ) RETURN
213 PHOTCT = CYIELD * STEPCR * STHETF
214 NSTEPC = PHOTCT / CERSIZ + 1
215 IF ( NSTEPC .LT. 1 ) RETURN
216 FSTEPI = 1. / REAL(NSTEPC)
217C CALCULATE INCREMENTS AND START VALUES FOR POSITION AND VELOCITY
218C LOOP 1000 RUNS FROM BOTTOM TO TOP OF STEP
219 DVCOR = VSTEP * FSTEPI
220 VCOR = -0.5 * DVCOR
221 DBETA = 2.D0*FSTEPI*EDEP*(PAMA(2)*1.D3)**2 / (ENEW**3*BETAF)
222 BETAE = BETAF - 0.5D0 * DBETA
223 ENDIF
224
225C LOOP OVER SUBSTEPS
226 DO 1000 ISTC = 1,NSTEPC
227 VCOR = VCOR + DVCOR
228 ZEMIS = H2 + VCOR * W(NP)
229 ETA1 = 1.D0 + ETADSN * RHOF(DBLE(ZEMIS))
230C VELOCITY IN THE MIDDLE OF SUBSTEP
231 BETAE = BETAE + DBETA
232 CTHETA = 1.D0 / (ETA1*BETAE)
233 STHETA = 1.D0 - CTHETA**2
234C PARTICLE IS AT ENERGY THRESHOLD IF EMISSION ANGLE BECOMES 0
235 IF ( STHETA .LE. 0.D0 ) RETURN
236C NUMBER OF EMITTED PHOTONS ON DISTANCE DVCOR
237 PHOTCM = CYIELD * STHETA * STEPCR * FSTEPI
238 STHETA = SQRT(STHETA)
239C ASSUME EMISSION POINT OF ALL PHOTONS IN THE MIDDLE OF THE STEP
240 XEMIS = X(NP) - VCOR * U(NP)
241 YEMIS = -Y(NP) + VCOR * V(NP)
242C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
243C GENERATE RANDOM WAVELENGTH FOR SINGLE C-PHOTON.
244 CALL RMMAR( RD,1,3 )
245 WAVELENGTH = 1. / (1/WAVLGL -
246 + RD(1)/(WAVLGL*WAVLGU/(WAVLGU-WAVLGL)))
247C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
248C CALCULATE PHOTON DIRECTION IN THE CORSIKA COORDINATE FRAME
249 CALL RMMAR( RD,1,3 )
250 PHICER = RD(1) * PI2
251 SINCR = SIN(PHICER)
252 COSCR = COS(PHICER)
253 A = U(NP)
254 B = -V(NP)
255 CC = W(NP)
256 SINPS2 = A**2 + B**2
257 IF ( SINPS2 .LT. 1.E-10 ) THEN
258 UEMIS = STHETA * COSCR
259 VEMIS = STHETA * SINCR
260 WEMIS = CTHETA * CC
261 ELSE
262 SINPSI = SQRT(SINPS2)
263 US = STHETA * COSCR
264 VS = STHETA * SINCR
265 SINDEL = B * (1./SINPSI)
266 COSDEL = A * (1./SINPSI)
267 UEMIS = CC * COSDEL * US - SINDEL * VS + A * CTHETA
268 VEMIS = CC * SINDEL * US + COSDEL * VS + B * CTHETA
269 WEMIS = -SINPSI * US + CC * CTHETA
270 ENDIF
271C EMISSION ANGLE WITHIN ZENITH ANGULAR CUT?
272 IF ( WEMIS .LT. C(29) ) GOTO 1000
273 RADINV = 1.5 - 0.5 * ( UEMIS**2 + VEMIS**2 + WEMIS**2 )
274 UEMIS2 = UEMIS * RADINV
275 VEMIS2 = VEMIS * RADINV
276 WEMIS = WEMIS * RADINV
277
278C CALCULATE DISTANCE FROM SHOWER AXIS AT THE DETECTOR LEVEL
279 PATHCR = ( ZEMIS - OBSLEV(NOBSLV) ) / WEMIS
280 XCER2 = XEMIS + PATHCR * UEMIS2 - XOFF(NOBSLV)
281 YCER2 = YEMIS + PATHCR * VEMIS2 - YOFF(NOBSLV)
282
283C ADD THE CERENKOV PHOTONS TO THE LONGITUDINAL DEVELOPMENT
284 IF ( LLONGI ) THEN
285C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
286 IF ( HLONG(NSTEP) .LE. ZEMIS ) THEN
287C FIND FIRST THE EQUIVALENT LEVELS
288 LPCT1 = LPCTE(NP)
289C ZEMIS IS ONLY LITTLE BELOW Z OLD, THEREFORE INCREMENTAL SEARCH
290C (REMEBER: LPCTE IS AT START OF ELECTRON STEP)
291 DO 6002 I1 = LPCT1,NSTEP
292 IF ( HLONG(I1) .LT. ZEMIS ) GOTO 6003
293 6002 CONTINUE
294 I1 = NSTEP + 1
295 6003 CONTINUE
296 DO 4862 I=I1,NSTEP
297 PLONG(I,9) = PLONG(I,9) + PHOTCM
298 4862 CONTINUE
299 ENDIF
300 ENDIF
301
302C TAKE INTO ACCOUNT A ROTATION OF ARRAY RELATIVE TO MAGNETIC NORD
303 XCER = XCER2 * COSANG + YCER2 * SINANG
304 YCER = YCER2 * COSANG - XCER2 * SINANG
305 UEMIS = UEMIS2 * COSANG + VEMIS2 * SINANG
306 VEMIS = VEMIS2 * COSANG - UEMIS2 * SINANG
307c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
308 CERELE = CERELE + PHOTCM
309 DO 7001 I=1,ICERML
310 DO 101 NCT=1,NCTELS
311c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
312c>> Modification to implement sphere algorithm >>>>>>>>>>>>>>>>>>>>>>>>
313c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
314c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
315c>>>>> this is the last (simple) check <<<<<
316C changes to ct frame
317c XG = XCER - CERXOS(I) - CTPARS(NCT,XCT)
318c YG = YCER - CERYOS(I) - CTPARS(NCT,YCT)
319c ZG = 0.0 - CTPARS(NCT,ZCT)
320c DIST2 = SQRT( XG**2 + YG**2 )
321c IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
322c>> New check >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
323 XG = XCER - CERXOS(I)
324 YG = YCER - CERYOS(I)
325 ZG = 0.0
326 DIST2 =
327 > SQRT((VEMIS*(-CTPARS(NCT,XCT) + XG) -
328 > UEMIS*(-CTPARS(NCT,YCT) + YG))**2 +
329 > (-(SQRT(1 - UEMIS**2 - VEMIS**2)*
330 > (-CTPARS(NCT,XCT) + XG)) +
331 > UEMIS*(-CTPARS(NCT,ZCT) + ZG))**2 +
332 > (SQRT(1 - UEMIS**2 - VEMIS**2)*
333 > (-CTPARS(NCT,YCT) + YG) -
334 > VEMIS*(-CTPARS(NCT,ZCT) + ZG))**2)
335 IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
336c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
337
338 101 CONTINUE
339
340 GOTO 7001
341
342C BUNCH FALLS ON A DETECTOR, CALCULATE ARRIVAL TIME (NSEC)
343
344 102 CARTIM = ((ETADSN*(THCKOB(NOBSLV)-THICK(DBLE(ZEMIS)))
345 * /WEMIS+PATHCR-VCOR/BETAE)/C(25)+TIME(NP))* 1.E9
346
347c CALL OUTPT2(IQ(NP),I)
348 CALL OUTPT2(NCT,I)
349C
350CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
351C
352C Modified by C. Bigongiari 2001 Jan 16
353C
354C
355Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
356C call jctime(cartim)
357Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
358C
359C
360CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
361C
362
363 GOTO 1000
364
365 7001 CONTINUE
366 1000 CONTINUE
367c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
368 RETURN
369 END
370
371C=======================================================================
372
373 SUBROUTINE CERENH( STEPCR,BETACR )
374
375C-----------------------------------------------------------------------
376C CEREN(KOV RADIATION FROM) H(ADRONS)
377C
378C CERENKOV RADIATION FROM HADRONS
379C CERENKOV RADIATION IS ONLY CALCULATED FOR LOWEST OBSERVATION LEVEL
380C THIS SUBROUTINE IS CALLED FROM UPDATE
381C ARGUMENTS:
382C STEPCR = STEP LENGTH FOR ELECTRON OR POSITRON
383C BETACR = VELOCITY OF PARTICLE IN UNITS OF SPEED OF LIGHT
384C
385C AUTHOR : M. ROZANSKA UNIVERSITY OF KRAKOW
386C S. MARTINEZ UNIVERSITY OF MADRID
387C F. ARQUEROS UNIVERSITY OF MADRID
388C CHANGES : D. HECK IK3 FZK KARLSRUHE
389C-----------------------------------------------------------------------
390
391c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
392 parameter (xct=1)
393 parameter (yct=2)
394 parameter (zct=3)
395 parameter (ctthet=4)
396 parameter (ctphi=5)
397 parameter (ctdiam=6)
398 parameter (ctfoc=7)
399c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
400
401c IMPLICIT NONE
402*KEEP,BUFFS.
403 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
404 INTEGER MAXBUF,MAXLEN
405 PARAMETER (MAXBUF=39*7)
406 PARAMETER (MAXLEN=12)
407 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
408 * RUNE(MAXBUF),DATAB(MAXBUF)
409 INTEGER LH
410 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
411 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
412 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
413*KEEP,CONST.
414 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
415 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
416*KEEP,EPCONT.
417 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
418 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
419 * TSCAT,IAUSFL
420 DOUBLE PRECISION EDEP,RATIO
421 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
422 * EKE,ELKE,BETA2,GLE,TSCAT
423 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
424*KEEP,LONGI.
425 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
426 * NSTEP,LLONGI,FLGFIT
427 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
428 * SPLONG(0:1040,9),THSTEP,THSTPI
429 INTEGER NSTEP
430 LOGICAL LLONGI,FLGFIT
431*KEEP,MAGANG.
432 COMMON /MAGANG/ ARRANG,ARRANR,COSANG,SINANG
433 DOUBLE PRECISION ARRANG,ARRANR,COSANG,SINANG
434*KEEP,OBSPAR.
435 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
436 * THETPR,PHIPR,NOBSLV
437 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
438 * THETAP,THETPR(2),PHIP,PHIPR(2)
439 INTEGER NOBSLV
440*KEEP,PAM.
441 COMMON /PAM/ PAMA,SIGNUM
442 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
443*KEEP,PARPAR.
444 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
445 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
446 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
447 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
448 INTEGER ITYPE,LEVL
449*KEEP,PARPAE.
450 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
451 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
452 * (CURPAR(4), PHI ), (CURPAR(5), H ),
453 * (CURPAR(6), T ), (CURPAR(7), X ),
454 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
455 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
456 * (CURPAR(12),ECM )
457*KEEP,RANDPA.
458 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
459 DOUBLE PRECISION FAC,U1,U2
460 REAL RD(3000)
461 INTEGER ISEED(103,10),NSEQ
462 LOGICAL KNOR
463*KEEP,RUNPAR.
464 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
465 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
466 * MONIOU,MDEBUG,NUCNUC,
467 * CETAPE,
468 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
469 * N1STTR,MDBASE,
470 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
471 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
472 * ,GHEISH,GHESIG
473 COMMON /RUNPAC/ DSN,HOST,USER
474 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
475 REAL STEPFC
476 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
477 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
478 * N1STTR,MDBASE
479 INTEGER CETAPE
480 CHARACTER*79 DSN
481 CHARACTER*20 HOST,USER
482
483 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
484 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
485 * ,GHEISH,GHESIG
486*KEEP,CEREN1.
487 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
488 * CERSIZ,LCERFI
489 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
490 REAL CERSIZ
491 LOGICAL LCERFI
492*KEEP,CEREN2.
493 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
494 * DCERX,DCERY,ACERX,ACERY,
495 * XCMAX,YCMAX,EPSX,EPSY,
496 * DCERXI,DCERYI,FCERX,FCERY,
497 * XSCATT,YSCATT,CERXOS,CERYOS,
498 * NCERX,NCERY,ICERML
499 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
500 * DCERX,DCERY,ACERX,ACERY,
501 * XCMAX,YCMAX,EPSX,EPSY,
502 * DCERXI,DCERYI,FCERX,FCERY,
503 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
504 INTEGER NCERX,NCERY,ICERML
505*KEEP,CERHDR.
506 COMMON/CERHDR/ TPART,UPART,VPART,WPART,XPART,YPART,ZPART
507 DOUBLE PRECISION TPART,UPART,VPART,WPART,XPART,YPART,ZPART
508c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
509*keep,certel.
510 common /certel/ cormxd,cord,coralp,ctpars,omega,
511 + photn,photnp,phpt,pht,vphot,
512 + vchi,veta,vzeta,vchim,vetam,vzetam,
513 + lambda,mu,nu,nctels,ncph
514 double precision cormxd,cord,coralp,ctpars(20,7),omega(20,3,3),
515 + photn(3),photnp(3),phpt(3),pht,vphot(3),
516 + vchi(3),veta(3),vzeta(3),vchim,vetam,vzetam,
517 + lambda,mu,nu
518 integer nctels,ncph(5)
519 double precision xg,yg,zg,xgp,ygp,zgp,up,vp,wp,xpcut,ypcut,zpcut
520 equivalence (photn(1) ,xg) ,(photn(2) ,yg) ,(photn(3) ,zg) ,
521 + (photnp(1),xgp) ,(photnp(2),ygp) ,(photnp(3),zgp),
522 + (phpt(1) ,xpcut),(phpt(2) ,ypcut),(phpt(3) ,zpcut),
523 + (vphot(1) ,up) ,(vphot(2) ,vp) ,(vphot(3) ,wp)
524 character *72 ctfile
525*keep,graal1.
526 common /graal1/ wavelength ! (nm)
527 real wavelength
528c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
529*KEND.
530
531 DOUBLE PRECISION BETACR,CINTEN,CTHETA,ETA,ETA1,HMID,PHOTCT,
532 * RHOF,STEPCR,STHETA,THICK
533 REAL A,B,CC,COSCR,COSDEL,DVCOR,DXXX,DYYY,FSTEPI,
534 * PATHCR,PHICER,RADINV,SINCR,SINDEL,SINPSI,SINPS2,
535 * UEMIS2,US,VEMIS2,VCOR,VS,WEMIS,XCER1,XCER2,XEMIS,
536 * XXX,YCER1,YCER2,YEMIS,YYY
537 INTEGER I,II,ISTC,I1,I2,NSTEPC
538 EXTERNAL RHOF,THICK
539C-----------------------------------------------------------------------
540
541c IF ( DEBUG ) WRITE(MDEBUG,*) 'CERENH: ZPART=',SNGL(ZPART),
542c * ' STEPCR=',SNGL(STEPCR),' BETACR=',SNGL(BETACR)
543
544C SKIP PARTICLE OUT OF ANGULAR ACCEPTANCE RANGE
545 IF ( WPART .LT. C(29) ) RETURN
546C CERENKOV INTENSITY FACTOR DEPENDS ON CHARGE STATE OF HEAVY IONS
547 CINTEN = CYIELD * ABS(SIGNUM(INT(CURPAR(1))))
548
549C REFRACTIVE INDEX PARAMETRISATION: N=1+ETA
550 HMID = ZPART + 0.5D0 * STEPCR * WPART
551 ETA1 = 1.D0 + ETADSN * RHOF(DBLE(HMID))
552 CTHETA = 1.D0 / ( ETA1 * BETACR )
553 STHETA = 1.D0 - CTHETA**2
554 IF ( STHETA .LE. 0.D0 ) RETURN
555
556 PHOTCT = CINTEN * STHETA * STEPCR
557 NSTEPC = PHOTCT / CERSIZ + 1
558 IF ( NSTEPC .LT. 1 ) RETURN
559 FSTEPI = 1. / REAL(NSTEPC)
560 VCOR = -0.5 * STEPCR * FSTEPI
561 DVCOR = -2. * VCOR
562C CERENKOV RADIATION IS ONLY CALCULATED FOR LOWEST OBSERVATION LEVEL
563 DO 1000 ISTC = 1,NSTEPC
564 VCOR = VCOR + DVCOR
565 ZEMIS = ZPART + VCOR * WPART
566 ETA = ETADSN * RHOF(DBLE(ZEMIS))
567 ETA1 = 1.D0 + ETA
568 CTHETA = 1.D0 / ( ETA1 * BETACR )
569 STHETA = 1.D0 - CTHETA**2
570 IF ( STHETA .LE. 0.D0 ) RETURN
571
572C NUMBER OF EMITTED PHOTONS ON DISTANCE STEPCR
573 PHOTCM = CINTEN * STHETA * STEPCR * FSTEPI
574 STHETA = SQRT(STHETA)
575
576C ASSUME EMISSION POINT OF ALL PHOTONS IN THE MIDDLE OF THE STEP
577C HAS TO BE CHECKED IF STEPS ARE NOT TOO LONG
578 XEMIS = XPART - VCOR * UPART
579 YEMIS = YPART - VCOR * VPART
580
581C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
582C GENERATE RANDOM WAVELENGTH FOR SINGLE C-PHOTON.
583 CALL RMMAR( RD,1,3 )
584 WAVELENGTH = 1. / (1/WAVLGL -
585 + RD(1)/(WAVLGL*WAVLGU/(WAVLGU-WAVLGL)))
586C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
587
588C CALCULATE PHOTON DIRECTION IN THE OVERALL COORDINATE FRAME
589 CALL RMMAR( RD,1,3 )
590 PHICER = RD(1) * PI2
591 SINCR = SIN(PHICER)
592 COSCR = COS(PHICER)
593 A = UPART
594 B = VPART
595 CC = WPART
596 SINPS2 = A**2 + B**2
597 IF ( SINPS2 .LT. 1.E-10 ) THEN
598 UEMIS = STHETA * COSCR
599 VEMIS = STHETA * SINCR
600 WEMIS = CTHETA * CC
601 ELSE
602 SINPSI = SQRT(SINPS2)
603 US = STHETA * COSCR
604 VS = STHETA * SINCR
605 SINDEL = B * (1./SINPSI)
606 COSDEL = A * (1./SINPSI)
607 UEMIS = CC * COSDEL * US - SINDEL * VS + A * CTHETA
608 VEMIS = CC * SINDEL * US + COSDEL * VS + B * CTHETA
609 WEMIS = -SINPSI * US + CC * CTHETA
610 ENDIF
611C EMISSION ANGLE WITHIN ZENITH ANGULAR CUT?
612 IF ( WEMIS .LT. C(29) ) GOTO 1000
613 RADINV = 1.5 - 0.5 * ( UEMIS**2 + VEMIS**2 + WEMIS**2 )
614 UEMIS2 = UEMIS * RADINV
615 VEMIS2 = VEMIS * RADINV
616 WEMIS = WEMIS * RADINV
617
618C CALCULATE DISTANCE FROM SHOWER AXIS AT THE DETECTOR LEVEL
619 PATHCR = ( ZEMIS - OBSLEV(NOBSLV) ) / WEMIS
620 XCER2 = XEMIS + PATHCR * UEMIS2 - XOFF(NOBSLV)
621 YCER2 = YEMIS + PATHCR * VEMIS2 - YOFF(NOBSLV)
622
623C ADD THE CERENKOV PHOTONS TO THE LONGITUDINAL DEVELOPMENT
624 IF ( LLONGI ) THEN
625C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
626 IF ( HLONG(NSTEP) .LE. ZEMIS ) THEN
627C FIND FIRST THE EQUIVALENT LEVELS
628 I1 = 0
629 I2 = NSTEP
630 6001 CONTINUE
631 II = (I1+I2)/2
632 IF ( HLONG(II) .LT. ZEMIS ) THEN
633 I2 = II
634 ELSE
635 I1 = II
636 ENDIF
637 IF ( I2-I1 .GT. 1 ) GOTO 6001
638 DO 4862 I=I2,NSTEP
639 PLONG(I,9) = PLONG(I,9) + PHOTCM
640 4862 CONTINUE
641 ENDIF
642 ENDIF
643
644C TAKE INTO ACCOUNT A ROTATION OF ARRAY RELATIVE TO MAGNETIC NORD
645 XCER = XCER2 * COSANG + YCER2 * SINANG
646 YCER = YCER2 * COSANG - XCER2 * SINANG
647 UEMIS = UEMIS2 * COSANG + VEMIS2 * SINANG
648 VEMIS = VEMIS2 * COSANG - UEMIS2 * SINANG
649c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
650 CERHAD = CERHAD + PHOTCM
651 DO 7001 I=1,ICERML
652 DO 101 NCT=1,NCTELS
653c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
654c>> Modification to implement sphere algorithm >>>>>>>>>>>>>>>>>>>>>>>>
655c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
656c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
657c>>>>> this is the last (simple) check <<<<<
658C changes to ct frame
659c XG = XCER - CERXOS(I) - CTPARS(NCT,XCT)
660c YG = YCER - CERYOS(I) - CTPARS(NCT,YCT)
661c ZG = 0.0 - CTPARS(NCT,ZCT)
662c DIST2 = SQRT( XG**2 + YG**2 )
663c IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
664c>> New check >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
665 XG = XCER - CERXOS(I)
666 YG = YCER - CERYOS(I)
667 ZG = 0.0
668 DIST2 =
669 > SQRT((VEMIS*(-CTPARS(NCT,XCT) + XG) -
670 > UEMIS*(-CTPARS(NCT,YCT) + YG))**2 +
671 > (-(SQRT(1 - UEMIS**2 - VEMIS**2)*
672 > (-CTPARS(NCT,XCT) + XG)) +
673 > UEMIS*(-CTPARS(NCT,ZCT) + ZG))**2 +
674 > (SQRT(1 - UEMIS**2 - VEMIS**2)*
675 > (-CTPARS(NCT,YCT) + YG) -
676 > VEMIS*(-CTPARS(NCT,ZCT) + ZG))**2)
677 IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
678c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
679 101 CONTINUE
680
681 GOTO 7001
682
683C BUNCH FALLS ON A DETECTOR, CALCULATE ARRIVAL TIME (NSEC)
684
685 102 CARTIM = ((ETADSN*(THCKOB(NOBSLV)-THICK(DBLE(ZEMIS)))
686 * /WEMIS+PATHCR-VCOR/BETACR)/C(25)+TPART)*1.E9
687
688c CALL OUTPT2(INT(CURPAR(1)),I)
689 CALL OUTPT2(NCT,I)
690
691C
692CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
693C
694C Modified by C. Bigongiari 2001 Jan 16
695C
696C
697Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
698C call jctime(cartim)
699Cc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
700C
701C
702CBC++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
703C
704
705 GOTO 1000
706
707 7001 CONTINUE
708 1000 CONTINUE
709c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
710 RETURN
711 END
712
713C=======================================================================
714
715 SUBROUTINE GETBUS( IPARTI,ENERGY,THETA,CERSZE )
716
717C-----------------------------------------------------------------------
718C GET BU(NCH) S(IZE)
719C
720C CALCULATES OPTIMAL BUNCH SIZE FOR CERENKOV PHOTONS. CERENKOV PHOTONS
721C ARE GROUPED IN BUNCHES IN ORDER TO ACCELERATE COMPUTING TIME.
722C HOWEVER, WE SET A MAXIMAL VALUE FOR THE GROUPING OF CERENKOV PHOTONS
723C SO THAT WE GET AT LEAST 100 BUNCHES/M**2 AT A CERENKOV FLUX OF 3000
724C PHOTONS/M**2. THIS IS THE MINIMUM CERENKOV FLUX WHICH CAN BE
725C DISTINGUISHED FROM THE NIGHT SKY LIGHT BACKGROUND IN THE HEGRA
726C EXPERIMENT AT THE ISLAND LA PALMA. SO THE PARAMETRIZATION OF THE
727C CERENKOV BUNCH AS CALCULATED IN THIS SUBROUTINE IS VALID FOR
728C OBSERVATION LEVELS SIMILAR TO THAT OF THE HEGRA EXPERIMENT.
729C FOR A GIVEN PRIMARY PARTICLE, INCIDENT ENERGY AND ANGLE, AN
730C OPTIMAL BUNCH SIZE IS CALCULATED BY INTERPOLATION IN A TABLE,
731C WHERE WE HAVE CHOSEN AN ENERGY RANGE UP TO 1000 TEV, INCIDENT
732C ANGLES 0 AND 40 DEGREES, AND 4 TYPES OF PRIMARIS: GAMMAS,
733C PROTONS, NITROGEN, AND IRON.
734C THIS SUBROUTINE IS CALLED FROM MAIN
735C ARGUMENTS:
736C IPARTI = TYPE OF PRIMARY PARTICLE
737C ENERGY (R4) = PARTICLES ENERGY IN GEV
738C THETA (R4) = ANGLE IN RAD
739C CERSZE (R4) = SIZE OF CERENKOV BUNCH
740C
741C AUTHORS : S. MARTINEZ UNIVERSITY OF MADRID
742C F. ARQUEROS UNIVERSITY OF MADRID
743C-----------------------------------------------------------------------
744
745 IMPLICIT NONE
746*KEEP,CONST.
747 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
748 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
749*KEEP,RUNPAR.
750 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
751 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
752 * MONIOU,MDEBUG,NUCNUC,
753 * CETAPE,
754 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
755 * N1STTR,MDBASE,
756 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
757 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
758 * ,GHEISH,GHESIG
759 COMMON /RUNPAC/ DSN,HOST,USER
760 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
761 REAL STEPFC
762 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
763 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
764 * N1STTR,MDBASE
765 INTEGER CETAPE
766 CHARACTER*79 DSN
767 CHARACTER*20 HOST,USER
768
769 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
770 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
771 * ,GHEISH,GHESIG
772*KEND.
773
774 REAL ANGLE(2),ENGAM(3),ENHAD(3),ENNIT(2),
775 * SIFE(3,2),SIGAM(3,2),SINIT(2),SIPRO(3,2)
776 REAL CERS1F,CERS1P,ENERGY,CERSZE,S1,S2,THETA
777 INTEGER I,IANFE,IANP,IATNUM,IPARTI,I1,I2
778
779 DATA ANGLE / 0., 40. /
780 DATA ENGAM / 100., 200., 500. /
781 DATA ENHAD / 100., 200., 1000. /
782 DATA ENNIT / 200., 1000. /
783 DATA ( SIFE (I,1),I=1,3 ) / 30., 30., 140. /
784 DATA ( SIFE (I,2),I=1,3 ) / 30., 30., 110. /
785 DATA ( SIGAM(I,1),I=1,3 ) / 30., 45., 100. /
786 DATA ( SIGAM(I,2),I=1,3 ) / 30., 40., 100. /
787 DATA SINIT / 30., 150. /
788 DATA ( SIPRO(I,1),I=1,3 ) / 30., 30., 120. /
789 DATA ( SIPRO(I,2),I=1,3 ) / 30., 30., 160. /
790 DATA IANP / 1 /, IANFE / 26 /
791C-----------------------------------------------------------------------
792
793c IF ( DEBUG ) WRITE(MDEBUG,100) IPARTI,SNGL(ENERGY),SNGL(THETA)
794c 100 FORMAT(' GETBUS: INPUT PARTICLE = ',I5,1P,2E10.3)
795
796C DEFAULT VALUE
797 CERSZE = 100.
798
799 ENERGY = ENERGY / 1000.
800 IF ( ENERGY .LE. 100. ) THEN
801 CERSZE = 30.
802 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
803 RETURN
804 ENDIF
805
806 THETA = THETA / PI * 180.
807
808C-----------------------------------------------------------------------
809C PHOTON, ELECTRON OR POSITRON AS PRIMARY PARTICLE
810 IF ( IPARTI .LE. 3 ) THEN
811C FIND ENERGY BIN FOR INTERPOLATION
812 IF ( ENERGY .LE. ENGAM(2) ) THEN
813 I1 = 1
814 I2 = 2
815 ELSE
816 I1 = 2
817 I2 = 3
818 ENDIF
819 S1 = SIGAM(I1,1) + (ENERGY - ENGAM(I1))
820 * / (ENGAM(I2) - ENGAM(I1))
821 * * (SIGAM(I2,1) - SIGAM(I1,1))
822 S2 = SIGAM(I1,2) + (ENERGY - ENGAM(I1))
823 * / (ENGAM(I2) - ENGAM(I1))
824 * * (SIGAM(I2,2) - SIGAM(I1,2))
825 CERSZE = S1 + (THETA-ANGLE(1))/(ANGLE(2)-ANGLE(1)) * (S2-S1)
826 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
827 RETURN
828 ENDIF
829
830C-----------------------------------------------------------------------
831C NITROGEN AS PRIMARY PARTICLE AND VERTICAL INCIDENCE
832CJOK WHY SPECIAL TREATMENT FOR NITROGEN ????
833CJOK WHY ONLY VERTICAL INCIDENCE ????
834 IF ( IPARTI .EQ. 1407 .AND. ABS(THETA) .LT. 1.E-1 ) THEN
835 IF ( ENERGY .LT. 200. ) THEN
836 CERSZE = 30.
837 ELSE
838 CERSZE = SINIT(1) + (ENERGY-ENNIT(1))
839 * / (ENNIT(2)-ENNIT(1)) * (SINIT(2)-SINIT(1))
840 ENDIF
841 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
842 RETURN
843 ENDIF
844
845C-----------------------------------------------------------------------
846C GET THE ATOMIC NUMBER OF THE NUCLEUS
847C Z IS 1, IF PROTON
848 IF ( IPARTI .EQ. 14 ) THEN
849 IATNUM = 1
850C REST OF POSSIBLE NUCLEI
851 ELSEIF ( IPARTI .GT. 100 ) THEN
852 IATNUM = MOD(IPARTI,100)
853 IF ( IATNUM .GT. 26 ) THEN
854 WRITE(MONIOU,*) 'GETBUS: UNEXPECTED PARTICLE CODE'
855 RETURN
856 ENDIF
857 ELSE
858 WRITE(MONIOU,*) 'GETBUS: UNEXPECTED PARTICLE CODE'
859 RETURN
860 ENDIF
861
862C FIND ENERGY BIN FOR INTERPOLATION IN CASE OF HADRONIC PRIMARY
863 IF ( ENERGY .LE. ENHAD(2) ) THEN
864 I1 = 1
865 I2 = 2
866 ELSE
867 I1 = 2
868 I2 = 3
869 ENDIF
870
871C INTERPOLATION FOR HADRONS
872 S1 = SIPRO(I1,1) + (ENERGY-ENHAD(I1))
873 * / (ENHAD(I2)-ENHAD(I1)) * (SIPRO(I2,1)-SIPRO(I1,1))
874 S2 = SIPRO(I1,2) + (ENERGY-ENHAD(I1))
875 * / (ENHAD(I2)-ENHAD(I1)) * (SIPRO(I2,2)-SIPRO(I1,2))
876 CERS1P = S1 + (THETA-ANGLE(1)) / (ANGLE(2)-ANGLE(1)) * (S2-S1)
877
878 S1 = SIFE(I1,1) + (ENERGY-ENHAD(I1)) / (ENHAD(I2)-ENHAD(I1))
879 * * (SIFE(I2,1)-SIFE(I1,1))
880 S2 = SIFE(I1,2) + (ENERGY-ENHAD(I1)) / (ENHAD(I2)-ENHAD(I1))
881 * * (SIFE(I2,2)-SIFE(I1,2))
882 CERS1F = S1 + (THETA-ANGLE(1)) / (ANGLE(2)-ANGLE(1)) * (S2-S1)
883
884 CERSZE = CERS1P + (IATNUM-IANP) * (CERS1F-CERS1P) / (IANFE-IANP)
885
886 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
887 101 FORMAT(' GETBUS: BUNCH SIZE = ',1P,1E10.3)
888
889 RETURN
890 END
891
892C=======================================================================
893
894 SUBROUTINE OUTND2
895
896C-----------------------------------------------------------------------
897C OUT(PUT AT E)ND (OF SHOWER)
898C
899C WRITE REST OF PARTICLES TO OUTPUT BUFFER
900C OUTND2 IS CALLED FROM MAIN
901C
902C AUTHORS : S. MARTINEZ, UNIVERSITY OF MADRID
903C F. ARQUEROS, UNIVERSITY OF MADRID
904C-----------------------------------------------------------------------
905
906 IMPLICIT NONE
907*KEEP,RUNPAR.
908 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
909 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
910 * MONIOU,MDEBUG,NUCNUC,
911 * CETAPE,
912 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
913 * N1STTR,MDBASE,
914 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
915 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
916 * ,GHEISH,GHESIG
917 COMMON /RUNPAC/ DSN,HOST,USER
918 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
919 REAL STEPFC
920 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
921 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
922 * N1STTR,MDBASE
923 INTEGER CETAPE
924 CHARACTER*79 DSN
925 CHARACTER*20 HOST,USER
926
927 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
928 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
929 * ,GHEISH,GHESIG
930*KEEP,CEREN2.
931 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
932 * DCERX,DCERY,ACERX,ACERY,
933 * XCMAX,YCMAX,EPSX,EPSY,
934 * DCERXI,DCERYI,FCERX,FCERY,
935 * XSCATT,YSCATT,CERXOS,CERYOS,
936 * NCERX,NCERY,ICERML
937 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
938 * DCERX,DCERY,ACERX,ACERY,
939 * XCMAX,YCMAX,EPSX,EPSY,
940 * DCERXI,DCERYI,FCERX,FCERY,
941 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
942 INTEGER NCERX,NCERY,ICERML
943*KEEP,CEREN3.
944 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
945 INTEGER MAXBF2
946 PARAMETER (MAXBF2 = 39 * 7)
947 DOUBLE PRECISION CERCNT
948 REAL DATAB2(MAXBF2)
949 INTEGER LHCER
950*KEND.
951
952 INTEGER I
953C-----------------------------------------------------------------------
954
955 IF ( LHCER .GT. 0 ) THEN
956 CALL TOBUFC( DATAB2,1 )
957 DO 2 I = 1,MAXBF2
958 DATAB2(I) = 0.
959 2 CONTINUE
960 ELSE
961 CALL TOBUFC( DATAB2,2 )
962 ENDIF
963
964 WRITE(MONIOU,*) 'CERCNT = ',SNGL( CERCNT )
965 CERCNT = 0.D0
966 LHCER = 0
967
968 RETURN
969 END
970
971C=======================================================================
972
973 SUBROUTINE OUTPT2(J,IMOV)
974
975C-----------------------------------------------------------------------
976C (WRITE CERENKOV RADIATION) OUTP(U)T
977C
978C OUTPUT ROUTINE FOR CERENKOV PHOTONS
979C THIS SUBROUTINE IS CALLED FROM CERENE AND CERENH
980C
981C AUTHORS : S. MARTINEZ, UNIVERSITY OF MADRID
982C F. ARQUEROS, UNIVERSITY OF MADRID
983C-----------------------------------------------------------------------
984
985 IMPLICIT NONE
986*KEEP,BUFFS.
987 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
988 INTEGER MAXBUF,MAXLEN
989 PARAMETER (MAXBUF=39*7)
990 PARAMETER (MAXLEN=12)
991 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
992 * RUNE(MAXBUF),DATAB(MAXBUF)
993 INTEGER LH
994 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
995 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
996 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
997*KEEP,RUNPAR.
998 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
999 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
1000 * MONIOU,MDEBUG,NUCNUC,
1001 * CETAPE,
1002 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1003 * N1STTR,MDBASE,
1004 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1005 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1006 * ,GHEISH,GHESIG
1007 COMMON /RUNPAC/ DSN,HOST,USER
1008 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
1009 REAL STEPFC
1010 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
1011 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1012 * N1STTR,MDBASE
1013 INTEGER CETAPE
1014 CHARACTER*79 DSN
1015 CHARACTER*20 HOST,USER
1016
1017 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1018 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1019 * ,GHEISH,GHESIG
1020*KEEP,CEREN1.
1021 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
1022 * CERSIZ,LCERFI
1023 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
1024 REAL CERSIZ
1025 LOGICAL LCERFI
1026*KEEP,CEREN2.
1027 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
1028 * DCERX,DCERY,ACERX,ACERY,
1029 * XCMAX,YCMAX,EPSX,EPSY,
1030 * DCERXI,DCERYI,FCERX,FCERY,
1031 * XSCATT,YSCATT,CERXOS,CERYOS,
1032 * NCERX,NCERY,ICERML
1033 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
1034 * DCERX,DCERY,ACERX,ACERY,
1035 * XCMAX,YCMAX,EPSX,EPSY,
1036 * DCERXI,DCERYI,FCERX,FCERY,
1037 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
1038 INTEGER NCERX,NCERY,ICERML
1039*KEEP,CEREN3.
1040 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
1041 INTEGER MAXBF2
1042 PARAMETER (MAXBF2 = 39 * 7)
1043 DOUBLE PRECISION CERCNT
1044 REAL DATAB2(MAXBF2)
1045 INTEGER LHCER
1046c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1047 COMMON /GRAAL1/ WAVELENGTH ! (NM)
1048 REAL WAVELENGTH
1049c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1050*KEND.
1051
1052 INTEGER I,J,IMOV
1053C-----------------------------------------------------------------------
1054
1055 IF(DEBUG)WRITE(MDEBUG,3)PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS
1056 3 FORMAT(' OUTPT2: ',1P,8E10.3)
1057C WRITE A BLOCK OF 39 PARTICLES TO THE CERENKOV OUTPUT BUFFER AND
1058C CLEAR FIELD
1059 IF ( LCERFI ) THEN
1060c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1061c DATAB2(LHCER+1) = PHOTCM
1062cc DATAB2(LHCER+1) = WAVELENGTH + J*1000.
1063 DATAB2(LHCER+1) = J*100000. + IMOV*1000. + WAVELENGTH
1064c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1065 DATAB2(LHCER+2) = XCER
1066 DATAB2(LHCER+3) = YCER
1067 DATAB2(LHCER+4) = UEMIS
1068 DATAB2(LHCER+5) = VEMIS
1069 DATAB2(LHCER+6) = CARTIM
1070 DATAB2(LHCER+7) = ZEMIS
1071c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1072c CERCNT = CERCNT + DBLE( PHOTCM )
1073 CERCNT = CERCNT + 1
1074c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1075 LHCER = LHCER + 7
1076 IF ( LHCER .GE. MAXBF2 ) THEN
1077 CALL TOBUFC( DATAB2,0 )
1078 DO 1 I = 1,MAXBF2
1079 DATAB2(I) = 0.
1080 1 CONTINUE
1081 LHCER = 0
1082 ENDIF
1083 ELSE
1084C WRITE A BLOCK OF 39 PARTICLES TO THE PARTICLE OUTPUT BUFFER AND
1085C CLEAR FIELD
1086 DATAB(LH+1) = 99.E5 + NINT(PHOTCM)*10. + 1.
1087 DATAB(LH+2) = XCER
1088 DATAB(LH+3) = YCER
1089 DATAB(LH+4) = UEMIS
1090 DATAB(LH+5) = VEMIS
1091 DATAB(LH+6) = CARTIM
1092 DATAB(LH+7) = ZEMIS
1093 LH = LH + 7
1094 NOPART = NOPART + 1
1095c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1096c CERCNT = CERCNT + DBLE( PHOTCM )
1097 CERCNT = CERCNT + 1
1098c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1099 IF ( LH .GE. MAXBUF ) THEN
1100 CALL TOBUF( DATAB,0 )
1101 DO 2 I = 1,MAXBUF
1102 DATAB(I) = 0.
1103 2 CONTINUE
1104 LH = 0
1105 ENDIF
1106 ENDIF
1107
1108 RETURN
1109 END
1110
1111C=======================================================================
1112
1113 SUBROUTINE TOBUFC( A,IFL )
1114
1115C-----------------------------------------------------------------------
1116C (WRITE) TO BUF(FER) C(ERENKOV DATA)
1117C
1118C COPY TO BUFFER CERENKOV DATA
1119C THIS SUBROUTINE IS CALLED FROM MAIN, INPRM, ELECTR, PHOTON, OUTND2,
1120C AND OUTPT2
1121C ARGUMENTS:
1122C A = ARRAY TO BE WRITTEN TO TAPE
1123C IFL = STARTING OF FINAL OUTPUT
1124C = 0 NORMAL BLOCK
1125C = 1 NORMAL BLOCK WITH END OF OUTPUT
1126C = 2 ONLY END OF OUTPUT
1127C-----------------------------------------------------------------------
1128
1129 IMPLICIT NONE
1130*KEEP,BUFFS.
1131 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
1132 INTEGER MAXBUF,MAXLEN
1133 PARAMETER (MAXBUF=39*7)
1134 PARAMETER (MAXLEN=12)
1135 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
1136 * RUNE(MAXBUF),DATAB(MAXBUF)
1137 INTEGER LH
1138 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
1139 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
1140 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
1141*KEEP,RECORD.
1142 COMMON /RECORD/ IRECOR
1143 INTEGER IRECOR
1144*KEEP,RUNPAR.
1145 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
1146 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
1147 * MONIOU,MDEBUG,NUCNUC,
1148 * CETAPE,
1149 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1150 * N1STTR,MDBASE,
1151 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1152 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1153 * ,GHEISH,GHESIG
1154 COMMON /RUNPAC/ DSN,HOST,USER
1155 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
1156 REAL STEPFC
1157 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
1158 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1159 * N1STTR,MDBASE
1160 INTEGER CETAPE
1161 CHARACTER*79 DSN
1162 CHARACTER*20 HOST,USER
1163
1164 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1165 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1166 * ,GHEISH,GHESIG
1167*KEEP,CEREN3.
1168 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
1169 INTEGER MAXBF2
1170 PARAMETER (MAXBF2 = 39 * 7)
1171 DOUBLE PRECISION CERCNT
1172 REAL DATAB2(MAXBF2)
1173 INTEGER LHCER
1174*KEEP,CEREN4.
1175 COMMON /CEREN4/ NRECER
1176 INTEGER NRECER
1177*KEND.
1178
1179 INTEGER NSUBBL
1180 PARAMETER (NSUBBL=21)
1181 REAL A(*)
1182C NSUBBL IS NUMBER OF SUBBLOCKS IN ONE OUTPUT RECORD
1183C (OUTPUT RECORD LENGTH = NSUBBL * 39 * 7 * 4 BYTES <= 22932 )
1184C IBLK2 IS COUNTER FOR SUBBLOCKS OF CERENKOV OUTPUT
1185C OUTPUT BUFFER FOR CERENKOV OUTPUT
1186 REAL OUTBF2(MAXBF2,NSUBBL)
1187 SAVE OUTBF2
1188 INTEGER I,IBLK2,IFL,K
1189 DATA IBLK2 / 0 /
1190C-----------------------------------------------------------------------
1191
1192 IF ( IFL .LE. 1 ) THEN
1193 IBLK2 = IBLK2 + 1
1194 DO 3 I = 1,MAXBF2
1195 OUTBF2(I,IBLK2) = A(I)
1196 3 CONTINUE
1197 ENDIF
1198
1199C WRITE TO TAPE IF BLOCK IS FULL OR IF IFL IS 1
1200 IF ( IFL .GE. 1 .OR. IBLK2 .EQ. NSUBBL ) THEN
1201 NRECER = NRECER + 1
1202
1203c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1204c WRITE(CETAPE) ((OUTBF2(I,K),I=1,MAXBF2),K=1,NSUBBL)
1205 call jccersave(outbf2)
1206c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1207 IBLK2 = 0
1208 DO 4 K = 1,NSUBBL
1209 DO 4 I = 1,MAXBF2
1210 OUTBF2(I,K) = 0.0
1211 4 CONTINUE
1212 ENDIF
1213
1214 RETURN
1215 END
Note: See TracBrowser for help on using the repository browser.