// $Id$ // // File: JEventProcessor_neutron.cc // Created: Wed Jun 1 13:54:40 EDT 2022 // Creator: staylor (on Linux ifarm1802.jlab.org 3.10.0-1160.11.1.el7.x86_64 x86_64) // #include "JEventProcessor_neutron.h" using namespace jana; // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_neutron()); } } // "C" //------------------ // JEventProcessor_neutron (Constructor) //------------------ JEventProcessor_neutron::JEventProcessor_neutron() { } //------------------ // ~JEventProcessor_neutron (Destructor) //------------------ JEventProcessor_neutron::~JEventProcessor_neutron() { } //------------------ // init //------------------ jerror_t JEventProcessor_neutron::init(void) { japp->CreateLock("myneutronlock"); PI0_CUT_VALUE=0.025; ETA_CUT_VALUE=0.05; ETAPRIME_CUT_VALUE=0.05; NUM_SIGMA_BG=2; gPARMS->SetDefaultParameter("SURVEY:NUM_SIGMA_BG",NUM_SIGMA_BG); CL_CUT=0.001; gPARMS->SetDefaultParameter("NEUTRON:CL_CUT",CL_CUT); BETA_CUT=0.9; gPARMS->SetDefaultParameter("NEUTRON:BETA_CUT",BETA_CUT); SPLIT_CUT=0.5; gPARMS->SetDefaultParameter("NEUTRON:SPLIT_CUT",SPLIT_CUT); FCAL_POS_CUT=10.0; gPARMS->SetDefaultParameter("NEUTRON:FCAL_POS_CUT",FCAL_POS_CUT); FCAL_RADIAL_CUT=104.0; gPARMS->SetDefaultParameter("NEUTRON:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT); BCAL_Z_CUT=384.0; gPARMS->SetDefaultParameter("NEUTRON:BCAL_Z_CUT",BCAL_Z_CUT); FCAL_THRESHOLD=0.1; gPARMS->SetDefaultParameter("NEUTRON:FCAL_THRESHOLD",FCAL_THRESHOLD); BCAL_THRESHOLD=0.05; gPARMS->SetDefaultParameter("NEUTRON:BCAL_THRESHOLD",BCAL_THRESHOLD); gDirectory->mkdir("neutronAnalysis")->cd(); gDirectory->mkdir("MC")->cd(); { MCM_vs_E=new TH2F("MCM_vs_E",";E_{#gamma} [GeV];M(#eta#pi^{+}) [GeV]",48,2.8,12.4,1600,0,4); MCM_vs_t=new TH2F("MCM_vs_t",";|t| [GeV^{2}];M(#eta#pi^{+}) [GeV]",200,0,2,1600,0,4); MC_EtaCosThetaGJ_n=new TH2F("MC_EtaCosThetaGJ_n",";M(#eta#pi^{+}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); { NeutralBeta=new TH1F("NeutralBeta",";#beta",100,0,1.5); GammaRFTimeDiff=new TH1F("GammaRFTimeDiff",";t(#gamma)-t(rf) [ns]", 1000,-10,10); BeamPhotonTimeDiff=new TH1F("BeamPhotonTimeDiff","; t(tag)-t(rf) [ns]", 1000,-50,50); MissingPionChiSq=new TH1F("MissingPionChiSq",";#chi^{2}/ndf",1000,0,100.); MissingPionCL=new TH1F("MissingPionCL",";CL",1000,0,1.); ExclusiveCL=new TH1F("ExclusiveCL",";CL",1000,0,1.); ExclusiveChiSq=new TH1F("ExclusiveChiSq",";#chi^{2}/ndf",1000,0,100.); UnusedEnergy=new TH1F("UnusedEnergy","unused shower energy",100,0,0.5); } gDirectory->cd("../"); gDirectory->mkdir("missingPip")->cd(); { TwoGammaMass_missingPip=new TH1F("TwoGammaMass_missingPip",";M(2#gamma) [GeV]",400,0,1.); DeltaPOverP_vs_theta=new TH2F("DeltaPOverP_vs_theta",";#theta [#circ];#deltap/p", 140,0,140,60,-0.3,0.3); DeltaPOverP_vs_P=new TH2F("DeltaPOverP_vs_P",";p [GeV/c];#deltap/p", 50,0,10,60,-0.3,0.3); DeltaTheta_vs_P=new TH2F("DeltaTheta_vs_P",";p [GeV/c];#delta#theta [#circ]", 50,0,10,100,-20.,20.); DeltaPhi_vs_P=new TH2F("DeltaPhi_vs_P",";p [GeV/c];#delta#phi [#circ]", 50,0,10,100,-20.,20.); } gDirectory->cd("../"); gDirectory->mkdir("n2gamma")->cd(); { TwoGammaMass_NPip=new TH1F("TwoGammaMass_NPip",";M(2#gamma) [GeV]",400,0,1.); PipPi0_NPi0_Dalitz=new TH2F("PipPi0_NPi0_Dalitz", ";M^{2}(#pi^{+}#pi^{0}) [GeV^{2};M^{2}(n#pi^{+}) [GeV^{2}]",1000,0,20,1000,0,20); PipEta_NEta_Dalitz=new TH2F("PipEta_NEta_Dalitz", ";M^{2}(#pi^{+}#eta}) [GeV^{2};M^{2}(n#pi^{+}) [GeV^{2}]",1000,0,20,1000,0,20); Pi0CosThetaGJ_n_bg=new TH2F("Pi0CosThetaGJ_n_bg",";M(#pi^{0}#pi^{+}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); Pi0CosThetaGJ_n=new TH2F("Pi0CosThetaGJ_n",";M(#pi^{#pi^{0}#pi^{+}}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); EtaCosThetaGJ_n_bg=new TH2F("EtaCosThetaGJ_n_bg",";M(#eta#pi^{+}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); EtaCosThetaGJ_n=new TH2F("EtaCosThetaGJ_n",";M(#eta#pi^{+}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); PipEtaMass_vs_E=new TH2F("PipEtaMass_vs_E",";E_{#gamma} [GeV];M(#eta#pi^{+}) [GeV]",48,2.8,12.4,1600,0,4); PipEtaMass_vs_t=new TH2F("PipEtaMass_vs_t",";|t| [GeV^{2}];M(#eta#pi^{+}) [GeV]",200,0,2,1600,0,4); PipPi0Mass_vs_E=new TH2F("PipPi0Mass_vs_E",";E_{#gamma} [GeV];M(#eta#pi^{+}) [GeV]",48,2.8,12.4,1600,0,4); PipPi0Mass_vs_t=new TH2F("PipPi0Mass_vs_t",";|t| [GeV^{2}];M(#eta#pi^{+}) [GeV]",200,0,2,1600,0,4); PipEtaPrimeMass_vs_E=new TH2F("PipEtaPrimeMass_vs_E",";E_{#gamma} [GeV];M(#eta#'pi^{+}) [GeV]",48,2.8,12.4,1600,0,4); PipEtaPrimeMass_vs_t=new TH2F("PipEtaPrimeMass_vs_t",";|t| [GeV^{2}];M(#eta'#pi^{+}) [GeV]",200,0,2,1600,0,4); EtaPrimeCosThetaGJ_n=new TH2F("EtaPrimeCosThetaGJ_n",";M(#eta'#pi^{+}) [GeV];cos(#theta_{GJ})",1600,0,4,100,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("n3gamma")->cd(); { ThreeGam_vs_NPip=new TH2F("ThreeGam_vs_NPip",";M(n#pi^{+}) [GeV];M(3#gamma) [GeV]",1600,0.9,4.9,1600,0.,4.); OmegaPip_n3g=new TH1F("OmegaPip_n3g",";M(#omega#pi^{+}) [GeV]",1600,0.5,4.5); } gDirectory->cd("../"); gDirectory->mkdir("n4gamma")->cd(); { FourGammaMass_NPip=new TH1F("FourGammaMass_NPip",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d_NPip=new TH2F("FourGamma2d_NPip", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); Pi0Pi0_NPip=new TH1F("Pi0Pi0_NPip",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_NPip=new TH1F("Pi0Eta_NPip",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_NPip=new TH1F("EtaEta_NPip",";M(2#eta) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("n6gamma")->cd(); { SixGammaMass_NPi=new TH1F("SixGammaMass_NPi",";M(6#gamma) [GeV]",1600,0,4.); SixGamma3D_NPi=new TH3F("SixGamma3D_NPi", ";M(2#gamma,pair1) [GeV];M(2#gamma,pair2) [GeV];M(2#gamma,pair3) [GeV]", 200,0,1,200,0,1,200,0,1); ThreePi0Mass_NPi=new TH1F("ThreePi0Mass_NPi",";M(3#pi^{0}) [GeV]",1600,0,4.); TwoPi0EtaMass_NPi=new TH1F("TwoPi0EtaMass_NPi",";M(2#pi^{0}#eta) [GeV]",1600,0,4.); ThreeEtaMass_NPi=new TH1F("ThreeEtaMass_NPi",";M(3#eta) [GeV]",1600,0,4.); TwoEtaPi0Mass_NPi=new TH1F("TwoEtaPi0Mass_NPi",";M(2#eta#pi^{0}) [GeV]",1600,0,4.); TwoEtaPi0Dalitz_NPi=new TH2F("TwoEtaPi0Dalitz_NPi", ";M^{2}(#eta_{1}#pi^{0}) [GeV^{2}];M^{2}(#eta_{2}#pi^{0}) [GeV^{2}]",100,0,10,100,0,10); Pip6G_N=new TH1F("Pip6G_N",";M(#pi^{+}6#gamma) [GeV]",1600,0,4.); PipEta_N6g=new TH1F("PipEta_N6g",";M(#pi^{+}#eta) [GeV]",1600,0,4.); Pi0gVsPi0g_NPi=new TH2F("Pi0gVsPi0g_NPi", ";M(#pi^{0}#gamma, group 1) [GeV];M(#pi^{0}#gamma, group 2) [GeV]",1600,0,4.,1600,0,4); } gDirectory->cd("../"); gDirectory->mkdir("n2PipPim")->cd(); { NeutronPipMass_2pip=new TH1F("NeutronPipMass_2pip",";M(n#pi^{+}) [GeV]",1600,0.9,4.9); PipPim_vs_NPip=new TH2F("PipPim_vs_NPip","M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); } gDirectory->cd("../"); gDirectory->mkdir("n2PipPim1G")->cd(); { NeutronPipMass_2pip1G=new TH1F("NeutronPipMass_2pip1G",";M(n#pi^{+}) [GeV]",1600,0.9,4.9); PipPim_vs_NPip_1G=new TH2F("PipPim_vs_NPip_1G","M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); PipPim1G_N=new TH1F("PipPim1G_N",";M(#pi^{+}#pi-#gamma) [GeV]",1600,0.,4.); } gDirectory->cd("../"); gDirectory->mkdir("n2PipPim2G")->cd(); { TwoGammaMass_n2pip=new TH1F("TwoGammaMass_n2pip",";M(2#gamma) [GeV]",400,0,1.); TwoGammaMass_n2pip_NU=new TH1F("TwoGammaMass_n2pip_NU",";M(2#gamma) [GeV]",400,0,1.); NeutronPipMass_2pip2g=new TH1F("NeutronPipMass_2pip2g",";M(n#pi^{+}) [GeV]",1600,0.9,4.9); TwoPipPimPi0_N=new TH1F("TwoPipPimPi0_N",";M(2#pi^{+}#pi^{-}#pi^{0}) [GeV]",1600,0.,4.); PipPimPi0_n=new TH1F("PipPimPi0_n",";M(pi^{+}#pi^{-}#pi^{0}) [GeV]",1600,0.,4.); EtaPip_N=new TH1F("EtaPip_N",";M(#eta#pi^{+}) [GeV]",1600,0.,4.); OmegaPip_N2g=new TH1F("OmegaPip_N2g",";M(#omega#pi^{+}) [GeV]",1600,0.,4.); NEta_vs_2PipPim=new TH2F("NEta_vs_2PipPim", ";M(2#pi^{+}#pi^{-}) [GeV];M(n#eta}) [GeV]", 1600,0.,4.,1600,0.9,4.9); NPi0_vs_2PipPim=new TH2F("NPi0_vs_2PipPim", ";M(2#pi^{+}#pi^{-}) [GeV];M(n#pi^{0}) [GeV]", 1600,0.,4.,1600,0.9,4.9); PipPim_vs_NPip_with_Pi0=new TH2F("PipPim_vs_NPip_with_Pi0", "M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); SigmaPlusPi0_SigmaPlusK0_Dalitz=new TH2F("SigmaPlusPi0_SigmaPlusK0_Dalitz", ";M^{2}(#Sigma^{+}#pi^{0}) [GeV^{2}];M^{2}(#Sigma^{+}K_{S}) [GeV^{2}]",100,1.,19.,100,1.,19.); TwoPipPimEta_N=new TH1F("TwoPipPimEta_N",";M(2#pi^{+}#pi^{-}#eta) [GeV]",1600,0.,4.); PipPimEta_n=new TH1F("PipPimEta_n",";M(pi^{+}#pi^{-}#eta) [GeV]",1600,0.,4.); EtaPrimePip_N=new TH1F("EtaPrimePip_N",";M(#eta#prime#pi^{+}) [GeV]",1600,0.,4.); EtaPrimePip_N_NU=new TH1F("EtaPrimePip_N_NU",";M(#eta#prime#pi^{+}) [GeV]",1600,0.,4.); PipPim_vs_NPip_with_Eta=new TH2F("PipPim_vs_NPip_with_Eta", "M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); SigmaPlusEta_SigmaPlusK0_Dalitz=new TH2F("SigmaPlusEta_SigmaPlusK0_Dalitz", ";M^{2}(#Sigma^{+}#eta) [GeV^{2}];M^{2}(#Sigma^{+}K_{S}) [GeV^{2}]",100,1.,19.,100,1.,19.); } gDirectory->cd("../"); gDirectory->mkdir("n2PipPim4G")->cd(); { FourGammaMass_N3Pi=new TH1F("FourGammaMass_N3Pi",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d_N3Pi=new TH2F("FourGamma2d_N3Pi", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); EtaPi0_N3Pi=new TH1F("EtaPi0_N3Pi",";M(#eta#pi^{0}) [GeV]",1600,0.,4.); TwoPi0_N3Pi=new TH1F("TwoPi0_N3Pi",";M(2#pi^{0}) [GeV]",1600,0.,4.); PipPim_vs_NPip_with_2Pi0=new TH2F("PipPim_vs_NPip_with_2Pi0", "M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); PipPim_vs_NPip_with_EtaPi0=new TH2F("PipPim_vs_NPip_with_EtaPi0", "M(#pi+#pi-) vs M(Npi+)", 1600,0.9,4.9,1600,0,4); PipPimPi0_N_with_Pi0=new TH1F("PipPimPi0_N_with_Pi0",";M(pi^{+}#pi^{-}#pi^{0}) [GeV]",1600,0.,4.); PipPimEta_N_with_Pi0=new TH1F("PipPimEta_N_with_Pi0",";M(pi^{+}#pi^{-}#eta}) [GeV]",1600,0.,4.); PipPimPi0_N_with_Eta=new TH1F("PipPimPi0_N_with_Eta",";M(pi^{+}#pi^{-}#pi^{0}) [GeV]",1600,0.,4.); EtaPi0_N_eta_3pi=new TH1F("EtaPi0_N_eta_3pi",";M(#eta#pi^{0}) [GeV]",1600,0.,4.); } gDirectory->cd("../"); gDirectory->mkdir("nKpPipPim")->cd(); { NPim_vs_KpPip=new TH2F("NPim_vs_KpPip", ";M(K^{+}#pi^{+}) [GeV];M(n#pi^{-}) [GeV]", 800,0.5,2.5,800,1.,3.); NPip_vs_KpPim=new TH2F("NPip_vs_KpPim", ";M(K^{+}#pi^{-}) [GeV];M(n#pi^{+}) [GeV]", 800,0.5,2.5,800,1.,3.); PipPim_vs_nKp=new TH2F("PipPim_vs_nKp", ";M(K^{+}n) [GeV];M(#pi^{+}#pi^{-}) [GeV]", 1200,1.3,4.3,800,0.2,2.2); NKs=new TH1F("NKs",";M(nK^{0}_{S}) [GeV]",1200,1.3,4.3); SigmapPim=new TH1F("SigmapPim",";M(#Sigma^{+}#pi^{-}) [GeV]",1200,1.,4.); SigmamPip=new TH1F("SigmamPip",";M(#Sigma^{-}#pi^{+}) [GeV]",1200,1.,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpPimPimN1G")->cd(); { NPipSqVsNPimSq_1G=new TH2F("NPipSqVSNPimSq_1G","M^2(n#pi+) vs M^2(n#pi-)", 240,1.,9,240,1.,9); SigmaPlusGamma= new TH1F("SigmaPlusGamma","#Sigma+#gamma mass",1600,0.9,4.9); SigmaMinusGamma= new TH1F("SigmaMinusGamma","#Sigma-#gamma mass",1600,0.9,4.9); } gDirectory->cd("../"); gDirectory->mkdir("KpPimPimN2G")->cd(); { NPipSqVsNPimSq_2G=new TH2F("NPipSqVSNPimSq_2G","M^2(n#pi+) vs M^2(n#pi-)", 240,1.,9,240,1.,9); TwoGamma_with_KNPipPim= new TH1F("TwoGamm_with_KNPipPim","2#gamma mass",1600,0,4.); SigmaPlusPi0_vs_KPim=new TH2F("SigmaPlusPi0_vs_KPim","M(#Sigma+#pi0) vs M(K+#pi-)",1600,0,4,1600,0,4); SigmaPlusPi0_SigmaPlusPim_Dalitz=new TH2F("SigmaPlusPi0_SigmaPlusPim_Dalitz","M^2(#Sigma+#pi-) vs M^2(#Sigma+#pi0)",1000,0,20,1000,0,20); SigmaPlusPim_vs_KPi0=new TH2F("SigmaPlusPim_vs_KPi0","M(#Sigma+#pi-) vs M(K+#pi0)",1600,0,4,1600,0,4); SigmaMinusPip_vs_KPi0=new TH2F("SigmaMinusPip_vs_KPi0","M(#Sigma-#pi+) vs M(K+#pi0)",1600,0,4,1600,0,4); SigmaMinusPi0_vs_KPip=new TH2F("SigmaMinusPi0_vs_KPip","M(#Sigma-#pi0) vs M(K+#pi+)",1600,0,4,1600,0,4); SigmaMinusPi0_SigmaMinusPip_Dalitz=new TH2F("SigmaMinusPi0_SigmaMinusPip_Dalitz","M^2(#Sigma-#pi+) vs M^2(#Sigma-#pi0)",1000,0,20,1000,0,20); SigmamPipPi0=new TH1F("SigmamPipPi0",";M(#Sigma-#pi+#pi0) [GeV]",1600,1.,5.); SigmapPimPi0=new TH1F("SigmapPimPi0",";M(#Sigma+#pi-#pi0) [GeV]",1600,1.,5.); SigmaPlusEta_vs_KPim=new TH2F("SigmaPlusEta_vs_KPim","M(#Sigma+#eta) vs M(K+#pi-)",1600,0,4,1600,0,4); SigmaPlusEta_SigmaPlusPim_Dalitz=new TH2F("SigmaPlusEta_SigmaPlusPim_Dalitz","M^2(#Sigma+#pi-) vs M^2(#Sigma+#eta)",1000,0,20,1000,0,20); SigmaPlusPim_vs_KEta=new TH2F("SigmaPlusPim_vs_KEta","M(#Sigma+#pi-) vs M(K+#eta)",1600,0,4,1600,0,4); SigmaMinusPip_vs_KEta=new TH2F("SigmaMinusPip_vs_KEta","M(#Sigma-#pi+) vs M(K+#eta)",1600,0,4,1600,0,4); SigmaMinusEta_vs_KPip=new TH2F("SigmaMinusEta_vs_KPip","M(#Sigma-#eta) vs M(K+#pi+)",1600,0,4,1600,0,4); SigmaMinusEta_SigmaMinusPip_Dalitz=new TH2F("SigmaMinusEta_SigmaMinusPip_Dalitz","M^2(#Sigma-#pi+) vs M^2(#Sigma-#eta)",1000,0,20,1000,0,20); SigmamPipEta=new TH1F("SigmamPipEta",";M(#Sigma-#pi+#eta) [GeV]",1600,1.,5.); SigmapPimEta=new TH1F("SigmapPimEta",";M(#Sigma+#pi-#eta) [GeV]",1600,1.,5.); } gDirectory->cd("../"); gDirectory->mkdir("Km2KpN2G")->cd(); { TwoGamma_2kpkmn= new TH1F("TwoGamma_2kpkmn",";M(2#gamma) [GeV]",1600,0,4.); NPi0_2kpkm= new TH1F("NPi0_2kpkm",";M(n#pi^{0})[GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpN2g")->cd(); { TwoGammaMass_NK= new TH1F("TwoGammaMass_NK","2#gamma mass",1600,0,4.); NPi0Mass= new TH1F("NPi0Mass","N#pi0 mass",1600,0,4.); NEtaMass= new TH1F("NEtaMass","N#eta mass",1600,0,4.); N2gMass= new TH1F("N2gMass","N2#gamma mass",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("p2n")->cd(); { TwoNeutronVertex=new TH2F("TwoNeutronVertex",";z[cm];r[cm]",140,30,100,20,0,2); TwoNeutronMomentum2D=new TH2F("TwoNeutronMomentum2D",";p(n_{1}) [GeV/c];p(n_2}) [GeV/c]",200,0,10,200,0,10); NeutronPVsTheta=new TH2F("NeutronPVsTheta",";#theta[^circ];p [GeV/c]",180,0,90,200,0,10); ProtonPVsThetaWithN=new TH2F("ProtonPVsThetaWithN",";#theta[^circ];p [GeV/c]",180,0,90,200,0,10); TwoNeutronMass= new TH1F("TwoNeutronMass",";M(2n) [GeV]",1600,1,5.); NantiN_vs_t=new TH2F("NantiN_vs_t",";-t [GeV^{2}];M(n#bar{n}) [GeV]", 1000,0,20.,1600,1.,5); NantiN_vs_E=new TH2F("NantiN_vs_E",";E_{beam} [GeV];M(n#bar{n}) [GeV]", 48,2.8,12.4,1600,1.,5); PNDalitz=new TH2F("PNDalitz",";M^{2}(pn_{1}) [GeV^{2}];M^{2}(pn_2}) [GeV^{2}]",1000,0,20,1000,0,20); } gDirectory->cd("../"); gDirectory->mkdir("p2nPipPim")->cd(); { NeutronPip_vs_NeutronPim= new TH2F("NeutronPip_vs_NeutronPim",";M(n#pi^{-}) [GeV];M(n#pi^{+}) [GeV]",1600,1,5.,1600,1.,5.); } gDirectory->cd("../"); gDirectory->mkdir("pnapPip")->cd(); { AntiPNPip= new TH1F("AntiPNPip",";M(n#bar{p}#pi^{+}) [GeV]",1600,1,5.); } gDirectory->cd("../"); gDirectory->mkdir("pnapPip2g")->cd(); { AntiPNPip_vs_2g= new TH2F("AntiPNPip_vs_2g",";M(2#gamma) [GeV];M(n#bar{p}#pi^{+}) [GeV]",1600,0,4.,1600,1.,5.); NPip_vs_AntiPPi0= new TH2F("NPip_vs_AntiPPi0",";M(#bar{p}#pi^{0|}) [GeV];M(n#pi^{+}) [GeV]",1600,0,4.,1600,0.,4.); NPi0_vs_AntiPPiP= new TH2F("NPi0_vs_AntiPPiP",";M(#bar{p}#pi^{+|}) [GeV];M(n#pi^{0}) [GeV]",1600,0,4.,1600,0.,4.); } gDirectory->cd("../"); gDirectory->mkdir("ppanPim2g")->cd(); { TwoGammaMass_with_2pnpim= new TH1F("TwoGammaMass_with_2pnpim",";M(2#gamma) [GeV]",1600,0,4.); PPim_vs_AntiNPi0= new TH2F("PPim_vs_AntiNPi0",";M(#bar{n}#pi^{0}) [GeV];M(p#pi^{-|}) [GeV]",1600,0,4.,1600,0.,4.); PPi0_vs_AntiNPim= new TH2F("PPi0_vs_AntiNPim",";M(bar{n}#pi^{-}) [GeV];M(p#pi^{-|}) [GeV]",1600,0,4.,1600,0.,4.); } gDirectory->cd("../"); gDirectory->cd("../"); // Make a container to keep track of doublets of photons vector>pairs; pairs.push_back(make_pair(0,1)); pairs.push_back(make_pair(2,3)); doublets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,2)); pairs.push_back(make_pair(1,3)); doublets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,3)); pairs.push_back(make_pair(2,1)); doublets.push_back(pairs); // .. and triplets of photons... MakeTriplets(); // set up maps for keeping track of pairs of mesons vector>my_pair_vector; particle_pairs.emplace(Pi0Pi0_,my_pair_vector); particle_pairs.emplace(Pi0Eta_,my_pair_vector); particle_pairs.emplace(Pi0EtaPrime_,my_pair_vector); particle_pairs.emplace(EtaEta_,my_pair_vector); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_neutron::brun(JEventLoop *eventLoop, int32_t runnumber) { // This is called whenever the run number changes DApplication *dapp=dynamic_cast(eventLoop->GetJApplication()); const DGeometry *geom=dapp->GetDGeometry(runnumber); double x0,y0,z0; geom->GetFCALPosition(x0,y0,z0); mFCALCenter.SetXYZ(x0,y0,z0); return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_neutron::evnt(JEventLoop *loop, uint64_t eventnumber) { vectorthrowns; loop->Get(throwns); if (throwns.size()>0){ const DTAGMGeometry *tagm_geom=NULL; loop->GetSingle(tagm_geom); if (tagm_geom==NULL) return RESOURCE_UNAVAILABLE; const DTAGHGeometry *tagh_geom=NULL; loop->GetSingle(tagh_geom); if (tagh_geom==NULL) return RESOURCE_UNAVAILABLE; const DMCReaction* reaction=NULL; loop->GetSingle(reaction); if (reaction==NULL) return RESOURCE_UNAVAILABLE; japp->RootWriteLock(); double mc_beam_E=reaction->beam.energy(),tagged_mc_beam_E=0.; unsigned int counter=0,column=0; if (tagm_geom->E_to_column(mc_beam_E,column)){ tagged_mc_beam_E=0.5*(tagm_geom->getEhigh(column) +tagm_geom->getElow(column)); } else if (tagh_geom->E_to_counter(mc_beam_E,counter)){ tagged_mc_beam_E=0.5*(tagh_geom->getEhigh(counter) +tagh_geom->getElow(counter)); } if (tagged_mc_beam_E>0.){ DLorentzVector target(0,0,0,ParticleMass(Proton)); DLorentzVector beam(0,0,tagged_mc_beam_E,tagged_mc_beam_E); vectormcpips,mcneutrons,mcetas,mcgammas,mcpi0s; for (unsigned int i=0;imech==0){ switch(throwns[i]->type){ case PiPlus: mcpips.push_back(throwns[i]->lorentzMomentum()); break; case Neutron: mcneutrons.push_back(throwns[i]->lorentzMomentum()); break; case Eta: mcetas.push_back(throwns[i]->lorentzMomentum()); break; case Gamma: mcgammas.push_back(throwns[i]->lorentzMomentum()); break; case Pi0: mcpi0s.push_back(throwns[i]->lorentzMomentum()); break; default: break; } } } if (mcneutrons.size()==1){ double t=(target-mcneutrons[0]).M2(); DLorentzVector missing=beam+target-mcneutrons[0]; double mass=missing.M(); MCM_vs_E->Fill(tagged_mc_beam_E,mass); MCM_vs_t->Fill(-t,mass); if (mcetas.size()==1){ FillGJHisto(1.,beam,missing,mcetas[0],MC_EtaCosThetaGJ_n); } } } japp->RootUnLock(); } vectorneutrals; loop->Get(neutrals); if (neutrals.size()==0) return RESOURCE_UNAVAILABLE; vectortracks; loop->Get(tracks); if (tracks.size()==0) return RESOURCE_UNAVAILABLE; vectorbeamphotons; loop->Get(beamphotons); if (beamphotons.size()==0) return RESOURCE_UNAVAILABLE; japp->WriteLock("myneutronlock"); // Assign charge track by PID map>chargedParticles; FillChargedParticleVectors(tracks,chargedParticles); // Find t0 (at "vertex") for event and assign list of neutral particles double t0_rf=tracks[0]->Get_BestTrackingFOM()->t0(); DVector3 vertex=tracks[0]->Get_BestTrackingFOM()->position(); // Assign neutrals by ID double unused_energy=0.; map>neutralParticles; FillNeutralParticleVectors(unused_energy,t0_rf,vertex,neutrals,neutralParticles); UnusedEnergy->Fill(unused_energy); if (neutralParticles[Neutron].size()>0){ vectorneutronHyps=neutralParticles[Neutron]; // Fill vector of beam photons to// Fill vector of beam photons to pass on the analysis along with weights // for in-time/ out-of-time pass on the analysis along with weights // for in-time/ out-of-time vector>beamPhotons; for (unsigned int i=0;itime()-t0_rf; BeamPhotonTimeDiff->Fill(dt_st); double weight=1.; bool got_beam_photon=false; if (fabs(dt_st)>6.012&&fabs(dt_st)<18.036){ weight=-1./6.; got_beam_photon=true; } if (fabs(dt_st)<2.004){// Fill vector of beam photons to pass on the analysis along with weights // for in-time/ out-of-time got_beam_photon=true; } if (got_beam_photon){ beamPhotons.push_back(make_pair(beamphotons[i],weight)); } } // Setup for performing kinematic fits DAnalysisUtilities *dAnalysisUtilities=new DAnalysisUtilities(loop); DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop); DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils); dKinFitUtils->Reset_NewEvent(); dKinFitter->Reset_NewEvent(); for (unsigned int i=0;iGet_Hypothesis(PiPlus); if (neutralParticles[Neutron].size()==1 &&neutralParticles[Gamma].size()==2 && tracks.size()==1&&pionHyp!=NULL ){ DoKinematicFit(t0_rf,vertex,beamPhotons[i].first,neutronHyps, chargedParticles,neutralParticles[Gamma],dKinFitUtils, dKinFitter,true); double chisq=dKinFitter->Get_ChiSq(); double ndf=dKinFitter->Get_NDF(); double CL=dKinFitter->Get_ConfidenceLevel(); MissingPionChiSq->Fill(chisq/ndf,weight); MissingPionCL->Fill(CL,weight); if (CL>0.1){ DLorentzVector beam_kf,missing_kf; map>final_kf; GetKF4vectors(dKinFitter,beam_kf,missing_kf,final_kf); const DTrackTimeBased *pion_track=pionHyp->Get_TrackTimeBased(); MissingPionAnalysis(pion_track,missing_kf,final_kf,weight); } } // Other topologies with a neutron in the final state DoKinematicFit(t0_rf,vertex,beamPhotons[i].first,neutronHyps, chargedParticles,neutralParticles[Gamma],dKinFitUtils, dKinFitter,false); double CL=dKinFitter->Get_ConfidenceLevel(); double chisq=dKinFitter->Get_ChiSq(); double ndf=dKinFitter->Get_NDF(); ExclusiveChiSq->Fill(chisq/ndf,weight); ExclusiveCL->Fill(CL,weight); if (CL>CL_CUT){ double weight=beamPhotons[i].second; DLorentzVector beam_kf,missing_kf; map>final_kf; GetKF4vectors(dKinFitter,beam_kf,missing_kf,final_kf); // Particle counts unsigned int numPiPlus=final_kf[PiPlus].size(); unsigned int numPiMinus=final_kf[PiMinus].size(); unsigned int numKPlus=final_kf[KPlus].size(); unsigned int numKMinus=final_kf[KMinus].size(); unsigned int numGamma=final_kf[Gamma].size(); unsigned int numNeutron=final_kf[Neutron].size(); unsigned int numProton=final_kf[Proton].size(); unsigned int numAntiProton=final_kf[AntiProton].size(); if (numNeutron==1){ if (numProton==0&&numAntiProton==0){ if (numPiPlus==1 && numPiMinus==0 && numKPlus==0 && numKMinus==0){ switch(numGamma){ case 2: NeutronTwoGammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 3: NeutronThreeGammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 4: NeutronFourGammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 6: NeutronSixGammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } if (numPiPlus==2 && numPiMinus==1 && numKPlus==0 && numKMinus==0){ switch(numGamma){ case 0: Neutron2PipPimAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 1: Neutron2PipPim1GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 2: Neutron2PipPim2GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 4: Neutron2PipPim4GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } if (numPiPlus==1 && numPiMinus==1 && numKPlus==1 && numKMinus==0){ switch(numGamma){ case 0: NeutronKpPipPimAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 1: NeutronKpPipPim1GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 2: NeutronKpPipPim2GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } if (numKPlus==2 && numPiMinus==0 && numKMinus==1 && numPiPlus==0){ switch(numGamma){ case 2: TwoKpKmN2GAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } if (numKPlus==1 && numPiMinus==0 && numKMinus==0 && numPiPlus==0){ switch(numGamma){ case 2: KpN2GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } } // nProton==0 && nAntiProton==0; else if (numProton==1 && numAntiProton==1 && numPiPlus==1){ switch(numGamma){ case 0: PAntiPNPipAnalysis(beam_kf,final_kf,weight); break; case 2: PAntiPNPip2GAnalysis(beam_kf,final_kf,weight); break; default: break; } } else if (numProton==2 && numAntiProton==0 && numPiMinus==1 && numPiPlus==0 && numKPlus==0 && numKMinus==0){ switch(numGamma){ case 2: TwoProtonPim2GAnalysis(beam_kf,final_kf,weight); break; default: break; } } } else if (numNeutron==2&&numProton==1&&numGamma==0){ if (numPiPlus==0 && numPiMinus==0 && numKPlus==0 && numKMinus==0){ TwoNeutronVertex->Fill(vertex.z(),vertex.Perp(),weight); TwoNeutronAnalysis(beam_kf,final_kf,weight); } if (numPiPlus==1 && numPiMinus==1 && numKPlus==0 && numKMinus==0){ TwoNeutronPipPimAnalysis(beam_kf,final_kf,weight); } } } // CL cut } // loop over beam photons delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; } // require 1 neutron japp->Unlock("myneutronlock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_neutron::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_neutron::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } // Run the kinematic fitter requiring energy and momentum conservation void JEventProcessor_neutron::DoKinematicFit(double t0_rf,const DVector3 &vertex, const DBeamPhoton *beamphoton, vector&neutronHyps, map >&chargedParticles, vector&gammas, DKinFitUtils_GlueX *dKinFitUtils, DKinFitter *dKinFitter, bool isMissing) const { set> InitialParticles, FinalParticles; dKinFitter->Reset_NewFit(); shared_ptrmyBeam=dKinFitUtils->Make_BeamParticle(beamphoton); shared_ptrmyTarget=dKinFitUtils->Make_TargetParticle(Proton); InitialParticles.insert(myBeam); InitialParticles.insert(myTarget); // Vertex info DLorentzVector vert4; vert4.SetVect(vertex); vert4.SetT(t0_rf); // final state neutron(s) for (unsigned int k=0;kmyNeutron=dKinFitUtils->Make_DetectedParticle(2112,0,ParticleMass(Neutron),vert4,neutronHyps[k]->momentum(),0.,neutronHyps[k]->errorMatrix()); FinalParticles.insert(myNeutron); } // Missing pion if (isMissing){ shared_ptrmyRecoil=dKinFitUtils->Make_MissingParticle(PiPlus); FinalParticles.insert(myRecoil); } else { // Charged particles Particle_t particle_list[6]={PiPlus,PiMinus,KPlus,KMinus,Proton,AntiProton}; for (unsigned int k=0;k<6;k++){ vectormyParticles=chargedParticles[particle_list[k]]; for (unsigned int m=0;mmyParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]); FinalParticles.insert(myParticle); } } } // Final state photons for (unsigned int k=0;kmyGamma=dKinFitUtils->Make_DetectedParticle(22,0,0.,vert4,gammaHyp->momentum(),0.,gammaHyp->errorMatrix()); FinalParticles.insert(myGamma); } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); } // Use the PIDFOM scheme to sort tracks according to particle type void JEventProcessor_neutron::FillChargedParticleVectors(vector&tracks,map>&chargedParticles) const { Particle_t particle_list[6]={KPlus,KMinus,PiPlus,PiMinus,Proton,AntiProton}; for (unsigned int j=0;j >probabilities; for (unsigned int i=0;i<6;i++){ if ((hyp=tracks[j]->Get_Hypothesis(particle_list[i]))!=NULL){ probabilities.push_back(make_pair(hyp->Get_FOM(),particle_list[i])); } } if (probabilities.size()>0){ sort(probabilities.begin(),probabilities.end(),SortParticleProbability); Particle_t myParticle=probabilities[0].second; hyp=tracks[j]->Get_Hypothesis(myParticle); chargedParticles[myParticle].push_back(hyp->Get_TrackTimeBased()); } } } // Sort neutral particles into lists of photons and neutrons void JEventProcessor_neutron::FillNeutralParticleVectors(double &unused_energy, double t0_rf, const DVector3 &vertex, vector&neutrals, map>&neutralParticles) const { for (unsigned int i=0;iGet_Hypothesis(Gamma); // Get shower info correspnding to this hypothesis const DNeutralShower *shower=gamHyp->Get_NeutralShower(); // Compute beta double ct=29.98*(gamHyp->t1()-t0_rf); DVector3 diff; if (shower->dDetectorSystem==SYS_FCAL){ const DFCALShower *fcalshower=dynamic_cast(shower->dBCALFCALShower); DVector3 pos=fcalshower->getPosition(); diff=pos-vertex; // exclude inner blocks where EM background is most likely to be seen pos-=mFCALCenter; if (fabs(pos.X())getEnergy(); continue; } } else if (shower->dDetectorSystem==SYS_BCAL){ const DBCALShower *bcalshower=dynamic_cast(shower->dBCALFCALShower); DVector3 pos(bcalshower->x,bcalshower->y,bcalshower->z); diff=pos-vertex; } double beta=diff.Mag()/ct; NeutralBeta->Fill(beta); if (betaGet_Hypothesis(Neutron); if (neutHyp!=NULL){ neutralParticles[Neutron].push_back(neutHyp); } } else{ double tdiff=gamHyp->time()-t0_rf; GammaRFTimeDiff->Fill(tdiff); if (fabs(tdiff)<1.){ bool got_photon=true; double E=gamHyp->lorentzMomentum().E(); if (shower->dDetectorSystem==SYS_FCAL){ if (shower->dQuality(shower->dBCALFCALShower); DVector3 pos=fcalshower->getPosition(); pos-=mFCALCenter; if (pos.Perp()>FCAL_RADIAL_CUT) got_photon=false; if (fabs(pos.X())dDetectorSystem==SYS_BCAL){ if (E(shower->dBCALFCALShower); if (bcalshower->z>BCAL_Z_CUT) got_photon=false; } if (got_photon){ neutralParticles[Gamma].push_back(gamHyp); } else unused_energy+=E; } } } } void JEventProcessor_neutron::GetKF4vectors(DKinFitter *dKinFitter, DLorentzVector &beam_kf, DLorentzVector &missing_kf, map>&final_kf) const{ set>myParticles=dKinFitter->Get_KinFitParticles(); set>::iterator locParticleIterator=myParticles.begin(); for(; locParticleIterator != myParticles.end(); ++locParticleIterator){ if ((*locParticleIterator)->Get_KinFitParticleType()==d_BeamParticle){ beam_kf=(*locParticleIterator)->Get_P4(); } else if ((*locParticleIterator)->Get_KinFitParticleType()==d_MissingParticle){ missing_kf=(*locParticleIterator)->Get_P4(); } else if ((*locParticleIterator)->Get_KinFitParticleType()==d_DetectedParticle){ Particle_t myPID=PDGtoPType((*locParticleIterator)->Get_PID()); final_kf[myPID].push_back((*locParticleIterator)->Get_P4()); } } } void JEventProcessor_neutron::FillGJHisto(double weight, const DLorentzVector &beam, const DLorentzVector &meson, const DLorentzVector &analyzer, TH2F *histo) const { DVector3 beta=(1./meson.E())*meson.Vect(); DLorentzVector mybeam=beam; mybeam.Boost(-beta); DLorentzVector myanalyzer=analyzer; myanalyzer.Boost(-beta); DVector3 beamdir=mybeam.Vect(); beamdir.SetMag(1.); DVector3 adir=myanalyzer.Vect(); adir.SetMag(1.); histo->Fill(meson.M(),adir.Dot(beamdir),weight); } void JEventProcessor_neutron::MakeMesonPairs(vector&gammas_kf, double weight,TH2F *histo){ // Clear the current set of vectors particle_pairs[Pi0Pi0_].clear(); particle_pairs[Pi0Eta_].clear(); particle_pairs[Pi0EtaPrime_].clear(); particle_pairs[EtaEta_].clear(); // Find pairs of mesons for (unsigned int i=0;iFill(pair1_mass,pair2_mass,weight); bool eta_1_cut=ETA_CUT(pair1_mass,ETA_CUT_VALUE); bool eta_2_cut=ETA_CUT(pair2_mass,ETA_CUT_VALUE); bool pi0_1_cut=PI0_CUT(pair1_mass,PI0_CUT_VALUE); bool pi0_2_cut=PI0_CUT(pair2_mass,PI0_CUT_VALUE); if (pi0_1_cut && pi0_2_cut){ particle_pairs[Pi0Pi0_].push_back(make_pair(pair1,pair2)); } else if (eta_1_cut && eta_2_cut){ particle_pairs[EtaEta_].push_back(make_pair(pair1,pair2)); } else if (eta_1_cut && pi0_2_cut){ particle_pairs[Pi0Eta_].push_back(make_pair(pair2,pair1)); } else if (eta_2_cut && pi0_1_cut){ particle_pairs[Pi0Eta_].push_back(make_pair(pair1,pair2)); } else if (pi0_1_cut && ETAPRIME_CUT(pair2_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair1,pair2)); } else if (pi0_2_cut && ETAPRIME_CUT(pair1_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair2,pair1)); } } } void JEventProcessor_neutron::MakeTriplets() { vector>pairs; pairs.push_back(make_pair(0,1)); pairs.push_back(make_pair(2,3)); pairs.push_back(make_pair(4,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,1)); pairs.push_back(make_pair(2,4)); pairs.push_back(make_pair(3,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,1)); pairs.push_back(make_pair(2,5)); pairs.push_back(make_pair(3,4)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,2)); pairs.push_back(make_pair(1,3)); pairs.push_back(make_pair(4,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,2)); pairs.push_back(make_pair(1,4)); pairs.push_back(make_pair(5,3)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,2)); pairs.push_back(make_pair(1,5)); pairs.push_back(make_pair(3,4)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,3)); pairs.push_back(make_pair(2,1)); pairs.push_back(make_pair(4,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,3)); pairs.push_back(make_pair(2,4)); pairs.push_back(make_pair(1,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,3)); pairs.push_back(make_pair(2,5)); pairs.push_back(make_pair(4,1)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,4)); pairs.push_back(make_pair(1,3)); pairs.push_back(make_pair(2,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,4)); pairs.push_back(make_pair(1,5)); pairs.push_back(make_pair(2,3)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,4)); pairs.push_back(make_pair(1,2)); pairs.push_back(make_pair(3,5)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,5)); pairs.push_back(make_pair(1,3)); pairs.push_back(make_pair(4,2)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,5)); pairs.push_back(make_pair(1,4)); pairs.push_back(make_pair(3,2)); triplets.push_back(pairs); pairs.clear(); pairs.push_back(make_pair(0,5)); pairs.push_back(make_pair(1,2)); pairs.push_back(make_pair(4,3)); triplets.push_back(pairs); } //------------------------------------------------------------------------------ // Neutron Pi+ Ngamma //------------------------------------------------------------------------------ void JEventProcessor_neutron::MissingPionAnalysis(const DTrackTimeBased *pion_track, const DLorentzVector &missing_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; double twogam_mass=twogam_kf.M(); TwoGammaMass_missingPip->Fill(twogam_mass,weight); if (PI0_CUT(twogam_mass,PI0_CUT_VALUE)){ const DVector3 measured_momentum=pion_track->momentum(); DVector3 missing_momentum=missing_kf.Vect(); double missing_p=missing_momentum.Mag(); double dP=measured_momentum.Mag()-missing_p; double dtheta=180./M_PI*(measured_momentum.Theta()-missing_momentum.Theta()); double dphi=measured_momentum.Phi()-missing_momentum.Phi(); if (dphi<-M_PI) dphi+=2.*M_PI; if (dphi>M_PI) dphi-=2.*M_PI; dphi*=180./M_PI; DeltaPOverP_vs_theta->Fill(180./M_PI*missing_momentum.Theta(), dP/missing_p,weight); DeltaPOverP_vs_P->Fill(missing_p,dP/missing_p,weight); DeltaTheta_vs_P->Fill(missing_p,dtheta,weight); DeltaPhi_vs_P->Fill(missing_p,dphi,weight); } } void JEventProcessor_neutron::NeutronTwoGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector pip2gam_kf=final_kf[PiPlus][0]+twogam_kf; DLorentzVector n2gam_kf=final_kf[Neutron][0]+twogam_kf; double twogam_mass=twogam_kf.M(); double t_mandelstam=(beam_kf-pip2gam_kf).M2(); TwoGammaMass_NPip->Fill(twogam_mass,weight); if (PI0_CUT(twogam_mass,PI0_CUT_VALUE)){ PipPi0_NPi0_Dalitz->Fill(pip2gam_kf.M2(),n2gam_kf.M2(),weight); PipPi0Mass_vs_t->Fill(-t_mandelstam,pip2gam_kf.M(),weight); PipPi0Mass_vs_E->Fill(beam_kf.E(),pip2gam_kf.M(),weight); FillGJHisto(weight,beam_kf,pip2gam_kf,twogam_kf,Pi0CosThetaGJ_n); } if (PI0_BG_CUT(twogam_mass,NUM_SIGMA_BG,PI0_CUT_VALUE) ||PI0_BG_CUT(twogam_mass,-NUM_SIGMA_BG,PI0_CUT_VALUE)){ FillGJHisto(weight,beam_kf,pip2gam_kf,twogam_kf,Pi0CosThetaGJ_n_bg); } if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE)){ PipEta_NEta_Dalitz->Fill(pip2gam_kf.M2(),n2gam_kf.M2(),weight); PipEtaMass_vs_t->Fill(-t_mandelstam,pip2gam_kf.M(),weight); PipEtaMass_vs_E->Fill(beam_kf.E(),pip2gam_kf.M(),weight); FillGJHisto(weight,beam_kf,pip2gam_kf,twogam_kf,EtaCosThetaGJ_n); } if (ETA_BG_CUT(twogam_mass,NUM_SIGMA_BG,ETA_CUT_VALUE) ||ETA_BG_CUT(twogam_mass,-NUM_SIGMA_BG,ETA_CUT_VALUE)){ FillGJHisto(weight,beam_kf,pip2gam_kf,twogam_kf,EtaCosThetaGJ_n_bg); } if (ETAPRIME_CUT(twogam_kf.M(),ETAPRIME_CUT_VALUE)){ PipEtaPrimeMass_vs_t->Fill(-t_mandelstam,pip2gam_kf.M(),weight); PipEtaPrimeMass_vs_E->Fill(beam_kf.E(),pip2gam_kf.M(),weight); FillGJHisto(weight,beam_kf,pip2gam_kf,twogam_kf,EtaPrimeCosThetaGJ_n); } } void JEventProcessor_neutron::NeutronThreeGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector threegam_kf=final_kf[Gamma][0]+final_kf[Gamma][1] +final_kf[Gamma][2]; double threegam_mass=threegam_kf.M(); DLorentzVector piplus_kf=final_kf[PiPlus][0]; DLorentzVector neutron_kf=final_kf[Neutron][0]; DLorentzVector npip=piplus_kf+neutron_kf; DLorentzVector pip3g=piplus_kf+threegam_kf; double pip3g_mass=pip3g.M(); double npip_mass=npip.M(); ThreeGam_vs_NPip->Fill(npip_mass,threegam_mass,weight); if (OMEGA_CUT(threegam_mass)){ OmegaPip_n3g->Fill(pip3g_mass,weight); } } void JEventProcessor_neutron::NeutronFourGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight){ DLorentzVector fourgam_kf=final_kf[Gamma][0]+final_kf[Gamma][1] +final_kf[Gamma][2]+final_kf[Gamma][3]; double mass=fourgam_kf.M(); FourGammaMass_NPip->Fill(mass,weight); MakeMesonPairs(final_kf[Gamma],weight,FourGamma2d_NPip); if (particle_pairs[Pi0Pi0_].size()>0){ Pi0Pi0_NPip->Fill(mass,weight); } if (particle_pairs[EtaEta_].size()>0){ EtaEta_NPip->Fill(mass,weight); } if (particle_pairs[Pi0Eta_].size()>0){ Pi0Eta_NPip->Fill(mass,weight); } } void JEventProcessor_neutron::NeutronSixGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector npip_kf=final_kf[Neutron][0]+final_kf[PiPlus][0]; DLorentzVector sixgam_kf=final_kf[Gamma][0]+final_kf[Gamma][1] +final_kf[Gamma][2]+final_kf[Gamma][3] +final_kf[Gamma][4]+final_kf[Gamma][5]; double sixgammamass=sixgam_kf.M(); DLorentzVector pip6g=final_kf[PiPlus][0]+sixgam_kf; SixGammaMass_NPi->Fill(sixgammamass,weight); Pip6G_N->Fill(pip6g.M(),weight); for (unsigned int i=0;iFill(pair1.M(),pair2.M(),pair3.M(),weight); vectorpi0s; vectoretas; vectorunused_pair; for (unsigned int j=0;jFill(sixgammamass,weight); if (ETA_CUT(sixgammamass,ETA_CUT_VALUE)){ PipEta_N6g->Fill(pip6g.M(),weight); } } if (pi0s.size()==2){ if (etas.size()==1){ TwoPi0EtaMass_NPi->Fill(sixgammamass,weight); } else{ DLorentzVector gamma1=final_kf[Gamma][triplets[i][unused_pair[0]].first]; DLorentzVector gamma2=final_kf[Gamma][triplets[i][unused_pair[0]].second]; DLorentzVector pi0_1_g_1=pi0s[0]+gamma1; DLorentzVector pi0_2_g_2=pi0s[1]+gamma2; Pi0gVsPi0g_NPi->Fill(pi0_1_g_1.M(),pi0_2_g_2.M(),weight); DLorentzVector pi0_1_g_2=pi0s[0]+gamma2; DLorentzVector pi0_2_g_1=pi0s[1]+gamma1; Pi0gVsPi0g_NPi->Fill(pi0_1_g_2.M(),pi0_2_g_1.M(),weight); } } if (pi0s.size()==1 && etas.size()==2){ DLorentzVector eta1pi0=etas[0]+pi0s[0]; DLorentzVector eta2pi0=etas[1]+pi0s[0]; TwoEtaPi0Dalitz_NPi->Fill(eta1pi0.M2(),eta2pi0.M2(),weight); TwoEtaPi0Mass_NPi->Fill(sixgammamass,weight); } if (etas.size()==3){ ThreeEtaMass_NPi->Fill(sixgammamass,weight); } } } //------------------------------------------------------------------------------ // Neutron 2pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_neutron::Neutron2PipPimAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ const DLorentzVector neutron=final_kf[Neutron][0]; const DLorentzVector pim=final_kf[PiMinus][0]; const DLorentzVector pip1=final_kf[PiPlus][0]; const DLorentzVector pip2=final_kf[PiPlus][1]; DLorentzVector npip1_kf=neutron+pip1; DLorentzVector npip2_kf=neutron+pip2; DLorentzVector pip1pim=pip1+pim; DLorentzVector pip2pim=pip2+pim; PipPim_vs_NPip->Fill(npip1_kf.M(),pip2pim.M(),weight); PipPim_vs_NPip->Fill(npip2_kf.M(),pip1pim.M(),weight); NeutronPipMass_2pip->Fill(npip1_kf.M(),weight); NeutronPipMass_2pip->Fill(npip2_kf.M(),weight); } void JEventProcessor_neutron::Neutron2PipPim1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ const DLorentzVector neutron=final_kf[Neutron][0]; const DLorentzVector pim=final_kf[PiMinus][0]; const DLorentzVector pip1=final_kf[PiPlus][0]; const DLorentzVector pip2=final_kf[PiPlus][1]; const DLorentzVector gamma=final_kf[Gamma][0]; DLorentzVector npip1_kf=neutron+pip1; DLorentzVector npip2_kf=neutron+pip2; DLorentzVector pip1pim=pip1+pim; DLorentzVector pip2pim=pip2+pim; DLorentzVector pip1pimg=pip1pim+gamma; DLorentzVector pip2pimg=pip2pim+gamma; PipPim_vs_NPip_1G->Fill(npip1_kf.M(),pip2pim.M(),weight); PipPim_vs_NPip_1G->Fill(npip2_kf.M(),pip1pim.M(),weight); NeutronPipMass_2pip1G->Fill(npip1_kf.M(),weight); NeutronPipMass_2pip1G->Fill(npip2_kf.M(),weight); PipPim1G_N->Fill(pip1pimg.M(),weight); PipPim1G_N->Fill(pip2pimg.M(),weight); } void JEventProcessor_neutron::Neutron2PipPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector npip1_kf=final_kf[Neutron][0]+final_kf[PiPlus][0]; DLorentzVector npip2_kf=final_kf[Neutron][0]+final_kf[PiPlus][1]; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector pimpip1=final_kf[PiMinus][0]+final_kf[PiPlus][0]; DLorentzVector pimpip2=final_kf[PiMinus][0]+final_kf[PiPlus][1]; DLorentzVector pippim2g_1=pimpip1+twogam; DLorentzVector pippim2g_2=pimpip2+twogam; DLorentzVector pim2pip=pimpip1+final_kf[PiPlus][1]; DLorentzVector pim2pip2g=pim2pip+twogam; DLorentzVector n2g=final_kf[Neutron][0]+twogam; if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ NPi0_vs_2PipPim->Fill(pim2pip.M(),n2g.M(),weight); TwoPipPimPi0_N->Fill(pim2pip2g.M(),weight); PipPimPi0_n->Fill(pippim2g_1.M(),weight); PipPimPi0_n->Fill(pippim2g_2.M(),weight); if (ETA_CUT(pippim2g_1.M(),ETA_CUT_VALUE)){ EtaPip_N->Fill(pim2pip2g.M(),weight); } if (ETA_CUT(pippim2g_2.M(),ETA_CUT_VALUE)){ EtaPip_N->Fill(pim2pip2g.M(),weight); } if (OMEGA_CUT(pippim2g_1.M())){ OmegaPip_N2g->Fill(pim2pip2g.M(),weight); } if (OMEGA_CUT(pippim2g_2.M())){ OmegaPip_N2g->Fill(pim2pip2g.M(),weight); } PipPim_vs_NPip_with_Pi0->Fill(npip1_kf.M(),pimpip2.M(),weight); PipPim_vs_NPip_with_Pi0->Fill(npip2_kf.M(),pimpip1.M(),weight); if (KSHORT_CUT(pimpip2.M())&&SIGMAP_CUT(npip1_kf.M())){ DLorentzVector SigmaPlusPi0=npip1_kf+twogam; DLorentzVector SigmaPlusKShort=pimpip2+npip1_kf; SigmaPlusPi0_SigmaPlusK0_Dalitz->Fill(SigmaPlusPi0.M2(), SigmaPlusKShort.M2(),weight); } if (KSHORT_CUT(pimpip1.M())&&SIGMAP_CUT(npip2_kf.M())){ DLorentzVector SigmaPlusPi0=npip2_kf+twogam; DLorentzVector SigmaPlusKShort=pimpip1+npip2_kf; SigmaPlusPi0_SigmaPlusK0_Dalitz->Fill(SigmaPlusPi0.M2(), SigmaPlusKShort.M2(),weight); } } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ NEta_vs_2PipPim->Fill(pim2pip.M(),n2g.M(),weight); TwoPipPimEta_N->Fill(pim2pip2g.M(),weight); PipPimEta_n->Fill(pippim2g_1.M(),weight); PipPimEta_n->Fill(pippim2g_2.M(),weight); if (ETAPRIME_CUT(pippim2g_1.M(),ETAPRIME_CUT_VALUE)){ EtaPrimePip_N->Fill(pim2pip2g.M(),weight); if (unused_energy<0.01){ EtaPrimePip_N_NU->Fill(pim2pip2g.M(),weight); } } if (ETAPRIME_CUT(pippim2g_2.M(),ETAPRIME_CUT_VALUE)){ EtaPrimePip_N->Fill(pim2pip2g.M(),weight); if (unused_energy<0.01){ EtaPrimePip_N_NU->Fill(pim2pip2g.M(),weight); } } PipPim_vs_NPip_with_Eta->Fill(npip1_kf.M(),pimpip2.M(),weight); PipPim_vs_NPip_with_Eta->Fill(npip2_kf.M(),pimpip1.M(),weight); if (KSHORT_CUT(pimpip2.M())&&SIGMAP_CUT(npip1_kf.M())){ DLorentzVector SigmaPlusEta=npip1_kf+twogam; DLorentzVector SigmaPlusKShort=pimpip2+npip1_kf; SigmaPlusEta_SigmaPlusK0_Dalitz->Fill(SigmaPlusEta.M2(), SigmaPlusKShort.M2(),weight); } if (KSHORT_CUT(pimpip1.M())&&SIGMAP_CUT(npip2_kf.M())){ DLorentzVector SigmaPlusEta=npip2_kf+twogam; DLorentzVector SigmaPlusKShort=pimpip1+npip2_kf; SigmaPlusEta_SigmaPlusK0_Dalitz->Fill(SigmaPlusEta.M2(), SigmaPlusKShort.M2(),weight); } } TwoGammaMass_n2pip->Fill(twogam.M(),weight); NeutronPipMass_2pip2g->Fill(npip1_kf.M(),weight); NeutronPipMass_2pip2g->Fill(npip2_kf.M(),weight); if (unused_energy<0.01){ TwoGammaMass_n2pip_NU->Fill(twogam.M(),weight); } } void JEventProcessor_neutron::Neutron2PipPim4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight){ DLorentzVector npip1_kf=final_kf[Neutron][0]+final_kf[PiPlus][0]; DLorentzVector npip2_kf=final_kf[Neutron][0]+final_kf[PiPlus][1]; DLorentzVector pimpip1=final_kf[PiMinus][0]+final_kf[PiPlus][0]; DLorentzVector pimpip2=final_kf[PiMinus][0]+final_kf[PiPlus][1]; DLorentzVector pim2pip=pimpip1+final_kf[PiPlus][1]; DLorentzVector fourgam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]+final_kf[Gamma][2]+final_kf[Gamma][3]; double fourgam_mass=fourgam_kf.M(); DLorentzVector pimpip1_4g=pimpip1+fourgam_kf; DLorentzVector pimpip2_4g=pimpip2+fourgam_kf; FourGammaMass_N3Pi->Fill(fourgam_mass,weight); MakeMesonPairs(final_kf[Gamma],weight,FourGamma2d_N3Pi); if (particle_pairs[Pi0Pi0_].size()>0){ PipPim_vs_NPip_with_2Pi0->Fill(npip1_kf.M(),pimpip2.M(),weight); PipPim_vs_NPip_with_2Pi0->Fill(npip2_kf.M(),pimpip1.M(),weight); TwoPi0_N3Pi->Fill(fourgam_mass,weight); for (unsigned int i=0;iFill(pimpip_1_pi0_1_mass,weight); if (ETA_CUT(pimpip_1_pi0_1_mass,ETA_CUT_VALUE)){ EtaPi0_N_eta_3pi->Fill(pimpip1_4g.M(),weight); } PipPimPi0_N_with_Pi0->Fill(pimpip_1_pi0_2_mass,weight); if (ETA_CUT(pimpip_1_pi0_2_mass,ETA_CUT_VALUE)){ EtaPi0_N_eta_3pi->Fill(pimpip1_4g.M(),weight); } PipPimPi0_N_with_Pi0->Fill(pimpip_2_pi0_2_mass,weight); if (ETA_CUT(pimpip_2_pi0_2_mass,ETA_CUT_VALUE)){ EtaPi0_N_eta_3pi->Fill(pimpip2_4g.M(),weight); } PipPimPi0_N_with_Pi0->Fill(pimpip_2_pi0_1_mass,weight); if (ETA_CUT(pimpip_2_pi0_1_mass,ETA_CUT_VALUE)){ EtaPi0_N_eta_3pi->Fill(pimpip2_4g.M(),weight); } } } if (particle_pairs[Pi0Eta_].size()>0){ PipPim_vs_NPip_with_EtaPi0->Fill(npip1_kf.M(),pimpip2.M(),weight); PipPim_vs_NPip_with_EtaPi0->Fill(npip2_kf.M(),pimpip1.M(),weight); EtaPi0_N3Pi->Fill(fourgam_mass,weight); for (unsigned int i=0;iFill(pimpip_1_pi0.M(),weight); DLorentzVector pimpip_2_pi0=pimpip2+pi0; PipPimPi0_N_with_Eta->Fill(pimpip_2_pi0.M(),weight); DLorentzVector pimpip_1_eta=pimpip1+eta; PipPimEta_N_with_Pi0->Fill(pimpip_1_eta.M(),weight); DLorentzVector pimpip_2_eta=pimpip2+eta; PipPimEta_N_with_Pi0->Fill(pimpip_2_eta.M(),weight); } } } //------------------------------------------------------------------------------ // Neutron K+ pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_neutron::NeutronKpPipPimAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector pippim=final_kf[PiPlus][0]+final_kf[PiMinus][0]; double pippim_mass=pippim.M(); DLorentzVector npippim=final_kf[Neutron][0]+pippim; double npippim_mass=npippim.M(); DLorentzVector kpippim=final_kf[KPlus][0]+pippim; DLorentzVector npip=final_kf[Neutron][0]+final_kf[PiPlus][0]; double npip_mass=npip.M(); DLorentzVector npim=final_kf[Neutron][0]+final_kf[PiMinus][0]; double npim_mass=npim.M(); DLorentzVector nKp=final_kf[Neutron][0]+final_kf[KPlus][0]; double nKp_mass=nKp.M(); DLorentzVector KpPip=final_kf[PiPlus][0]+final_kf[KPlus][0]; double KpPip_mass=KpPip.M(); DLorentzVector KpPim=final_kf[PiMinus][0]+final_kf[KPlus][0]; double KpPim_mass=KpPim.M(); PipPim_vs_nKp->Fill(nKp_mass,pippim_mass,weight); if (KSHORT_CUT(pippim_mass)){ NKs->Fill(npippim_mass,weight); } NPim_vs_KpPip->Fill(KpPip_mass,npim_mass,weight); if (SIGMAM_CUT(npim_mass)){ SigmamPip->Fill(npippim_mass,weight); } NPip_vs_KpPim->Fill(KpPim_mass,npip_mass,weight); if (SIGMAP_CUT(npip_mass)){ SigmapPim->Fill(npippim_mass,weight); } } void JEventProcessor_neutron::NeutronKpPipPim1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector neutron_pim=final_kf[Neutron][0]+final_kf[PiMinus][0]; DLorentzVector neutron_pip=final_kf[Neutron][0]+final_kf[PiPlus][0]; NPipSqVsNPimSq_1G->Fill(neutron_pim.M2(),neutron_pip.M2(),weight); if (SIGMAM_CUT(neutron_pim.M())){ DLorentzVector sigm_gam=neutron_pim+final_kf[Gamma][0]; SigmaMinusGamma->Fill(sigm_gam.M(),weight); } if (SIGMAP_CUT(neutron_pip.M())){ DLorentzVector sigp_gam=neutron_pip+final_kf[Gamma][0]; SigmaPlusGamma->Fill(sigp_gam.M(),weight); } } void JEventProcessor_neutron::NeutronKpPipPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector neutron_pim=final_kf[Neutron][0]+final_kf[PiMinus][0]; DLorentzVector neutron_pip=final_kf[Neutron][0]+final_kf[PiPlus][0]; DLorentzVector neutron_pippim=neutron_pip+final_kf[PiMinus][0]; DLorentzVector k2g=twogam+final_kf[KPlus][0]; DLorentzVector npippim2g=neutron_pippim+twogam; NPipSqVsNPimSq_2G->Fill(neutron_pim.M2(),neutron_pip.M2(),weight); TwoGamma_with_KNPipPim->Fill(twogam.M(),weight); if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ if (SIGMAM_CUT(neutron_pim.M())){ DLorentzVector sigm_pi0=neutron_pim+twogam; DLorentzVector kp_pip=final_kf[KPlus][0]+final_kf[PiPlus][0]; SigmaMinusPip_vs_KPi0->Fill(k2g.M(),neutron_pippim.M(),weight); SigmaMinusPi0_vs_KPip->Fill(kp_pip.M(),sigm_pi0.M(),weight); SigmaMinusPi0_SigmaMinusPip_Dalitz->Fill(sigm_pi0.M2(),neutron_pippim.M2(), weight); SigmamPipPi0->Fill(npippim2g.M(),weight); } if (SIGMAP_CUT(neutron_pip.M())){ DLorentzVector sigp_pi0=neutron_pip+twogam; DLorentzVector kp_pim=final_kf[KPlus][0]+final_kf[PiMinus][0]; SigmaPlusPi0_vs_KPim->Fill(kp_pim.M(),sigp_pi0.M(),weight); SigmaPlusPim_vs_KPi0->Fill(k2g.M(),neutron_pippim.M(),weight); SigmaPlusPi0_SigmaPlusPim_Dalitz->Fill(sigp_pi0.M2(),neutron_pippim.M2(), weight); SigmapPimPi0->Fill(npippim2g.M(),weight); } } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ if (SIGMAM_CUT(neutron_pim.M())){ DLorentzVector sigm_eta=neutron_pim+twogam; DLorentzVector kp_pip=final_kf[KPlus][0]+final_kf[PiPlus][0]; SigmaMinusPip_vs_KEta->Fill(k2g.M(),neutron_pippim.M(),weight); SigmaMinusEta_vs_KPip->Fill(kp_pip.M(),sigm_eta.M(),weight); SigmaMinusEta_SigmaMinusPip_Dalitz->Fill(sigm_eta.M2(),neutron_pippim.M2(), weight); SigmamPipEta->Fill(npippim2g.M(),weight); } if (SIGMAP_CUT(neutron_pip.M())){ DLorentzVector sigp_eta=neutron_pip+twogam; DLorentzVector kp_pim=final_kf[KPlus][0]+final_kf[PiMinus][0]; SigmaPlusEta_vs_KPim->Fill(kp_pim.M(),sigp_eta.M(),weight); SigmaPlusPim_vs_KEta->Fill(k2g.M(),neutron_pippim.M(),weight); SigmaPlusEta_SigmaPlusPim_Dalitz->Fill(sigp_eta.M2(),neutron_pippim.M2(), weight); SigmapPimEta->Fill(npippim2g.M(),weight); } } } //--------------------------------------------------------------------------- // K+ n 2gamma //--------------------------------------------------------------------------- void JEventProcessor_neutron::KpN2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf,double weight) const{ DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector neutron_pi0=final_kf[Neutron][0]+twogam; TwoGammaMass_NK->Fill(twogam.M(),weight); N2gMass->Fill(neutron_pi0.M(),weight); if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ NPi0Mass->Fill(neutron_pi0.M(),weight); } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ NEtaMass->Fill(neutron_pi0.M(),weight); } } //--------------------------------------------------------------------------- // Events with 2K+ //--------------------------------------------------------------------------- void JEventProcessor_neutron::TwoKpKmN2GAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector n2g=twogam+final_kf[Neutron][0]; TwoGamma_2kpkmn->Fill(twogam.M(),weight); if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ NPi0_2kpkm->Fill(n2g.M(),weight); } } //------------------------------------------------------------------------------ // 2 neutron events //------------------------------------------------------------------------------ void JEventProcessor_neutron::TwoNeutronAnalysis(const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ const DLorentzVector neutron1=final_kf[Neutron][0]; const DLorentzVector neutron2=final_kf[Neutron][1]; const DLorentzVector nn=neutron1+neutron2; double t_mandelstam=(beam_kf-nn).M2(); const DLorentzVector proton=final_kf[Proton][0]; const DLorentzVector pn1=proton+neutron1; const DLorentzVector pn2=proton+neutron2; PNDalitz->Fill(pn1.M2(),pn2.M2(),weight); double p_p=proton.Vect().Mag(); double p_theta=proton.Vect().Theta()*180/M_PI; double n1_p=neutron1.Vect().Mag(); double n2_p=neutron2.Vect().Mag(); double n1_theta=neutron1.Vect().Theta()*180/M_PI; double n2_theta=neutron2.Vect().Theta()*180/M_PI; TwoNeutronMomentum2D->Fill(n1_p,n2_p,weight); NeutronPVsTheta->Fill(n1_theta,n1_p,weight); NeutronPVsTheta->Fill(n2_theta,n2_p,weight); ProtonPVsThetaWithN->Fill(p_theta,p_p,weight); TwoNeutronMass->Fill(nn.M(),weight); NantiN_vs_t->Fill(-t_mandelstam,nn.M(),weight); NantiN_vs_E->Fill(beam_kf.E(),nn.M(),weight); } void JEventProcessor_neutron::TwoNeutronPipPimAnalysis(const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector neutron1=final_kf[Neutron][0]; DLorentzVector neutron2=final_kf[Neutron][1]; DLorentzVector pim=final_kf[PiMinus][0]; DLorentzVector pip=final_kf[PiPlus][0]; DLorentzVector n1pip=neutron1+pip; DLorentzVector n2pip=neutron2+pip; DLorentzVector n1pim=neutron1+pim; DLorentzVector n2pim=neutron2+pim; DLorentzVector pippim=pip+pim; NeutronPip_vs_NeutronPim->Fill(n1pim.M(),n2pip.M(),weight); NeutronPip_vs_NeutronPim->Fill(n2pim.M(),n1pip.M(),weight); } //------------------------------------------------------------------------------ // p pbar n pi+ Ngamma events //----------------------------------------------------------------------------- void JEventProcessor_neutron::PAntiPNPipAnalysis(const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector neutron=final_kf[Neutron][0]; DLorentzVector proton=final_kf[Proton][0]; DLorentzVector antiproton=final_kf[AntiProton][0]; DLorentzVector pip=final_kf[PiPlus][0]; DLorentzVector pbar_n_pip=antiproton+neutron+pip; AntiPNPip->Fill(pbar_n_pip.M(),weight); } void JEventProcessor_neutron::PAntiPNPip2GAnalysis(const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector neutron=final_kf[Neutron][0]; DLorentzVector proton=final_kf[Proton][0]; DLorentzVector antiproton=final_kf[AntiProton][0]; DLorentzVector pip=final_kf[PiPlus][0]; DLorentzVector pbar_n_pip=antiproton+neutron+pip; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; AntiPNPip_vs_2g->Fill(twogam.M(),pbar_n_pip.M(),weight); if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ DLorentzVector antip_pi0=antiproton+twogam; DLorentzVector neutron_pip=neutron+pip; NPip_vs_AntiPPi0->Fill(antip_pi0.M(),neutron_pip.M(),weight); DLorentzVector antip_pip=antiproton+pip; DLorentzVector neutron_pi0=neutron+twogam; NPi0_vs_AntiPPiP->Fill(antip_pip.M(),neutron_pi0.M(),weight); } } void JEventProcessor_neutron::TwoProtonPim2GAnalysis(const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector neutron=final_kf[Neutron][0]; DLorentzVector proton1=final_kf[Proton][0]; DLorentzVector proton2=final_kf[Proton][1]; DLorentzVector pim=final_kf[PiMinus][0]; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; double twog_mass=twogam.M(); TwoGammaMass_with_2pnpim->Fill(twogam.M(),weight); DLorentzVector npim=neutron+pim; DLorentzVector n2g=neutron+twogam; DLorentzVector p1pim=proton1+pim; DLorentzVector p2pim=proton2+pim; DLorentzVector p1_2g=proton1+twogam; DLorentzVector p2_2g=proton2+twogam; if (PI0_CUT(twog_mass,PI0_CUT_VALUE)){ PPim_vs_AntiNPi0->Fill(p1pim.M(),n2g.M(),weight); PPim_vs_AntiNPi0->Fill(p2pim.M(),n2g.M(),weight); PPi0_vs_AntiNPim->Fill(p1_2g.M(),npim.M(),weight); PPi0_vs_AntiNPim->Fill(p2_2g.M(),npim.M(),weight); } }