subroutine digi_hodo2 c_begin_doc c Documentation for subroutine digi_ec1 c $Id: digi_hodo2.f,v 1.1.1.1 2007/12/04 20:12:15 jones Exp $ c Purpose: c -------- c performs ec1 digitizations c Author: c ------- c Elliott Wolin, College of William and Mary, 13-apr-95 c Major revisions: C Created: 23-Dec-1993 S.Boiarinov C Modified: 22-Mar-1995 H.Avakian C Purpose and methods : Digitization of the LAC CC C HITS: numsv(1) - number of layer (1-33) C hitsh(1) - X/Y-coordinate in strip system (cm) C hitsh(2) - Z-coordinate in strip system (cm) C hitsh(3) - lenght of light in strip (cm) C (from negative end) C hitsh(4) - energy lost (MeV) C hitsh(5) - TOF (nsec) C C DIGI: numvs(1) - number of sector (1-2) C numvs(2) - number of block (1-4) C numvs(3) - number of strip (1-24 or 1-40) C ldigi(1) - ADC- (count) C ldigi(2) - TDC- (count) C ldigi(3) - ADC+ (count) C ldigi(4) - TDC+ (count) c ---------------- c_end_doc c include files: #include #include #include include 'hodo2geom.inc' c#include "ffpar_ec1.inc" c -------------- c_end_inc c local variables: c ---------------- INTEGER lattlen(256) REAL laclumi(256,2) COMMON /LACATTN/ lattlen,laclumi integer i,j,k,l,nbeg,nend,kod,nsub,nabs integer nset,iset,js,ndet,idet,idig integer numvs(3),numsv(3,10000),itrack(10000),nhit integer ldigi(4),itemp,ndigi(2) real*4 hitsh(5,10000),hitshs(7,10000),lstrip,attleng,temp(2) real*4 Nel_to_charge,sigma_x5,sigma_x7,sigma5,sigma7 real*4 ggauss, pe_mevl,pe_mevr,sigma_pel,sigma_per,lstack parameter (Nel_to_charge=0.605) ! Gain*el_charge (G=4.5E+6)[PC/electron] C External references: external ggauss cc executable code: cc ---------------- cc if (jset.le.0) go to 999 cc nset = iq(jset-1) cc if (nset.le.0) go to 999 cc call glook('LUC2',iq(jset+1),nset,iset) cc if (iset.eq.0) go to 999 cc js = lq(jset-iset) cc ndet = iq(js-1) cc if (ndet.le.0) go to 999 cc call glook('HODO',iq(js+1),ndet,idet) cc if (idet.eq.0) go to 999 cC cC-- retrieve hits: cC c numvs(1)=0 c numvs(2)=0 c numvs(3)=0 c call GFHITS('LUC2','HODO',1,5,10000,0,numvs, c & itrack,numsv,hitsh,nhit) cC 10000 maximum number of hits numsv(3,10000),hitsh(7,10000) cC 0 all tracks are taken cC hitsh(7,10000) contain NHIT hits cC c if(nhit.le.0) goto 999 cC cC attenuation, Nlayer-->Nblock cC c write(*,*)hitsh c do i=1,nhit cC-- copy hitsh --> hitshs ( split energy between 2 ends ): c hitshs(1,i)=hitsh(1,i) c hitshs(2,i)=hitsh(2,i) c hitshs(3,i)=hitsh(3,i) c hitshs(4,i)=hitsh(4,i)/2. c hitshs(5,i)=hitsh(5,i) c hitshs(6,i)=hitsh(4,i)/2. c hitshs(7,i)=hitsh(5,i) cC-- index for attenuation array: c j=numsv(1,i) c k=numsv(2,i) c l=numsv(3,i) cC-- strip length: c nsub=(k-1)-((k-1)/2)*2 c if(nsub.eq.0) then c lstrip=(HS1+DS*(k-1))*2. ! short strip c else c lstrip=(HB1+DB*k)*2. ! long strips c endif cC-- layer number ---> block number: c itemp=numsv(2,i) c numsv(2,i)=numsv(2,i)-((numsv(2,i)-1)/2)*2 cC c if(itemp.gt.17) numsv(2,i)=numsv(2,i)+2 cc the block number 1-4 c k=numsv(2,i) cC-- atten. correction of amplitude (take 1/2 energy for each PM): c if(iatten.ne.0) then c if(iatten.lt.0) then c nabs=128*(j-1)+8*mod(k+1,2)+32*(k-1)+l c if((nabs.ge.1).and.(nabs.le.256)) then cC-- get the corresponding attenuation length c attleng=float(lattlen(nabs)) c else cc generate some warning message c print *,'nnnnnnnnnnnnaaaaaaaaaaaaaaabbbbbbbbbbbssss',nabs c attleng=-float(iatten) c endif c elseif(iatten.gt.0) then c attleng=float(iatten) c endif c hitshs(4,i)=hitshs(4,i)*(exp(-hitshs(3,i)/attleng) c 1 +reflect*exp((-2.0*lstrip+hitshs(3,i))/attleng)) c hitshs(6,i)=hitshs(6,i)*(exp(-(lstrip-hitshs(3,i))/ c 1 attleng)+reflect*exp((-lstrip-hitshs(3,i))/attleng)) c endif ! non 0 attlengths cccc print *,'first: ijk,nabs, attlen, lstrip',i,j,k,nabs,attleng,lstrip c enddo cC cC find the first cC c if(nhit.gt.0) then c nbeg=1 c nend=nhit c i=nbeg c kod=numsv(1,1)*10000+numsv(2,1)*100+numsv(3,1) c do while(nbeg.LT.nend) c i=i+1 c k=numsv(1,i)*10000+numsv(2,i)*100+numsv(3,i) c do while(k.eq.kod.and.i.le.nend) c hitshs(4,nbeg)=hitshs(4,nbeg)+hitshs(4,i) c hitshs(6,nbeg)=hitshs(6,nbeg)+hitshs(6,i) c IF(hitshs(5,i).LT.hitshs(5,nbeg)) hitshs(5,nbeg)=hitshs(5,i) c if(hitshs(7,i).lt.hitshs(7,nbeg)) hitshs(7,nbeg)=hitshs(7,i) c numsv(1,i)=numsv(1,nend) c numsv(2,i)=numsv(2,nend) c numsv(3,i)=numsv(3,nend) c hitshs(3,i)=hitshs(3,nend) c hitshs(4,i)=hitshs(4,nend) c hitshs(5,i)=hitshs(5,nend) c hitshs(6,i)=hitshs(6,nend) c hitshs(7,i)=hitshs(7,nend) c k=numsv(1,i)*10000+numsv(2,i)*100+numsv(3,i) c nend=nend-1 c enddo c if(i.eq.nend.or.i.eq.nend+1) then c nbeg=nbeg+1 c kod=numsv(1,nbeg)*10000+numsv(2,nbeg)*100+numsv(3,nbeg) c i=nbeg c endif c enddo c endif c do i=1,nend cC-- tof: cC cC Compute sigma for tof distribution cC cC-- index for strip length: c j=numsv(1,i) c k=numsv(2,i) c l=numsv(3,i) cC-- strip length: cc nsub=(k-1)-((k-1)/2)*2 c if(k.eq.1) then c lstack=(HS1+(HS2-HS1)/2.)*2. ! short stack inner c else if(k.eq.2) then c lstack=(HB1+(HB2-HB1)/2.)*2. ! long stack inner c else if(k.eq.3) then c lstack=(HS2)*2. ! short stack outer c else if(k.eq.4) then c lstack=(HB2)*2. ! long stack outer c endif cc--------------- c if(iatten.lt.0) then c nabs=128*(j-1)+8*mod(k+1,2)+32*(k-1)+l c if((nabs.ge.1).and.(nabs.le.256)) then cC-- get the corresponding N p.e. c pe_mevl=laclumi(nabs,1) c pe_mevr=laclumi(nabs,2) c else cc generate some warning message c print *,'nnnnnnnnnnnnaaaaaaaaaaaaaaabbbbbbbbbbbssss',nabs c pe_mevl=photoel c pe_mevr=photoel c endif c elseif(iatten.gt.0) then c pe_mevl=photoel c pe_mevr=photoel c endif cc c sigma_pel = min(1., c & 0.14+7./(4.*hitshs(4,i)*pe_mevl*Nel_to_charge)) ! nsec c sigma_per = min(1., c & 0.14+7./(4.*hitshs(6,i)*pe_mevr*Nel_to_charge)) ! nsec c sigma_x5 = 0.53*exp(hitshs(3,i)/590.-1.) c sigma_x7 = 0.53*exp((lstack-hitshs(3,i))/590.-1.) c sigma5 = sqrt(sigma_pel*sigma_pel + sigma_x5*sigma_x5) c sigma7 = sqrt(sigma_per*sigma_per + sigma_x7*sigma_x7) c hitshs(5,i)=hitshs(5,i)+delay*hitshs(3,i)/100. c hitshs(7,i)=hitshs(7,i)+delay*(lstack-hitshs(3,i))/100. c hitshs(5,i)=ggauss(hitshs(5,i),sigma5) c hitshs(7,i)=ggauss(hitshs(7,i),sigma7) c cC cC Poisson, store DIGI cC cC-- energy : MeV --> N p.m. c hitshs(4,i)=hitshs(4,i)*pe_mevl c hitshs(6,i)=hitshs(6,i)*pe_mevr cC-- Poisson c temp(1)=ggauss(0.,1.) c temp(2)=ggauss(0.,1.) c temp(1)=max(0.,hitshs(4,i)+temp(1)*sqrt(hitshs(4,i))) c temp(2)=max(0.,hitshs(6,i)+temp(2)*sqrt(hitshs(6,i))) cC-- N p.e. --> ADC counts c ndigi(1)=temp(1)*Nel_to_charge*charge_to_adc c ndigi(2)=temp(2)*Nel_to_charge*charge_to_adc cC cC store digi cC c if(ndigi(1).ge.ithreshold) then c ldigi(1)=ndigi(1) c else c ldigi(1)=0.0 c endif c if(ndigi(2).ge.ithreshold) then c ldigi(3)=ndigi(2) c else c ldigi(3)=0.0 c endif c if(ldigi(1).gt.0.0.or.ldigi(3).gt.0.0) then cc-- tof: c ldigi(2)=hitshs(5,i)*20. ! TDC counts c ldigi(4)=hitshs(7,i)*20. cC-- sector number c numvs(1)=numsv(1,i) cC-- block number c numvs(2)=numsv(2,i) cC-- strip number c numvs(3)=numsv(3,i) c call GSDIGI(iset,idet,itrack,1,numvs,ldigi,idig) c endif c enddo C********************************************************************* 1200 format ('?ec1digi error : digitization for LAEC ' 1 ,'no. ',i3,' in sector ',i2,' could not be stored.') C 999 continue return end c---------------------------------------------------------------------------------