c **************************************************************** PROGRAM MAIN c **************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c ADDRESS: Departamento de Fisica Atomica, Molecular y Nuclear. c Facultad de Fisicas. Universidad Complutense. c Madrid. E-28040. c e-MAIL: martinez@eucmdx.gae.ucm.es c LAST UPDATE: 13 / Dec / 95 c c Program to simulate the AIROBICC detector. c Input files: CER00000# (from CORSIKA); airobicc.in c Output files: - airqt000# with TDC and qADC signals in each C_hut. c - airsta00# with info about average properties of c the C-light disk on the ground (before detector). c - airobicc.hbook with histograms for several showers. c - airobicc.out with log information c **************************************************************** c ==================================== c SHOWER INFORMATION c ==================================== c C_header contains global info of the shower, and is partially filled c in read_Cfile when reading the 1st record of the file CER# COMMON / C__header / C_header(20), corsika_version equivalence (C_header(1),PartType), (C_header(2),Energy) equivalence (C_header(3),Theta), (C_header(4),Phi) equivalence (C_header(5),height_1stint) equivalence (C_header(6),x_core), (C_header(7),y_core) equivalence (C_header(8),Huts_X), (C_header(9),ArrayGrid) equivalence (C_header(10),HutSide) equivalence (C_header(11),Angle_X_North) equivalence (C_header(12),cersiz) equivalence (C_header(13),vernum) c C_header(1): Incident particle code c (2): Incident energy (GeV) c (3): incident theta (rad) c (4): incident phi (rad) c (5): height of 1st interaction (cm) (not always available) c (6): x coordinate of core position (m) c (7): y coordinate of core position (m) c (8): Number of huts in X (Y) direction c (9): Grid spacing of huts (m) c (10): Side Length of hut (m) c (11): Angle between X-direction of array and North (rad) c (12): CERSIZ = photons/bunch (not always available) c (13): VERNUM = corsika version (not always available) common /EASfront/ X_delay1 , Y_delay1, X_delay2, Y_delay2 c C_hut_count(i,j) = No C-photons in square (i,j) c C_hut_arrTime(i,j) = of C-photons in square (i,j) c _At -> with the effect of atmospheric scattering DIMENSION C_hut_count(50,50) , C_hut_arrTime(50,50) DIMENSION C_hut_count_At(50,50) , C_hut_arrTime_At(50,50) c ==================================== c Array Characteristics c ==================================== common /HEGRA_level/ H_LAPALMA c ArraySide: Distance between square (1,1) and (1,N), where N=Number of c squares in extended array COMMON /C_Array/ C_ArrayHuts , C_ArrayGrid , C_ArraySide, * imiddle , jmiddle, C_HutSide INTEGER C_ArrayHuts c ============================ c Hut Geometry c ============================ c r_cone: radii of conical sections of the mirror COMMON /Chutgeom/ r_cone(0:3),zu(0:3),ta(3),ca(3),nz(3),zd(3) REAL nz REAL LightSpeed c ============================================ c Cerenkov Transmission (up to PM) c =========================================== c curphot: current C-bunch being simulated DIMENSION curphot(7) PARAMETER ( nlambdas = 4 ) COMMON /cerenkov/ C_spectrum(nlambdas) , lambda(nlambdas) REAL lambda COMMON /cerenkov_2/ extended_spectrum_factor COMMON /TRAtmos/ TR_Atmos(nlambdas) , iAtmos_absor common /TRCounter/ TR_Plexiglas(nlambdas), * Reflex_mirror(nlambdas), * TR_Filter(nlambdas), * QE_Photocath(nlambdas) DIMENSION TR_total(nlambdas) c ================================================================ c ============ ELECTRONICS =========== c ================================================================ c NSL: Night Sky Light COMMON /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM c Transit Time Curve of PM COMMON /TrTimePM/ nTrTime , TrTime_PM , Dist_TrTime_PM DIMENSION TrTime_PM(100) , Dist_TrTime_PM(100) COMMON /histog/ nC_bins , C_TimeWindow , C_time_shift COMMON /histog2/ nbins_PM , TimeWindow_PM , time_shift_PM COMMON /histog3/ nbins_Amp1_2 , TimeWindow_Amp1_2 COMMON /histog4/ Window_shift PARAMETER ( nbins = 2000 ) DIMENSION PMcurrent_C(nbins) COMMON / anodePM / PMcurrent(nbins), PMcurrent_NSL(nbins) common /Amp1_2_current/ currentAmp1_2( nbins ) COMMON / CFD / CFD_delay1 , CFD_delay2 , U_Threshold, Resist_CFD PARAMETER ( nhuts = 2000 ) c 1 -> after 1st amplifier c 1_2 -> after 2nd amplifier COMMON /ADC_TDC / q_ADC1( nhuts ), q_ADC1_2( nhuts ) * , TDC1_2( nhuts ) , Gate_ADC1_2 DIMENSION TDC_qADC( nhuts , 3 ) c ---------- Conversion factor between charge (pC) c ---------- and C-photons at the hut entrance dimension Conversion_qPhot_HG( 50, 50 ) dimension Conversion_qPhot_LG( 50, 50 ) c c ==================================== c I/O variables c ==================================== COMMON /datac1/ InpDirectory , OutDirectory , dsn0 CHARACTER*60 InpDirectory , OutDirectory CHARACTER*9 dsn0 , dsn COMMON /datac2/ i_C_file1 , n_C_files COMMON /LUNs/ Lun_Inp , Lun_Out , Lun_CORS_Cfile CHARACTER*70 C_corsika_file , Hbook_file , C_stat_file CHARACTER*70 C_QT_file INTEGER QT_RECL COMMON /Cfile/ cerbuf(273,30), irecor, n_C_blks, isb, ieof, lh c --------------- CHARACTER*40 HTitlePM CHARACTER*40 HTitleAmp1_2 c --------------- double precision C_phot_total c ===================================================================== c ASSIGN VALUES TO VARIABLES c DATA Lun_Inp / 5 / , Lun_Out / 6 / , Lun_CORS_Cfile / 3 / , * Lun_HBOOK / 4 / , Lun_C_STAT / 7 / , Lun_C_QT / 8 / c DATA nData_C_hut / 3 / DATA nBytes_C_QT_Data / 4 / c DATA LightSpeed / 30. / ! cm/ns DATA pi / 3.141593 / c ///////////////////////////////////////////////////////////////////// c ///////////////////////////////////////////////////////////////////// c c Initialize variables , arrays and histograms, and read input file c with steering cards CALL Init_AIROBICC raddeg = 180 / pi degrad = pi / 180 c ---------------------------------------- c Open output (log) unit OPEN ( unit=Lun_Out , file='airobicc.out',form='formatted', * access='sequential',status='unknown') c c ----------------------------------------- c Open unit for output of histograms ind_I = index( InpDirectory , ' ' ) ind_O = index( OutDirectory , ' ' ) DO i=1,70 Hbook_file(i:i) = ' ' ENDDO Hbook_file = OutDirectory WRITE(Hbook_file(ind_O:ind_O+13) , '(a14)') 'airobicc.hbook' icode = system('rm -f '//Hbook_file) OPEN ( unit=Lun_HBOOK ,file=Hbook_file,form='unformatted', * access='direct',status='new',recl=1024) c CALL hrfile( Lun_HBOOK ,'hbook','N') c c ========================* SHOWER LOOP ========================* c DO 2000 iEvent = i_C_file1 , i_C_file1 + n_C_files - 1 c c open more I/O units c DO i=1,70 C_corsika_file(i:i) = ' ' C_stat_file(i:i) = ' ' C_QT_file(i:i) = ' ' ENDDO c ==================================== c OPEN CERENKOV FILE FROM CORSIKA dsn = dsn0 WRITE(dsn(4:9),'(i6)') iEvent DO i=4,9 IF ( dsn(i:i) .eq. ' ' ) dsn(i:i) = '0' ENDDO c C_corsika_file = InpDirectory WRITE(C_corsika_file(ind_I:ind_I+8) , '(a9)') dsn c lrec = record length (words) in CER-file lrec = 273 * n_C_blks OPEN ( unit=Lun_CORS_Cfile , file=C_corsika_file , form= * 'unformatted', access='direct',status='old',recl=lrec * , err=2000, readonly) c =============================================================== c OPEN OUTPUT FILE FOR STATISTICS OF C-PHOTONS BEFORE C-DETECTOR dsn = 'airsta000' WRITE(dsn(7:9),'(i3)') iEvent DO i=7,9 IF ( dsn(i:i) .eq. ' ' ) dsn(i:i) = '0' ENDDO C_stat_file = OutDirectory WRITE(C_stat_file(ind_O:ind_O+8) , '(a9)') dsn c icode = system('rm -f '//C_stat_file) OPEN ( unit=Lun_C_STAT ,file=C_stat_file,form='formatted', * access='sequential',status='new') c ---------------------------------------- c Initialize some variables for information of single events c CALL vzero( C_hut_count , 2500 ) CALL vzero( C_hut_arrTime , 2500 ) c CALL vzero( C_hut_count_At , 2500 ) CALL vzero( C_hut_arrTime_At , 2500 ) c C_phot_total = 0.d0 c n_Bad_uv = 0 c c Initialize variables for reading C-files c irecor = 0 ieof = 0 isb = n_C_blks c c LOOP FOR READING IN DATA OF C-PHOTONS c 20 isb = isb + 1 IF ( isb . gt . n_C_blks ) THEN CALL Read_Cfile IF ( irecor .eq. 1 ) THEN call Init_Histos c --------------------------------------- WRITE( Lun_Out , '(/)' ) if (iEvent .eq. i_C_file1) then write( Lun_Out ,*) 'Atmospheric absorption ->', * iAtmos_absor endif WRITE( Lun_Out ,*) ' EVENT NUMBER = ' , iEvent WRITE( Lun_Out , '(a,f5.0)') ' PARTICLE = ', PartType WRITE( Lun_Out , '(a,f6.1,a)') ' ENERGY = ', * Energy/1e3, ' TeV' WRITE( Lun_Out , '(a,f5.1,a)') ' THETA = ', & Theta*raddeg, ' deg' WRITE( Lun_Out , '(a,f5.1,a)') ' PHI = ', & Phi*raddeg, ' deg' WRITE( Lun_Out , '(/)' ) c c --------------------------------------- WRITE( Lun_C_STAT , * ) ' SHOWER #', iEvent WRITE(Lun_C_STAT, '(a,i5)' ) ' PARTICLE = ', int(PartType) WRITE( Lun_C_STAT , 7002 ) ' ENERGY = ', * Energy/1e3, ' TeV ; ' , ' THETA = ', Theta*raddeg, * ' deg',' PHI = ', Phi*raddeg, ' deg' 7002 FORMAT( a , f6.1 , a , a , f5.1 , a , a , f5.1 , a) write( Lun_C_STAT ,*) 'Atmospheric absorption ->', * iAtmos_absor c --------------------------------------- c c RedFac: reduction factor to get correct intensity at the huts c RedFac = pi * (r_cone(3)/100.) ** 2 / C_HutSide ** 2 GOTO 30 ENDIF !END OF if ( irecor .eq. 1 ) ENDIF !END OF if ( isb . gt . n_C_blks ) c c End of C_corsika_file ? c IF ( ieof .eq. 1 ) THEN irecor = irecor - 1 WRITE( Lun_Out ,*) 'Early end of Cfile at record ', irecor GOTO 1000 ! Process signals in each hut ENDIF if ( cerbuf(1, isb) .eq. 'CELO' ) then write( Lun_Out ,*) 'Mark CELO found' GOTO 1000 ENDIF c lh+1: C-bunch counter in current isb-block (it does not mean Luft Hansa!) c lh range: 0-38 lh = 0 30 CONTINUE DO 40 i = 1,7 curphot(i) = cerbuf( 7 * lh + i , isb ) 40 CONTINUE c c End of C_corsika_file ? c IF ( (curphot(1) .eq. 0.) .and. (curphot(7) .eq. 0.) ) THEN WRITE( Lun_Out ,*) ' Normal end of C-file' GOTO 1000 c Maybe the C-bunch has null content (?) ELSEIF ( curphot(1) .eq. 0. ) THEN GOTO 50 ENDIF c ============================================================= c ======= GET PROPERTIES OF C-BUNCH ========== c ============================================================= c C_photons_0: number of C-photons in bunch C_photons_0 = curphot(1) c c Correct for change in C-spectral range: CORSIKA::300-450 nm -> C AIROBICC::300-500 nm c C_photons_0 = C_photons_0 * extended_spectrum_factor c c C_x, C_y (meters): coordinates respect to shower axis if (corsika_version .eq. 3.02) then C_x = curphot(2) ! m C_y = curphot(3) ! m else C_x = curphot(2) / 100. ! m C_y = curphot(3) / 100. ! m endif c c C_u_emis,C_v_emis: direction cosines in x,y axis C_u_emis = curphot(4) C_v_emis = curphot(5) c C_arrTime: arrival time (ns) ,respect to axis arrival time. if (corsika_version .gt. 4.06) then time_1stint = - (height_1stint - H_LAPALMA) / & cos( Theta ) / LightSpeed else time_1stint = 0. endif C_ArrTime = curphot(6) - time_1stint c C_z_emis = height (a.s.l, cm) of C-bunch emission C_z_emis = curphot(7) C_phot_total = C_phot_total + C_photons_0 C_w_emis = 1. - C_u_emis ** 2 - C_v_emis ** 2 IF ( C_w_emis .LT. 0. ) THEN n_Bad_uv = n_Bad_uv + 1 GOTO 50 ENDIF C_w_emis = sqrt( C_w_emis ) c x,y: Coordinates respect to lower left corner of extended array x = C_x + C_ArraySide / 2. y = C_y + C_ArraySide / 2. c --------------------------------------------- c check if bunch inside extended C_array c IF ( (x .lt. -C_HutSide/2) .or. (x .gt. (C_ArraySide + * C_HutSide/2) ) ) then write( Lun_Out,*) 'C-bunch out of array (x)' GOTO 50 endif IF ( (y .lt. -C_HutSide/2) .or. (y .gt. (C_ArraySide + * C_HutSide/2) ) ) then write( Lun_Out,*) 'C-bunch out of array (y) ' GOTO 50 endif c ------------------------------------------------- c calculate the index (ihut,jhut) of the C-hut: >0 c jhut -> along X-axis c ihut -> along Y-axis ihut = ( y + C_HutSide/2 ) / C_ArrayGrid + 1 jhut = ( x + C_HutSide/2 ) / C_ArrayGrid + 1 c Get ID of hut for filling histograms id_hut = C_ArrayHuts * (ihut-1) + jhut c ------------------------------------------------- c calculate relative coordinates to the center of the square c x0 = x - C_ArrayGrid * (jhut-1) y0 = y - C_ArrayGrid * (ihut-1) if ( (abs(x0) .gt. C_HutSide/2) .or. & (abs(y0) .gt. C_HutSide/2) ) then write( 6,*) 'C-bunch out of the square (ihut,jhut)' goto 50 endif c ************************************************************* c ****** C-BUNCH HAS HITTED A HUT OF EXTENDED ARRAY ******** c ************************************************************* c ==================================================================== c get Shower Properties in squares of extended matrix c ==================================================================== C_hut_count(ihut, jhut) = C_hut_count(ihut, jhut) + C_photons_0 C_hut_arrTime( ihut , jhut ) = * C_hut_arrTime( ihut , jhut ) + C_photons_0 * C_ArrTime c----------------------------------------------------------------------* IF ( iAtmos_absor .eq. 1 ) THEN CALL Atmos_absorption( C_z_emis , C_w_emis ) C_phot_Atmos = 0. DO i = 1 , nlambdas C_phot_Atmos = C_phot_Atmos + C_spectrum(i) * * TR_Atmos(i) * C_photons_0 ENDDO C_hut_count_At( ihut , jhut ) = * C_hut_count_At( ihut, jhut) + C_phot_Atmos C_hut_arrTime_At( ihut , jhut ) = * C_hut_arrTime_At(ihut, jhut) + C_phot_Atmos * C_ArrTime ENDIF c ========================================================== c ======= SIMULATE C-BUNCH IN AIROBICC COUNTER ========= c ========================================================== c --------------------------------------------------- c get a random position at the entrance of the C-hut c 10 rndm01 = rndm( dummy ) c c the number of C-bunches between r , r+dr is c proportional to r c rndm02 = rndm( dummy ) IF ( rndm02 . gt . rndm01 ) GOTO 10 rndm03 = rndm( dummy ) r_rndm = r_cone(3) * rndm01 phi_rndm = 2 * pi * rndm03 x_rndm = r_rndm * cos( phi_rndm ) y_rndm = r_rndm * sin( phi_rndm ) c ----------------------------------------------- c Correct arrival time after shifting impact point, due to c shower front delay c delta_x = x0 * 100 - x_rndm ! cm delta_y = y0 * 100 - y_rndm ! cm time_corr = ( delta_x * C_u_emis + delta_y * C_v_emis ) * / LightSpeed C_ArrTime = C_ArrTime - time_corr c c ----------------------------------------------------- c tracking of the C-photons down to the photocathode c tracklength: length of trayectory inside counter c nreflex: number of reflexions in mirror CALL Ctracking( x_rndm , y_rndm , C_u_emis , C_v_emis , * tracklength , nreflex ) c c Check if C-bunch reflected out of counter IF ( tracklength . le . 0.0 ) GOTO 50 c ------------------------------------------------------------------ c Get the number of photoelectrons at the photocathode : Photoelect_PM c c Correct for ratio between areas of square and counter C_photons_0 = C_photons_0 * RedFac Photoelect_PM = 0. DO i=1,nlambdas TR_total(i) = TR_Atmos(i) * TR_Plexiglas(i) * * Reflex_mirror(i) ** nreflex * * TR_Filter(i) * QE_Photocath(i) Photoelect_PM = Photoelect_PM + C_photons_0 * * C_spectrum(i) * TR_total(i) ENDDO c --------------------------------------------------------- c Fill histograms: get the current at the anode of the PM due to C-bunch c c delay_hut: Time of flight inside counter delay_hut = tracklength / LightSpeed DO i=1,nTrTime time = C_ArrTime + delay_hut + TrTime_PM(i) current_PM = Photoelect_PM * Gain_PM * * current_factor * Dist_TrTime_PM(i) CALL hfill(id_hut,time,0.,current_PM) ENDDO c 50 lh = lh + 1 IF ( lh .eq. 39 ) GOTO 20 GOTO 30 c =================================================== c END OF LOOP OVER CERENKOV PHOTONS c =================================================== c 1000 CONTINUE WRITE( Lun_Out ,'(a,a,f10.0)') 'Total number of Cerenkov ', * 'photons in extended C-Array = ', sngl( C_phot_total ) WRITE( Lun_Out, '(a,a,i)') 'Number of Cerenkov bunches with ', * 'ilegal values of direction cosines u, v = ', n_Bad_uv WRITE( Lun_Out ,*) 'irecor = ', irecor WRITE( Lun_Out ,*) 'isb = ', isb WRITE( Lun_Out ,*) 'lh = ', lh CLOSE( Lun_CORS_Cfile ) c C =====================================================================* c ======* Get the ADC and TDC SIGNALS for each HUT ============ C =====================================================================* c c WARNING: Now the central hut (imiddle, jmiddle) is at (0,0) HTitlePM = 'AnodeCurrent x = 111111 m, y = 111111 m' HTitleAmp1_2 = 'Amp-Current x = 111111 m , y = 111111 m' DO 60 j=1,C_ArrayHuts jhut = j X_tmin = X_delay1 * ( jhut - jmiddle ) + X_delay2 xhut = C_ArrayGrid * ( jhut - jmiddle ) WRITE(HTitlePM(18:23),'(f6.1)') xhut WRITE(HTitleAmp1_2(17:22),'(f6.1)') xhut DO 70 i=1,C_ArrayHuts ihut = i Y_tmin = Y_delay1 * ( ihut - imiddle ) + Y_delay2 tmin_PM = X_tmin + Y_tmin - time_shift_PM tmax_PM = tmin_PM + TimeWindow_PM tmin_Amp1_2 = tmin_PM + Window_shift tmax_Amp1_2 = tmin_Amp1_2 + TimeWindow_Amp1_2 yhut = C_ArrayGrid * ( ihut - imiddle ) WRITE(HTitlePM(32:37),'(f6.1)') yhut WRITE(HTitleAmp1_2(32:37),'(f6.1)') yhut c ----------------------------------------------------------- c Book 2 histograms for each counter: 1) Current at PM anode c (id = 2000 + id_hut) c 2) Current after 2nd Amplifier c (id = 4000 + id_hut) id_hut = C_ArrayHuts * ( ihut - 1 ) + jhut CALL hbook1(2000+id_hut, HTitlePM, nbins_PM, tmin_PM, & tmax_PM, 0.) CALL hbook1(4000+id_hut, HTitleAmp1_2 , * nbins_Amp1_2 , tmin_Amp1_2 , tmax_Amp1_2 , 0.) c c ------------------------------------- c Fill vector NSLNoise (anode current due to NSL background) CALL Get_NSLNoise do ibin=1, nbins_PM PMcurrent(ibin) = PMcurrent_NSL(ibin) enddo c -------------------------------------- c Get anode current due to C-light CALL hunpak( id_hut , PMcurrent_C , ' ' , 0 ) c -------------------------------------- c Add Anode Currents from NSL & C-bunches ibin1 = 2 * (time_shift_PM - C_time_shift) + 1 ibin2 = ibin1 + nC_bins - 1 DO ibin = ibin1, ibin2 PMcurrent(ibin) = PMcurrent(ibin) + * PMcurrent_C(ibin - ibin1 + 1) ENDDO c ---------------------------------------------------- c I) Get integrated current after 1st amplifier (LOW-GAIN BRANCH) c CALL Get_qADC1(id_hut) c ------------------------------------------- c II) Get current after 2nd amplifier c (HIGH-GAIN BRANCH) C CALL Get_current_Amp1_2 c -------------------------------------------------- c III) Get the arrival time of the C-light front and integrate the c current of the pulse using a timewindow of fixed duration c CALL Get_TDC_qADC1_2_DATA(id_hut) c ----------------------------- c Output histograms for 1 showers and specific counters (those spaced c 3 * C_ArrayGrid from the central hut) if (iEvent .eq. i_C_file1) then krest_x = mod(iabs(jhut-jmiddle), 3) krest_y = mod(iabs(ihut-imiddle), 3) if ( (krest_x .eq. 0) .and. (krest_y .eq. 0) ) then CALL hpak( id_hut+2000 , PMcurrent ) CALL hpak( id_hut+4000 , currentAmp1_2 ) CALL hcdir( '//hbook', ' ' ) call hrout( id_hut+2000, icycle, ' ' ) call hrout( id_hut+4000, icycle, ' ' ) CALL hcdir( '//PAWC', ' ' ) endif endif CALL hdelet(2000+id_hut) CALL hdelet(4000+id_hut) 70 CONTINUE 60 CONTINUE c ============================================================== c OPEN OUTPUT FILE FOR THE FINAL DATA ( TIME & CHARGE ) c ============================================================== dsn = 'airqt0000' WRITE(dsn(7:9),'(i3)') iEvent DO i=7,9 IF ( dsn(i:i) .eq. ' ' ) dsn(i:i) = '0' ENDDO C_QT_file = OutDirectory WRITE(C_QT_file(ind_O:ind_O+8) , '(a9)') dsn c c QT_RECL = record length (words) of file airqt# QT_RECL = C_ArrayHuts * nData_C_hut * nBytes_C_QT_Data / 4 icode = system('rm -f '//C_QT_file) OPEN ( unit=Lun_C_QT ,file=C_QT_file,form='unformatted', * access='direct',status='new',recl = QT_RECL) c c Read out C_header to 1st record of output data file c C_header(14) = sngl(C_phot_total) ! in extended array C_header(15) = U_Threshold ! mV (<0) WRITE( Lun_C_QT , REC= 1 ) ( C_header(i) , i=1, 20 ) c c write to disk information about TDC, q_ADC1_2 , and q_ADC1 c DO id_hut = 1 , C_ArrayHuts ** 2 TDC_qADC( id_hut , 1 ) = TDC1_2( id_hut ) TDC_qADC( id_hut , 2 ) = q_ADC1_2( id_hut ) TDC_qADC( id_hut , 3 ) = q_ADC1( id_hut ) ENDDO c DO ihut= 1 , C_ArrayHuts id_hut1 = (ihut-1) * C_ArrayHuts + 1 id_hut2 = ihut * C_ArrayHuts WRITE( Lun_C_QT , REC= ihut + 1 ) * ( ( TDC_qADC( id_hut , j ) , j=1,3 ) * , id_hut = id_hut1 , id_hut2 ) ENDDO CLOSE( Lun_C_QT ) c =================================================================== c Get conversion factor from charge (pC) to C-density (phots/m2) for c every counter of the extended AIROBICC array, both for the low-gain c and high-gain branches. c =================================================================== DO 1010 i = 1 , C_ArrayHuts DO 1010 j = 1 , C_ArrayHuts Conversion_qPhot_HG(i,j) = 0. Conversion_qPhot_LG(i,j) = 0. IF ( C_hut_count(i,j) . gt . 0. ) then C_hut_arrTime(i,j) = C_hut_arrTime(i,j) / * C_hut_count(i,j) C_hut_count(i,j) = C_hut_count(i,j) / C_HutSide**2 IF ( iAtmos_absor .eq. 1 ) C_hut_arrTime_At(i,j) = * C_hut_arrTime_At(i,j) / C_hut_count_At(i,j) IF ( iAtmos_absor .eq. 1 ) C_hut_count_At(i,j) = & C_hut_count_At(i,j) / C_HutSide**2 else goto 1010 ENDIF id_hut = (i-1) * C_ArrayHuts + j if ( TDC_qADC(id_hut,2) * iAtmos_absor .gt. 0.) then Conversion_qPhot_HG(i,j) = C_hut_count_At(i,j) / * TDC_qADC(id_hut,2) Conversion_qPhot_LG(i,j) = C_hut_count_At(i,j) / * TDC_qADC(id_hut,3) elseif (TDC_qADC(id_hut,2) .gt. 0.) then Conversion_qPhot_HG(i,j) = C_hut_count(i,j) / * TDC_qADC(id_hut,2) Conversion_qPhot_LG(i,j) = C_hut_count(i,j) / * TDC_qADC(id_hut,3) endif 1010 CONTINUE c ====================================================================== c write statistics for Cerenkov photons to disk c ====================================================================== c WRITE( Lun_C_STAT ,'(/)') WRITE( Lun_C_STAT ,*) ' NUMBER OF C-PHOTONS IN C-HUTS /M**2' WRITE( Lun_C_STAT ,*) ' NOT INCLUDING ATMOSPHERIC ABSORPTION ' WRITE( Lun_C_STAT ,'()') kmax = (C_ArrayHuts - imiddle) / 3 DO 1020 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7010) & ( int( C_hut_count(i,j) ), * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1020 CONTINUE 7010 FORMAT(40i8) C WRITE(Lun_C_STAT,'(/)') WRITE(Lun_C_STAT,'(a,a)') ' MEAN ARRIVAL TIME OF CERENKOV ', * 'PHOTONS IN SQUARES (NS)' WRITE( Lun_C_STAT ,*) ' NOT INCLUDING ATMOSPHERIC ABSORPTION ' WRITE( Lun_C_STAT ,'()') DO 1030 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7020) ( C_hut_arrTime(i,j) , * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1030 CONTINUE c 7020 FORMAT(40f8.1,/) IF ( iAtmos_absor .eq. 1 ) THEN c WRITE( Lun_C_STAT ,'(//)') WRITE(Lun_C_STAT,*) ' NUMBER OF C-PHOTONS IN C-HUTS /M**2' write(Lun_C_STAT,*) ' INCLUDING ATMOSPHERIC ABSORPTION' write( Lun_C_STAT ,'()') DO 1040 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7010) * ( int( C_hut_count_At(i,j) ), * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1040 CONTINUE C WRITE( Lun_C_STAT ,'(/)') WRITE(Lun_C_STAT,'(a,a)') ' MEAN ARRIVAL TIME OF CERENKOV', * ' PHOTONS IN SQUARES (NS)' WRITE( Lun_C_STAT ,*) ' INCLUDING ATMOSPHERIC ABSORPTION' write( Lun_C_STAT ,'()') DO 1050 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7020) * ( C_hut_arrTime_At(i,j), * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1050 CONTINUE ENDIF !END OF if ( iAtmos_absor .eq. 1 ) c ====================================================================== c Write out conversion factors (pC -> C-density) c for several counters c ====================================================================== WRITE( Lun_C_STAT ,'(//)') WRITE( Lun_C_STAT ,*) ' CONVERSION FACTOR FOR', * ' CHARGE (pC) --> PHOTONS/M2' write( Lun_C_STAT ,'(/)') WRITE( Lun_C_STAT ,*) ' 1) HIGH GAIN BRANCH' write( Lun_C_STAT ,'()') DO 1060 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7020) * ( Conversion_qPhot_HG(i,j), * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1060 CONTINUE write( Lun_C_STAT ,'(/)') WRITE( Lun_C_STAT ,*) ' 2) LOW GAIN BRANCH' write( Lun_C_STAT ,'()') DO 1070 i = imiddle+3*kmax, imiddle-3*kmax, -3 WRITE( Lun_C_STAT , 7020) * ( Conversion_qPhot_LG(i,j), * j = jmiddle-3*kmax, jmiddle+3*kmax, 3 ) 1070 CONTINUE c ====================================================================== c Get center of gravity of the C-light disk (exclude central hut): c (clight_xcore, clight_ycore) c ====================================================================== clight_xcore = 0.0 clight_ycore = 0.0 clight_photons = 0.0 DO i=1, C_ArrayHuts yhut = (i - imiddle) * C_ArrayGrid DO j=1, C_ArrayHuts xhut = (j - jmiddle) * C_ArrayGrid if (iAtmos_absor .eq. 1) then weight = C_hut_count_At(i,j) else weight = C_hut_count(i,j) endif clight_xcore = clight_xcore + xhut * weight clight_ycore = clight_ycore + yhut * weight clight_photons = clight_photons + weight ENDDO ENDDO c c Remove the central hut (since it has a large weight) c if (iAtmos_absor .eq. 1) then clight_photons = clight_photons - * C_hut_count_At( imiddle, jmiddle ) else clight_photons = clight_photons - * C_hut_count( imiddle, jmiddle ) endif c clight_xcore = clight_xcore / clight_photons clight_ycore = clight_ycore / clight_photons c ===================================================================== c Get center of gravity of the AIROBICC signals (exclude central hut) in: c the extended array: (qADC_xcore, qADC_ycore) c ===================================================================== qADC_xcore = 0.0 qADC_ycore = 0.0 qADC_signal = 0.0 DO i=1, C_ArrayHuts yhut = (i - imiddle) * C_ArrayGrid DO j=1, C_ArrayHuts xhut = (j - jmiddle) * C_ArrayGrid id_hut = C_ArrayHuts * (i-1) + j weight = TDC_qADC( id_hut , 2) qADC_xcore = qADC_xcore + xhut * weight qADC_ycore = qADC_ycore + yhut * weight qADC_signal = qADC_signal + weight ENDDO ENDDO c ----------------------------------------------------- c Remove the central hut (since it has a large weight) c qADC_signal = qADC_signal - * TDC_qADC( C_ArrayHuts**2 / 2 + 1 , 2) if (qADC_signal .gt. 0.0) then qADC_xcore = qADC_xcore / qADC_signal qADC_ycore = qADC_ycore / qADC_signal else qADC_xcore = -10000. qADC_ycore = -10000. endif c ====================================================================== write(Lun_C_STAT, '(//)') WRITE(Lun_C_STAT, 7025) x_core, y_core WRITE(Lun_C_STAT, 7030) clight_xcore, clight_ycore WRITE(Lun_C_STAT, 7040) qADC_xcore, qADC_ycore 7025 FORMAT('true xcore (m)= ',f10.1,5x,'true ycore (m)= ',f10.1) 7030 FORMAT('clight_xcore (m)= ',f10.1,5x,'clight_ycore (m)= ',f10.1) 7040 FORMAT('qADC_xcore (m)= ',f10.1,5x,'qADC_ycore (m)= ',f10.1) CLOSE ( Lun_C_STAT ) c ====================================================================== CALL hfill(6001, clight_xcore - x_core, 0.,1.) CALL hfill(6002, clight_ycore - y_core, 0.,1.) CALL hfill(6003, qADC_xcore - x_core, 0., 1.) CALL hfill(6004, qADC_ycore - y_core, 0., 1.) c 2000 CONTINUE c c ================== END OF SHOWER LOOP ===================== c CLOSE( Lun_Out ) CALL hcdir( '//hbook', ' ' ) DO id=6001, 6004 CALL hrout( id, icycle, ' ' ) ENDDO CALL hrend( 'hbook' ) CLOSE( Lun_HBOOK ) END c c ********************************************************************** SUBROUTINE Init_AIROBICC c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 c c Reads inputs, initializes variables and histograms c ********************************************************************** C ============================ c HBOOK parameter ( mempaw = 300 000 ) common /pawc/ hpaw(mempaw) common /histog/ nC_bins , C_TimeWindow , C_time_shift common /histog2/ nbins_PM , TimeWindow_PM , time_shift_PM common /histog3/ nbins_Amp1_2 , TimeWindow_Amp1_2 common /histog4/ Window_shift C ============================ c EXTENDED ARRAY common /HEGRA_level/ H_LAPALMA common /C_Array/ C_ArrayHuts , C_ArrayGrid , C_ArraySide , * imiddle , jmiddle, C_HutSide integer C_ArrayHuts C ============================ C HUT GEOMETRY common /Chutgeom/ r_cone(0:3),zu(0:3),ta(3),ca(3),nz(3),zd(3) real nz common /sphere_PM/ r_sph,zu_sph,sph1,sph2 C ============================ C CERENKOV BUNCHES IN AIROBICC COUNTER parameter ( nlambdas = 4 ) common /cerenkov/ C_spectrum(nlambdas) , lambda(nlambdas) real lambda COMMON /cerenkov_2/ extended_spectrum_factor COMMON /TRAtmos/ TR_Atmos(nlambdas) , iAtmos_absor common /TRCounter/ TR_Plexiglas(nlambdas), & Reflex_mirror(nlambdas), & TR_Filter(nlambdas), & QE_Photocath(nlambdas) C ============================ C PM common /poiscalls/ npois_calls common /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM common /TrTimePM/ nTrTime , TrTime_PM , Dist_TrTime_PM dimension TrTime_PM(100) , Dist_TrTime_PM(100) C ============================ C ELECTRONICS common / Amplifiers / Gain_Amp1 , Gain_Amp1_2 , FWHM_CATS * , CATS_Dist( 200 ) , nbins_CATS common / CFD / CFD_delay1 , CFD_delay2 , U_Threshold, Resist_CFD parameter ( nhuts = 2000 ) common /ADC_TDC / q_ADC1( nhuts ), q_ADC1_2( nhuts ) * , TDC1_2( nhuts ) , Gate_qADC1_2 c ==================================================================== c data H_LAPALMA / 2.2e5 / ! cm data r_cone / 10. , 12.35 , 15. , 20. / data zu / 0. , 2.6 , 10.1 , 46.7 / data r_sph / 11. / , zu_sph / 6.418 / c Values taken from figure (3.6) in Albrecht's Thesis (p. 44) c data lambda / 325. , 375. , 425. , 475. / data TR_Plexiglas / 0.8 , 0.88 , 0.9 , 0.9 / DATA Reflex_mirror / 0.86 , 0.9 , 0.9 , 0.92 / DATA TR_Filter / 0.89 , 0.91 , 0.85 , 0.5 / DATA QE_Photocath / 0.19 , 0.26 , 0.22 , 0.15 / c data C_TimeWindow / 100. / data electron_charge / 1.6e-19 / data nTrTime / 30 / data ( Dist_TrTime_PM(i) , i = 1 , 5 ) / 0.0 , 0.0 , 0.0 , 0.0 * , 0.0 / data ( Dist_TrTime_PM(i) , i = 6 , 10 ) / 0.4 , 1.0 , 2.1 , 4.9 * , 7.9 / data ( Dist_TrTime_PM(i) , i = 11 , 15 ) / 9.0 , 9.7 , 10.0 * , 9.9 , 9.6 / data ( Dist_TrTime_PM(i) , i = 16 , 20 ) / 9.0 , 7.6 , 5.5 * , 4.0 , 2.9 / data ( Dist_TrTime_PM(i) , i = 21 , 25 ) / 2.1 , 1.5 , 1.0 * , 0.7 , 0.5 / data ( Dist_TrTime_PM(i) , i = 26 , 30 ) / 0.3 , 0.2 , 0.1 * , 0.0 , 0.0 / data TimeWindow_PM / 800. / data Transm_cable / 0.6 / data TimeWindow_Amp1_2 / 600. / ! duration of linear gate for ADCs data FWHM_CATS / 7. / data Gain_Amp1 / 25. / , Gain_Amp1_2 / 500. / ! Gain in Charge data Resist_CFD / 50. / ! Ohm data Gate_qADC1_2 / 30. / c ////////////////////////////////////////////////////////////////////// c ================================================================== c Read input variables specified by users from LUN = Lun_Inp c ================================================================== call Read_DataCards c ------------------------------------------------ c Normalize the Cerenkov spectrum sum_lambda_2 = 0. do i=1,nlambdas sum_lambda_2 = sum_lambda_2 + 1. / lambda(i) ** 2 enddo do i=1,nlambdas C_spectrum(i) = 1. / sum_lambda_2 / lambda(i) ** 2 if (iAtmos_absor .eq. 0) TR_Atmos(i) = 1 enddo c ------------------------------------------------------- c Correct for change in C-spectral range: 300-450 nm -> 300-500 nm c extended_spectrum_factor = 1. / (1. - C_spectrum(4)) c ----------------------------------------------------- c calculate geometrical parameters of the C-hut mirrors c do i = 1,3 a = r_cone(i-1) * zu(i) - r_cone(i) * zu(i-1) b = r_cone(i) - r_cone(i-1) zd(i) = a / b ta(i) = r_cone(i) / ( zu(i) + zd(i) ) ca(i) = 1 / sqrt ( 1 + ta(i) ** 2 ) nz(i) = ta(i) * ca(i) enddo c ------------------------------------------------------ c calculate geometrical parameters for the Photocathode c zd_sph = r_sph - zu_sph sph1 = 2 * zd_sph sph2 = zd_sph ** 2 - r_sph ** 2 c ------------------------------------------------ c Normalize the TrTime distribution c sum_DistTrTime = 0. do i=1,nTrTime TrTime_PM(i) = i * 0.5 - 0.25 sum_DistTrTime = sum_DistTrTime + Dist_TrTime_PM(i) enddo do i=1,nTrTime Dist_TrTime_PM(i) = Dist_TrTime_PM(i) / sum_DistTrTime enddo c ------------------------------------- c Warm up poisson & rndm generators c do i=1, npois_calls call poissn( 1000. , n , ierror ) xxx = rndm() enddo c -------------------------------------------------- c factor to convert from No of electrons/bin to anode current c in microAmpere (bin width = 0.5 ns) current_factor = electron_charge * 2 * 1e9 * 1e6 c ----------------------------------------------------- c Get Average No. of photoelectrons/bin at photocathode, due to NSL Av_NSLNoise = Av_NSLCurrent / current_factor / Gain_PM c --------------------------------------- Gain_Amp1 = Gain_Amp1 * Transm_cable * 0.3 Gain_Amp1_2 = Gain_Amp1_2 * Transm_cable * 0.6 c --------------------------------------------------- c Get the Cables+Amplifiers Time Spread (CATS) distribution c Each bin is 1/2 ns wide nbins_CATS = FWHM_CATS * 2 * 3 sigma_CATS = FWHM_CATS / 2 / sqrt( alog(4.) ) sum_CATS_Dist = 0. do i=1, nbins_CATS xx = (i - nbins_CATS/2)/2 /sigma_CATS xx = xx ** 2 / 2 CATS_Dist(i) = exp( -xx ) sum_CATS_Dist = sum_CATS_Dist + CATS_Dist(i) enddo c c Normalize the CATS distribution do i=1, nbins_CATS CATS_Dist(i) = CATS_Dist(i) / sum_CATS_Dist enddo c ------------------------------------------- c Init HBOOK memory call hlimit(mempaw) c --------------------------------------------------------------- c C_TimeWindow: time width of histograms for anode current due to c C-pulses (ns) nC_bins = C_TimeWindow * 2 c Time window for recording the anode current due to C-light is opened c `C_time_shift' ns before arrival of shower tangent plane to the c corresponding counter. C_time_shift = 10. c ------------------------------------------------------------- c TimeWindow_PM: time width of histogram for storing anode current c due to C-pulse (ns) + NSL contribution nbins_PM = TimeWindow_PM * 2. c Time window of histogram containing total anode current starts c `time_shift_PM' ns before arrival of shower front to the c corresponding counter. time_shift_PM = TimeWindow_PM / 2. c -------------------------------------------------------------- c TimeWindow_Amp1_2: time width of histogram for total current c after 2nd amplifier nbins_Amp1_2 = TimeWindow_Amp1_2 * 2 c Delay between the opening of the time window for recording anode c current and time window to integrate pulse in ADCs Window_shift = (TimeWindow_PM - TimeWindow_Amp1_2) / 2 c -------------------------------------------------------- c Book histograms for studying core reconstruction CALL hbook1(6001, 'clight_dx_core(m)', 200, -40., 40., 0.) CALL hbook1(6002, 'clight_dy_core(m)', 200, -40., 40., 0.) CALL hbook1(6003, 'qADC_dx_core(m)', 200, -40., 40., 0.) CALL hbook1(6004, 'qADC_dy_core(m)', 200, -40., 40., 0.) DO id = 6001, 6004 CALL hidopt(id, 'STAT') ENDDO c RETURN END c c ********************************************************************** SUBROUTINE Read_DataCards c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 c c Read input file airobicc.in c ********************************************************************** COMMON / C__header / C_header(20), corsika_version COMMON /datac1/ InpDirectory , OutDirectory , dsn0 CHARACTER*60 InpDirectory , OutDirectory CHARACTER*9 dsn0 COMMON /datac2/ i_C_file1 , n_C_files COMMON /LUNs/ Lun_Inp , Lun_Out , Lun_CORS_Cfile COMMON/Cfile/ cerbuf(273,30) , irecor , n_C_blks , isb , ieof, lh COMMON/IncPart_Inputs/ InpPartType, InpEnergy, InpTheta, InpPhi REAL InpPartType , InpEnergy(2) , InpTheta(2) , InpPhi(2) COMMON /check/ iCheck_PartType, iCheck_Energy, iCheck_Theta, & iCheck_Phi common /Array_Inputs/ InpArraySide, InpArrayGrid, InpHutSide real InpArraySide, InpArrayGrid, InpHutSide COMMON /check_2/ iCheck_Array, iCheck_HutSide COMMON /C_Array/ C_ArrayHuts , C_ArrayGrid , C_ArraySide , * imiddle , jmiddle, C_HutSide INTEGER C_ArrayHuts PARAMETER ( nlambdas = 4 ) COMMON /TRAtmos/ TR_Atmos(nlambdas) , iAtmos_absor COMMON /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM COMMON /poiscalls/ npois_calls COMMON / CFD / CFD_delay1 , CFD_delay2 , U_Threshold, Resist_CFD CHARACTER*10 KEYW CHARACTER*80 LINE c ///////////////////////////////////////////////////////////////////////// corsika_version = 3.02 DO i=1,60 InpDirectory(i:i) = ' ' OutDirectory(i:i) = ' ' ENDDO dsn0 = 'CER ' i_C_file1 = 1 n_C_files = 1 n_C_blks = 30 InpPartType = 1. InpEnergy(1) = 1e2 InpEnergy(2) = 1e2 InpTheta(1) = 10. InpTheta(2) = 10. InpPhi(1) = 0. InpPhi(2) = 360. iCheck_PartType = 0 iCheck_Energy = 0 iCheck_Theta = 0 iCheck_Phi = 0 InpHutSide = 1. InpArrayGrid = 15. InpArraySide = 390. iCheck_Array = 0 iCheck_HutSide = 0 iAtmos_absor = 1 Av_NSLCurrent = 32. Gain_PM = 8e3 CFD_delay1 = 5. CFD_delay2 = 10. U_Threshold = -1. c ///////////////////////////////////////////////////////////////////// OPEN ( unit=Lun_Inp , file='airobicc.in', form='FORMATTED', * access='SEQUENTIAL',status='OLD') c C GET A NEW INPUT LINE AND PRINT IT c 1 CONTINUE READ ( Lun_Inp , 500 , END=1000) LINE 500 FORMAT(A80) c C GET TYPE OF DATACARD c READ ( LINE(1:10) , 100) KEYW 100 FORMAT(A10) C----------------------------------------------------------------------* C INTERPRET DATA CARD AND READ PARAMETERS C END OF DATA CARD INPUT IF ( KEYW .EQ. 'EXIT ' ) THEN CLOSE( Lun_Inp ) ind_O = index( OutDirectory , ' ' ) IF (ind_O .eq. 1) OutDirectory = InpDirectory RETURN C GET Input Directory ELSEIF ( KEYW .EQ. 'INPUT_DIR ' ) THEN READ(LINE(1:71),101) KEYW , InpDirectory C GET Output Directory ELSEIF ( KEYW .EQ. 'OUTPUT_DIR' ) THEN READ(LINE(1:71),101) KEYW , OutDirectory C GET leading characters of Cerenkov files ELSEIF ( KEYW .EQ. 'DSN0 ' ) THEN READ(LINE(1:19),104) KEYW , dsn0 C GET TYPE OF PRIMARY PARTICLE ELSEIF ( KEYW .EQ. 'INC_PARTYP' ) THEN READ(LINE(1:20),102) KEYW , InpPartType C GET ENERGY OF PRIMARY PARTICLE ELSEIF ( KEYW .EQ. 'INC_ENERGY' ) THEN READ(LINE(1:30),107) KEYW , InpEnergy C GET THETA OF PRIMARY PARTICLE ELSEIF ( KEYW .EQ. 'INC_THETA ' ) THEN READ(LINE(1:30),108) KEYW , InpTheta C GET PHI OF PRIMARY PARTICLE ELSEIF ( KEYW .EQ. 'INC_PHI ' ) THEN READ(LINE(1:30),108) KEYW , InpPhi C GET first event number ELSEIF ( KEYW .EQ. 'FIRSTEVENT' ) THEN READ(LINE(1:20),103) KEYW , i_C_file1 C GET number of events being processed ELSEIF ( KEYW .EQ. 'NUM_EVENTS' ) THEN READ(LINE(1:20),103) KEYW , n_C_files C GET number of blocks of the Cerenkov_file from CORSIKA ELSEIF ( KEYW .EQ. 'NUM_C_BLKS' ) THEN READ(LINE(1:20),103) KEYW , n_C_blks C GET length of side for extended C_Array ELSEIF ( KEYW .EQ. 'CARRAYSIDE' ) THEN READ(LINE(1:20),102) KEYW , InpArraySide C GET the grid spacing of the extended C_Array ELSEIF ( KEYW .EQ. 'CARRAYGRID' ) THEN READ(LINE(1:20),102) KEYW , InpArrayGrid C GET the side length of Chut assumed in CORSIKA simul. ELSEIF ( KEYW .EQ. 'C_HUT_SIDE' ) THEN READ(LINE(1:20),102) KEYW , InpHutSide C GET flag for use of atmospheric absorption ELSEIF ( KEYW .EQ. 'ATMOSABSOR' ) THEN READ(LINE(1:20),103) KEYW , iAtmos_absor C GET flag for checking header & incident particle properties c of Cerenkov file ELSEIF ( KEYW .EQ. 'CHECK_PART' ) THEN READ(LINE(1:50),106) KEYW , iCheck_PartType, iCheck_Energy, * iCheck_Theta, iCheck_Phi C GET flag for checking C-array parameters ELSEIF ( KEYW .EQ. 'CHECKARRAY' ) THEN READ(LINE(1:50),109) KEYW , iCheck_Array, iCheck_HutSide C GET number of calls to warm up poisson generator ELSEIF ( KEYW .EQ. 'POIS_CALLS' ) THEN READ(LINE(1:20),103) KEYW , npois_calls C GET average current at anode of PM from NSL ELSEIF ( KEYW .EQ. 'NSLCURRENT' ) THEN READ(LINE(1:20),102) KEYW , Av_NSLCurrent C GET gain of PM ELSEIF ( KEYW .EQ. 'GAIN_PM ' ) THEN READ(LINE(1:20),105) KEYW , Gain_PM C GET 1st delay factor of CFD ELSEIF ( KEYW .EQ. 'CFD_DELAY1' ) THEN READ(LINE(1:20),102) KEYW , CFD_delay1 C GET 2nd delay factor of CFD ELSEIF ( KEYW .EQ. 'CFD_DELAY2' ) THEN READ(LINE(1:20),102) KEYW , CFD_delay2 C GET U_threshold of CFD (in miliVolts) ELSEIF ( KEYW .EQ. 'UTHRESHOLD' ) THEN READ(LINE(1:20),102) KEYW , U_Threshold C GET CORSIKA version ELSEIF ( KEYW .EQ. 'CORS_VERS ' ) THEN READ(LINE(1:20),102) KEYW , corsika_version C DUMMY LINE ELSEIF ( KEYW .EQ. ' ' ) THEN C ILLEGAL KEYWORD ELSE WRITE(Lun_Out, *) 'Read_DataCards: UNKNOWN KEYWORD :', KEYW ENDIF GOTO 1 C----------------------------------------------------------------------* 1000 CONTINUE 101 FORMAT(A10,1X,A60) 102 FORMAT(a10,f10.1) 103 FORMAT(a10,i10) 104 FORMAT(A10,A9) 105 FORMAT(a10,e10.2) 106 FORMAT(a10,4i10) 107 FORMAT(a10,2e10.2) 108 FORMAT(a10,2f10.2) 109 FORMAT(a10,2i10) RETURN END c c ********************************************************************** SUBROUTINE Read_Cfile c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 c c Read one record of the file CER# c ********************************************************************** COMMON/Cfile/ cerbuf(273,30) , irecor , n_C_blks , isb , ieof, lh COMMON /check/ iCheck_PartType, iCheck_Energy, iCheck_Theta, & iCheck_Phi COMMON /LUNs/ Lun_Inp , Lun_Out , Lun_CORS_Cfile c C_header contains global info of the shower COMMON / C__header / C_header(20), corsika_version equivalence (C_header(1),PartType), (C_header(2),Energy) equivalence (C_header(3),Theta), (C_header(4),Phi) equivalence (C_header(5),height_1stint) equivalence (C_header(6),x_core), (C_header(7),y_core) equivalence (C_header(8),Huts_X), (C_header(9),ArrayGrid) equivalence (C_header(10),HutSide) equivalence (C_header(11),Angle_X_North) equivalence (C_header(12),cersiz) equivalence (C_header(13),vernum) data pi / 3.141593 / * ///////////////////////////////////////////////////////////////// degrad = pi / 180. irecor = irecor + 1 READ ( Lun_CORS_Cfile , rec = irecor , err = 10) * ((cerbuf(i,isb),i=1,273),isb=1,n_C_blks) isb = 1 IF ( irecor .EQ. 1 ) THEN c Get header for output data file c if (corsika_version .le. 4.11) then if (corsika_version .eq. 3.02) then write(Lun_Out,*) 'WARNING: In this version of corsika' write(Lun_Out,*) 'the shower core is located in a hut' write(Lun_Out,*) 'and you will get unrealistic results' write(Lun_Out,*) 'when reconstructing the shower ' write(Lun_Out,*) 'parameters.' c n_C_blks := 30 PartType = cerbuf(1,1) Energy = cerbuf(2,1) Theta = cerbuf(3,1) * degrad Phi = cerbuf(4,1) * degrad height_1stint = -1. x_core = 0. y_core = 0. Huts_X = 27. ArrayGrid = 15. HutSide = 1. Angle_X_North = 0. cersiz = cerbuf(5,1) vernum = cerbuf(6,1) / 100. isb = 1 lh = 1 elseif (corsika_version .eq. 4.06) then PartType = cerbuf(3,1) Energy = cerbuf(4,1) Theta = cerbuf(11,1) Phi = cerbuf(12,1) height_1stint = cerbuf(7,1) x_core = cerbuf(44,1) / 100. y_core = cerbuf(45,1) / 100. ArrayGrid = cerbuf(11,2) Huts_X = 2 * int(cerbuf(10,2) / 2 / ArrayGrid) + 1 HutSide = cerbuf(12,2) Angle_X_North = cerbuf(46,1) ! rad cersiz = cerbuf(85,1) vernum = cerbuf(5,2) c c Skip 1st block (Corsika EVTH) and 14 words in 2nd block c isb = 2 lh = 2 elseif (corsika_version .eq. 4.11) then write(Lun_Out,*) 'WARNING: In this version of corsika' write(Lun_Out,*) 'the shower core is located in a hut' write(Lun_Out,*) 'and you will get unrealistic results' write(Lun_Out,*) 'when reconstructing the shower ' write(Lun_Out,*) 'parameters.' PartType = cerbuf(3,1) Energy = cerbuf(4,1) Theta = cerbuf(11,1) Phi = cerbuf(12,1) height_1stint = cerbuf(7,1) x_core = 0. y_core = 0. Huts_X = cerbuf(86,1) ArrayGrid = cerbuf(88,1) / 100 HutSide = cerbuf(90,1) / 100 Angle_X_North = 0. ! IT IS NOT RECORDED IN CORSIKA !!! cersiz = cerbuf(85,1) vernum = 4.11 c c Skip 1st block (Corsika EVTH) c isb = 2 lh = 0 else write(Lun_Out,'(a,a)') 'THIS CORSIKA VERSION IS NOT ', * 'IMPLEMENTED (YET).' stop endif else write(Lun_Out,*) 'THIS CORSIKA VERSION IS NOT ' write(Lun_Out,*) 'IMPLEMENTED (YET).' stop endif !END OF if (corsika_version if (Phi .gt. 2*pi) Phi = Phi - 2*pi if (Phi .lt. 0.) Phi = Phi + 2*pi CALL Check_Inputs ENDIF !END OF if ( (irecor .eq. 1) .and. c GOTO 20 10 ieof = 1 20 RETURN END c c ************************************************************* SUBROUTINE Check_Inputs c ************************************************************* c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 C C Check energy, theta and type of incident particle c Also, get and check array parameters c ************************************************************* C COMMON/IncPart_Inputs/ InpPartType, InpEnergy, InpTheta, InpPhi REAL InpPartType , InpEnergy(2) , InpTheta(2) , InpPhi(2) COMMON /check/ iCheck_PartType, iCheck_Energy, iCheck_Theta, & iCheck_Phi common /Array_Inputs/ InpArraySide, InpArrayGrid, InpHutSide real InpArraySide, InpArrayGrid, InpHutSide COMMON /check_2/ iCheck_Array, iCheck_HutSide COMMON /LUNs/ Lun_Inp , Lun_Out , Lun_CORS_Cfile c C_header contains global info of the shower, and is partially filled c in read_Cfile when reading the 1st record of the file CER# COMMON / C__header / C_header(20), corsika_version equivalence (C_header(1),PartType), (C_header(2),Energy) equivalence (C_header(3),Theta), (C_header(4),Phi) equivalence (C_header(5),height_1stint) equivalence (C_header(6),x_core), (C_header(7),y_core) equivalence (C_header(8),Huts_X), (C_header(9),ArrayGrid) equivalence (C_header(10),HutSide) equivalence (C_header(11),Angle_X_North) equivalence (C_header(12),cersiz) equivalence (C_header(13),vernum) DATA pi / 3.141593 / c ++++++++++++++++++++++++++++++++++++ raddeg = 180 / pi degrad = pi / 180 c Convert input angle from degrees to radians; input energy from TeV to GeV do i=1, 2 InpEnergy(i) = InpEnergy(i) * 1e3 InpTheta(i) = InpTheta(i) * degrad InpPhi(i) = InpPhi(i) * degrad enddo if ( iCheck_PartType .ne. 1 ) goto 5 CALL compare( PartType, InpPartType, InpPartType, JFLAG ) IF ( JFLAG . ne . 0 ) THEN WRITE( Lun_Out ,'(a,a)')' the incident particle type is not', * ' what was expected' write( Lun_Out ,*)' Instead , it is ', PartType STOP endif 5 continue if ( iCheck_Energy .ne. 1 ) goto 10 call compare( Energy, InpEnergy(1), InpEnergy(2), JFLAG ) if (JFLAG. ne. 0) then write( Lun_Out ,'(a,a)')'the incident energy is not in the ', * 'expected range' write( Lun_Out ,*)' It is ', Energy, ' GeV' STOP endif 10 continue if ( iCheck_Theta .ne. 1 ) goto 20 call compare( Theta, InpTheta(1), InpTheta(2), JFLAG ) if (JFLAG. ne. 0) then write(Lun_Out,'(a,a)') 'the incident angle theta is not in ', * ' the expected range' write( Lun_Out ,*)' Instead , it is ', Theta*raddeg, ' deg' STOP endif 20 continue if ( iCheck_Phi .ne. 1 ) goto 30 call compare( Phi, InpPhi(1), InpPhi(2), JFLAG ) if (JFLAG. ne. 0) then write( Lun_Out,'(a,a)') 'the incident angle phi is not in ', * 'the expected range' write( Lun_Out ,*)' Instead , it is ', Phi*raddeg, ' deg' STOP endif 30 continue if ( iCheck_Array .ne. 1 ) goto 40 ArraySide = (Huts_X - 1) * ArrayGrid call compare(ArraySide, InpArraySide, InpArraySide, JFLAG) if (JFLAG. ne. 0) then write( Lun_Out , '(a,a)') & ' The side length of the C-array used in CORSIKA', & ' is not the expected one' write( Lun_Out ,*)' Instead , it is ',ArraySide, ' m' STOP endif call compare(ArrayGrid, InpArrayGrid, InpArrayGrid, JFLAG) if (JFLAG. ne. 0) then write( Lun_Out , '(a,a)') & ' The grid spacing of the C-array used in CORSIKA', & ' is not the expected one' write( Lun_Out ,*)' Instead , it is ',ArrayGrid, ' m' STOP endif 40 continue if ( iCheck_HutSide .ne. 1 ) goto 50 call compare( HutSide, InpHutSide, InpHutSide, JFLAG ) if (JFLAG. ne. 0) then write( Lun_Out , '(a,a)') & ' The side length of the Chut used in CORSIKA', & ' is not the expected one' write( Lun_Out ,*)' Instead , it is ', HutSide, ' m' STOP endif 50 continue return end c c ********************************************************************** SUBROUTINE compare(a,b1,b2,JFLAG) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 C C COMPARES THE VALUES a AND b AND DECIDES IF THEY ARE EQUAL IF THEY C DIFFER IN LESS THAN SOME PER CENT ( JFLAG = 0 ) c ********************************************************************** data RelativeError / 1.e-4 / C JFLAG = 0 c1 = b1 * (1.0 - RelativeError) c2 = b2 * (1.0 + RelativeError) if ( (a .lt. c1) .or. (a .gt. c2) ) JFLAG = 1 return end c c ********************************************************************** SUBROUTINE Init_Histos c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 c c Book histos for the anode current due to C-pulse c ********************************************************************** c C_header contains global info of the shower, and is partially filled c in read_Cfile when reading the 1st record of the file CER# COMMON / C__header / C_header(20), corsika_version equivalence (C_header(1),PartType), (C_header(2),Energy) equivalence (C_header(3),Theta), (C_header(4),Phi) equivalence (C_header(5),height_1stint) equivalence (C_header(6),x_core), (C_header(7),y_core) equivalence (C_header(8),Huts_X), (C_header(9),ArrayGrid) equivalence (C_header(10),HutSide) equivalence (C_header(11),Angle_X_North) equivalence (C_header(12),cersiz) equivalence (C_header(13),vernum) common /C_Array/ C_ArrayHuts , C_ArrayGrid , C_ArraySide , * imiddle , jmiddle, C_HutSide integer C_ArrayHuts common /EASfront/ X_delay1 , Y_delay1, X_delay2, Y_delay2 common /histog/ nC_bins , C_TimeWindow , C_time_shift character*40 HistTitle logical hexist real LightSpeed c ---------------------------- data pi / 3.141593 / data LightSpeed / 0.3 / ! m/ns c //////////////////////////////////////////////////////////////////// c Variables to characterize extended array ArraySide = (int(Huts_X) - 1) * ArrayGrid C_ArraySide = ArraySide C_ArrayGrid = ArrayGrid C_ArrayHuts = 1 + int(C_ArraySide / C_ArrayGrid) imiddle = ( C_ArrayHuts + 1 ) / 2 jmiddle = imiddle C_HutSide = HutSide c c Book histograms for C_pulse timing c Theta_Inc = Theta Phi_Inc = Phi + pi c u_Inc = sin( Theta_Inc ) * cos( Phi_Inc ) X_delay1 = - C_ArrayGrid * u_Inc / LightSpeed X_delay2 = x_core * u_Inc / LightSpeed v_Inc = sin( Theta_Inc ) * sin( Phi_Inc ) Y_delay1 = - C_ArrayGrid * v_Inc / LightSpeed Y_delay2 = y_core * v_Inc / LightSpeed c HistTitle = ' PM-Current x = 111111 m , y = 111111 m$' c DO 10 j=1,C_ArrayHuts jhut = j X_tmin = X_delay1 * ( jhut - jmiddle ) + X_delay2 xhut = C_ArrayGrid * ( jhut - jmiddle ) WRITE(HistTitle(17:22),'(f6.1)') xhut DO 10 i=1,C_ArrayHuts ihut = i Y_tmin = Y_delay1 * ( ihut - imiddle ) + Y_delay2 tmin = X_tmin + Y_tmin - C_time_shift tmax = tmin + C_TimeWindow id_hut = C_ArrayHuts * ( ihut - 1 ) + jhut yhut = C_ArrayGrid * ( ihut - imiddle ) WRITE(HistTitle(32:37),'(f6.1)') yhut if ( hexist(id_hut) ) call hdelet( id_hut) CALL hbook1( id_hut , HistTitle , nC_bins , * tmin , tmax , 0. ) 10 CONTINUE c RETURN END c c ********************************************************************** SUBROUTINE Atmos_absorption( height , cos_Theta ) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Calculate atmospheric absorption (Rayleigh+aerosol) for each interval c of the C-photon spectrum c ********************************************************************** parameter ( nlambdas = 4 ) common /TRAtmos/ TR_Atmos(nlambdas) , iAtmos_absor common /cerenkov/ C_spectrum(nlambdas) , lambda(nlambdas) real lambda , L_Mie common /HEGRA_level/ H_LAPALMA c ///////////////////////////////////////////////// data H_Rayleigh / 7.1e5 / data Rho_0 / 1.29e-3 / data XR_400nm / 2970. / c data L_Mie / 14e5 / data H_Mie / 1.2e5 / c ------------------------------------------------- c calculate absorption due to scattering Rayleigh c do i=1,nlambdas f1 = Rho_0 / XR_400nm * ( 400. / lambda(i) ) ** 4 * * H_Rayleigh / cos_Theta TR_Atmos(i) = exp( f1 * ( exp( -height / H_Rayleigh ) - * exp( -H_LAPALMA / H_Rayleigh ) ) ) enddo c ------------------------------------------------ c calculate absorption due to aerosol scattering c f2 = H_Mie / L_Mie / cos_Theta TRAtmos_Mie = exp( f2 * ( exp( -height / H_Mie ) - * exp( -H_LAPALMA / H_Mie ) ) ) c do i=1,nlambdas TR_Atmos(i) = TR_Atmos(i) * TRAtmos_Mie enddo c return end c c ********************************************************************** SUBROUTINE Ctracking(x,y,u,v,track_length,nreflex) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Subroutine to track a Cerenkov photon in the Winston cone down to c the PM c ********************************************************************** common /Chutgeom/ r(0:3),zu(0:3),ta(3),ca(3),nz(3),zd(3) real nx,ny,nz,ni0 common /sphere_PM/ r_sph,zu_sph,sph1,sph2 common /st_line/ a2,b2,c2 double precision a2 , b2 , c2 double precision a , b , c , discrim double precision ax , bx , ay , by double precision u0 , v0 , w0 , u1 , v1 , w1 double precision norm c data epsilon / 1e-3 / radius = sqrt( x ** 2 + y ** 2 ) c if ( radius . gt . r(3) ) then track_length = 0. return endif c uv2 = u ** 2 + v ** 2 if ( uv2 . gt . 1. ) then track_length = -1. return endif c tl = 0. ! tracklength nreflex = 0 ic = 3 ! number of conical sections x0 = x y0 = y u0 = u v0 = v z0 = zu(3) w0 = - dsqrt( 1. - u0 ** 2 - v0 ** 2 ) c 10 ax = x0 - u0 * z0 / w0 ay = y0 - v0 * z0 / w0 bx = u0 / w0 by = v0 / w0 a2 = bx ** 2 + by ** 2 b2 = 2 * ( ax * bx + ay * by ) c2 = ax ** 2 + ay ** 2 call hit_sphere(zu_sph,zd_sph) if ( zu_sph . ge . 0. ) then if ( z0 . ge . zu_sph ) then z1 = zu_sph else z1 = zd_sph endif x1 = ax + bx * z1 y1 = ay + by * z1 dx = x1 - x0 dy = y1 - y0 dz = z1 - z0 step = sqrt ( dx ** 2 + dy ** 2 + dz ** 2 ) track_length = tl + step return endif c do i = ic , 1 ,-1 a = (1d0*ta(i)) ** 2 - a2 b = 2 * ( (1d0*ta(i)) ** 2 ) * zd(i) - b2 c = ( 1d0 * ta(i) * zd(i) ) ** 2 - c2 if ( a . ne . 0. ) then discrim = dsqrt ( b ** 2 - 4 * a * c ) z2 = ( - b + discrim ) / ( 2 * a ) z3 = ( - b - discrim ) / ( 2 * a ) z1 = amin1 ( z2 , z3 ) if ( z1 . lt . - zd(i) ) z1 = amax1 ( z2 , z3 ) else z1 = - c / b endif if ( (i . gt . 1) . and . (z1 . gt . zu(i-1)) ) then ic = i goto 20 endif enddo c ic = 1 20 ratio = abs ( z1 - z0 ) / r(ic-1) c c check if C-photon is moving out of the hut c if ( (ratio . lt . epsilon) . and . (nreflex . gt . 0) ) then track_length = 0. return endif c x1 = ax + bx * z1 y1 = ay + by * z1 r1 = sqrt ( x1 ** 2 + y1 ** 2 ) dx = x1 - x0 dy = y1 - y0 dz = z1 - z0 step = sqrt ( dx ** 2 + dy ** 2 + dz ** 2 ) tl = tl + step nx = ( - x1 / r1 ) * ca(ic) ny = ( - y1 / r1 ) * ca(ic) ni0 = nx * u0 + ny * v0 + nz(ic) * w0 u1 = 2 * ni0 * nx - u0 v1 = 2 * ni0 * ny - v0 w1 = 2 * ni0 * nz(ic) - w0 x0 = x1 y0 = y1 z0 = z1 u0 = u1 v0 = v1 w0 = w1 norm = dsqrt ( u0 ** 2 + v0 ** 2 + w0 ** 2 ) u0 = u0 / norm v0 = v0 / norm w0 = w0 / norm if ( dabs ( w0 ) . lt . 1d-5 ) then track_length = 0. return endif nreflex = nreflex + 1 goto 10 end c c ********************************************************************** SUBROUTINE hit_sphere(z_u,z_d) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Subroutine to calculate the meeting point of the straight path c of the C-photon with the photocathode of the (hemispherical) PM c ********************************************************************** common /sphere_PM/ r_sph,zu_sph,sph1,sph2 common /st_line/ a2,b2,c2 double precision a2 , b2 , c2 double precision a , b , c a = 1. + a2 b = sph1 + b2 c = sph2 + c2 disc = b ** 2 - 4 * a * c if ( disc . ge . 0. ) then z_u = ( - b + sqrt ( disc ) ) / ( 2 * a ) z_d = ( - b - sqrt ( disc ) ) / ( 2 * a ) else z_u = -1. z_d = -1. endif return end c c ********************************************************************** SUBROUTINE Get_NSLNoise c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 13 / Dec / 95 C c ********************************************************************** c common /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM common /histog2/ nbins_PM , TimeWindow_PM , time_shift_PM common /TrTimePM/ nTrTime , TrTime_PM , Dist_TrTime_PM dimension TrTime_PM(100) , Dist_TrTime_PM(100) parameter ( nbins = 2000 ) COMMON / anodePM / PMcurrent(nbins), PMcurrent_NSL(nbins) c call vzero( PMcurrent_NSL , nbins ) do 20 iNoise = 1 , nbins_PM + nTrTime - 1 jchan1 = iNoise - nTrTime + 1 jchan1 = max0( 1 , jchan1 ) jchan2 = iNoise jchan2 = min0( nbins_PM , jchan2 ) call poissn( Av_NSLNoise , n , ierror ) do 30 jchan = jchan1 , jchan2 jTr = nTrTime - iNoise + jchan if ( (jTr .lt. 1) .or. (jTr .gt. nTrTime) ) goto 30 PMcurrent_NSL(jchan) = PMcurrent_NSL(jchan) + * n * Gain_PM * current_factor * * Dist_TrTime_PM(jTr) 30 continue 20 continue return end c c ********************************************************************** SUBROUTINE Get_qADC1(id_hut) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Get charge (pC) in ADC of low-gain branch c ********************************************************************** COMMON /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM PARAMETER ( nbins = 2000 ) COMMON / anodePM / PMcurrent(nbins), PMcurrent_NSL(nbins) COMMON /histog2/ nbins_PM , TimeWindow_PM , time_shift_PM COMMON / Amplifiers / Gain_Amp1 , Gain_Amp1_2 , FWHM_CATS * , CATS_Dist( 200 ) , nbins_CATS COMMON /histog3/ nbins_Amp1_2 , TimeWindow_Amp1_2 PARAMETER ( nhuts = 2000 ) COMMON /ADC_TDC / q_ADC1( nhuts ), q_ADC1_2( nhuts ) * , TDC1_2( nhuts ) , Gate_qADC1_2 c ================================================= ibin1_PM = (nbins_PM - nbins_Amp1_2)/2 + 1 ibin2_PM = ibin1_PM + nbins_Amp1_2 - 1 c c ----------------------------------- c Get integrated charge into q_ADC1 c q_ADC1(id_hut) = 0.0 DO 20 ibin_PM = ibin1_PM , ibin2_PM ! a 600-ns window ! centered in the C-pulse q_ADC1(id_hut) = q_ADC1(id_hut) + * (PMcurrent(ibin_PM) - Av_NSLCurrent) 20 CONTINUE q_ADC1(id_hut) = q_ADC1(id_hut) * Gain_Amp1 * / 2 / 1e15 * 1e12 c RETURN END c c ********************************************************************** SUBROUTINE Get_current_Amp1_2 c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Gets the output current after the 2nd amplifier into histogram # 2000 c ********************************************************************** c COMMON /NSL/ Av_NSLNoise , Av_NSLCurrent , current_factor , * Gain_PM COMMON /histog2/ nbins_PM , TimeWindow_PM , time_shift_PM PARAMETER ( nbins = 2000 ) COMMON / anodePM / PMcurrent( nbins ), PMcurrent_NSL(nbins) COMMON / Amplifiers / Gain_Amp1 , Gain_Amp1_2 , FWHM_CATS * , CATS_Dist( 200 ) , nbins_CATS COMMON /histog3/ nbins_Amp1_2 , TimeWindow_Amp1_2 common /Amp1_2_current/ currentAmp1_2( nbins ) c ///////////////////////////////////////////////////////////// ibin1_PM = (nbins_PM - nbins_Amp1_2)/2 - nbins_CATS/2 ibin2_PM = nbins_PM - ibin1_PM c CALL vzero( currentAmp1_2 , nbins ) c DO 20 ibin_PM = ibin1_PM + 1, ibin2_PM + 1 jbin1_Amp = (ibin_PM - ibin1_PM) - nbins_CATS + 1 jbin1_Amp = max0( 1 , jbin1_Amp ) jbin2_Amp = (ibin_PM - ibin1_PM) jbin2_Amp = min0( nbins_Amp1_2 , jbin2_Amp ) DO 30 jbin_Amp = jbin1_Amp , jbin2_Amp jCATS = nbins_CATS - (ibin_PM - ibin1_PM) + jbin_Amp IF ( (jCATS .LT. 1) .OR. (jCATS .GT. nbins_CATS) ) * GOTO 30 currentAmp1_2(jbin_Amp) = currentAmp1_2(jbin_Amp) + * (PMcurrent(ibin_PM) - Av_NSLCurrent) * * Gain_Amp1_2 * CATS_Dist(jCATS) 30 CONTINUE 20 CONTINUE RETURN END c c ********************************************************************** SUBROUTINE Get_TDC_qADC1_2_DATA(id_hut) c ********************************************************************** c AUTHOR: Sergio Martinez (Madrid-HEGRA group) c LAST UPDATE: 10 / Mar / 95 c c Get Time t_0 of CFD trigger & charge (pC) in ADC of high-gain branch, c integrated between t_0-5ns and t_0-5ns + 30ns c ********************************************************************** COMMON /histog3/ nbins_Amp1_2 , TimeWindow_Amp1_2 PARAMETER ( nbins = 2000 ) common /Amp1_2_current/ currentAmp1_2( nbins ) COMMON / CFD / CFD_delay1 , CFD_delay2 , U_Threshold, Resist_CFD PARAMETER ( nhuts = 2000 ) COMMON /ADC_TDC / q_ADC1( nhuts ), q_ADC1_2( nhuts ) * , TDC1_2( nhuts ) , Gate_qADC1_2 c /////////////////////////////////////////////////////////////////////// data t0 / 5. / c ------------------------------ idel1 = 2 * CFD_delay1 idel2 = 2 * CFD_delay2 Current_THRES = U_Threshold / Resist_CFD * 1e3 ! U_threshold<0 c TDC1_2( id_hut ) = 1 000 000. q_ADC1_2( id_hut ) = 0. ibin1_Amp = max0( 0 , idel2 ) ibin_Amp_max = min0(nbins_Amp1_2 , nbins_Amp1_2 + idel1) c ---------------------------------------------------------------------* DO 20 ibin_Amp = ibin1_Amp + 1 , ibin_Amp_max IF ( (-0.5*currentAmp1_2( ibin_Amp ) ) .GT. * Current_THRES ) GOTO 20 c S_CFD = 3. * currentAmp1_2( ibin_Amp - idel2 ) - * currentAmp1_2( ibin_Amp - idel1 ) IF (S_CFD . lt . 0.) THEN GOTO 20 ELSE ibin0 = (ibin_Amp - idel2) - 2 * t0 CALL hix( id_hut+4000 , ibin0 , TDC1_2_temp ) TDC1_2( id_hut ) = TDC1_2_temp GOTO 30 ENDIF 20 CONTINUE 30 CONTINUE c c Get the integrated charge into q_ADC1_2 c IF ( TDC1_2( id_hut ) . ne . 1 000 000. ) THEN ibin1 = ibin0 ibin2 = ibin1 + 2 * Gate_qADC1_2 IF (ibin2 .GT. nbins_Amp1_2) ibin2 = nbins_Amp1_2 DO ibin = ibin1, ibin2 q_ADC1_2( id_hut ) = q_ADC1_2( id_hut ) * + currentAmp1_2(ibin) * 0.45 * /2 / 1e15 * 1e12 ENDDO ENDIF RETURN END