program calorimeter c include 'clrm.inc' integer nwpawc integer lrecl integer istat parameter (ngbank=1000000,nwpawc=1000000) parameter (lrecl=1024) common/gcbank/q(ngbank) common/pawc/h(nwpawc) c This routine is needed to book job time (in seconds) since the default c value is 999 seconds. call TIMEST(10000.) c Initialises HBOOK and GEANT memory call gzebra(ngbank) c init_hbook call hlimit(-nwpawc) c Geant Initialisation call uginit c Start event processing call grun c call GRFILE (15,'examp.rz','NO') c End of run call uglast END c********************************************* c********************************************************* subroutine uginit c********************************************************* c Initialise GEANT include 'clrm.inc' parameter (lrecl=1024) common/withme/n_reg,n_lost,thresh call ginit c-- File for hbooking energy lost data --------------------- call hropen(1,'myfile','clrm.dst', + 'N',lrecl,istat) call hbookn(2000,'ntuple',6,' ',1024,chtags) call hbook1(100,'X shower profile',24,1.,25.,0.) call hbook1(200,'Y shower profile',24,1.,25.,0.) call hbook1(3000,'en. dep. per plane',60,1.,61.,0.) call hbook1(10,'cum. en. dep. per plane',60,1.,61.,0.) c do i=1,60 c do j=1,24 c call hbook1(i*100+j,' ',1024,0.,1024.,0.) c enddo c enddo c------------------------------------------------------------- c call GRFILE(20,'minical.rz','N') n_reg=0 n_lost=0 c write(*,*) 'Input threshold for triger (MeV)' c read *, thresh thresh=0.3 c read(*,*) numeropart,momentpart,nevent c read(*,*) momentpart momentpart=2.0 numeropart=9 nevent=1000 c Esplane - energy deposition per scint plane c Efplane - energy deposition per Fe plane c Ealplane - energy deposition per Al plane do i=1,60 Esplane(i)=0. Efplane(i)=0. Ealplane(i)=0. enddo c Estrip2 - energy deposition per strip for X and Y shower profile calc. do i=1,60 do j=1,24 Estrip2(i,j)=0. enddo enddo c cutele=0.00001 c cutgam=0.00001 call gzinit c call gdinit call gpart c Description of materials call gmate c call GROUT ('MATE',20,' ') c Geometry description call ugeom call gphysi c call gprint('volu',0) return end c*********************************************************************** subroutine ugeom c******************************************************** include 'clrm.inc' dimension volair(3) data volair /55.,55.,200./ dimension volst(3) data volst /50.,50.,1.27/ dimension volsc(3) data volsc /49.2,2.05,0.5/ dimension vols2(3) data vols2 /2.05,49.2,0.5/ dimension volal(3) data volal /49.2,49.2,0.05/ dimension AARY(2) dimension ZARY(2) dimension WMAT(2) C PURE POLYSTYRENE : C6H5CH=CH2 NPAR = -2 ! -n means numbers by atomic proportions (not weights) AARY(1) = 1.00794 AARY(2) = 12.011 ZARY(1) = 1.0 ZARY(2) = 6.0 WMAT(1) = 8. WMAT(2) = 8. DENS = 1.032 CALL GSMIXT(21,'POLY',AARY,ZARY,DENS,NPAR,WMAT) c Automatic calculation of the tracking parameters call gstmed (1,'IRON',10,0,0,0.,0.,STEMAX, + DEEMAX,EPSIL,STMIN,0,0) CALL GSTMED (2,'POLY',21,1,0,0.,0.,STEMAX, & DEEMAX,EPSIL,STMIN,0,0) c CALL GSTMED (3,'VACM',16,0,0,0.,0.,STEMAX, c + DEEMAX,EPSIL,STMIN,0,0) call gstmed (4,'ALUM',9,0,0,0.,0.,STEMAX, + DEEMAX,EPSIL,STMIN,0,0) call gstmed (5,'AIRY',15,0,0,0.,0.,STEMAX, + DEEMAX,EPSIL,STMIN,0,0) c iargo=1 call gsvolu('AIRY','BOX ',5,volair,3,iair) call gsvolu('IRON','BOX ',1,volst,3,iiron) call gsvolu('SCNT','BOX ',2,volsc,3,iscnt) call gsvolu('SNT2','BOX ',2,vols2,3,isnt2) call gsvolu('ALUM','BOX ',4,volal,3,ialu) c call gsvolu('ALU2','BOX ',4,volal2,3,ialu2) c call gsvolu('VACM','BOX ',3,volvac,3,ivac) c copy number key: xxxx are scintillator strips in the main calorimeter c 70xx are iron plates in main calorimeter c 80xx an d 81xx are aluminium plates in main cal c call gspos('VACM',9000,'AIRY',0.,0.,22.1,0,'ONLY') do i=1,60 call gspos('IRON',7000+i,'AIRY',0.,0.,184.73-(6.0*i),0, + 'ONLY') do j=1,24 if (mod(i,2).eq.0) then call gspos('SCNT',(i*100)+j,'AIRY',0.,51.25-(j*4.1), + 181.73-(6.0*i),0,'ONLY') else call gspos('SNT2',(i*100)+j,'AIRY',51.25-(j*4.1),0., + 181.73-(6.0*i),0,'ONLY') endif enddo call gspos('ALUM',8100+i,'AIRY',0.,0.,182.28-(6.0*i),0, + 'ONLY') call gspos('ALUM',8000+i,'AIRY',0.,0.,181.18-(6.0*i),0, + 'ONLY') enddo c call GROUT ('VOLU',20,' ') call gpvolu(0) call gptmed(0) call ggclos return end c******************************************************************* subroutine gukine c************************************************************************ include 'clrm.inc' dimension plab(3) c data plab/0.,0.,-2.0/ c plab(1)=0. c plab(2)=0. c plab(3)=-momentpart sigres=0.01*momentpart call norran(xres) p1=xres*sigres+momentpart c Angular distribution phi1=2.*pi*rndm(d) c phi1=7.*pi/4. c phi1=0. ctet1=-0.98-0.02*rndm(d) c ctet1=-0.978 c ctet1=-1. stet1=sqrt(1.-ctet1*ctet1) plab(1)=p1*stet1*cos(phi1) plab(2)=p1*stet1*sin(phi1) plab(3)=p1*ctet1 c nevent=3000 vert(1)=0. ! ..cm from 12strip vert(2)=0. ! ..cm from 12strip vert(3)=181. call gsvert (vert,0,0,0,0,nvtx) call gskine (plab,numeropart,nvtx,0,0,nt) return end c****************************************************************** subroutine gutrev c****************************************************************** include 'clrm.inc' do i=1,6 n_tuple(i)=0. enddo do i=1,60 do j=1,24 Estrip(i,j)=0. Emip(i,j)=0. enddo enddo call gtreve return end c************************************************************************* subroutine gutrak c**************************************************** include 'clrm.inc' call gtrack return end c***************************************************************** subroutine gustep c************************************************************ include 'clrm.inc' c write(*,*) numed,destep if(numed.eq.1) n_tuple(1)=n_tuple(1)+destep ! Iron if(numed.eq.2) n_tuple(2)=n_tuple(2)+destep ! Scintillator total if(numed.eq.4) n_tuple(3)=n_tuple(3)+destep ! Al c write(*,*) number(2) if (number(2).gt.100.and.number(2).lt.7000) then nplane = number(2)/100 ! scintillator plane number istrip=number(2)-(nplane*100) ! scin.strip ind at each plane c nstrip=(nplane-1)*24+nstrip_ind ! scint. strip number c write(*,*) nplane,nstrip Esplane(nplane) = Esplane(nplane)+destep Estrip(nplane,istrip) = Estrip(nplane,istrip)+destep Estrip2(nplane,istrip) = Estrip2(nplane,istrip)+destep elseif(number(2).gt.7000.and.number(2).lt.8000) then nplane=number(2)-7000 Efplane(nplane)=Efplane(nplane)+destep elseif(number(2).gt.8000.and.number(2).lt.8100) then nplane=number(2)-8000 Ealplane(nplane)=Ealplane(nplane)+destep elseif(number(2).gt.8100) then nplane=number(2)-8100 Ealplane(nplane)=Ealplane(nplane)+destep endif cc-- Including secondary particles ----------------- if(ngkine.ne.0) then do i=1,ngkine call GSKING(i) enddo endif cc--------------------------------------------------------- return end c********************************************************************** subroutine guout c******************************************************************** include 'clrm.inc' logical hicheck,hexist parameter (lrecl=1024) common/withme/n_reg,n_lost,thresh nat=0 c n_tuple(1) ! Fe in GeV n_tuple(2) = n_tuple(2)*1000./1.94 ! scint in mip's n_tuple(3) = n_tuple(3)*1000. ! Al in MeV c n_tuple(4) = vect(1) n_tuple(5) = vect(2) n_tuple(6) = vect(3) DO nplane=1,60 do istrip=1,24 Estrip(nplane,istrip)=Estrip(nplane,istrip)*1000. if(Estrip(nplane,istrip).ge.thresh) then nat=nat+1 indhist=nplane*100+istrip indhist2=nplane*10000+istrip hicheck = hexist(indhist) if(.NOT.hicheck) then c 12 bit ADC call hbook1(indhist,' ',4096,0.,4096.,0.) call hbook1(indhist2,' ',6000,0.,300.,0.) endif call hf1(indhist2,Estrip(nplane,istrip)/1.94,1.) CALL PEOUT(nplane,istrip) endif n_tuple(4)=n_tuple(4)+Emip(nplane,istrip) enddo ENDDO c call PEOUT(36) c write(*,*) number(2) c write(*,*) vect(3) c write(*,*) lvolum(2) call HFN(2000,n_tuple) if(nat.ge.1) then n_reg=n_reg+1 else n_lost=n_lost+1 endif return end c*********************************************************************** subroutine uglast c************************************************************** include 'clrm.inc' logical hicheck,hexist dimension Ecum(60) ! Cumulative energy dep per plane real Eprofile_x(24),Eprofile_y(24) parameter (lrecl=1024) common/withme/n_reg,n_lost,thresh ped_wid=3.5 ped_lev=33.4 DO nplane=1,60 do istrip=1,24 indhist=nplane*100+istrip hicheck = hexist(indhist) if(hicheck) then ihevent=hstati(indhist,3,' ',0) if(ihevent.lt.nevent) then do i=1,(nevent-ihevent) call norran(x) ped=x*ped_wid+ped_lev call hf1(indhist,ped,1.) enddo endif endif enddo ENDDO do iplane=1,60 Ecum(iplane)=Esplane(iplane)+Efplane(iplane)+ + Ealplane(iplane) Esplane(iplane)=Esplane(iplane)*1000./1.94 call hf1(3000,float(iplane),Esplane(iplane)/real(n_reg)) enddo do iplane=1,60 if(iplane.gt.1) Ecum(iplane)=Ecum(iplane)+Ecum(iplane-1) call hf1(10,float(iplane),Ecum(iplane)/real(n_reg)) enddo c ----------- X and Y shower profiles ---------------------- do istrip=1,24 do iplane=1,60 Estrip2(iplane,istrip)=Estrip2(iplane,istrip)*1000./1.94 ! in mip if(mod(iplane,2).ne.0) then Eprofile_x(istrip)=Eprofile_x(istrip)+ + Estrip2(iplane,istrip) else Eprofile_y(istrip)=Eprofile_y(istrip)+ + Estrip2(iplane,istrip) endif enddo call hf1(100,float(istrip),Eprofile_x(istrip)/real(n_reg)) call hf1(200,float(istrip),Eprofile_y(istrip)/real(n_reg)) c write(*,*) Eprofile_x(istrip),Eprofile_y(istrip) enddo c -------------------------------------------------------------- call hrout(0,icycle,' ') call hrend('myfile') close(1) c call GREND(15) c call GREND(20) call glast write(*,*), 'lost events ----', n_lost write(*,*), 'registrated events ----', n_reg EF=real(n_reg)/real(nevent) write(*,*), 'EFFICIENCY --- ', EF return END c*********************************************************** subroutine peout(nplane,istrip) include 'clrm.inc' INTEGER i, npe REAL x, ped_level, ped_width, chargeperadc, electroncharge REAL zeros, xmean, xsig, hstati, meanpe, diff, gainfactor REAL prob_missed REAL pedestal, extra, chain(12), se(12), amplifier, prob_noise REAL npetot, xmeancorr, xsigcorr, semission_meas, vd(13) REAl total, high_voltage, epsilon, dsebydv, eps,collection real MeV2pe real fraction ! fraction of last stage gain (space charge effect) INTEGER npe3,npe4,n2,n3,n4 INTEGER nzero, nentries, ierr, j, adc, npe2, n, m INTEGER ndynode, nd, npeini DATA ped_level/33.4/, ped_width/3.5/, chargeperadc/0.25e-12/ DATA electroncharge/1.6e-19/, high_voltage/941.0/ DATA ndynode/12/, epsilon/1.0/, fraction /0.5/ DATA dsebydv/0.04/, amplifier/3.0/ DATA prob_noise/0.0/, prob_missed/0.0/ c 1.94 MeV/mip and 2 p.e./mip data MeV2pe/1.030928/ data chain/3,2,2,1,1,1,1,1,1,1,2,5/ !hamamastu c data chain/3.55,2.2,2.2,1,1,1,1,1,1,1,1,1/ !oxford c data chain/3,1,1,1,1,1,1,1,1,1,1,1/ ! proba total = 0.0 DO i = 1, 12 total = total + chain(i) END DO DO i = 1, ndynode eps = epsilon vd(i) = chain(i)/total*high_voltage c se(i) = (0.259*vd(i)**0.7)*eps se(i) = (0.163*vd(i)**0.74)*eps c write(*,*)vd(i),se(i) END DO c DO istrip = 1,1440 meanpe = Estrip(nplane,istrip)*MeV2pe nzero = 0 nentries = 0 c... initial photoelectrons from photocathode: meanpe INCLUDES collection c -------------------------------------------------------------- call rnpssn(meanpe,npeini,ierr) call rnpssn(npeini*se(1),npe2,ierr) call rnpssn(npe2*se(2),npe3,ierr) call rnpssn(npe3*se(3),npe4,ierr) call rnpssn(npe4*se(4),npe5,ierr) call rnpssn(npe5*se(5),npe6,ierr) call rnpssn(npe6*se(6),npe7,ierr) call rnpssn(npe7*se(8),npe8,ierr) call rnpssn(npe8*se(9),npe9,ierr) call rnpssn(npe9*se(10),npe10,ierr) call rnpssn(npe10*se(11),npe11,ierr) call rnpssn(fraction*npe11*se(12),npe12,ierr) npetot = npe12 c --------------------------------------------------------- adc = (npetot*electroncharge/chargeperadc)*amplifier gainfactor = (electroncharge/chargeperadc)*amplifier call norran(x) extra = x*ped_width pedestal = ped_level+extra Emip(nplane,istrip) = 0.5*real(adc)/17.137823 adc = (adc+pedestal) indhist=nplane*100+istrip call hf1(indhist,float(adc),1.0) c call hfill(istrip+500,float(adc)-ped_level,0.0,1.0) c ENDDO RETURN END c************************************************************