*********************************************************************** subroutine uhbook * * book histograms and ntuples for detector components * * last modified on: 05/12/2008 by I. Koenig *********************************************************************** implicit none #include "geantdef.h" #include "user.inc" #include "rztunits.inc" #include "kinetups.inc" #include "pmc.inc" #include "hcons.inc" integer i, jpart real T1, T2, frac, b, A2, A4 logical vdm, hexist data vdm /.true./ real E, mass, weight, angle, cost real vdmform, qedform, blast character*20 cdum1 real dum2, dum3, dum4, dum5, dum6 call hcdir('//PAWC',' ') ! make sure we are in PAWC if(.not.hexist(9995)) then call hbook1(9995,'[p]^0! Dalitz',2000,0.0,200.0,0.0) do i=1,2000 ! set up pi0 Dalitz e+e- mass distribution call hix(9995,i,mass) mass = mass + 0.05 if(vdm) then ! use VDM form factor weight = 1.0e5*vdmform(mass,134.9766,1) else ! use QED form factor weight = 1.0e5*qedform(mass,134.9766,1) end if call hf1(9995,mass,weight) end do end if if(.not.hexist(9996)) then call hbook1(9996,'[c] Dalitz',6000,0.0,600.0,0.0) do i=1,6000 ! set up eta Dalitz e+e- mass distribution call hix(9996,i,mass) mass = mass + 0.05 if(vdm) then weight = 1.0e5*vdmform(mass,547.75,1) ! PDG 2004 value else weight = 1.0e5*qedform(mass,547.75,1) end if call hf1(9996,mass,weight) end do end if if(.not.hexist(9997)) then call hbook1(9997,'[w] Dalitz',8000,0.0,800.0,0.0) do i=1,8000 ! set up omega Dalitz e+e- mass distribution call hix(9997,i,mass) mass = mass + 0.05 if(vdm) then weight = 1.0e5*vdmform(mass,782.59,2) ! PDG 2004 value else weight = 1.0e5*qedform(mass,782.59,2) end if call hf1(9997,mass,weight) end do end if if(pmc2(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc2(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9902)) then T1 = pmc2(3) T2 = pmc2(4) frac = pmc2(5) b = pmc2(6) call hbook1(9902,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9902,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9902,E,weight) end do end if if(.not.hexist(9912)) then A2 = pmc2(7) A4 = pmc2(8) call hbook1(9912,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9912,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9912,angle,weight) end do end if end if if(pmc3(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc3(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9903)) then T1 = pmc3(3) T2 = pmc3(4) frac = pmc3(5) b = pmc3(6) call hbook1(9903,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9903,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9903,E,weight) end do end if if(.not.hexist(9913)) then A2 = pmc3(7) A4 = pmc3(8) call hbook1(9913,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9913,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9913,angle,weight) end do end if end if if(pmc4(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc4(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9904)) then T1 = pmc4(3) T2 = pmc4(4) frac = pmc4(5) b = pmc4(6) call hbook1(9904,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9904,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9901,E,weight) end do end if if(.not.hexist(9914)) then A2 = pmc4(7) A4 = pmc4(8) call hbook1(9914,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9914,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9914,angle,weight) end do end if end if if(pmc5(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc5(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9905)) then T1 = pmc5(3) T2 = pmc5(4) frac = pmc5(5) b = pmc5(6) call hbook1(9905,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9905,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9905,E,weight) end do end if if(.not.hexist(9915)) then A2 = pmc5(7) A4 = pmc5(8) call hbook1(9915,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9915,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9915,angle,weight) end do end if end if if(pmc6(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc6(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9906)) then T1 = pmc6(3) T2 = pmc6(4) frac = pmc6(5) b = pmc6(6) call hbook1(9906,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9906,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9906,E,weight) end do end if if(.not.hexist(9916)) then A2 = pmc6(7) A4 = pmc6(8) call hbook1(9916,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9916,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9916,angle,weight) end do end if end if if(pmc7(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc7(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9907)) then T1 = pmc7(3) T2 = pmc7(4) frac = pmc7(5) b = pmc7(6) call hbook1(9907,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9901,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9907,E,weight) end do end if if(.not.hexist(9917)) then A2 = pmc7(7) A4 = pmc7(8) call hbook1(9917,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9917,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9917,angle,weight) end do end if end if if(pmc8(2).gt.0.) then ! prepare thermal source generator jpart = ifix(pmc8(2) + 0.01) ! get Geant particle and its mass call gfpart(jpart,cdum1,dum2,mass,dum3,dum4,dum5,dum6) mass = mass*1000. if(.not.hexist(9908)) then T1 = pmc8(3) T2 = pmc8(4) frac = pmc8(5) b = pmc8(6) call hbook1(9908,'Thermal source',1000,mass,mass+2000.,0.0) do i=1,1000 call hix(9908,i,E) if(T1.gt.0.) weight = frac*blast(E,T1,b,mass) if(T2.gt.1.) weight = weight+(1.-frac)*blast(E,T2,b,mass) weight = 1.0e5*weight call hf1(9908,E,weight) end do end if if(.not.hexist(9918)) then A2 = pmc8(7) A4 = pmc8(8) call hbook1(9918,'Thermal source dN/dTheta',180,0.,180.,0.) do i=1,181 ! set up polar distribution call hix(9918,i,angle) cost = cos(dtor*angle) weight = sin(dtor*angle)*(1.0 + A2*cost**2 + A4*cost**4) call hf1(9918,angle,weight) end do end if end if if (ikin.ge.1) then ! kinematics ntuple(s) if (iroot.eq.0) then c maxtrk = maxtrk1 call hcdir('//KINETUP',' ') call kinebook else i = 0 #ifdef WITHROOT call kinebookr #endif endif endif if (istart.ge.1) then ! START branch i = 0 #ifdef WITHROOT call startbookr #endif endif if(irich.ge.1) then ! RICH ntuple(s) if (iroot.eq.0) then call hcdir('//RICHTUP',' ') call richbook else i = 0 #ifdef WITHROOT call richbookr #endif endif endif if (imdc.ge.1) then ! book MDC ntuple(s) if (iroot.eq.0) then call hcdir('//MDCTUP',' ') call mdcbook else i = 0 #ifdef WITHROOT call mdcbookr #endif endif endif if(ishow.ge.1) then ! SHOWER detector ntuple(s) if (iroot.eq.0) then call hcdir('//SHOWTUP',' ') call showbook else i = 0 #ifdef WITHROOT call showbookr #endif endif endif if (itof.ge.1) then ! TOF wall ntuples(s) if (iroot.eq.0) then call hcdir('//TOFTUP',' ') call tofbook else i = 0 #ifdef WITHROOT call tofbookr #endif endif endif if (irpc.ge.1) then ! RPC branch i = 0 #ifdef WITHROOT call rpcbookr #endif endif if (iwall.ge.1) then ! Forward Wall branch i = 0 #ifdef WITHROOT call wallbookr #endif endif if (iemc.ge.1) then ! EMC branch i = 0 #ifdef WITHROOT call emcbookr #endif endif if (ists.ge.1) then ! Forward Dector STS branch i = 0 #ifdef WITHROOT call stsbookr #endif endif if (ifrpc.ge.1) then ! Forward Dector RPC branch i = 0 #ifdef WITHROOT call frpcbookr #endif endif if(ireft.gt.0) call refbook ! ntuple(s) for reference trajectories #ifdef WITHROOT if(iroot.eq.1 .and. itup.eq.1) call booktree(nsplit) ! set up ROOT tree #endif return *********************************************************************** entry uhbook_last ! cleanup after simulation is done call hcdir('//PAWC',' ') if(irich.ge.1) then if (iroot.eq.0) then call hcdir('//RICHTUP',' ') call richbooklast call hrend('RICHTUP') endif endif if (imdc.ge.1) then if (iroot.eq.0) then call hcdir('//MDCTUP',' ') call mdcbooklast call hrend('MDCTUP') endif endif if(ishow.ge.1) then if (iroot.eq.0) then call hcdir('//SHOWTUP',' ') call showbooklast call hrend('SHOWTUP') endif endif if (itof.ge.1)then if (iroot.eq.0) then call hcdir('//TOFTUP',' ') call tofbooklast call hrend('TOFTUP') endif endif if(irpc.ge.1) then endif if(iwall.ge.1) then endif if(ikin.ge.1) then if (iroot.eq.0) then call hcdir('//KINETUP',' ') call kinebooklast call hrend('KINETUP') endif endif #ifdef GEANT_REF write(6,401)ireft 401 format(' enter reflast at 401: ireft = ',i4) call reflast write(6,402)ireft 402 format(' exit reflast at 402: ireft = ',i4) #endif if(iroot.eq.0) then do i=1,nrztfiles ! close all open rzt-files close(ifileun(i)) enddo else i = 0 #ifdef WITHROOT if(itup.eq.1) call closetree ! or close ROOT tree/file #endif endif call hcdir('//PAWC',' ') return end c c Thermal source + blast according to Siemens and Rasmussen c c c real function blast(E,T,beta,mass) implicit none real E, T, beta, mass real alpha, gamma, p, p2 blast = 0.0 if(E.le.0.0 .or. T.le.0.0 .or. mass.le.0.0) return p2 = E*E - mass*mass if(p2.le.0.0) return p = sqrt(p2) if(beta.le.0.01) then ! no blast, just Boltzmann blast = E*p*exp(-E/T) return else gamma = 1./sqrt(1.-beta*beta) alpha = beta*gamma*p/T blast = (gamma+T/E)*sinh(alpha)/alpha - T/E*cosh(alpha) blast = blast * E*p*exp(-gamma*E/T) return end if end