*======================================================================= subroutine s_cal_calib(mode) *======================================================================= c SOS calorimeter calibration with electrons. c c Input paramater mode = 0 means collect data for calibration, c otherwise calibrate. * implicit none * integer mode include 'sos_data_structures.cmn' include 'sos_calorimeter.cmn' include 'sos_tracking.cmn' include 'gen_run_info.cmn' include 'gen_event_info.cmn' * integer ihit integer nblk real adc_pos,adc_neg integer nct_hit_blk(66),ipmt !sos logical write_out c common/scal_calib/nct_hit_blk,ncall,spare_id integer ncall data ncall/0/ real thr_lo,thr_hi !thresholds on sammed raw calorimeter signal. c integer spare_id logical ABORT character*80 err c ncall=ncall+1 if (ncall .eq. 1) then call g_IO_control(spare_id,'ANY',ABORT,err) !get IO channel open(spare_id,file='s_cal_calib.raw_data') do ipmt=1,66 nct_hit_blk(ipmt)=0 enddo endif if(mode.eq.0) then !collect data. c Choose clean single electron tracks within SOS momentum acceptance. if( (sntracks_fp.eq.1).and. & (snclusters_cal.eq.1).and. & (sntracks_cal.eq.1).and. & (scer_npe_sum.gt.4).and. & (abs(sdelta_tar(1)).lt.20.).and. & (abs(sbeta(1)-1.).lt.0.1).and. & spare_id .ne. 0 ) then *** & (sbeta_chisq(1).ge.0.).and.(sbeta_chisq(1).lt.1.) ) then c write_out = .false. do ihit=1,scal_num_hits nblk=(scal_cols(ihit)-1)*smax_cal_rows+scal_rows(ihit) nct_hit_blk(nblk) = nct_hit_blk(nblk) + 1 if (nct_hit_blk(nblk) .lt. 4000) write_out = .true. enddo c if (write_out) then c write(spare_id,'(i2,1x,f7.4,2(1x,f5.1,1x,f9.6))') & scal_num_hits,sp_tar(1), & strack_xc(1),sxp_fp(1),strack_yc(1),syp_fp(1) do ihit=1,scal_num_hits if(scal_cols(ihit).le.scal_num_neg_columns) then adc_neg=scal_adcs_neg(ihit) else adc_neg=0. end if adc_pos=scal_adcs_pos(ihit) nblk=(scal_cols(ihit)-1)*smax_cal_rows+scal_rows(ihit) write(spare_id,'(2(f9.3,1x),i2)'), & adc_pos,adc_neg,nblk end do endif ! if write_out c end if !electron in acceptance else !mode<>0, calibrate. close(spare_id) print*,'==========================================================' print*,'Calibrating SOS Calorimeter at event #',gen_event_id_number call scal_raw_thr(spare_id,thr_lo,thr_hi) print*,'lo & hi thresholds:', thr_lo,thr_hi call scal_clb_det(spare_id,gen_run_number,thr_lo,thr_hi) print*,'==========================================================' end if !mode=0 end *======================================================================= subroutine scal_raw_thr(lun,thr_lo,thr_hi) implicit none integer lun real thr_lo,thr_hi c Get thresholds around electron peak in summed raw calorimeter signal. integer*4 num_negs parameter (num_negs=22) !sos integer*4 nhit real*4 adc_pos,adc_neg integer*4 nh integer*4 nb real*8 eb c integer*4 nrow parameter (nrow=11) !sos real*4 zbl parameter (zbl=10.) real*4 x,xp,y,yp real*4 xh,yh integer*4 nc real*4 sig,avr,t real*4 qdc integer nev real s_correct_cal_neg, s_correct_cal_pos, s_correct_cal * * Get thresholds on total_signal/p_tar. * open(lun,file='s_cal_calib.raw_data',err=989) avr=0. sig=0. nev=0 do while(.true.) read(lun,*,end=3) nhit,eb,x,xp,y,yp qdc=0. do nh=1,nhit read(lun,*,end=3) adc_pos,adc_neg,nb nc=(nb-1)/nrow+1 xh=x+xp*(nc-0.5)*zbl yh=y+yp*(nc-0.5)*zbl if(nb.le.num_negs) then qdc=qdc+adc_pos*s_correct_cal_pos(xh,yh)*0.5 qdc=qdc+adc_neg*s_correct_cal_neg(xh,yh)*0.5 else qdc=qdc+adc_pos*s_correct_cal(xh,yh) end if enddo eb=eb*1000. t=qdc/eb c write(lun,*) t c write(lun,*) t,nhit,eb,x,xp,y,yp,nev avr=avr+t sig=sig+t*t nev=nev+1 c print*,eb,qdc,nev end do 3 close(lun) c print*,avr,sig,nev avr=avr/nev sig=sqrt(sig/nev-avr*avr) thr_lo=avr-3.*sig thr_hi=avr+3.*sig c write(*,*) 'thr_lo=',thr_lo,' thr_hi=',thr_hi return 989 write(*,*) ' error opening file s_cal_calib.raw_data, channel',lun, * ' in scal_raw_thr.f' c end *======================================================================= subroutine scal_clb_det(lun,nrun,thr_lo,thr_hi) implicit none c include 'sos_data_structures.cmn' include 'sos_calorimeter.cmn' c integer lun,nrun real thr_lo,thr_hi integer npmts parameter (npmts=66) !sos integer npmts2 parameter (npmts2=npmts*npmts) integer nrow parameter (nrow=11) !sos real*8 q0(npmts) real*8 qm(npmts,npmts) real*8 qe(npmts) real*8 q(npmts) real*8 eb real*8 e0 real*8 ac(npmts) real*8 au(npmts) c real*8 t real*8 s integer nev logical*1 eod integer i,j integer nf(npmts) integer minf parameter (minf=200) ! minimum number to hit pmt before including pmt in calib integer nums(npmts) integer numsel real*8 q0s(npmts) real*8 qes(npmts) integer nsi,nsj real*8 acs(npmts) real*8 aus(npmts) real*8 aux(npmts2) integer jp integer spare_id logical ABORT character*80 err character*40 fn real xh,yh open(lun,file='s_cal_calib.raw_data') do i=1,npmts q0(i)=0. qe(i)=0. do j=1,npmts qm(i,j)=0. end do au(i)=0. ac(i)=0. nf(i)=0 end do e0=0. c nev=0 eod=.false. do while(.not.eod) call s_get_data(lun,eb,q,xh,yh,eod,thr_lo,thr_hi) if(.not.eod) then do i=1,npmts if(q(i).gt.0.) then q0(i)=q0(i)+q(i) qe(i)=qe(i)+eb*q(i) do j=1,npmts qm(i,j)=qm(i,j)+q(i)*q(j) end do nf(i)=nf(i)+1 end if end do e0=e0+eb nev=nev+1 c if(nev/1000*1000.eq.nev) write(*,'(e10.3,i7)') e0,nev end if end do close(lun) do i=1,npmts q0(i)=q0(i)/nev qe(i)=qe(i)/nev do j=1,npmts qm(i,j)=qm(i,j)/nev end do end do e0=e0/nev numsel=0 do i=1,npmts if(nf(i).ge.minf) then numsel=numsel+1 nums(numsel)=i c print*,nums(numsel),numsel,nf(i) else write(*,*) ' PMT ',i,' only ',nf(i),' events. Will not to be calibrated. Gain is set to 0.' end if end do c print*,'numsel =',numsel write(*,'(''Number of events for each PMT for calib for run '',i7,'', '', 1 i6,'' events processed'')') nrun,nev write(*,*) ' PMT with less than', minf,' events are not included in calibration.' write(*,*) write(*,11) 'scal_pos_gain_cor=',(nf(i),i= 1, nrow) write(*,11) ' ',(nf(i),i= nrow+1,2*nrow) write(*,11) ' ',(nf(i),i=2*nrow+1,3*nrow) write(*,11) ' ',(nf(i),i=3*nrow+1,4*nrow) write(*,11) 'scal_neg_gain_cor=',(nf(i),i=4*nrow+1,5*nrow) write(*,11) ' ',(nf(i),i=5*nrow+1,6*nrow) write(*,11) ' ',(0., i=6*nrow+1,7*nrow) write(*,11) ' ',(0., i=7*nrow+1,8*nrow) c do i=1,numsel nsi=nums(i) q0s(i)=q0(nsi) qes(i)=qe(nsi) do j=1,numsel nsj=nums(j) jp=j+(i-1)*numsel aux(jp)=qm(nsj,nsi) c write(65,'(e12.5)') aux(jp) end do end do call calib(e0,q0s,qes,aux,numsel,numsel*numsel,aus,acs) do i=1,numsel nsi=nums(i) au(nsi)=aus(i) ac(nsi)=acs(i) end do c write(*,'(2e10.3,i5)') (ac(i),au(i),i,i=1,npmts) write(fn,'(a17,i5.5)') 'PARAM/scal.param.',nrun call g_IO_control(spare_id,'ANY',ABORT,err) !get IO channel open(spare_id,file=fn) write(spare_id,'(''; Calibration constants for run '',i7,'', '', 1 i6,'' events processed'')') nrun,nev write(spare_id,*) write(spare_id,10) 'scal_pos_gain_cor=',(ac(i)*1.D+3,i= 1, nrow) write(spare_id,10) ' ',(ac(i)*1.D+3,i= nrow+1,2*nrow) write(spare_id,10) ' ',(ac(i)*1.D+3,i=2*nrow+1,3*nrow) write(spare_id,10) ' ',(ac(i)*1.D+3,i=3*nrow+1,4*nrow) write(spare_id,10) 'scal_neg_gain_cor=',(ac(i)*1.D+3,i=4*nrow+1,5*nrow) write(spare_id,10) ' ',(ac(i)*1.D+3,i=5*nrow+1,6*nrow) write(spare_id,10) ' ',(0., i=6*nrow+1,7*nrow) write(spare_id,10) ' ',(0., i=7*nrow+1,8*nrow) * write(spare_id,10) 'scal_pos_gain_cor=',(ac(i)*2.D+3,i= 1, nrow) * write(spare_id,10) ' ',(ac(i)*2.D+3,i= nrow+1,2*nrow) * write(spare_id,10) ' ',(ac(i)*1.D+3,i=2*nrow+1,3*nrow) * write(spare_id,10) ' ',(ac(i)*1.D+3,i=3*nrow+1,4*nrow) * write(spare_id,10) 'scal_neg_gain_cor=',(ac(i)*2.D+3,i=4*nrow+1,5*nrow) * write(spare_id,10) ' ',(ac(i)*2.D+3,i=5*nrow+1,6*nrow) * write(spare_id,10) ' ',(0., i=6*nrow+1,7*nrow) * write(spare_id,10) ' ',(0., i=7*nrow+1,8*nrow) close(spare_id) call G_IO_control(spare_ID,'FREE',ABORT,err) !free up IO channel write(*,*) write(*,'(''Calibration constants for run '',i7,'', '', 1 i6,'' events processed'')') nrun,nev write(*,*) write(*,*) ' constants written to ',fn write(*,*) write(*,10) 'scal_pos_gain_cor=',(ac(i)*1.D+3,i= 1, nrow) write(*,10) ' ',(ac(i)*1.D+3,i= nrow+1,2*nrow) write(*,10) ' ',(ac(i)*1.D+3,i=2*nrow+1,3*nrow) write(*,10) ' ',(ac(i)*1.D+3,i=3*nrow+1,4*nrow) write(*,10) 'scal_neg_gain_cor=',(ac(i)*1.D+3,i=4*nrow+1,5*nrow) write(*,10) ' ',(ac(i)*1.D+3,i=5*nrow+1,6*nrow) write(*,10) ' ',(0., i=6*nrow+1,7*nrow) write(*,10) ' ',(0., i=7*nrow+1,8*nrow) * write(*,10) 'scal_pos_gain_cor=',(ac(i)*2.D+3,i= 1, nrow) * write(*,10) ' ',(ac(i)*2.D+3,i= nrow+1,2*nrow) * write(*,10) ' ',(ac(i)*1.D+3,i=2*nrow+1,3*nrow) * write(*,10) ' ',(ac(i)*1.D+3,i=3*nrow+1,4*nrow) * write(*,10) 'scal_neg_gain_cor=',(ac(i)*2.D+3,i=4*nrow+1,5*nrow) * write(*,10) ' ',(ac(i)*2.D+3,i=5*nrow+1,6*nrow) * write(*,10) ' ',(0., i=6*nrow+1,7*nrow) * write(*,10) ' ',(0., i=7*nrow+1,8*nrow) 10 format(a18,13(f6.3,',')) 11 format(a18,13(i5,',')) open(lun,file='s_cal_calib.raw_data') call g_IO_control(spare_id,'ANY',ABORT,err) !get IO channel open(spare_id,file='s_cal_calib.cal_data') write(*,*) 'In sos shower cal creating file s_cal_calib.cal_data, ', * 'channel ',spare_id nev=0 eod=.false. do while(.not.eod) call s_get_data(lun,eb,q,xh,yh,eod,0.,1.E+8) if(.not.eod) then s=0. * t=0. do i=1,npmts s=s+q(i)*ac(i) * t=t+q(i)*au(i) end do write(spare_id,*) s,eb,xh,yh end if end do close(lun) close(spare_id) call G_IO_control(spare_ID,'FREE',ABORT,err) !free up IO channel end *======================================================================= c c... Done already in h_cal_calib.f. c c subroutine calib(e0,q0,qe,aux,npmts,npmts2,au,ac) c implicit none c integer npmts,npmts2 c real*8 e0 c real*8 q0(npmts) c real*8 qe(npmts) c real*8 aux(npmts2) c real*8 ac(npmts) c real*8 au(npmts) c real*8 qm(npmts,npmts) c real*8 t c real*8 s c integer ifail c integer i,j c integer jp c c do i=1,npmts c do j=1,npmts c jp=j+(i-1)*npmts c qm(j,i)=aux(jp) cD write(66,'(e12.5)') qm(j,i) c end do c end do c c print*,'Calib: npmts =',npmts c print*,' ' c c print*,'Inversing the Matrix...' c call smxinv(qm,npmts,ifail) c if(ifail.ne.0) then c stop '*** Singular Matrix ***' c else c print*,' ...done.' c end if c c do i=1,npmts c au(i)=0. c do j=1,npmts c au(i)=au(i)+qm(i,j)*qe(j) c end do c end do c c s=0. c do i=1,npmts c t=0. c do j=1,npmts c t=t+qm(i,j)*q0(j) c end do c s=s+q0(i)*t c end do c c t=0. c do i=1,npmts c t=t+au(i)*q0(i) c end do c s=(e0-t)/s c c do i=1,npmts c t=0. c do j=1,npmts c t=t+qm(i,j)*q0(j) c end do c ac(i)=s*t+au(i) c end do c c end *----------------------------------------------------------------------- subroutine s_get_data(lun,eb,q,xh,yh,eod,thr_lo,thr_hi) implicit none c integer lun real*8 eb integer*4 num_blocks,num_negs,num_pmts parameter (num_blocks=44,num_negs=22,num_pmts=66) !sos. real*8 q(num_pmts) logical*1 eod integer*4 nhit real*4 adc_pos,adc_neg integer*4 nh integer*4 nb c integer*4 nrow parameter (nrow=11) !sos real*4 zbl parameter (zbl=10.) real*4 x,xp,y,yp real*4 xh,yh integer*4 nc real*4 s_correct_cal real*4 s_correct_cal_pos,s_correct_cal_neg real*4 thr_lo,thr_hi logical*1 good_ev real*4 qnet good_ev=.false. do while(.not.good_ev) eb=0.d0 do nb=1,num_pmts q(nb)=0.d0 end do qnet=0. eod=.true. read(lun,*,end=5) nhit,eb,x,xp,y,yp do nh=1,nhit read(lun,*,end=5) adc_pos,adc_neg,nb nc=(nb-1)/nrow+1 xh=x+xp*(nc-0.5)*zbl yh=y+yp*(nc-0.5)*zbl if(nb.le.num_negs) then q(nb)=adc_pos*s_correct_cal_pos(xh,yh) q(num_blocks+nb)=adc_neg*s_correct_cal_neg(xh,yh) qnet=qnet+0.5*(q(nb)+q(num_blocks+nb)) else q(nb)=adc_pos*s_correct_cal(xh,yh) qnet=qnet+q(nb) end if enddo eod=.false. qnet=qnet/(eb*1000.) good_ev=(qnet.gt.thr_lo).and.(qnet.lt.thr_hi) c write(99,*) qnet end do !.not.good_ev 5 continue end *----------------------------------------------------------------------- c c... Done already in h_cal_calib.f. c c SUBROUTINE SMXINV (A,NDIM,IFAIL) C C CERN PROGLIB# F107 SMXINV .VERSION KERNFOR 1.0 720503 C ORIG. 03/05/72 CL C c REAL*8 A(*),RI(100) c INTEGER*4 INDEX(100) C c DATA TOL / 1.D-14/ C c IFAIL=0 c N=NDIM c NP1=N+1 c DO 10 I=1,N c 10 INDEX(I)=1 C c DO 80 I=1,N C C-- FIND PIVOT c PIVOT=0.0D0 c JJ=1 c DO 20 J=1,N c IF (INDEX(J).EQ.0) GO TO 19 c ELM=DABS (A(JJ)) c IF (ELM.LE.PIVOT) GO TO 19 c PIVOT=ELM c K=J c KK=JJ c 19 JJ=JJ+NP1 c 20 CONTINUE c IF (PIVOT/DABS(A(1)).LT.TOL) GO TO 100 c INDEX(K)=0 c PIVOT=-A(KK) C C-- ELIMINATION c KJ=K c NP=N C c DO 70 J=1,N c IF (J-K) 34,30,34 C c 30 A(KJ)=1.0D0/PIVOT c RI(J)=0.0D0 c NP=1 c GO TO 70 C c 34 ELM=-A(KJ) c 40 RI(J)=ELM/PIVOT c IF (ELM.EQ.0.0D0) GO TO 50 C c JL=J c DO 45 L=1,J c A(JL)=A(JL)+ELM*RI(L) c 45 JL=JL+N C c 50 A(KJ)=RI(J) C c 70 KJ=KJ+NP C c 80 CONTINUE C C-- CHANGE THE SIGN AND PROVISIONAL FILL-UP c IJ0=1 c JI0=1 c DO 95 I=1,N c IJ=IJ0 c JI=JI0 C c DO 90 J=1,I c A(IJ)=-A(IJ) c A(JI)=A(IJ) c IJ=IJ+N c JI=JI+1 c 90 CONTINUE C c IJ0=IJ0+1 c JI0=JI0+N c 95 CONTINUE c RETURN C C-- FAILURE RETURN c 100 IFAIL=1 c RETURN c END *=======================================================================