Return to results_write.f CVS log | Up to [HallC] / simc_semi |
File: [HallC] / simc_semi / results_write.f
(download)
Revision: 1.4, Thu Jan 31 19:01:00 2008 UTC (16 years, 7 months ago) by jones Branch: MAIN CVS Tags: HEAD Changes since 1.3: +2 -2 lines modified for calorimeter, mainly add electron = 7 (HMS side) and 8 (SOS side) |
subroutine results_ntu_write(main,vertex,orig,recon,success) implicit none include 'radc.inc' include 'hbook.inc' include 'simulate.inc' real*4 ntu(80) record /event_main/ main record /event/ vertex, orig, recon !local (e,e'p) calculations: real*8 poftheta !p as calculated from theta, assuming elastic. real*8 corrsing !'corrected singles' for elastic ! real*8 mm,mm2 !missing mass (assuming struck nucleon) ! real*8 mmA,mmA2 !missing mass (assuming struck nucleus) real*8 Pm_Heepx,Pm_Heepy,Pm_Heepz !Pm components for Heepcheck. !local (e,e'pi/K) calculations: ! real*8 t !t real*8 dummy logical success if (debug(2)) write(6,*)'r_ntu_write: entering' !If event did not make it thru spectrometers, return. Later, we will want to ! add option to write event even if it failed when doing_phsp. if (.not.success) return * if (phot1.gt.lim1) write(6,*) 'phot1,lim1=',phot1,lim1 * if (phot2.gt.lim2) then * write(6,*) 'phot2,lim2=',phot2,lim2 * write(6,*) egamma_tot_max-egamma_used(1) * write(6,*) vertex.e.E - edge.e.E.min * endif * if (phot3.gt.lim3) write(6,*) 'phot3,lim3=',phot3,lim3 !Calculate some proton/pion/kaon specific stuff: ! ! if (doing_pion .or. doing_kaon) then ! mm2 = recon.Em**2 - recon.Pm**2 ! mm = sqrt(abs(mm2)) * abs(mm2)/mm2 ! mmA2= (recon.nu + targ.M - recon.p.E)**2 - recon.Pm**2 ! mmA = sqrt(abs(mmA2)) * abs(mmA2)/mmA2 ! t = recon.Q2 - Mh2 ! > + 2*(recon.nu*recon.p.E - recon.p.P*recon.q*cos(recon.theta_pq)) ! endif if (doing_hyd_elast .or. doing_deuterium .or. doing_heavy) then poftheta = Mp*Ebeam / (2*ebeam*sin(recon.e.theta/2.)**2 + Mp) corrsing = recon.e.P - poftheta Pm_Heepz = -(recon.Pmy*recon.uq.y+recon.Pmz*recon.uq.z) > / sqrt(recon.uq.y**2+recon.uq.z**2) Pm_Heepy = (recon.Pmz*recon.uq.y-recon.Pmy*recon.uq.z) > / sqrt(recon.uq.y**2+recon.uq.z**2) Pm_Heepx = -recon.Pmx endif if(electron_arm.eq.1 .or. electron_arm.eq.3.or. electron_arm.eq.7)then !electron = right side. ntu(1) = recon.e.delta ntu(2) = recon.e.yptar !mr ntu(3) = recon.e.xptar !mr ntu(4) = recon.e.z ntu(5) = main.FP.e.x ntu(6) = main.FP.e.dx !mr ntu(7) = main.FP.e.y ntu(8) = main.FP.e.dy !mr ntu(9) = orig.e.delta ntu(10) = orig.e.yptar !mr ntu(11) = orig.e.xptar !mr ntu(12) = main.target.z*spec.e.sin_th ntu(13) = recon.p.delta ntu(14) = recon.p.yptar !mr ntu(15) = recon.p.xptar !mr ntu(16) = recon.p.z ntu(17) = main.FP.p.x ntu(18) = main.FP.p.dx !mr ntu(19) = main.FP.p.y ntu(20) = main.FP.p.dy !mr ntu(21) = orig.p.delta ntu(22) = orig.p.yptar !mr ntu(23) = orig.p.xptar !mr ntu(24) = -main.target.z*spec.p.sin_th else if (electron_arm.eq.2 .or. electron_arm.eq.4 .or. > electron_arm.eq.5 .or. electron_arm.eq.6.or. electron_arm.eq.8) then !e- = left. ntu(1) = recon.p.delta ntu(2) = recon.p.yptar !mr ntu(3) = recon.p.xptar !mr ntu(4) = recon.p.z ntu(5) = main.FP.p.x ntu(6) = main.FP.p.dx !mr ntu(7) = main.FP.p.y ntu(8) = main.FP.p.dy !mr ntu(9) = vertex.p.delta ntu(10) = vertex.p.yptar !mr ntu(11) = vertex.p.xptar !mr ntu(12) = main.target.z*spec.p.sin_th ntu(13) = recon.e.delta ntu(14) = recon.e.yptar !mr ntu(15) = recon.e.xptar !mr ntu(16) = recon.e.z ntu(17) = main.FP.e.x ntu(18) = main.FP.e.dx !mr ntu(19) = main.FP.e.y ntu(20) = main.FP.e.dy !mr ntu(21) = orig.e.delta ntu(22) = orig.e.yptar !mr ntu(23) = orig.e.xptar !mr ntu(24) = -main.target.z*spec.e.sin_th else write (6,*) 'results_write not yet set up for your spectrometers.' endif ntu(25) = recon.q/1000. !q - GeV/c ntu(26) = recon.nu/1000. !nu - GeV ntu(27) = recon.Q2/1.d6 !Q^2 - (GeV/c)^2 ntu(28) = recon.W/1000. !W - GeV/c ntu(29) = recon.epsilon !epsilon ntu(30) = recon.Em/1000. !GeV ntu(31) = recon.Pm/1000. !GeV/c ntu(32) = recon.theta_pq !theta_pq - radians ntu(33) = recon.phi_pq !phi_pq - radians if (doing_pion .or. doing_kaon .or. doing_delta) then ntu(34) = ntup.mm/1000. !missmass (nucleon) ntu(35) = ntup.mmA/1000. !missmass (nucleus) ntu(36) = recon.p.P/1000. !ppi - GeV/c ntu(37) = ntup.t/1.d6 !t - GeV^2 ntu(38) = recon.PmPar/1000. ntu(39) = recon.PmPer/1000. ntu(40) = recon.PmOop/1000. ntu(41) = -main.target.rastery !fry - cm ntu(42) = ntup.radphot/1000. !radphot - GeV dummy = pferx*vertex.uq.x + pfery*vertex.uq.y + pferz*vertex.uq.z if (dummy.eq.0) dummy=1.d-20 ntu(43) = pfer/1000.*abs(dummy)/dummy !p_fermi - GeV/c ntu(44) = main.sigcc !d5sig ntu(45) = ntup.sigcm !pion sig_cm ntu(46) = main.weight ntu(47) = decdist !decay distance (cm) ntu(48) = sqrt(Mh2_final) ntu(49) = pfer/1000.*dummy !p_fermi along q. ntu(50) = vertex.Q2/1.d6 ntu(51) = main.w/1.d3 ntu(52) = main.t/1.d6 ntu(53) = main.phi_pq if(using_tgt_field) then ntu(54) = recon.theta_tarq ntu(55) = recon.phi_targ ntu(56) = recon.beta ntu(57) = recon.phi_s ntu(58) = recon.phi_c ntu(59) = main.beta ntu(60) = vertex.phi_s ntu(61) = vertex.phi_c if (doing_kaon) then ntu(62) = ntup.sigcm1 !sigcm - saghai model ntu(63) = ntup.sigcm2 !sigcm - factorized. endif else if (doing_kaon) then ntu(54) = ntup.sigcm1 !sigcm - saghai model ntu(55) = ntup.sigcm2 !sigcm - factorized. endif endif else if (doing_semi.or.doing_rho) then ntu(34) = ntup.mm/1000. !missmass (nucleon) ntu(35) = recon.p.P/1000. !ppi - GeV/c ntu(36) = ntup.t/1.d6 !t - GeV^2 ntu(37) = -main.target.rastery !fry - cm ntu(38) = ntup.radphot/1000. !radphot - GeV ntu(39) = main.sigcc !d5sig ntu(40) = main.sigcent !d5sig - central kin. ntu(41) = main.weight ntu(42) = decdist !decay distance (cm) ntu(43) = sqrt(Mh2_final) ntu(44) = recon.zhad ntu(45) = vertex.zhad ntu(46) = recon.pt2/1.d06 ntu(47) = vertex.pt2/1.d06 ntu(48) = recon.xbj ntu(49) = vertex.xbj ntu(50) = acos(vertex.uq.z) ntu(51) = ntup.sigcm ntu(52) = main.davejac ntu(53) = main.johnjac dummy = pferx*vertex.uq.x + pfery*vertex.uq.y + pferz*vertex.uq.z ntu(54) = pfer/1000.*abs(dummy)/dummy !p_fermi - GeV/c ntu(55) = ntup.xfermi ntu(56) = main.phi_pq if(using_tgt_field) then ntu(57) = recon.theta_tarq ntu(58) = recon.phi_targ ntu(59) = recon.beta ntu(60) = recon.phi_s ntu(61) = recon.phi_c ntu(62) = main.beta ntu(63) = vertex.phi_s ntu(64) = vertex.phi_c if(doing_rho) then ntu(65) = ntup.rhomass ntu(66) = ntup.rhotheta endif else if(doing_rho) then ntu(57) = ntup.rhomass ntu(58) = ntup.rhotheta endif endif else if (doing_hyd_elast .or. doing_deuterium .or. doing_heavy) then ntu(34) = corrsing/1000. ntu(35) = Pm_Heepx/1000. ntu(36) = Pm_Heepy/1000. ntu(37) = Pm_Heepz/1000. ntu(38) = recon.PmPar/1000. ntu(39) = recon.PmPer/1000. ntu(40) = recon.PmOop/1000. ntu(41) = -main.target.rastery !fry - cm ntu(42) = ntup.radphot/1000. !radphot - GeV ntu(43) = main.sigcc ntu(44) = main.weight endif call HFN(NtupleID,ntu) if (debug(2)) write(6,*)'r_ntu_write: ending' return end subroutine results_ntu_write1(vertex,recon,main,success) implicit none include 'hbook.inc' include 'simulate.inc' integer*4 nentries parameter (nentries = 9) real*8 ntu(nentries) logical success record /event_main/ main record /event/ vertex, recon if (debug(2)) write(6,*)'r_ntu_write: entering' ntu(1) = vertex.e.delta ntu(2) = vertex.e.yptar ntu(3) = -vertex.e.xptar ntu(4) = main.SP.e.z if(success)then ntu(5) = recon.e.delta ntu(6) = recon.e.yptar ntu(7) = recon.e.xptar ntu(8) = recon.e.z else ntu(5) = 30. ntu(6) = 0.1 ntu(7) = 0.1 ntu(8) = 4. endif ntu(9) = main.weight call HFN(NtupleID,ntu) if (debug(2)) write(6,*)'r_ntu_write: ending' return end
Analyzer/Replay: Mark Jones, Documents: Stephen Wood |
Powered by ViewCVS 0.9.2-cvsgraph-1.4.0 |