source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/cerenkov.f@ 10099

Last change on this file since 10099 was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 45.6 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)
349
350c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
351 call jctime(cartim)
352c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
353
354 GOTO 1000
355
356 7001 CONTINUE
357 1000 CONTINUE
358c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
359 RETURN
360 END
361
362C=======================================================================
363
364 SUBROUTINE CERENH( STEPCR,BETACR )
365
366C-----------------------------------------------------------------------
367C CEREN(KOV RADIATION FROM) H(ADRONS)
368C
369C CERENKOV RADIATION FROM HADRONS
370C CERENKOV RADIATION IS ONLY CALCULATED FOR LOWEST OBSERVATION LEVEL
371C THIS SUBROUTINE IS CALLED FROM UPDATE
372C ARGUMENTS:
373C STEPCR = STEP LENGTH FOR ELECTRON OR POSITRON
374C BETACR = VELOCITY OF PARTICLE IN UNITS OF SPEED OF LIGHT
375C
376C AUTHOR : M. ROZANSKA UNIVERSITY OF KRAKOW
377C S. MARTINEZ UNIVERSITY OF MADRID
378C F. ARQUEROS UNIVERSITY OF MADRID
379C CHANGES : D. HECK IK3 FZK KARLSRUHE
380C-----------------------------------------------------------------------
381
382c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
383 parameter (xct=1)
384 parameter (yct=2)
385 parameter (zct=3)
386 parameter (ctthet=4)
387 parameter (ctphi=5)
388 parameter (ctdiam=6)
389 parameter (ctfoc=7)
390c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
391
392c IMPLICIT NONE
393*KEEP,BUFFS.
394 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
395 INTEGER MAXBUF,MAXLEN
396 PARAMETER (MAXBUF=39*7)
397 PARAMETER (MAXLEN=12)
398 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
399 * RUNE(MAXBUF),DATAB(MAXBUF)
400 INTEGER LH
401 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
402 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
403 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
404*KEEP,CONST.
405 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
406 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
407*KEEP,EPCONT.
408 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
409 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
410 * TSCAT,IAUSFL
411 DOUBLE PRECISION EDEP,RATIO
412 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
413 * EKE,ELKE,BETA2,GLE,TSCAT
414 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
415*KEEP,LONGI.
416 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
417 * NSTEP,LLONGI,FLGFIT
418 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
419 * SPLONG(0:1040,9),THSTEP,THSTPI
420 INTEGER NSTEP
421 LOGICAL LLONGI,FLGFIT
422*KEEP,MAGANG.
423 COMMON /MAGANG/ ARRANG,ARRANR,COSANG,SINANG
424 DOUBLE PRECISION ARRANG,ARRANR,COSANG,SINANG
425*KEEP,OBSPAR.
426 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
427 * THETPR,PHIPR,NOBSLV
428 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
429 * THETAP,THETPR(2),PHIP,PHIPR(2)
430 INTEGER NOBSLV
431*KEEP,PAM.
432 COMMON /PAM/ PAMA,SIGNUM
433 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
434*KEEP,PARPAR.
435 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
436 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
437 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
438 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
439 INTEGER ITYPE,LEVL
440*KEEP,PARPAE.
441 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
442 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
443 * (CURPAR(4), PHI ), (CURPAR(5), H ),
444 * (CURPAR(6), T ), (CURPAR(7), X ),
445 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
446 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
447 * (CURPAR(12),ECM )
448*KEEP,RANDPA.
449 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
450 DOUBLE PRECISION FAC,U1,U2
451 REAL RD(3000)
452 INTEGER ISEED(103,10),NSEQ
453 LOGICAL KNOR
454*KEEP,RUNPAR.
455 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
456 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
457 * MONIOU,MDEBUG,NUCNUC,
458 * CETAPE,
459 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
460 * N1STTR,MDBASE,
461 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
462 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
463 * ,GHEISH,GHESIG
464 COMMON /RUNPAC/ DSN,HOST,USER
465 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
466 REAL STEPFC
467 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
468 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
469 * N1STTR,MDBASE
470 INTEGER CETAPE
471 CHARACTER*79 DSN
472 CHARACTER*20 HOST,USER
473
474 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
475 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
476 * ,GHEISH,GHESIG
477*KEEP,CEREN1.
478 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
479 * CERSIZ,LCERFI
480 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
481 REAL CERSIZ
482 LOGICAL LCERFI
483*KEEP,CEREN2.
484 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
485 * DCERX,DCERY,ACERX,ACERY,
486 * XCMAX,YCMAX,EPSX,EPSY,
487 * DCERXI,DCERYI,FCERX,FCERY,
488 * XSCATT,YSCATT,CERXOS,CERYOS,
489 * NCERX,NCERY,ICERML
490 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
491 * DCERX,DCERY,ACERX,ACERY,
492 * XCMAX,YCMAX,EPSX,EPSY,
493 * DCERXI,DCERYI,FCERX,FCERY,
494 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
495 INTEGER NCERX,NCERY,ICERML
496*KEEP,CERHDR.
497 COMMON/CERHDR/ TPART,UPART,VPART,WPART,XPART,YPART,ZPART
498 DOUBLE PRECISION TPART,UPART,VPART,WPART,XPART,YPART,ZPART
499c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
500*keep,certel.
501 common /certel/ cormxd,cord,coralp,ctpars,omega,
502 + photn,photnp,phpt,pht,vphot,
503 + vchi,veta,vzeta,vchim,vetam,vzetam,
504 + lambda,mu,nu,nctels,ncph
505 double precision cormxd,cord,coralp,ctpars(20,7),omega(20,3,3),
506 + photn(3),photnp(3),phpt(3),pht,vphot(3),
507 + vchi(3),veta(3),vzeta(3),vchim,vetam,vzetam,
508 + lambda,mu,nu
509 integer nctels,ncph(5)
510 double precision xg,yg,zg,xgp,ygp,zgp,up,vp,wp,xpcut,ypcut,zpcut
511 equivalence (photn(1) ,xg) ,(photn(2) ,yg) ,(photn(3) ,zg) ,
512 + (photnp(1),xgp) ,(photnp(2),ygp) ,(photnp(3),zgp),
513 + (phpt(1) ,xpcut),(phpt(2) ,ypcut),(phpt(3) ,zpcut),
514 + (vphot(1) ,up) ,(vphot(2) ,vp) ,(vphot(3) ,wp)
515 character *72 ctfile
516*keep,graal1.
517 common /graal1/ wavelength ! (nm)
518 real wavelength
519c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
520*KEND.
521
522 DOUBLE PRECISION BETACR,CINTEN,CTHETA,ETA,ETA1,HMID,PHOTCT,
523 * RHOF,STEPCR,STHETA,THICK
524 REAL A,B,CC,COSCR,COSDEL,DVCOR,DXXX,DYYY,FSTEPI,
525 * PATHCR,PHICER,RADINV,SINCR,SINDEL,SINPSI,SINPS2,
526 * UEMIS2,US,VEMIS2,VCOR,VS,WEMIS,XCER1,XCER2,XEMIS,
527 * XXX,YCER1,YCER2,YEMIS,YYY
528 INTEGER I,II,ISTC,I1,I2,NSTEPC
529 EXTERNAL RHOF,THICK
530C-----------------------------------------------------------------------
531
532c IF ( DEBUG ) WRITE(MDEBUG,*) 'CERENH: ZPART=',SNGL(ZPART),
533c * ' STEPCR=',SNGL(STEPCR),' BETACR=',SNGL(BETACR)
534
535C SKIP PARTICLE OUT OF ANGULAR ACCEPTANCE RANGE
536 IF ( WPART .LT. C(29) ) RETURN
537C CERENKOV INTENSITY FACTOR DEPENDS ON CHARGE STATE OF HEAVY IONS
538 CINTEN = CYIELD * ABS(SIGNUM(INT(CURPAR(1))))
539
540C REFRACTIVE INDEX PARAMETRISATION: N=1+ETA
541 HMID = ZPART + 0.5D0 * STEPCR * WPART
542 ETA1 = 1.D0 + ETADSN * RHOF(DBLE(HMID))
543 CTHETA = 1.D0 / ( ETA1 * BETACR )
544 STHETA = 1.D0 - CTHETA**2
545 IF ( STHETA .LE. 0.D0 ) RETURN
546
547 PHOTCT = CINTEN * STHETA * STEPCR
548 NSTEPC = PHOTCT / CERSIZ + 1
549 IF ( NSTEPC .LT. 1 ) RETURN
550 FSTEPI = 1. / REAL(NSTEPC)
551 VCOR = -0.5 * STEPCR * FSTEPI
552 DVCOR = -2. * VCOR
553C CERENKOV RADIATION IS ONLY CALCULATED FOR LOWEST OBSERVATION LEVEL
554 DO 1000 ISTC = 1,NSTEPC
555 VCOR = VCOR + DVCOR
556 ZEMIS = ZPART + VCOR * WPART
557 ETA = ETADSN * RHOF(DBLE(ZEMIS))
558 ETA1 = 1.D0 + ETA
559 CTHETA = 1.D0 / ( ETA1 * BETACR )
560 STHETA = 1.D0 - CTHETA**2
561 IF ( STHETA .LE. 0.D0 ) RETURN
562
563C NUMBER OF EMITTED PHOTONS ON DISTANCE STEPCR
564 PHOTCM = CINTEN * STHETA * STEPCR * FSTEPI
565 STHETA = SQRT(STHETA)
566
567C ASSUME EMISSION POINT OF ALL PHOTONS IN THE MIDDLE OF THE STEP
568C HAS TO BE CHECKED IF STEPS ARE NOT TOO LONG
569 XEMIS = XPART - VCOR * UPART
570 YEMIS = YPART - VCOR * VPART
571
572C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
573C GENERATE RANDOM WAVELENGTH FOR SINGLE C-PHOTON.
574 CALL RMMAR( RD,1,3 )
575 WAVELENGTH = 1. / (1/WAVLGL -
576 + RD(1)/(WAVLGL*WAVLGU/(WAVLGU-WAVLGL)))
577C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
578
579C CALCULATE PHOTON DIRECTION IN THE OVERALL COORDINATE FRAME
580 CALL RMMAR( RD,1,3 )
581 PHICER = RD(1) * PI2
582 SINCR = SIN(PHICER)
583 COSCR = COS(PHICER)
584 A = UPART
585 B = VPART
586 CC = WPART
587 SINPS2 = A**2 + B**2
588 IF ( SINPS2 .LT. 1.E-10 ) THEN
589 UEMIS = STHETA * COSCR
590 VEMIS = STHETA * SINCR
591 WEMIS = CTHETA * CC
592 ELSE
593 SINPSI = SQRT(SINPS2)
594 US = STHETA * COSCR
595 VS = STHETA * SINCR
596 SINDEL = B * (1./SINPSI)
597 COSDEL = A * (1./SINPSI)
598 UEMIS = CC * COSDEL * US - SINDEL * VS + A * CTHETA
599 VEMIS = CC * SINDEL * US + COSDEL * VS + B * CTHETA
600 WEMIS = -SINPSI * US + CC * CTHETA
601 ENDIF
602C EMISSION ANGLE WITHIN ZENITH ANGULAR CUT?
603 IF ( WEMIS .LT. C(29) ) GOTO 1000
604 RADINV = 1.5 - 0.5 * ( UEMIS**2 + VEMIS**2 + WEMIS**2 )
605 UEMIS2 = UEMIS * RADINV
606 VEMIS2 = VEMIS * RADINV
607 WEMIS = WEMIS * RADINV
608
609C CALCULATE DISTANCE FROM SHOWER AXIS AT THE DETECTOR LEVEL
610 PATHCR = ( ZEMIS - OBSLEV(NOBSLV) ) / WEMIS
611 XCER2 = XEMIS + PATHCR * UEMIS2 - XOFF(NOBSLV)
612 YCER2 = YEMIS + PATHCR * VEMIS2 - YOFF(NOBSLV)
613
614C ADD THE CERENKOV PHOTONS TO THE LONGITUDINAL DEVELOPMENT
615 IF ( LLONGI ) THEN
616C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
617 IF ( HLONG(NSTEP) .LE. ZEMIS ) THEN
618C FIND FIRST THE EQUIVALENT LEVELS
619 I1 = 0
620 I2 = NSTEP
621 6001 CONTINUE
622 II = (I1+I2)/2
623 IF ( HLONG(II) .LT. ZEMIS ) THEN
624 I2 = II
625 ELSE
626 I1 = II
627 ENDIF
628 IF ( I2-I1 .GT. 1 ) GOTO 6001
629 DO 4862 I=I2,NSTEP
630 PLONG(I,9) = PLONG(I,9) + PHOTCM
631 4862 CONTINUE
632 ENDIF
633 ENDIF
634
635C TAKE INTO ACCOUNT A ROTATION OF ARRAY RELATIVE TO MAGNETIC NORD
636 XCER = XCER2 * COSANG + YCER2 * SINANG
637 YCER = YCER2 * COSANG - XCER2 * SINANG
638 UEMIS = UEMIS2 * COSANG + VEMIS2 * SINANG
639 VEMIS = VEMIS2 * COSANG - UEMIS2 * SINANG
640c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
641 CERHAD = CERHAD + PHOTCM
642 DO 7001 I=1,ICERML
643 DO 101 NCT=1,NCTELS
644c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
645c>> Modification to implement sphere algorithm >>>>>>>>>>>>>>>>>>>>>>>>
646c>> JCG Wed Sep 21 10:49:14 MET DST 1998 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
647c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
648c>>>>> this is the last (simple) check <<<<<
649C changes to ct frame
650c XG = XCER - CERXOS(I) - CTPARS(NCT,XCT)
651c YG = YCER - CERYOS(I) - CTPARS(NCT,YCT)
652c ZG = 0.0 - CTPARS(NCT,ZCT)
653c DIST2 = SQRT( XG**2 + YG**2 )
654c IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
655c>> New check >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
656 XG = XCER - CERXOS(I)
657 YG = YCER - CERYOS(I)
658 ZG = 0.0
659 DIST2 =
660 > SQRT((VEMIS*(-CTPARS(NCT,XCT) + XG) -
661 > UEMIS*(-CTPARS(NCT,YCT) + YG))**2 +
662 > (-(SQRT(1 - UEMIS**2 - VEMIS**2)*
663 > (-CTPARS(NCT,XCT) + XG)) +
664 > UEMIS*(-CTPARS(NCT,ZCT) + ZG))**2 +
665 > (SQRT(1 - UEMIS**2 - VEMIS**2)*
666 > (-CTPARS(NCT,YCT) + YG) -
667 > VEMIS*(-CTPARS(NCT,ZCT) + ZG))**2)
668 IF ( DIST2 .LT. (CTPARS(NCT,CTDIAM)/2.) ) GOTO 102
669c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
670 101 CONTINUE
671
672 GOTO 7001
673
674C BUNCH FALLS ON A DETECTOR, CALCULATE ARRIVAL TIME (NSEC)
675
676 102 CARTIM = ((ETADSN*(THCKOB(NOBSLV)-THICK(DBLE(ZEMIS)))
677 * /WEMIS+PATHCR-VCOR/BETACR)/C(25)+TPART)*1.E9
678
679c CALL OUTPT2(INT(CURPAR(1)),I)
680 CALL OUTPT2(NCT,I)
681
682c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
683 call jctime(cartim)
684c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
685
686 GOTO 1000
687
688 7001 CONTINUE
689 1000 CONTINUE
690c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
691 RETURN
692 END
693
694C=======================================================================
695
696 SUBROUTINE GETBUS( IPARTI,ENERGY,THETA,CERSZE )
697
698C-----------------------------------------------------------------------
699C GET BU(NCH) S(IZE)
700C
701C CALCULATES OPTIMAL BUNCH SIZE FOR CERENKOV PHOTONS. CERENKOV PHOTONS
702C ARE GROUPED IN BUNCHES IN ORDER TO ACCELERATE COMPUTING TIME.
703C HOWEVER, WE SET A MAXIMAL VALUE FOR THE GROUPING OF CERENKOV PHOTONS
704C SO THAT WE GET AT LEAST 100 BUNCHES/M**2 AT A CERENKOV FLUX OF 3000
705C PHOTONS/M**2. THIS IS THE MINIMUM CERENKOV FLUX WHICH CAN BE
706C DISTINGUISHED FROM THE NIGHT SKY LIGHT BACKGROUND IN THE HEGRA
707C EXPERIMENT AT THE ISLAND LA PALMA. SO THE PARAMETRIZATION OF THE
708C CERENKOV BUNCH AS CALCULATED IN THIS SUBROUTINE IS VALID FOR
709C OBSERVATION LEVELS SIMILAR TO THAT OF THE HEGRA EXPERIMENT.
710C FOR A GIVEN PRIMARY PARTICLE, INCIDENT ENERGY AND ANGLE, AN
711C OPTIMAL BUNCH SIZE IS CALCULATED BY INTERPOLATION IN A TABLE,
712C WHERE WE HAVE CHOSEN AN ENERGY RANGE UP TO 1000 TEV, INCIDENT
713C ANGLES 0 AND 40 DEGREES, AND 4 TYPES OF PRIMARIS: GAMMAS,
714C PROTONS, NITROGEN, AND IRON.
715C THIS SUBROUTINE IS CALLED FROM MAIN
716C ARGUMENTS:
717C IPARTI = TYPE OF PRIMARY PARTICLE
718C ENERGY (R4) = PARTICLES ENERGY IN GEV
719C THETA (R4) = ANGLE IN RAD
720C CERSZE (R4) = SIZE OF CERENKOV BUNCH
721C
722C AUTHORS : S. MARTINEZ UNIVERSITY OF MADRID
723C F. ARQUEROS UNIVERSITY OF MADRID
724C-----------------------------------------------------------------------
725
726 IMPLICIT NONE
727*KEEP,CONST.
728 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
729 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
730*KEEP,RUNPAR.
731 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
732 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
733 * MONIOU,MDEBUG,NUCNUC,
734 * CETAPE,
735 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
736 * N1STTR,MDBASE,
737 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
738 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
739 * ,GHEISH,GHESIG
740 COMMON /RUNPAC/ DSN,HOST,USER
741 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
742 REAL STEPFC
743 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
744 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
745 * N1STTR,MDBASE
746 INTEGER CETAPE
747 CHARACTER*79 DSN
748 CHARACTER*20 HOST,USER
749
750 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
751 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
752 * ,GHEISH,GHESIG
753*KEND.
754
755 REAL ANGLE(2),ENGAM(3),ENHAD(3),ENNIT(2),
756 * SIFE(3,2),SIGAM(3,2),SINIT(2),SIPRO(3,2)
757 REAL CERS1F,CERS1P,ENERGY,CERSZE,S1,S2,THETA
758 INTEGER I,IANFE,IANP,IATNUM,IPARTI,I1,I2
759
760 DATA ANGLE / 0., 40. /
761 DATA ENGAM / 100., 200., 500. /
762 DATA ENHAD / 100., 200., 1000. /
763 DATA ENNIT / 200., 1000. /
764 DATA ( SIFE (I,1),I=1,3 ) / 30., 30., 140. /
765 DATA ( SIFE (I,2),I=1,3 ) / 30., 30., 110. /
766 DATA ( SIGAM(I,1),I=1,3 ) / 30., 45., 100. /
767 DATA ( SIGAM(I,2),I=1,3 ) / 30., 40., 100. /
768 DATA SINIT / 30., 150. /
769 DATA ( SIPRO(I,1),I=1,3 ) / 30., 30., 120. /
770 DATA ( SIPRO(I,2),I=1,3 ) / 30., 30., 160. /
771 DATA IANP / 1 /, IANFE / 26 /
772C-----------------------------------------------------------------------
773
774c IF ( DEBUG ) WRITE(MDEBUG,100) IPARTI,SNGL(ENERGY),SNGL(THETA)
775c 100 FORMAT(' GETBUS: INPUT PARTICLE = ',I5,1P,2E10.3)
776
777C DEFAULT VALUE
778 CERSZE = 100.
779
780 ENERGY = ENERGY / 1000.
781 IF ( ENERGY .LE. 100. ) THEN
782 CERSZE = 30.
783 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
784 RETURN
785 ENDIF
786
787 THETA = THETA / PI * 180.
788
789C-----------------------------------------------------------------------
790C PHOTON, ELECTRON OR POSITRON AS PRIMARY PARTICLE
791 IF ( IPARTI .LE. 3 ) THEN
792C FIND ENERGY BIN FOR INTERPOLATION
793 IF ( ENERGY .LE. ENGAM(2) ) THEN
794 I1 = 1
795 I2 = 2
796 ELSE
797 I1 = 2
798 I2 = 3
799 ENDIF
800 S1 = SIGAM(I1,1) + (ENERGY - ENGAM(I1))
801 * / (ENGAM(I2) - ENGAM(I1))
802 * * (SIGAM(I2,1) - SIGAM(I1,1))
803 S2 = SIGAM(I1,2) + (ENERGY - ENGAM(I1))
804 * / (ENGAM(I2) - ENGAM(I1))
805 * * (SIGAM(I2,2) - SIGAM(I1,2))
806 CERSZE = S1 + (THETA-ANGLE(1))/(ANGLE(2)-ANGLE(1)) * (S2-S1)
807 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
808 RETURN
809 ENDIF
810
811C-----------------------------------------------------------------------
812C NITROGEN AS PRIMARY PARTICLE AND VERTICAL INCIDENCE
813CJOK WHY SPECIAL TREATMENT FOR NITROGEN ????
814CJOK WHY ONLY VERTICAL INCIDENCE ????
815 IF ( IPARTI .EQ. 1407 .AND. ABS(THETA) .LT. 1.E-1 ) THEN
816 IF ( ENERGY .LT. 200. ) THEN
817 CERSZE = 30.
818 ELSE
819 CERSZE = SINIT(1) + (ENERGY-ENNIT(1))
820 * / (ENNIT(2)-ENNIT(1)) * (SINIT(2)-SINIT(1))
821 ENDIF
822 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
823 RETURN
824 ENDIF
825
826C-----------------------------------------------------------------------
827C GET THE ATOMIC NUMBER OF THE NUCLEUS
828C Z IS 1, IF PROTON
829 IF ( IPARTI .EQ. 14 ) THEN
830 IATNUM = 1
831C REST OF POSSIBLE NUCLEI
832 ELSEIF ( IPARTI .GT. 100 ) THEN
833 IATNUM = MOD(IPARTI,100)
834 IF ( IATNUM .GT. 26 ) THEN
835 WRITE(MONIOU,*) 'GETBUS: UNEXPECTED PARTICLE CODE'
836 RETURN
837 ENDIF
838 ELSE
839 WRITE(MONIOU,*) 'GETBUS: UNEXPECTED PARTICLE CODE'
840 RETURN
841 ENDIF
842
843C FIND ENERGY BIN FOR INTERPOLATION IN CASE OF HADRONIC PRIMARY
844 IF ( ENERGY .LE. ENHAD(2) ) THEN
845 I1 = 1
846 I2 = 2
847 ELSE
848 I1 = 2
849 I2 = 3
850 ENDIF
851
852C INTERPOLATION FOR HADRONS
853 S1 = SIPRO(I1,1) + (ENERGY-ENHAD(I1))
854 * / (ENHAD(I2)-ENHAD(I1)) * (SIPRO(I2,1)-SIPRO(I1,1))
855 S2 = SIPRO(I1,2) + (ENERGY-ENHAD(I1))
856 * / (ENHAD(I2)-ENHAD(I1)) * (SIPRO(I2,2)-SIPRO(I1,2))
857 CERS1P = S1 + (THETA-ANGLE(1)) / (ANGLE(2)-ANGLE(1)) * (S2-S1)
858
859 S1 = SIFE(I1,1) + (ENERGY-ENHAD(I1)) / (ENHAD(I2)-ENHAD(I1))
860 * * (SIFE(I2,1)-SIFE(I1,1))
861 S2 = SIFE(I1,2) + (ENERGY-ENHAD(I1)) / (ENHAD(I2)-ENHAD(I1))
862 * * (SIFE(I2,2)-SIFE(I1,2))
863 CERS1F = S1 + (THETA-ANGLE(1)) / (ANGLE(2)-ANGLE(1)) * (S2-S1)
864
865 CERSZE = CERS1P + (IATNUM-IANP) * (CERS1F-CERS1P) / (IANFE-IANP)
866
867 IF ( DEBUG ) WRITE(MDEBUG,101) CERSZE
868 101 FORMAT(' GETBUS: BUNCH SIZE = ',1P,1E10.3)
869
870 RETURN
871 END
872
873C=======================================================================
874
875 SUBROUTINE OUTND2
876
877C-----------------------------------------------------------------------
878C OUT(PUT AT E)ND (OF SHOWER)
879C
880C WRITE REST OF PARTICLES TO OUTPUT BUFFER
881C OUTND2 IS CALLED FROM MAIN
882C
883C AUTHORS : S. MARTINEZ, UNIVERSITY OF MADRID
884C F. ARQUEROS, UNIVERSITY OF MADRID
885C-----------------------------------------------------------------------
886
887 IMPLICIT NONE
888*KEEP,RUNPAR.
889 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
890 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
891 * MONIOU,MDEBUG,NUCNUC,
892 * CETAPE,
893 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
894 * N1STTR,MDBASE,
895 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
896 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
897 * ,GHEISH,GHESIG
898 COMMON /RUNPAC/ DSN,HOST,USER
899 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
900 REAL STEPFC
901 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
902 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
903 * N1STTR,MDBASE
904 INTEGER CETAPE
905 CHARACTER*79 DSN
906 CHARACTER*20 HOST,USER
907
908 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
909 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
910 * ,GHEISH,GHESIG
911*KEEP,CEREN2.
912 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
913 * DCERX,DCERY,ACERX,ACERY,
914 * XCMAX,YCMAX,EPSX,EPSY,
915 * DCERXI,DCERYI,FCERX,FCERY,
916 * XSCATT,YSCATT,CERXOS,CERYOS,
917 * NCERX,NCERY,ICERML
918 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
919 * DCERX,DCERY,ACERX,ACERY,
920 * XCMAX,YCMAX,EPSX,EPSY,
921 * DCERXI,DCERYI,FCERX,FCERY,
922 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
923 INTEGER NCERX,NCERY,ICERML
924*KEEP,CEREN3.
925 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
926 INTEGER MAXBF2
927 PARAMETER (MAXBF2 = 39 * 7)
928 DOUBLE PRECISION CERCNT
929 REAL DATAB2(MAXBF2)
930 INTEGER LHCER
931*KEND.
932
933 INTEGER I
934C-----------------------------------------------------------------------
935
936 IF ( LHCER .GT. 0 ) THEN
937 CALL TOBUFC( DATAB2,1 )
938 DO 2 I = 1,MAXBF2
939 DATAB2(I) = 0.
940 2 CONTINUE
941 ELSE
942 CALL TOBUFC( DATAB2,2 )
943 ENDIF
944
945 WRITE(MONIOU,*) 'CERCNT = ',SNGL( CERCNT )
946 CERCNT = 0.D0
947 LHCER = 0
948
949 RETURN
950 END
951
952C=======================================================================
953
954 SUBROUTINE OUTPT2(J,IMOV)
955
956C-----------------------------------------------------------------------
957C (WRITE CERENKOV RADIATION) OUTP(U)T
958C
959C OUTPUT ROUTINE FOR CERENKOV PHOTONS
960C THIS SUBROUTINE IS CALLED FROM CERENE AND CERENH
961C
962C AUTHORS : S. MARTINEZ, UNIVERSITY OF MADRID
963C F. ARQUEROS, UNIVERSITY OF MADRID
964C-----------------------------------------------------------------------
965
966 IMPLICIT NONE
967*KEEP,BUFFS.
968 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
969 INTEGER MAXBUF,MAXLEN
970 PARAMETER (MAXBUF=39*7)
971 PARAMETER (MAXLEN=12)
972 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
973 * RUNE(MAXBUF),DATAB(MAXBUF)
974 INTEGER LH
975 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
976 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
977 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
978*KEEP,RUNPAR.
979 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
980 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
981 * MONIOU,MDEBUG,NUCNUC,
982 * CETAPE,
983 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
984 * N1STTR,MDBASE,
985 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
986 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
987 * ,GHEISH,GHESIG
988 COMMON /RUNPAC/ DSN,HOST,USER
989 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
990 REAL STEPFC
991 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
992 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
993 * N1STTR,MDBASE
994 INTEGER CETAPE
995 CHARACTER*79 DSN
996 CHARACTER*20 HOST,USER
997
998 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
999 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1000 * ,GHEISH,GHESIG
1001*KEEP,CEREN1.
1002 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
1003 * CERSIZ,LCERFI
1004 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
1005 REAL CERSIZ
1006 LOGICAL LCERFI
1007*KEEP,CEREN2.
1008 COMMON /CEREN2/ PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
1009 * DCERX,DCERY,ACERX,ACERY,
1010 * XCMAX,YCMAX,EPSX,EPSY,
1011 * DCERXI,DCERYI,FCERX,FCERY,
1012 * XSCATT,YSCATT,CERXOS,CERYOS,
1013 * NCERX,NCERY,ICERML
1014 REAL PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS,
1015 * DCERX,DCERY,ACERX,ACERY,
1016 * XCMAX,YCMAX,EPSX,EPSY,
1017 * DCERXI,DCERYI,FCERX,FCERY,
1018 * XSCATT,YSCATT,CERXOS(20),CERYOS(20)
1019 INTEGER NCERX,NCERY,ICERML
1020*KEEP,CEREN3.
1021 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
1022 INTEGER MAXBF2
1023 PARAMETER (MAXBF2 = 39 * 7)
1024 DOUBLE PRECISION CERCNT
1025 REAL DATAB2(MAXBF2)
1026 INTEGER LHCER
1027c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1028 COMMON /GRAAL1/ WAVELENGTH ! (NM)
1029 REAL WAVELENGTH
1030c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1031*KEND.
1032
1033 INTEGER I,J,IMOV
1034C-----------------------------------------------------------------------
1035
1036 IF(DEBUG)WRITE(MDEBUG,3)PHOTCM,XCER,YCER,UEMIS,VEMIS,CARTIM,ZEMIS
1037 3 FORMAT(' OUTPT2: ',1P,8E10.3)
1038C WRITE A BLOCK OF 39 PARTICLES TO THE CERENKOV OUTPUT BUFFER AND
1039C CLEAR FIELD
1040 IF ( LCERFI ) THEN
1041c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1042c DATAB2(LHCER+1) = PHOTCM
1043cc DATAB2(LHCER+1) = WAVELENGTH + J*1000.
1044 DATAB2(LHCER+1) = J*100000. + IMOV*1000. + WAVELENGTH
1045c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1046 DATAB2(LHCER+2) = XCER
1047 DATAB2(LHCER+3) = YCER
1048 DATAB2(LHCER+4) = UEMIS
1049 DATAB2(LHCER+5) = VEMIS
1050 DATAB2(LHCER+6) = CARTIM
1051 DATAB2(LHCER+7) = ZEMIS
1052c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1053c CERCNT = CERCNT + DBLE( PHOTCM )
1054 CERCNT = CERCNT + 1
1055c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1056 LHCER = LHCER + 7
1057 IF ( LHCER .GE. MAXBF2 ) THEN
1058 CALL TOBUFC( DATAB2,0 )
1059 DO 1 I = 1,MAXBF2
1060 DATAB2(I) = 0.
1061 1 CONTINUE
1062 LHCER = 0
1063 ENDIF
1064 ELSE
1065C WRITE A BLOCK OF 39 PARTICLES TO THE PARTICLE OUTPUT BUFFER AND
1066C CLEAR FIELD
1067 DATAB(LH+1) = 99.E5 + NINT(PHOTCM)*10. + 1.
1068 DATAB(LH+2) = XCER
1069 DATAB(LH+3) = YCER
1070 DATAB(LH+4) = UEMIS
1071 DATAB(LH+5) = VEMIS
1072 DATAB(LH+6) = CARTIM
1073 DATAB(LH+7) = ZEMIS
1074 LH = LH + 7
1075 NOPART = NOPART + 1
1076c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1077c CERCNT = CERCNT + DBLE( PHOTCM )
1078 CERCNT = CERCNT + 1
1079c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1080 IF ( LH .GE. MAXBUF ) THEN
1081 CALL TOBUF( DATAB,0 )
1082 DO 2 I = 1,MAXBUF
1083 DATAB(I) = 0.
1084 2 CONTINUE
1085 LH = 0
1086 ENDIF
1087 ENDIF
1088
1089 RETURN
1090 END
1091
1092C=======================================================================
1093
1094 SUBROUTINE TOBUFC( A,IFL )
1095
1096C-----------------------------------------------------------------------
1097C (WRITE) TO BUF(FER) C(ERENKOV DATA)
1098C
1099C COPY TO BUFFER CERENKOV DATA
1100C THIS SUBROUTINE IS CALLED FROM MAIN, INPRM, ELECTR, PHOTON, OUTND2,
1101C AND OUTPT2
1102C ARGUMENTS:
1103C A = ARRAY TO BE WRITTEN TO TAPE
1104C IFL = STARTING OF FINAL OUTPUT
1105C = 0 NORMAL BLOCK
1106C = 1 NORMAL BLOCK WITH END OF OUTPUT
1107C = 2 ONLY END OF OUTPUT
1108C-----------------------------------------------------------------------
1109
1110 IMPLICIT NONE
1111*KEEP,BUFFS.
1112 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
1113 INTEGER MAXBUF,MAXLEN
1114 PARAMETER (MAXBUF=39*7)
1115 PARAMETER (MAXLEN=12)
1116 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
1117 * RUNE(MAXBUF),DATAB(MAXBUF)
1118 INTEGER LH
1119 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
1120 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
1121 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
1122*KEEP,RECORD.
1123 COMMON /RECORD/ IRECOR
1124 INTEGER IRECOR
1125*KEEP,RUNPAR.
1126 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
1127 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
1128 * MONIOU,MDEBUG,NUCNUC,
1129 * CETAPE,
1130 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1131 * N1STTR,MDBASE,
1132 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1133 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1134 * ,GHEISH,GHESIG
1135 COMMON /RUNPAC/ DSN,HOST,USER
1136 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
1137 REAL STEPFC
1138 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
1139 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
1140 * N1STTR,MDBASE
1141 INTEGER CETAPE
1142 CHARACTER*79 DSN
1143 CHARACTER*20 HOST,USER
1144
1145 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
1146 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
1147 * ,GHEISH,GHESIG
1148*KEEP,CEREN3.
1149 COMMON /CEREN3/ CERCNT,DATAB2,LHCER
1150 INTEGER MAXBF2
1151 PARAMETER (MAXBF2 = 39 * 7)
1152 DOUBLE PRECISION CERCNT
1153 REAL DATAB2(MAXBF2)
1154 INTEGER LHCER
1155*KEEP,CEREN4.
1156 COMMON /CEREN4/ NRECER
1157 INTEGER NRECER
1158*KEND.
1159
1160 INTEGER NSUBBL
1161 PARAMETER (NSUBBL=21)
1162 REAL A(*)
1163C NSUBBL IS NUMBER OF SUBBLOCKS IN ONE OUTPUT RECORD
1164C (OUTPUT RECORD LENGTH = NSUBBL * 39 * 7 * 4 BYTES <= 22932 )
1165C IBLK2 IS COUNTER FOR SUBBLOCKS OF CERENKOV OUTPUT
1166C OUTPUT BUFFER FOR CERENKOV OUTPUT
1167 REAL OUTBF2(MAXBF2,NSUBBL)
1168 SAVE OUTBF2
1169 INTEGER I,IBLK2,IFL,K
1170 DATA IBLK2 / 0 /
1171C-----------------------------------------------------------------------
1172
1173 IF ( IFL .LE. 1 ) THEN
1174 IBLK2 = IBLK2 + 1
1175 DO 3 I = 1,MAXBF2
1176 OUTBF2(I,IBLK2) = A(I)
1177 3 CONTINUE
1178 ENDIF
1179
1180C WRITE TO TAPE IF BLOCK IS FULL OR IF IFL IS 1
1181 IF ( IFL .GE. 1 .OR. IBLK2 .EQ. NSUBBL ) THEN
1182 NRECER = NRECER + 1
1183
1184c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1185c WRITE(CETAPE) ((OUTBF2(I,K),I=1,MAXBF2),K=1,NSUBBL)
1186 call jccersave(outbf2)
1187c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1188 IBLK2 = 0
1189 DO 4 K = 1,NSUBBL
1190 DO 4 I = 1,MAXBF2
1191 OUTBF2(I,K) = 0.0
1192 4 CONTINUE
1193 ENDIF
1194
1195 RETURN
1196 END
Note: See TracBrowser for help on using the repository browser.