// $Id$ // // File: JEventProcessor_survey.cc // Created: Mon Dec 6 14:28:07 EST 2021 // Creator: staylor (on Linux ifarm1901.jlab.org 3.10.0-1160.11.1.el7.x86_64 x86_64) // #include "JEventProcessor_survey.h" using namespace jana; // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_survey()); } } // "C" //------------------ // JEventProcessor_survey (Constructor) //------------------ JEventProcessor_survey::JEventProcessor_survey() { } //------------------ // ~JEventProcessor_survey (Destructor) //------------------ JEventProcessor_survey::~JEventProcessor_survey() { } //------------------ // init //------------------ jerror_t JEventProcessor_survey::init(void) { japp->CreateLock("mykinlock"); CL_CUT=0.01; gPARMS->SetDefaultParameter("SURVEY:CL_CUT",CL_CUT); 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); MIN_BEAM_E=0.0; gPARMS->SetDefaultParameter("SURVEY:MIN_BEAM_E",MIN_BEAM_E); MAX_BEAM_E=12.0; gPARMS->SetDefaultParameter("SURVEY:MAX_BEAM_E",MAX_BEAM_E); BETA_CUT=0.9; gPARMS->SetDefaultParameter("SURVEY:BETA_CUT",BETA_CUT); SPLIT_CUT=0.5; gPARMS->SetDefaultParameter("SURVEY:SPLIT_CUT",SPLIT_CUT); FCAL_POS_CUT=10.0; gPARMS->SetDefaultParameter("SURVEY:FCAL_POS_CUT",FCAL_POS_CUT); FCAL_RADIAL_CUT=104.0; gPARMS->SetDefaultParameter("SURVEY:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT); BCAL_Z_CUT=384.0; gPARMS->SetDefaultParameter("SURVEY:BCAL_Z_CUT",BCAL_Z_CUT); FCAL_THRESHOLD=0.1; gPARMS->SetDefaultParameter("SURVEY:FCAL_THRESHOLD",FCAL_THRESHOLD); BCAL_THRESHOLD=0.05; gPARMS->SetDefaultParameter("SURVEY:BCAL_THRESHOLD",BCAL_THRESHOLD); TOF_FCAL_CUT=10.; gPARMS->SetDefaultParameter("SURVEY:TOF_FCAL_CUT",TOF_FCAL_CUT); gDirectory->mkdir("MC")->cd(); { MC_beam_E_all = new TH1F("MC_beam_E_all","Photon energy",48,2.8,12.4); MC_beam_E = new TH1F("MC_beam_E","Photon energy",48,2.8,12.4); MC_meson_mass =new TH1F("MC_meson_mass","mX",1600,0,4.); MC_t=new TH1F("MC_t","t distribution",1000,0,2.); } gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); { NeutralBeta=new TH1F("NeutralBeta","#beta for neutrals",100,0,1.5); GammaRFTimeDiff=new TH1F("GammaRFTimeDiff","t(#gamma)-t(rf) at vertex", 1000,-10,10); BeamPhotonTimeDiff=new TH1F("BeamPhotonTimeDiff","t(tag)-t(rf) at vertex", 1000,-50,50); ExclusiveCL=new TH1F("ExclusiveCL",";CL",1000,0,1.); MissingProtonCL=new TH1F("MissingProtonCL",";CL",1000,0,1.); NFDCPseudos=new TH1F("NFDCPsuedos","Number of FDC pseudopoints", 1000,0.5,1000.5); TofFcaldt=new TH1F("TofFcaldt",";#Deltat [ns]",400,-20,20); TofFcaldxdy=new TH2F("TofFcaldxdy",";#Deltax[cm];#Deltay[cm]", 400,-20,20,400,-20,20); TOFdE_with_FCAL=new TH2F("TOFdE_with_FCAL",";#DeltaE1 [MeV];#DeltaE2 [MeV]", 300,0,30,300,0,30); } gDirectory->cd("../"); gDirectory->mkdir("survey")->cd(); gDirectory->mkdir("2gamma")->cd(); { TwoGammaMass=new TH1F("TwoGammaMass",";M(2#gamma) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("4gamma")->cd(); { FourGammaMass=new TH1F("FourGammaMass",";M(4#gamma) [GeV]",1600,0,4.); FourGamma_Insertcut=new TH1F("FourGamma_Insertcut",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d=new TH2F("FourGamma2d", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); FourGamma2d_NoMesonPairs=new TH2F("FourGamma2d_NoMesonPairs", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); FourGammaNoPairs=new TH1F("FourGammaNoPairs",";M(4#gamma) [GeV]",1600,0,4.); FourGammaNoPairs_FDCcut=new TH1F("FourGammaNoPairs_FDCcut",";M(4#gamma) [GeV]",1600,0,4.); FourGammaNoPairs_TOFcut=new TH1F("FourGammaNoPairs_TOFcut",";M(4#gamma) [GeV]",1600,0,4.); FourGammaNoPairs_Insertcut=new TH1F("FourGammaNoPairs_Insertcut",";M(4#gamma) [GeV]",1600,0,4.); Pi0Pi0=new TH1F("Pi0Pi0",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg=new TH1F("Pi0Pi0bg",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta=new TH1F("Pi0Eta",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg=new TH1F("Pi0Etabg",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta=new TH1F("EtaEta",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg=new TH1F("EtaEtabg",";M(2#eta) [GeV]",1600,0,4.); Pi0CosThetaGF_bg=new TH2F("Pi0CosThetaGF_bg",";M(2#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); Pi0CosThetaGF=new TH2F("Pi0CosThetaGF",";M(2#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); EtaCosThetaGF_bg=new TH2F("EtaCosThetaGF_bg",";M(#eta#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); EtaCosThetaGF=new TH2F("EtaCosThetaGF",";M(#eta#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); Pi02g_vs_Pi0g=new TH2F("Pi02g_vs_Pi0g", ";M(#pi^{0}#gamma) [GeV];M(#pi^{0}#gamma#gamma) [GeV]", 800,0,2.,1600,0,4.); Pi02g_vs_2g=new TH2F("Pi02g_vs_2g", ";M(2#gamma) [GeV];M(#pi^{0}#gamma#gamma) [GeV]", 400,0,1.,1600,0,4.); } 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); Pi0CosThetaGF_n_bg=new TH2F("Pi0CosThetaGF_n_bg",";M(2#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); Pi0CosThetaGF_n=new TH2F("Pi0CosThetaGF_n",";M(2#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); EtaCosThetaGF_n_bg=new TH2F("EtaCosThetaGF_n_bg",";M(#eta#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); EtaCosThetaGF_n=new TH2F("EtaCosThetaGF_n",";M(#eta#pi^{0}) [GeV];cos(#theta_{GF})",1600,0,4,100,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("6gamma")->cd(); { SixGammaMass=new TH1F("SixGammaMass",";M(6#gamma) [GeV]",1600,0,4.); } 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.); Pi0Pi0bg_NPip=new TH1F("Pi0Pi0bg_NPip",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_NPip=new TH1F("Pi0Eta_NPip",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_NPip=new TH1F("Pi0Etabg_NPip",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_NPip=new TH1F("EtaEta_NPip",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_NPip=new TH1F("EtaEtabg_NPip",";M(2#eta) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("PipPim2gamma")->cd(); { TwoGammaMass_PipPim=new TH1F("TwoGammaMass_PipPim",";M(2#gamma) [GeV]",400,0,1.); PipPimPi0=new TH1F("PipPimPi0",";M(#pi^{=}#pi^{-}#pi^{0}) [GeV]",1600,0,4.); PipPimEta=new TH1F("PipPimEta",";M(#pi^{=}#pi^{-}#eta) [GeV]",1600,0.5,4.5); } gDirectory->cd("../"); gDirectory->mkdir("PipPim4gamma")->cd(); { PipPimMass_4g=new TH1F("PipPim_4g",";M(#pi^{+}#pi^{-}) [GeV]",1600,0.,4.); FourGammaMass_PipPim=new TH1F("FourGammaMass_PipPim",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d_PipPim=new TH2F("FourGamma2d_PipPim", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); Pi0Pi0_PipPim=new TH1F("Pi0Pi0_PipPim",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg_PipPim=new TH1F("Pi0Pi0bg_PipPim",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_PipPim=new TH1F("Pi0Eta_PipPim",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_PipPim=new TH1F("Pi0Etabg_PipPim",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_PipPim=new TH1F("EtaEta_PipPim",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_PipPim=new TH1F("EtaEtabg_PipPim",";M(2#eta) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpKm2gamma")->cd(); { TwoGamma_vs_KpKm=new TH2F("TwoGamma_vs_KpKm", ";M(K^{+}K^{-}) [GeV];M(2#gamma) [GeV]", 1600,0.9,4.9,400,0,1.); EtaPhiDalitz=new TH2F("EtaPhiDalitz", ";M^{2}(p#eta) [GeV^{2}];M^{2}(#phi#eta) [GeV^{2}]", 1000,2.,20.,1000,2.,20.); } gDirectory->cd("../"); gDirectory->mkdir("KpKm4gamma")->cd(); { KpKmMass_4g=new TH1F("KpKm_4g",";M(K^{+}K^{-}) [GeV]",1600,0.9,4.9); FourGammaMass_KpKm=new TH1F("FourGammaMass_KpKm",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d_KpKm=new TH2F("FourGamma2d_KpKm", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); Pi0Pi0_KpKm=new TH1F("Pi0Pi0_KpKm",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg_KpKm=new TH1F("Pi0Pi0bg_KpKm",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_KpKm=new TH1F("Pi0Eta_KpKm",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_KpKm=new TH1F("Pi0Etabg_KpKm",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_KpKm=new TH1F("EtaEta_KpKm",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_KpKm=new TH1F("EtaEtabg_KpKm",";M(2#eta) [GeV]",1600,0,4.); Pi0Pi0_Phi=new TH1F("Pi0Pi0_Phi",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg_Phi=new TH1F("Pi0Pi0bg_Phi",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_Phi=new TH1F("Pi0Eta_Phi",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_Phi=new TH1F("Pi0Etabg_Phi",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_Phi=new TH1F("EtaEta_Phi",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_Phi=new TH1F("EtaEtabg_Phi",";M(2#eta) [GeV]",1600,0,4.); Pi0Pi0_Phi_bg=new TH1F("Pi0Pi0_Phi_bg",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg_Phi_bg=new TH1F("Pi0Pi0bg_Phi_bg",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_Phi_bg=new TH1F("Pi0Eta_Phi_bg",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_Phi_bg=new TH1F("Pi0Etabg_Phi_bg",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_Phi_bg=new TH1F("EtaEta_Phi_bg",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_Phi_bg=new TH1F("EtaEtabg_Phi_bg",";M(2#eta) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("TwoKpKm")->cd(); { } gDirectory->cd("../"); gDirectory->mkdir("TwoKpKm2gamma")->cd(); { } gDirectory->cd("../"); gDirectory->mkdir("KpKmPipPim2gamma")->cd(); { TwoGamma_vs_KpKm_with_PipPim=new TH2F("TwoGamma_vs_KpKm_with_PipPim", ";M(K^{+}K^{-}) [GeV];M(2#gamma) [GeV]", 1600,0.9,4.9,400,0,1.); PipPimPi0_KpKm=new TH1F("PipPimPi0_KpKm",";M(#pi^{=}#pi^{-}#pi^{0}) [GeV]",1600,0,4.); PipPimEta_KpKm=new TH1F("PipPimEta_KpKm",";M(#pi^{=}#pi^{-}#eta) [GeV]",1600,0.5,4.5); EtaCosThetaGF_with_Phi=new TH2F("EtaCosThetaGF_with_Phi",";M(#eta#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); OmegaCosThetaGF_with_Phi=new TH2F("OmegaCosThetaGF_with_Phi",";M(#omega#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); EtaPrimeCosThetaGF_with_Phi=new TH2F("EtaPrimeCosThetaGF_with_Phi",";M(#eta'#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); PipPimPi0_KpKm_bg=new TH1F("PipPimPi0_KpKm_bg",";M(#pi^{=}#pi^{-}#pi^{0}) [GeV]",1600,0,4.); PipPimEta_KpKm_bg=new TH1F("PipPimEta_KpKm_bg",";M(#pi^{=}#pi^{-}#eta) [GeV]",1600,0.5,4.5); EtaCosThetaGF_with_Phi_bg=new TH2F("EtaCosThetaGF_with_Phi_bg",";M(#eta#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); OmegaCosThetaGF_with_Phi_bg=new TH2F("OmegaCosThetaGF_with_Phi_bg",";M(#omega#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); EtaPrimeCosThetaGF_with_Phi_bg=new TH2F("EtaPrimeCosThetaGF_with_Phi_bg",";M(#eta'#phi}) [GeV];cos(#theta_{GF})",1600,1,5,100,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("KpKmPipPim4gamma")->cd(); { KpKmMass_4gPipPim=new TH1F("KpKm_4gPipPim",";M(K^{+}K^{-}) [GeV]",1600,0.9,4.9); FourGammaMass_KpKmPipPim=new TH1F("FourGammaMass_KpKmPipPim",";M(4#gamma) [GeV]",1600,0,4.); FourGamma2d_KpKmPipPim=new TH2F("FourGamma2d_KpKmPipPim", ";M(2#gamma,pair1)[GeV]; M(2#gamma,pair2) [GeV]", 500,0,1,500,0,1); Pi0Pi0_KpKmPipPim=new TH1F("Pi0Pi0_KpKmPipPim",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Pi0bg_KpKmPipPim=new TH1F("Pi0Pi0bg_KpKmPipPim",";M(2#pi^{0}) [GeV]",1600,0,4.); Pi0Eta_KpKmPipPim=new TH1F("Pi0Eta_KpKmPipPim",";M(#eta#pi^{0}) [GeV]",1600,0,4.); Pi0Etabg_KpKmPipPim=new TH1F("Pi0Etabg_KpKmPipPim",";M(#eta#pi^{0}) [GeV]",1600,0,4.); EtaEta_KpKmPipPim=new TH1F("EtaEta_KpKmPipPim",";M(2#eta) [GeV]",1600,0,4.); EtaEtabg_KpKmPipPim=new TH1F("EtaEtabg_KpKmPipPim",";M(2#eta) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("EpEm")->cd(); { EpEmMass=new TH1F("EpEmMass",";M(e^{+}e^{-}) [GeV]",1600,0.,4.); } gDirectory->cd("../"); gDirectory->mkdir("EpEm2gamma")->cd(); { TwoGamma_vs_EpEm=new TH2F("TwoGamma_vs_EpEm", ";M(e^{+}e^{-}) [GeV];M(2#gamma) [GeV]", 1600,0.,4.,400,0,1.); } gDirectory->cd("../"); gDirectory->mkdir("aP2PPipPim2gamma")->cd(); { PPim_vs_aPPip=new TH2F("PPim_vs_aPPip", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,1.9,400,0.9,1.9); AntiLambdaGamma_vs_LambdaGamma=new TH2F("AntiLambdaGamma_vs_LambdaGamma", ";M(#Lambda#gamma) [GeV];M(#bar{#Lambda}#gamma) [GeV]", 400,0.9,1.9,400,0.9,1.9); Sigma0AntiSigma0=new TH1F("Sigma0AntiSigma0", ";M(#Sigma^{0}#bar{#Sigma}^{0}) [GeV]", 800,2.,4.); CosThetaAntiLambdaRest=new TH1F("CosThetaAntiLambdaRest","cos(#theta(anti-proton #gamma))",20,-1,1); CosThetaLambdaRest=new TH1F("CosThetaLambdaRest","cos(#theta(proton #gamma))",20,-1,1); } 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->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); // 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(EtaEtaPrime_,my_pair_vector); particle_pairs.emplace(EtaEta_,my_pair_vector); bg_particle_pairs.emplace(Pi0Pi0_,my_pair_vector); bg_particle_pairs.emplace(Pi0Eta_,my_pair_vector); bg_particle_pairs.emplace(Pi0EtaPrime_,my_pair_vector); bg_particle_pairs.emplace(EtaEtaPrime_,my_pair_vector); bg_particle_pairs.emplace(EtaEta_,my_pair_vector); vectormy_weight_vector; bg_weights.emplace(Pi0Pi0_,my_weight_vector); bg_weights.emplace(Pi0Eta_,my_weight_vector); bg_weights.emplace(Pi0EtaPrime_,my_weight_vector); bg_weights.emplace(EtaEtaPrime_,my_weight_vector); bg_weights.emplace(EtaEta_,my_weight_vector); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_survey::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_survey::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.; MC_beam_E_all->Fill(mc_beam_E); 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); vectormcprotons,mcneutrons,mcetas,mcgammas,mcpi0s; for (unsigned int i=0;imech==0){ switch(throwns[i]->type){ case Proton: mcprotons.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; } } } MC_beam_E->Fill(tagged_mc_beam_E); if (mcprotons.size()==1){ double t=(target-mcprotons[0]).M2(); DLorentzVector missing=beam+target-mcprotons[0]; double mass=missing.M(); MC_meson_mass->Fill(mass); MC_t->Fill(-t); } if (mcneutrons.size()==1){ } } japp->RootUnLock(); } vectortracks; loop->Get(tracks); if (tracks.size()==0) return RESOURCE_UNAVAILABLE; vectorbeamphotons; loop->Get(beamphotons); if (beamphotons.size()==0) return RESOURCE_UNAVAILABLE; vectorneutrals; loop->Get(neutrals); const DEventHitStatistics *stats=NULL; loop->GetSingle(stats); vectortofs; loop->Get(tofs); japp->WriteLock("mykinlock"); int num_fdc_pseudos=0; if (stats!=NULL){ NFDCPseudos->Fill(stats->fdc_pseudos); num_fdc_pseudos=stats->fdc_pseudos; } // 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 int num_in_insert=0,num_tof_match=0; map>neutralParticles; FillNeutralParticleVectors(t0_rf,vertex,tofs,neutrals,neutralParticles, num_in_insert,num_tof_match); // Fill vector of beam photons to pass on the analysis along with weights // for in-time/ out-of-time vector>beamPhotons; for (unsigned int i=0;ilorentzMomentum().E()lorentzMomentum().E()>MAX_BEAM_E) continue; double dt_st=beamphotons[i]->time()-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){ 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_ConfidenceLevel(); ExclusiveCL->Fill(CL); 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 numElectron=final_kf[Electron].size(); unsigned int numPositron=final_kf[Positron].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(); // Proton e+ e- Ngamma if (numElectron==1 && numPositron==1 && numPiPlus==0 && numPiMinus==0 && numKPlus==0 && numKMinus==0 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 0: EpEmAnalysis(beam_kf,final_kf,weight); break; case 2: EpEmTwoGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // 2 Proton AntiProton Pi+ Pi- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==1 && numPiMinus==1 && numKPlus==0 && numKMinus==0 && numProton==2 && numAntiProton==1 && numNeutron==0){ switch(numGamma){ case 2: AntiProtonPipPim2GammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Proton Pi+ Pi- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==1 && numPiMinus==1 && numKPlus==0 && numKMinus==0 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 2: PipPimTwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: PipPimFourGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Proton K+ K- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==0 && numPiMinus==0 && numKPlus==1 && numKMinus==1 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 2: KpKmTwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: KpKmFourGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Proton 2K+ 2K- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==0 && numPiMinus==0 && numKPlus==2 && numKMinus==2 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 0: // TwoKpKmAnalysis(beam_kf,final_kf,weight); break; case 2: //TwpKpKmTwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: //TwoKpKmFourGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Proton K+ K- pi+ pi- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==1 && numPiMinus==1 && numKPlus==1 && numKMinus==1 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 2: KpKmPipPimTwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: KpKmPipPimFourGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Neutron K+ pi+ pi- Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==1 && numPiMinus==1 && numKPlus==1 && numKMinus==0 && numProton==0 && numAntiProton==0 && numNeutron==1){ switch(numGamma){ case 0: NeutronKpPipPimAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Neutron Pi+ Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==1 && numPiMinus==0 && numKPlus==0 && numKMinus==0 && numProton==0 && numAntiProton==0 && numNeutron==1){ switch(numGamma){ case 2: NeutronTwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: NeutronFourGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // Proton Ngamma if (numElectron==0 && numPositron==0 && numPiPlus==0 && numPiMinus==0 && numKPlus==0 && numKMinus==0 && numProton==1 && numAntiProton==0 && numNeutron==0){ switch(numGamma){ case 2: TwoGammaAnalysis(beam_kf,final_kf,weight); break; case 4: FourGammaAnalysis(num_fdc_pseudos,num_in_insert,num_tof_match,beam_kf,final_kf,weight); break; case 6: SixGammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } } } delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; japp->Unlock("mykinlock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_survey::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_survey::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } // Use the PIDFOM scheme to sort tracks according to particle type void JEventProcessor_survey::FillChargedParticleVectors(vector&tracks,map>&chargedParticles) const { Particle_t particle_list[8]={Proton,AntiProton,KPlus,KMinus,PiPlus,PiMinus, Positron,Electron}; for (unsigned int j=0;j >probabilities; for (unsigned int i=0;i<8;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()); } } } void JEventProcessor_survey::FillNeutralParticleVectors(double t0_rf, const DVector3 &vertex, vector&tofs, vector&neutrals, map>&neutralParticles, int &num_in_insert, int &num_tof_match) 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; } 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.){ 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) continue; if (fabs(pos.X())t-fcalshower->getTime(); TofFcaldt->Fill(dt); DVector3 diff2=tofs[j]->pos-pos; TofFcaldxdy->Fill(diff2.x(),diff2.y()); double dr=diff2.Perp(); if (drFill(1000.*tofs[jmatch]->dE1,1000.*tofs[jmatch]->dE2); if (1000.*tofs[jmatch]->dE1>4. && 1000*tofs[jmatch]->dE2>4.) num_tof_match++; } } else if (shower->dDetectorSystem==SYS_BCAL){ if (E(shower->dBCALFCALShower); if (bcalshower->z>BCAL_Z_CUT) continue; } neutralParticles[Gamma].push_back(gamHyp); } } } } // Run the kinematic fitter requiring energy and momentum conservation void JEventProcessor_survey::DoKinematicFit(const DBeamPhoton *beamphoton, map >&chargedParticles, map>&neutralParticles, 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); // Charged particles Particle_t particle_list[8]={PiPlus,PiMinus,KPlus,KMinus,Proton, AntiProton,Electron,Positron}; TLorentzVector vert; bool got_vertex_track=false; for (unsigned int k=0;k<8;k++){ vectormyParticles=chargedParticles[particle_list[k]]; for (unsigned int m=0;mmyParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]); FinalParticles.insert(myParticle); if (got_vertex_track==false){ vert.SetVect(myParticles[m]->position()); vert.SetT(myParticles[m]->t0()); got_vertex_track=true; } } } // Final state photons for (unsigned int k=0;kmyGamma=dKinFitUtils->Make_DetectedParticle(22,0,0.,vert,gammaHyp->momentum(),0.,gammaHyp->errorMatrix()); FinalParticles.insert(myGamma); } // Missing proton if (isMissing){ shared_ptrmyRecoil =dKinFitUtils->Make_MissingParticle(Proton); FinalParticles.insert(myRecoil); } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); } void JEventProcessor_survey::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()); } } } bool JEventProcessor_survey::MakeMesonPairs(vector&paired_doublets, 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[EtaEtaPrime_].clear(); particle_pairs[EtaEta_].clear(); bg_particle_pairs[Pi0Pi0_].clear(); bg_particle_pairs[Pi0Eta_].clear(); bg_particle_pairs[Pi0EtaPrime_].clear(); bg_particle_pairs[EtaEtaPrime_].clear(); bg_particle_pairs[EtaEta_].clear(); bg_weights[Pi0Pi0_].clear(); bg_weights[Pi0Eta_].clear(); bg_weights[Pi0EtaPrime_].clear(); bg_weights[EtaEtaPrime_].clear(); bg_weights[EtaEta_].clear(); bool got_meson_pairs=false; // 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); bool pi0_bg2_hi_cut=PI0_BG_CUT(pair2_mass,NUM_SIGMA_BG,PI0_CUT_VALUE); bool pi0_bg2_lo_cut=PI0_BG_CUT(pair2_mass,-NUM_SIGMA_BG,PI0_CUT_VALUE); bool pi0_bg1_hi_cut=PI0_BG_CUT(pair1_mass,NUM_SIGMA_BG,PI0_CUT_VALUE); bool pi0_bg1_lo_cut=PI0_BG_CUT(pair1_mass,-NUM_SIGMA_BG,PI0_CUT_VALUE); bool eta_bg2_hi_cut=ETA_BG_CUT(pair2_mass,NUM_SIGMA_BG,ETA_CUT_VALUE); bool eta_bg2_lo_cut=ETA_BG_CUT(pair2_mass,-NUM_SIGMA_BG,ETA_CUT_VALUE); bool eta_bg1_hi_cut=ETA_BG_CUT(pair1_mass,NUM_SIGMA_BG,ETA_CUT_VALUE); bool eta_bg1_lo_cut=ETA_BG_CUT(pair1_mass,-NUM_SIGMA_BG,ETA_CUT_VALUE); double bgscale=0.; if (mass>2.*(ParticleMass(Pi0)-PI0_CUT_VALUE)){ if (pi0_1_cut && pi0_bg2_hi_cut) bgscale=1.; if (pi0_1_cut && pi0_bg2_lo_cut) bgscale=1.; if (pi0_2_cut && pi0_bg1_hi_cut) bgscale=1.; if (pi0_2_cut && pi0_bg1_lo_cut) bgscale=1.; if (pi0_bg1_hi_cut && pi0_bg2_hi_cut) bgscale=-1.; if (pi0_bg1_lo_cut && pi0_bg2_lo_cut) bgscale=-1.; if (pi0_bg2_lo_cut && pi0_bg1_hi_cut) bgscale=-1.; if (pi0_bg2_hi_cut && pi0_bg1_lo_cut) bgscale=-1.; } if (fabs(bgscale)>0.){ bg_particle_pairs[Pi0Pi0_].push_back(make_pair(pair1,pair2)); bg_weights[Pi0Pi0_].push_back(weight*bgscale); } double bgscale_2eta=0.; if (mass>2.*(ParticleMass(Eta)-ETA_CUT_VALUE)){ if (eta_1_cut && eta_bg2_hi_cut) bgscale_2eta=1.; if (eta_1_cut && eta_bg2_lo_cut) bgscale_2eta=1.; if (eta_2_cut && eta_bg1_hi_cut) bgscale_2eta=1.; if (eta_2_cut && eta_bg1_lo_cut) bgscale_2eta=1.; if (eta_bg1_hi_cut && eta_bg2_hi_cut) bgscale_2eta=-1.; if (eta_bg1_lo_cut && eta_bg2_lo_cut) bgscale_2eta=-1.; if (eta_bg2_lo_cut && eta_bg1_hi_cut) bgscale_2eta=-1.; if (eta_bg2_hi_cut && eta_bg1_lo_cut) bgscale_2eta=-1.; } if (fabs(bgscale_2eta)>0.){ bg_particle_pairs[EtaEta_].push_back(make_pair(pair1,pair2)); bg_weights[EtaEta_].push_back(weight*bgscale_2eta); } double bgscale_pi0eta_1=0.; if (mass>ParticleMass(Pi0)+ParticleMass(Eta)-ETA_CUT_VALUE-PI0_CUT_VALUE){ if (pi0_1_cut && eta_bg2_hi_cut) bgscale_pi0eta_1=1.; if (pi0_1_cut && eta_bg2_lo_cut) bgscale_pi0eta_1=1.; if (eta_2_cut && pi0_bg1_hi_cut) bgscale_pi0eta_1=1.; if (eta_2_cut && pi0_bg1_lo_cut) bgscale_pi0eta_1=1.; if (pi0_bg1_hi_cut && eta_bg2_hi_cut) bgscale_pi0eta_1=-1.; if (pi0_bg1_lo_cut && eta_bg2_lo_cut) bgscale_pi0eta_1=-1.; if (eta_bg2_lo_cut && pi0_bg1_hi_cut) bgscale_pi0eta_1=-1.; if (eta_bg2_hi_cut && pi0_bg1_lo_cut) bgscale_pi0eta_1=-1.; } if (fabs(bgscale_pi0eta_1)>0.){ double bgweight=weight*bgscale_pi0eta_1; bg_particle_pairs[Pi0Eta_].push_back(make_pair(pair1,pair2)); bg_weights[Pi0Eta_].push_back(bgweight); } double bgscale_pi0eta_2=0.; if (mass>ParticleMass(Pi0)+ParticleMass(Eta)-ETA_CUT_VALUE-PI0_CUT_VALUE){ if (eta_1_cut && pi0_bg2_hi_cut) bgscale_pi0eta_2=1.; if (eta_1_cut && pi0_bg2_lo_cut) bgscale_pi0eta_2=1.; if (pi0_2_cut && eta_bg1_hi_cut) bgscale_pi0eta_2=1.; if (pi0_2_cut && eta_bg1_lo_cut) bgscale_pi0eta_2=1.; if (eta_bg1_hi_cut && pi0_bg2_hi_cut) bgscale_pi0eta_2=-1.; if (eta_bg1_lo_cut && pi0_bg2_lo_cut) bgscale_pi0eta_2=-1.; if (pi0_bg2_lo_cut && eta_bg1_hi_cut) bgscale_pi0eta_2=-1.; if (pi0_bg2_hi_cut && eta_bg1_lo_cut) bgscale_pi0eta_2=-1.; } if (fabs(bgscale_pi0eta_2)>0.){ double bgweight=weight*bgscale_pi0eta_2; bg_particle_pairs[Pi0Eta_].push_back(make_pair(pair2,pair1)); bg_weights[Pi0Eta_].push_back(bgweight); } if (pi0_1_cut && pi0_2_cut){ particle_pairs[Pi0Pi0_].push_back(make_pair(pair1,pair2)); got_meson_pairs=true; } if (eta_1_cut && eta_2_cut){ particle_pairs[EtaEta_].push_back(make_pair(pair1,pair2)); got_meson_pairs=true; } if (eta_1_cut && pi0_2_cut){ particle_pairs[Pi0Eta_].push_back(make_pair(pair2,pair1)); got_meson_pairs=true; } if (eta_2_cut && pi0_1_cut){ particle_pairs[Pi0Eta_].push_back(make_pair(pair1,pair2)); got_meson_pairs=true; } if (pi0_1_cut && ETAPRIME_CUT(pair2_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair1,pair2)); got_meson_pairs=true; } if (pi0_2_cut && ETAPRIME_CUT(pair1_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair2,pair1)); got_meson_pairs=true; } if (eta_1_cut && ETAPRIME_CUT(pair2_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[EtaEtaPrime_].push_back(make_pair(pair1,pair2)); got_meson_pairs=true; } if (eta_2_cut && ETAPRIME_CUT(pair1_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[EtaEtaPrime_].push_back(make_pair(pair2,pair1)); got_meson_pairs=true; } } return got_meson_pairs; } void JEventProcessor_survey::FillPairedMesonMassHistos(double weight,double mass, TH1F *myPi0Pi0,TH1F *myPi0Pi0bg, TH1F *myPi0Eta,TH1F *myPi0Etabg, TH1F *myEtaEta,TH1F *myEtaEtabg ){ if (particle_pairs[Pi0Pi0_].size()>0){ myPi0Pi0->Fill(mass,weight); } unsigned int num_pi0pi0_bg=bg_particle_pairs[Pi0Pi0_].size(); if (num_pi0pi0_bg>0){ double bg_scale=1./double(num_pi0pi0_bg); for (unsigned int i=0;iFill(mass,bg_weight); } } if (particle_pairs[EtaEta_].size()>0){ myEtaEta->Fill(mass,weight); } unsigned int num_etaeta_bg=bg_particle_pairs[EtaEta_].size(); if (num_etaeta_bg>0){ double bg_scale=1./double(num_etaeta_bg); for (unsigned int i=0;iFill(mass,bg_weight); } } if (particle_pairs[Pi0Eta_].size()>0){ myPi0Eta->Fill(mass,weight); } unsigned int num_pi0eta_bg=bg_particle_pairs[Pi0Eta_].size(); if (num_pi0eta_bg>0){ double bg_scale=1./double(num_pi0eta_bg); for (unsigned int i=0;iFill(mass,bg_weight); } } } void JEventProcessor_survey::FillGFHisto(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); } //------------------------------------------------------------------------------ // Proton Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::TwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; TwoGammaMass->Fill(twogam_kf.M(),weight); } void JEventProcessor_survey::FourGammaAnalysis(int num_fdcs, int num_in_insert, int num_tof_match, 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->Fill(mass,weight); if (num_in_insert>=3){ FourGamma_Insertcut->Fill(mass,weight); } vectorpaired_doublets(doublets.size()); bool got_meson_pairs=MakeMesonPairs(paired_doublets,final_kf[Gamma],weight, FourGamma2d); if (got_meson_pairs){ FillPairedMesonMassHistos(weight,mass,Pi0Pi0,Pi0Pi0bg,Pi0Eta,Pi0Etabg, EtaEta,EtaEtabg); unsigned int num_pi0pi0=particle_pairs[Pi0Pi0_].size(); if (num_pi0pi0>0){ double scaled_weight=weight/double(num_pi0pi0); for (unsigned int i=0;i0){ double bg_scale=1./double(num_pi0pi0_bg); for (unsigned int i=0;i0){ double scaled_weight=weight/double(num_pi0eta); for (unsigned int i=0;i0){ double bg_scale=1./double(num_pi0eta_bg); for (unsigned int i=0;iFill(mass,weight); if (num_fdcs==0){ FourGammaNoPairs_FDCcut->Fill(mass,weight); } if (num_tof_match==0){ FourGammaNoPairs_TOFcut->Fill(mass,weight); } if (num_in_insert>=3){ FourGammaNoPairs_Insertcut->Fill(mass,weight); } for (unsigned int i=0;iFill(pair1_mass,pair2_mass,weight); if (PI0_CUT(pair1_mass,PI0_CUT_VALUE)){ DLorentzVector pi0g_1=pair1+g3; DLorentzVector pi0g_2=pair1+g4; Pi02g_vs_Pi0g->Fill(pi0g_1.M(),mass,weight); Pi02g_vs_Pi0g->Fill(pi0g_2.M(),mass,weight); Pi02g_vs_2g->Fill(pair2.M(),mass,weight); } if (PI0_CUT(pair2_mass,PI0_CUT_VALUE)){ DLorentzVector pi0g_1=pair2+g1; DLorentzVector pi0g_2=pair2+g2; Pi02g_vs_Pi0g->Fill(pi0g_1.M(),mass,weight); Pi02g_vs_Pi0g->Fill(pi0g_2.M(),mass,weight); Pi02g_vs_2g->Fill(pair1.M(),mass,weight); } } } } void JEventProcessor_survey::SixGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector sixgam_kf; for (unsigned int i=0;i<6;i++) sixgam_kf+=final_kf[Gamma][i]; SixGammaMass->Fill(sixgam_kf.M(),weight); } //------------------------------------------------------------------------------ // Neutron Pi+ Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::NeutronTwoGammaAnalysis(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(); 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); FillGFHisto(weight,beam_kf,pip2gam_kf,twogam_kf,Pi0CosThetaGF_n); } if (PI0_BG_CUT(twogam_mass,NUM_SIGMA_BG,PI0_CUT_VALUE) ||PI0_BG_CUT(twogam_mass,-NUM_SIGMA_BG,PI0_CUT_VALUE)){ FillGFHisto(weight,beam_kf,pip2gam_kf,twogam_kf,Pi0CosThetaGF_n_bg); } if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE)){ PipEta_NEta_Dalitz->Fill(pip2gam_kf.M2(),n2gam_kf.M2(),weight); FillGFHisto(weight,beam_kf,pip2gam_kf,twogam_kf,EtaCosThetaGF_n); } if (ETA_BG_CUT(twogam_mass,NUM_SIGMA_BG,ETA_CUT_VALUE) ||ETA_BG_CUT(twogam_mass,-NUM_SIGMA_BG,ETA_CUT_VALUE)){ FillGFHisto(weight,beam_kf,pip2gam_kf,twogam_kf,EtaCosThetaGF_n_bg); } } void JEventProcessor_survey::NeutronFourGammaAnalysis(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); vectorpaired_doublets(doublets.size()); MakeMesonPairs(paired_doublets,final_kf[Gamma],weight,FourGamma2d_NPip); FillPairedMesonMassHistos(weight,mass,Pi0Pi0_NPip,Pi0Pi0bg_NPip, Pi0Eta_NPip,Pi0Etabg_NPip,EtaEta_NPip, EtaEtabg_NPip); } //------------------------------------------------------------------------------ // Proton Pi+ Pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::PipPimTwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector pippim_kf=final_kf[PiPlus][0]+final_kf[PiMinus][0]; DLorentzVector pippim2g_kf=pippim_kf+twogam_kf; TwoGammaMass_PipPim->Fill(twogam_kf.M(),weight); if (PI0_CUT(twogam_kf.M(),PI0_CUT_VALUE)){ PipPimPi0->Fill(pippim2g_kf.M(),weight); } if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE)){ PipPimEta->Fill(pippim2g_kf.M(),weight); } } void JEventProcessor_survey::PipPimFourGammaAnalysis(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 mass4g=fourgam_kf.M(); DLorentzVector pippim_kf=final_kf[PiPlus][0]+final_kf[PiMinus][0]; FourGammaMass_PipPim->Fill(mass4g,weight); PipPimMass_4g->Fill(pippim_kf.M(),weight); vectorpaired_doublets(doublets.size()); MakeMesonPairs(paired_doublets,final_kf[Gamma],weight,FourGamma2d_PipPim); FillPairedMesonMassHistos(weight,mass4g,Pi0Pi0_PipPim,Pi0Pi0bg_PipPim, Pi0Eta_PipPim,Pi0Etabg_PipPim,EtaEta_PipPim, EtaEtabg_PipPim); } //------------------------------------------------------------------------------ // Proton e+ e- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::EpEmAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; EpEmMass->Fill(epem_kf.M(),weight); } void JEventProcessor_survey::EpEmTwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; TwoGamma_vs_EpEm->Fill(epem_kf.M(),twogam_kf.M(),weight); } //------------------------------------------------------------------------------ // Proton K+ K- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::KpKmTwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector kpkm_kf=final_kf[KPlus][0]+final_kf[KMinus][0]; DLorentzVector kpkm2g=twogam_kf+kpkm_kf; DLorentzVector p2g=final_kf[Proton][0]+twogam_kf; TwoGamma_vs_KpKm->Fill(kpkm_kf.M(),twogam_kf.M(),weight); if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE) && PHI_CUT(kpkm_kf.M())){ EtaPhiDalitz->Fill(p2g.M2(),kpkm2g.M2(),weight); } } void JEventProcessor_survey::KpKmFourGammaAnalysis(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 mass4g=fourgam_kf.M(); DLorentzVector kpkm_kf=final_kf[KPlus][0]+final_kf[KMinus][0]; FourGammaMass_KpKm->Fill(mass4g,weight); KpKmMass_4g->Fill(kpkm_kf.M(),weight); vectorpaired_doublets(doublets.size()); MakeMesonPairs(paired_doublets,final_kf[Gamma],weight,FourGamma2d_KpKm); FillPairedMesonMassHistos(weight,mass4g,Pi0Pi0_KpKm,Pi0Pi0bg_KpKm, Pi0Eta_KpKm,Pi0Etabg_KpKm,EtaEta_KpKm, EtaEtabg_KpKm); if (PHI_CUT(kpkm_kf.M())){ FillPairedMesonMassHistos(weight,mass4g,Pi0Pi0_Phi,Pi0Pi0bg_Phi, Pi0Eta_Phi,Pi0Etabg_Phi,EtaEta_Phi, EtaEtabg_Phi); } if (PHI_BG_CUT(kpkm_kf.M(),NUM_SIGMA_BG)){ FillPairedMesonMassHistos(weight,mass4g,Pi0Pi0_Phi_bg,Pi0Pi0bg_Phi_bg, Pi0Eta_Phi_bg,Pi0Etabg_Phi_bg,EtaEta_Phi_bg, EtaEtabg_Phi_bg); } } //------------------------------------------------------------------------------ // Proton K+ K- pi+ pi- Ngamma //------------------------------------------------------------------------------ /*void JEventProcessor_survey::TwoKpKmTwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; TwoGamMass_2kpkm->Fill(twogam_kf.M(),weight); DLorentzVector kpkm_1=final_kf[KPlus][0]+final_kf[KMinus][0]; DLorentzVector kpkm_2=final_kf[KPlus][1]+final_kf[KMinus][1]; double kpkm_mass_1=kpkm_1.M(); double kpkm_mass_2=kpkm_2.M(); KpKm_vs_KpKm_2g->Fill(kpkm_mass_1,kpkm_mass_2,weight); DLorentzVector kpkm_3=final_kf[KPlus][0]+final_kf[KMinus][1]; DLorentzVector kpkm_4=final_kf[KPlus][1]+final_kf[KMinus][0]; double kpkm_mass_3=kpkm_3.M(); double kpkm_mass_4=kpkm_4.M(); KpKm_vs_KpKm_2g->Fill(kpkm_mass_3,kpkm_mass_4,weight); } */ //------------------------------------------------------------------------------ // Proton K+ K- pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::KpKmPipPimTwoGammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector twogam_kf=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector kpkm_kf=final_kf[KPlus][0]+final_kf[KMinus][0]; double kpkm_mass=kpkm_kf.M(); DLorentzVector pippim_kf=final_kf[PiPlus][0]+final_kf[PiMinus][0]; DLorentzVector kpkm2g=twogam_kf+kpkm_kf; DLorentzVector p2g=final_kf[Proton][0]+twogam_kf; DLorentzVector pippim2g=pippim_kf+twogam_kf; double pippim2g_mass=pippim2g.M(); DLorentzVector kpkmpippim2g=kpkm_kf+pippim2g; TwoGamma_vs_KpKm_with_PipPim->Fill(kpkm_mass,twogam_kf.M(),weight); if (PHI_CUT(kpkm_mass)){ if (PI0_CUT(twogam_kf.M(),PI0_CUT_VALUE)){ PipPimPi0_KpKm->Fill(pippim2g_mass,weight); if (ETA_CUT(pippim2g_mass,ETA_CUT_VALUE)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, EtaCosThetaGF_with_Phi); } if (OMEGA_CUT(pippim2g_mass)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, OmegaCosThetaGF_with_Phi); } } if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE)){ PipPimEta_KpKm->Fill(pippim2g_mass,weight); if (ETAPRIME_CUT(pippim2g_mass,ETAPRIME_CUT_VALUE)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, EtaPrimeCosThetaGF_with_Phi); } } } if (PHI_BG_CUT(kpkm_mass,NUM_SIGMA_BG) || PHI_BG_CUT(kpkm_mass,-NUM_SIGMA_BG)){ if (PI0_CUT(twogam_kf.M(),PI0_CUT_VALUE)){ PipPimPi0_KpKm_bg->Fill(pippim2g_mass,weight); if (ETA_CUT(pippim2g_mass,ETA_CUT_VALUE)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, EtaCosThetaGF_with_Phi_bg); } if (OMEGA_CUT(pippim2g_mass)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, OmegaCosThetaGF_with_Phi_bg); } } if (ETA_CUT(twogam_kf.M(),ETA_CUT_VALUE)){ PipPimEta_KpKm_bg->Fill(pippim2g_mass,weight); if (ETAPRIME_CUT(pippim2g_mass,ETAPRIME_CUT_VALUE)){ FillGFHisto(weight,beam_kf,kpkmpippim2g,pippim2g, EtaPrimeCosThetaGF_with_Phi_bg); } } } } void JEventProcessor_survey::KpKmPipPimFourGammaAnalysis(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 mass4g=fourgam_kf.M(); DLorentzVector kpkm_kf=final_kf[KPlus][0]+final_kf[KMinus][0]; FourGammaMass_KpKmPipPim->Fill(mass4g,weight); KpKmMass_4gPipPim->Fill(kpkm_kf.M(),weight); vectorpaired_doublets(doublets.size()); MakeMesonPairs(paired_doublets,final_kf[Gamma],weight,FourGamma2d_KpKmPipPim); FillPairedMesonMassHistos(weight,mass4g,Pi0Pi0_KpKmPipPim,Pi0Pi0bg_KpKmPipPim, Pi0Eta_KpKmPipPim,Pi0Etabg_KpKmPipPim,EtaEta_KpKmPipPim, EtaEtabg_KpKmPipPim); } //------------------------------------------------------------------------------ // AntiProton 2Proton pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::FillSigma0Histos(const DLorentzVector &lambda, const DLorentzVector &lambda_bar, DLorentzVector gamma1, DLorentzVector gamma2, DLorentzVector proton, DLorentzVector antiproton, double weight) const { DLorentzVector lg=lambda+gamma1; double lg_mass=lg.M(); DLorentzVector alg=lambda_bar+gamma2; double alg_mass=alg.M(); AntiLambdaGamma_vs_LambdaGamma->Fill(lg_mass,alg_mass,weight); if (SIGMA0_CUT(alg_mass) && SIGMA0_CUT(lg_mass)){ DLorentzVector lg_alg=lg+alg; double lg_alg_mass=lg_alg.M(); Sigma0AntiSigma0->Fill(lg_alg_mass,weight); DVector3 beta_al=(1./lambda_bar.E())*lambda_bar.Vect(); antiproton.Boost(-beta_al); gamma2.Boost(-beta_al); DVector3 apdir=antiproton.Vect(); apdir.SetMag(1.); DVector3 g2dir=gamma2.Vect(); g2dir.SetMag(1.); double cosTheta_apg=apdir.Dot(g2dir); CosThetaAntiLambdaRest->Fill(cosTheta_apg,weight); DVector3 beta_l=(1./lambda.E())*lambda.Vect(); proton.Boost(-beta_l); gamma1.Boost(-beta_l); DVector3 pdir=proton.Vect(); pdir.SetMag(1.); DVector3 g1dir=gamma1.Vect(); g1dir.SetMag(1.); double cosTheta_pg=pdir.Dot(g1dir); CosThetaLambdaRest->Fill(cosTheta_pg,weight); } } void JEventProcessor_survey::AntiProtonPipPim2GammaAnalysis(DLorentzVector beam_kf, map>&final_kf, double weight) const{ DLorentzVector antiproton=final_kf[AntiProton][0]; DLorentzVector appip=antiproton+final_kf[PiPlus][0]; double appip_mass=appip.M(); DLorentzVector gamma1=final_kf[Gamma][0]; DLorentzVector gamma2=final_kf[Gamma][1]; for (unsigned int i=0;i<2;i++){ DLorentzVector proton=final_kf[Proton][i]; DLorentzVector ppim=proton+final_kf[PiMinus][0]; double ppim_mass=ppim.M(); PPim_vs_aPPip->Fill(appip_mass,ppim_mass,weight); if (LAMBDA_CUT(appip_mass)&&LAMBDA_CUT(ppim_mass)){ FillSigma0Histos(ppim,appip,gamma1,gamma2,proton,antiproton,weight); FillSigma0Histos(ppim,appip,gamma2,gamma1,proton,antiproton,weight); } } } //------------------------------------------------------------------------------ // Neutron K+ pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_survey::NeutronKpPipPimAnalysis(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); } }