// $Id$ // // File: JEventProcessor_hyperon.cc // Created: Wed Aug 30 09:03:49 EDT 2017 // Creator: staylor (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64) // #include "JEventProcessor_hyperon.h" using namespace jana; #include #include #include #include "../include/mass_cuts.h" // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_hyperon()); } } // "C" //------------------ // JEventProcessor_hyperon (Constructor) //------------------ JEventProcessor_hyperon::JEventProcessor_hyperon() { } //------------------ // ~JEventProcessor_hyperon (Destructor) //------------------ JEventProcessor_hyperon::~JEventProcessor_hyperon() { } //------------------ // init //------------------ jerror_t JEventProcessor_hyperon::init(void) { japp->CreateLock("mykpilock"); EMULATE_TRIGGER=true; gPARMS->SetDefaultParameter("HYPERON:EMULATE_TRIGGER",EMULATE_TRIGGER); DROP_ONE_CELL_BCAL_SHOWERS=false; gPARMS->SetDefaultParameter("HYPERON:DROP_ONE_CELL_BCAL_SHOWERS",DROP_ONE_CELL_BCAL_SHOWERS); KAON_GAMMA_DT_CUT=2.; gPARMS->SetDefaultParameter("HYPERON:KAON_GAMMA_DT_CUT", KAON_GAMMA_DT_CUT); PIM_CL_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:PIM_CL_CUT",PIM_CL_CUT); PIP_CL_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:PIP_CL_CUT",PIM_CL_CUT); PROTON_CL_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:PROTON_CL_CUT",PROTON_CL_CUT); KP_CL_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:KP_CL_CUT",KP_CL_CUT); PI0_CUT_VALUE=0.025; ETA_CUT_VALUE=0.05; ETAPRIME_CUT_VALUE=0.05; PI0_VETO_CUT=0.042; // 3 sigma? gPARMS->SetDefaultParameter("HYPERON:PI0_VETO_CUT",PI0_VETO_CUT); SPLIT_CUT=0.5; gPARMS->SetDefaultParameter("HYPERON:SPLIT_CUT",SPLIT_CUT); FCAL_RADIAL_CUT=104.0; gPARMS->SetDefaultParameter("HYPERON:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT); BCAL_Z_CUT=384.0; gPARMS->SetDefaultParameter("HYPERON:BCAL_Z_CUT",BCAL_Z_CUT); NUM_SIGMA_BG=2; gPARMS->SetDefaultParameter("HYPERON:NUM_SIGMA_BG",NUM_SIGMA_BG); BETA_CUT=0.9; // This is called once at program startup. gDirectory->mkdir("hyperon")->cd(); gDirectory->mkdir("MC")->cd(); MCBeam=new TH1F("MCBeam","Thrown beam energy",440,2.8,11.6); MCt= new TH1F("MCt","; |t| [GeV^{2}]",2000,0.,20.); gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); { KpBeamPhotonTimeDiff=new TH2F("KpBeamPhotonTimeDiff","t(tag)-t(rf) at vertex", 48,2.8,12.4,360,-18,18); RFGammaTimeDiff=new TH1F("RFGammaTimeDiff","#deltat between rf time and #gamma at vertex",200,-10,10); RFGammaTimeDiffNoNeutron=new TH1F("RFGammaTimeDiffNoNeutron","#deltat between rf time and #gamma at vertex",200,-10,10); NeutralBeta=new TH1F("NeutralBeta","#beta for neutrals",100,0,1.5); KinFitCL=new TH1F("KinFitCL","fit probability",100,0,1); KinFitChiSq=new TH1F("KinFitChiSq","#chi^2/ndf",1000,0,100); NeutrondTvsZ=new TH2F("NeutrondTvsZ","#delta t vs z at cal",650,0,650,100,-10,10); NeutrondXvsZ=new TH2F("NeutrondXvsZ","#delta x vs z at cal",650,0,650,100,-20,20); NeutrondYvsZ=new TH2F("NeutrondYvsZ","#delta y vs z at cal",650,0,650,100,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("KpPipPimX")->cd(); { KpPipPimNCL=new TH1F("KpPipPimNCL","fit probability",100,0,1); KpPipPimMissingMass=new TH1F("KpPipPimMissingMass","missing mass off K+pi+pi-",1000,0,2); KpPipPimNMissingMassSq=new TH1F("KpPipPimNMissingMassSq","missing mass sq. off K+pi+pi-n",1000,-1,1); NPipSqVsNPimSq_kf=new TH2F("NPipSqVSNPimSq_kf","M^2(n#pi+) vs M^2(n#pi-)", 240,1.,9,240,1.,9); NPipPimMass= new TH1F("NPipPimMass","n#pi+#pi- mass",1600,0.9,4.9); NPip_vs_Kpim=new TH2F("NPip_vs_Kpim","M(n#pi+) vs M(K#pi-)", 800,0.,2,1600,1.,5); NKp_vs_pippim=new TH2F("NKp_vs_pippim","M(nK+) vs M(#pi+#pi-)", 800,0.,2,1600,1.,5); SigmaPlusPimMass= new TH1F("SigmaPlusPimMass","n#pi+#pi- mass",1600,0.9,4.9); SigmaMinusPipMass= new TH1F("SigmaMinusPipMass","n#pi+#pi- mass",1600,0.9,4.9); KpK0sqVsNK0sq=new TH2F("KpK0sqVsNK0sq","M^2(K+K0) vs M^2(nK0)",1000,1.,11., 1000,0.5,10.5); } gDirectory->cd("../"); gDirectory->mkdir("KpPim")->cd(); { MissingPChiSq=new TH1F("MissingPChiSq","#chi^2/ndf",1000,0,100); MissingMassOffKaon=new TH1F("MissingMassOffKaon","K missing mass",1600,0,4.); MissingMassOffKaonCut=new TH1F("MissingMassOffKaonCut","K missing mass",1600,0,4.); MissingMassOffKaonKF=new TH1F("MissingMassOffKaonKF","K missing mass",1600,0,4.); MissingMassOffKPi=new TH1F("MissingMassOffKPi","K#pi missing mass",1600,0,4.); MissingProton=new TH2F("MissingProton","Predicted p vs #theta",120,0,120,100,0,10); MissingProtonRecon=new TH2F("MissingProtonRecon","Recon p vs #theta",120,0,120,100,0,10); MissingProtonReconCut=new TH2F("MissingProtonReconCut","Recon p vs #theta",120,0,120,100,0,10); ReconMissingPDiff=new TH1F("ReconMissingPDiff",";#Deltap [GeV/c]",100,-2.5,2.5); ReconMissingThetaDiff=new TH1F("ReconMissingThetaDiff",";#Delta#theta [radians]",100,-2.5,2.5); ReconMissingPhiDiff=new TH1F("ReconMissingPhiDiff",";#Delta#phi [radians]",100,-2.5,-2.5); PPimMass_with_K= new TH1F("PPimMass_with_K","p#pi- mass",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpPimGamma")->cd(); { PPimGammaMass_with_K= new TH1F("PPimGammaMass_with_K","p#pi-#gamma mass",1600,0,4.); PPimMass_with_Kg= new TH1F("PPimMass_with_Kg","p#pi- mass",1600,0,4.); LambdaGammaMass= new TH1F("LambdaGammaMass","p#pi-#gamma mass",1600,0,4.); LambdaGammaThetaVsMass= new TH2F("LambdaGammaThetaVsMass",";#Lambda#gamma mass [GeV];cos(#theta_{H})",800,0,4.,50,-1,1); LambdaGammaMassVsE= new TH2F("LambdaGammaMassVsE","p#pi-#gamma mass vs E", 48,2.8,12.4,1600,0,4.); LambdaGammaMassVst= new TH2F("LambdaGammaMassVst","p#pi-#gamma mass vs t", 2000,0.,20.,1600,0,4.); K_mom_vs_theta_1g=new TH2F("K_mom_vs_theta_1g","K+ p vs #theta",180,0,180,100,0,10); } gDirectory->cd("../"); gDirectory->mkdir("KpPim2Gamma")->cd(); { PPimGammaMass_with_Kg= new TH1F("PPimGammaMass_with_Kg","p#pi-#gamma mass",1600,0,4.); PPimMass_with_K2g= new TH1F("PPimMass_with_K2g","p#pi- mass",1600,0,4.); LambdaGammaMass_K2g= new TH1F("LambdaGammaMass_K2g","p#pi-#gamma mass",1600,0,4.); TwoGammaMass_Kppim= new TH1F("TwoGammaMass_Kppim","2#gamma mass",1600,0,4.); LambdaPi0Mass= new TH1F("LambdaPi0Mass","p#pi-#pi0 mass",1600,0,4.); PPi0_with_K= new TH1F("PPi0_with_K","p#pi0 mass",1600,0,4.); LambdaEtaMass= new TH1F("LambdaEtaMass","p#pi-#eta mass",1600,0,4.); Sigma0GammaMass= new TH1F("Sigma0GammaMass","p#pi-2#gamma mass",1600,0,4.); SigmaPlusPiMinusMass= new TH1F("SigmaPlusPiMinusMass","p#pi-2#gamma mass",1600,0,4.); Sigma0GammaThetaVsMass= new TH2F("Sigma0GammaThetaVsMass",";#Sigma^{0}#gamma mass [GeV];cos(#theta_{H})",800,0,4.,50,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("KpPim3Gamma")->cd(); { TwoGammaMass_kppim3g= new TH1F("TwoGammaMass_kppim3g","2#gamma mass",1600,0,4.); PPimMass_kppim3g= new TH1F("PPimMass_kppim3g","p#pi- mass",1600,0,4.); LambdaGammaMass_kppim3g= new TH1F("LambdaGammaMass_kppim3g","p#pi-#gamma mass",1600,0,4.); LambdaPi0Mass_kppim3g= new TH1F("LambdaPi0Mass_kppim3g","p#pi-#pi0 mass",1600,0,4.); Sigma0Pi0Mass_kppim3g= new TH1F("Sigma0Pi0Mass_kppim3g","p#pi-#pi0 mass",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpPim4Gamma")->cd(); { PPim4GammaMass_kppim4g=new TH1F("PPim4GammaMass_kppim4g","p#pi-4#gamma invariant mass",1600,1,5.); PPim4GammaMass_kppim4g->SetXTitle("M(p#pi^{-}4#gamma) [GeV]"); PPim4GammaMass_kppim4g->SetYTitle("counts / 0.0025 GeV"); PPimMass_kppim4g= new TH1F("PPimMass_kppim4g","p#pi- mass",1600,0,4.); FourGamma2d_kppim4g=new TH2F("FourGamma2d_kppim4g","m(2#gamma,pair2) vs m(2#gamma,pair1)",100,0,1,100,0,1); LambdaPi0Pi0Dalitz=new TH2F("LambdaPi0Pi0Dalitz","M^2(#Lambda#pi0_1) vs M^2(#Lambda#pi02)",300,1,16,300,1,16); LambdaEtaPi0Dalitz=new TH2F("LambdaEtaPi0Dalitz","M^2(#Lambda#eta) vs M^2(#Lambda#pi0)",300,1,16,300,1,16); FourGammaMass_kppim4g=new TH1F("FourGammaMass_kppim4g","M(4#gamma), p#pi#pi events",1600,0,4.); FourGamma2d_kppim4g->SetYTitle("M(2#gamma,pair 2) [GeV]"); FourGamma2d_kppim4g->SetXTitle("M(2#gamma,pair 1) [GeV]"); LambdaPi0Pi0Mass= new TH1F("LambdaPi0Pi0Mass","p#pi-2#pi0 mass",1600,1,5.); LambdaPi0Pi0Dalitz_bg=new TH2F("LambdaPi0Pi0Dalitz_bg","M^2(#Lambda#pi0_1) vs M^2(#Lambda#pi02)",300,1,16,300,1,16); LambdaEtaPi0Dalitz_bg=new TH2F("LambdaEtaPi0Dalitz_bg","M^2(#Lambda#eta) vs M^2(#Lambda#pi0)",300,1,16,300,1,16); LambdaPi0Pi0Mass_bg= new TH1F("LambdaPi0Pi0Mass_bg","p#pi-2#pi0 mass",1600,1,5.); FourGamma2d_bg=new TH2F("FourGamma2d_kppim4g_bg","m(2#gamma,pair2) vs m(2#gamma,pair1)",100,0,1,100,0,1); } gDirectory->cd("../"); gDirectory->mkdir("2Kp2Pim")->cd(); { TwoKpTwoPimMass= new TH1F("TwoKpTwoPimMass","2K+2#pi- mass",1600,0,4.); KpPim_vs_KpPim=new TH2F("KpPim_vs_KpPim","M(K+#pi- pair 1) vs M(K+#pi- pair 2)",200,0,2,200,0,2); P2PimMass_2kp= new TH1F("P2PimMass_2kp","p2#pi- mass",1600,0,4.); PPimMass_2kp= new TH1F("PPimMass_2kp","p#pi- mass",400,1,3.); LambdaPimMass= new TH1F("LambdaPimMass","p#pi-#pi- mass",400,1,3.); } gDirectory->cd("../"); gDirectory->mkdir("2Kp2Pim2g")->cd(); { TwoGammaMass_2kp2g= new TH1F("TwoGammaMass_2kp2g","2#gamma mass",1600,0,4.); PPimMass_2kp2g= new TH1F("PPimMass_2kp2g","p#pi- mass",400,1,3.); LambdaPimMass_2g= new TH1F("LambdaPimMass_2g","p#pi-#pi- mass",400,1,3.); LambdaPimPi0Mass= new TH1F("LambdaPimPi0Mass_2g","p#pi-#pi-#pi0 mass",300,1,4.); XiMinusPi0Mass= new TH1F("XiMinusPi0Mass_2g","p#pi-#pi-#pi0 mass",300,1,4.); } gDirectory->cd("../"); gDirectory->mkdir("K2Pip3Pim")->cd(); { PPim_k2pip3pim= new TH1F("PPim_k2pip3pim","; p#pi^{-} mass [GeV]",400,1,3.); PipPim2D_k2pip3pim=new TH2F("PipPim2D_k2pip3pim",";M(#pi^{+}#pi^{-}, group 1) [GEV];M(#pi^{+}#pi^{-}, group 2) [GeV]",200,0,2,200,0,2); PipPim2D_k2pip3pim_LambdaCut=new TH2F("PipPim2D_k2pip3pim_LambdaCut",";M(#pi^{+}#pi^{-}, group 1) [GEV];M(#pi^{+}#pi^{-}, group 2) [GeV]",200,0,2,200,0,2); LambdaPip_k2pip3pim_Rcut = new TH1F("LambdaPip_k2pip3pim_Rcut",";#Lambda#pi^{+} mass [GeV]",1600,0.9,4.9); PipPim2D_k2pip3pim_Rcut=new TH2F("PipPim2D_k2pip3pim_Rcut",";M(#pi^{+}#pi^{-} pair 1) [GeV]; M(#pi^{+}#pi^{-} pair 2) [GeV]",1000,0,2,1000,0,2); PPimDOCA_k2pip3pim=new TH1F("PPimDOCA_k2pip3pim",";DOCA [cm]",1000,0,10); PPimMassVsR_k2pip3pim=new TH2F("PPimMassVsR_k2pip3pim",";R [cm];M(p#pi^{-}) [GeV]",130,0,65, 1000,0.9,3.9); LambdaPipPimDalitz_Rcut=new TH2F("LambdaPipPimDalitz_Rcut",";M^{2}(#Lambda#pi^{-}) [GeV^{2}];M^{2}(#Lambda#pi^{+}) [GeV^{2}]",800,1.,17.,800,1.,17.); } gDirectory->cd("../"); gDirectory->mkdir("K2PimPip")->cd(); { PPimDOCA_k2pimpip=new TH1F("PPimDOCA_k2pimpip",";DOCA [cm]",1000,0,10); PPimMassVsR_k2pimpip=new TH2F("PPimMassVsR_k2pimpip",";R [cm];M(p#pi^{-}) [GeV]",130,0,65, 1000,0.9,3.9); PPim_vs_Kpim_with_pip=new TH2F("PPim_vs_Kpim_with_pip",";M(K^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); PPim_vs_pippim_with_Kp=new TH2F("PPim_vs_pippim_with_Kp",";M(#pi^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); LambdaPip_vs_Kpim=new TH2F("LambdaPip_vs_Kpim",";M(K^{+}#pi^{-}) [GeV];M(#Lambda#pi^{+}) [GeV]",1600,0.,4.,1600,0.9,4.9); } gDirectory->cd("../"); gDirectory->mkdir("K2PimPip2g")->cd(); { PPimDOCA_k2pimpip2g=new TH1F("PPimDOCA_k2pimpip2g",";DOCA [cm]",1000,0,10); PPimMassVsR_k2pimpip2g=new TH2F("PPimMassVsR_k2pimpip2g",";R [cm];M(p#pi^{-}) [GeV]",130,0,65, 1000,0.9,3.9); PPim_vs_Kpim_with_pip2g=new TH2F("PPim_vs_Kpim_with_pip2g",";M(K^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); PPim_vs_pippim_with_Kp2g=new TH2F("PPim_vs_pippim_with_Kp2g",";M(#pi^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); LambdaPip_vs_Kpim_with_2g=new TH2F("LambdaPip_vs_Kpim_with_2g",";M(K^{+}#pi^{-}) [GeV];M(#Lambda#pi^{+}) [GeV]",1600,0.,4.,1600,0.9,4.9); TwoGamma_k2pimpip=new TH1F("TwoGamma_k2pimpip",";M(2#gamma) [GeV]",1600,0.,4.); LambdaPip_vs_KPimPi0=new TH2F("LambdaPip_vs_KpimPi0",";M(K^{+}#pi^{-}#pi^{0}) [GeV];M(#Lambda#pi^{+}) [GeV]",1600,0.,4.,1600,0.9,4.9); LambdaPipPi0_vs_KPim=new TH2F("LambdaPipPi0_vs_Kpim",";M(K^{+}#pi^{-}) [GeV];M(#Lambda#pi^{+}#pi^{0}) [GeV]",1600,0.,4.,1600,0.9,4.9); } gDirectory->cd("../"); gDirectory->mkdir("K2PimPip4g")->cd(); { PPimDOCA_k2pimpip4g=new TH1F("PPimDOCA_k2pimpip4g",";DOCA [cm]",1000,0,10); PPimMassVsR_k2pimpip4g=new TH2F("PPimMassVsR_k2pimpip4g",";R [cm];M(p#pi^{-}) [GeV]",130,0,65, 1000,0.9,3.9); PPim_vs_Kpim_with_pip4g=new TH2F("PPim_vs_Kpim_with_pip4g",";M(K^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); PPim_vs_pippim_with_Kp4g=new TH2F("PPim_vs_pippim_with_Kp4g",";M(#pi^{+}#pi^{-}) [GeV];M(p#pi^{-}) [GeV]",1600,0.,4.,1600,0.9,4.9); LambdaPip_vs_Kpim_with_4g=new TH2F("LambdaPip_vs_Kpim_with_4g",";M(K^{+}#pi^{-}) [GeV];M(#Lambda#pi^{+}) [GeV]",1600,0.,4.,1600,0.9,4.9); FourGamma_k2pimpip=new TH1F("FourGamma_k2pimpip",";M(4#gamma) [GeV]",1600,0.,4.); FourGamma2d_k2pimpip=new TH2F("FourGamma2d_k2pimpip","m(2#gamma,pair2) vs m(2#gamma,pair1)",100,0,1,100,0,1); LambdaPipPi0_vs_KPimPi0=new TH2F("LambdaPipPi0_vs_KPimPi0",";M(K^{+}#pi^{-}#pi^{0}) [GeV];M(#Lambda#pi^{+}#pi^{0}) [GeV]",1600,0.,4.,1600,0.9,4.9); LambdaPipPi0_vs_KPimPi0_bg=new TH2F("LambdaPipPi0_vs_KPimPi0_bg",";M(K^{+}#pi^{-}#pi^{0}) [GeV];M(#Lambda#pi^{+}#pi^{0}) [GeV]",1600,0.,4.,1600,0.9,4.9); } gDirectory->cd("../"); gDirectory->cd("../"); MIN_BEAM_E=0.0; gPARMS->SetDefaultParameter("HYPERON:MIN_BEAM_E",MIN_BEAM_E); MAX_BEAM_E=12.0; gPARMS->SetDefaultParameter("HYPERON:MAX_BEAM_E",MAX_BEAM_E); CL_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:CL_CUT",CL_CUT); TRACK_CL_CUT=0.; gPARMS->SetDefaultParameter("HYPERON:TRACK_CL_CUT",TRACK_CL_CUT); FCAL_POS_CUT=10.0; gPARMS->SetDefaultParameter("HYPERON:FCAL_POS_CUT",FCAL_POS_CUT); BCAL_THRESHOLD=0.05; gPARMS->SetDefaultParameter("HYPERON:BCAL_THRESHOLD",BCAL_THRESHOLD); FCAL_THRESHOLD=0.1; gPARMS->SetDefaultParameter("HYPERON:FCAL_THRESHOLD",FCAL_THRESHOLD); UNUSED_ENERGY_CUT=0.01; gPARMS->SetDefaultParameter("HYPERON:UNUSED_ENERGY_CUT",UNUSED_ENERGY_CUT); // Make a container to keep track of doublets of photons MakeDoublets(); // 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_hyperon::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_hyperon::evnt(JEventLoop *loop, uint64_t eventnumber) { // This is called for every event. Use of common resources like writing // to a file or filling a histogram should be mutex protected. Using // loop->Get(...) to get reconstructed objects (and thereby activating the // reconstruction algorithm) should be done outside of any mutex lock // since multiple threads may call this method at the same time. // Here's an example: // // vector mydataclasses; // loop->Get(mydataclasses); // // japp->RootFillLock(this); // ... fill historgrams or trees ... // japp->RootFillUnLock(this); vectorbeamphotons; loop->Get(beamphotons); vectortracks; loop->Get(tracks); vectorneutrals; loop->Get(neutrals); vectorthrowns; loop->Get(throwns); vectortagm_geom; loop->Get(tagm_geom); vectortagh_geom; loop->Get(tagh_geom); const DMCReaction* reaction=NULL; vector triggers; if (throwns.size()>0){ loop->GetSingle(reaction); //loop->Get(triggers); } /* if (EMULATE_TRIGGER) { loop->Get(triggers); if (triggers.size()>0){ int trig_mask=triggers[0]->Get_L1TriggerBits(); //cout << "trigger: " << trig_mask <WriteLock("mykpilock"); // Handle thrown events if (reaction!=NULL&& tagm_geom.size()>0 && tagh_geom.size()>0){ FillThrownHistos(reaction,throwns,tagm_geom,tagh_geom); } if (tracks.size()==0){ japp->Unlock("mykpilock"); return NOERROR; } // Assign charge track by PID map>chargedParticles; FillChargedParticleVectors(tracks,chargedParticles); if (chargedParticles[KPlus].size()==0){ japp->Unlock("mykpilock"); return RESOURCE_UNAVAILABLE; } // 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(); // Get final state photon candidates vectorgammaParticles; FillGammaParticleVector(t0_rf,vertex,neutrals,gammaParticles); // 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(); unsigned int NDoF=dKinFitter->Get_NDF(); double ReducedChiSq=ChiSq/double(NDoF); KinFitCL->Fill(CL); KinFitCL->Fill(locConfidenceLevel); if (CL>CL_CUT){ DLorentzVector beam_kf,missing_kf; map>final_kf; GetKF4vectors(dKinFitter,beam_kf,missing_kf,final_kf); // Particle counts unsigned int numKPlus=final_kf[KPlus].size(); unsigned int numKMinus=final_kf[KMinus].size(); unsigned int numPiPlus=final_kf[PiPlus].size(); unsigned int numPiMinus=final_kf[PiMinus].size(); unsigned int numGamma=final_kf[Gamma].size(); unsigned int numProton=final_kf[Proton].size(); unsigned int numAntiProton=final_kf[AntiProton].size(); if (numKPlus==1 && numPiMinus==2 && numKMinus==0 && numProton==1 && numPiPlus==1){ switch(numGamma){ case 0: Kp2PimPipAnalysis(beam_kf,final_kf,weight); break; case 2: Kp2PimPip2GammaAnalysis(beam_kf,final_kf,weight); break; case 4: Kp2PimPip4GammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } if (numKPlus==1 && numPiMinus==3 && numKMinus==0 && numProton==1 && numPiPlus==2 && numGamma==0){ Kp3Pim2PipAnalysis(beam_kf,final_kf,weight); } if (numKPlus==1 && numPiMinus==1 && numKMinus==0 && numProton==1 && numPiPlus==0){ switch(numGamma){ case 0: KpPimAnalysis(beam_kf,final_kf,weight); break; case 1: KpPimGammaAnalysis(beam_kf,final_kf,weight); break; case 2: KpPim2GammaAnalysis(beam_kf,final_kf,weight); break; case 3: KpPim3GammaAnalysis(beam_kf,final_kf,weight); break; case 4: KpPim4GammaAnalysis(beam_kf,final_kf,weight); break; case 6: // KpPim6GammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } if (numKPlus==2 && numPiMinus==2 && numPiPlus==0 && numKMinus==0 && numProton==1){ switch(numGamma){ case 0: TwoKpTwoPimAnalysis(beam_kf,final_kf,weight); break; case 2: TwoKpTwoPim2GammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } } //delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; } japp->Unlock("mykpilock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_hyperon::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_hyperon::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } void JEventProcessor_hyperon::KpPimAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf,kp_kf,pim_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -211: // pi- pim_kf=(*locParticleIterator)->Get_P4(); break; default: break; } } } DLorentzVector ppim_kf=proton_kf+pim_kf; PPimMass_with_K->Fill(ppim_kf.M(),weight); } void JEventProcessor_hyperon::KpPimGammaAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf,gamma_kf; DLorentzVector proton_kf,kp_kf,pim_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -211: // pi- pim_kf=(*locParticleIterator)->Get_P4(); break; case 1: // Gamma gamma_kf=(*locParticleIterator)->Get_P4(); break; default: break; } } } DLorentzVector ppim_kf=proton_kf+pim_kf; DLorentzVector ppimgam=ppim_kf+gamma_kf; double t_mandelstam=(beam_kf-kp_kf).M2(); PPimMass_with_Kg->Fill(ppim_kf.M(),weight); PPimGammaMass_with_K->Fill(ppimgam.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ // Angular distribution in helicity frame DVector3 beta_lab=(1./(ppimgam.E()))*ppimgam.Vect(); gamma_kf.Boost(-beta_lab); double cosTheta=(1./gamma_kf.E())*gamma_kf.Vect().Dot(beta_lab); LambdaGammaThetaVsMass->Fill(ppimgam.M(),cosTheta,weight); LambdaGammaMass->Fill(ppimgam.M(),weight); LambdaGammaMassVsE->Fill(beam_kf.E(),ppimgam.M(),weight); LambdaGammaMassVst->Fill(-t_mandelstam,ppimgam.M(),weight); if (SIGMA0_CUT(ppimgam.M())){ double theta=180./M_PI*kp_kf.Theta(); double p=kp_kf.P(); K_mom_vs_theta_1g->Fill(theta,p,weight); } } } void JEventProcessor_hyperon::KpPim2GammaAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf,kp_kf,pim_kf; vectorgammas_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -211: // pi- pim_kf=(*locParticleIterator)->Get_P4(); break; case 1: // Gamma gammas_kf.push_back((*locParticleIterator)->Get_P4()); break; default: break; } } } DLorentzVector ppim_kf=proton_kf+pim_kf; DLorentzVector ppimgam1=ppim_kf+gammas_kf[0]; DLorentzVector ppimgam2=ppim_kf+gammas_kf[1]; DLorentzVector twogam=gammas_kf[0]+gammas_kf[1]; DLorentzVector ppim2g=ppim_kf+twogam; DLorentzVector ppi0=proton_kf+twogam; PPimMass_with_K2g->Fill(ppim_kf.M(),weight); PPimGammaMass_with_Kg->Fill(ppimgam1.M(),weight); PPimGammaMass_with_Kg->Fill(ppimgam2.M(),weight); TwoGammaMass_Kppim->Fill(twogam.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ if (PI0_ANTI_CUT(twogam.M(),PI0_VETO_CUT) && ETA_ANTI_CUT(twogam.M())){ LambdaGammaMass_K2g->Fill(ppimgam1.M(),weight); DVector3 beta_lab=(1./(ppim2g.E()))*ppim2g.Vect(); if (SIGMA0_CUT(ppimgam1.M())){ Sigma0GammaMass->Fill(ppim2g.M(),weight); // Angular distribution in helicity frame gammas_kf[1].Boost(-beta_lab); double cosTheta=(1./gammas_kf[1].E())*gammas_kf[1].Vect().Dot(beta_lab); Sigma0GammaThetaVsMass->Fill(ppim2g.M(),cosTheta,weight); } LambdaGammaMass_K2g->Fill(ppimgam2.M(),weight); if (SIGMA0_CUT(ppimgam2.M())){ Sigma0GammaMass->Fill(ppim2g.M(),weight); // Angular distribution in helicity frame gammas_kf[0].Boost(-beta_lab); double cosTheta=(1./gammas_kf[0].E())*gammas_kf[0].Vect().Dot(beta_lab); Sigma0GammaThetaVsMass->Fill(ppim2g.M(),cosTheta,weight); } } if (PI0_CUT(twogam.M())){ LambdaPi0Mass->Fill(ppim2g.M(),weight); } if (ETA_CUT(twogam.M())){ LambdaEtaMass->Fill(ppim2g.M(),weight); } } if (PI0_CUT(twogam.M())){ PPi0_with_K->Fill(ppi0.M(),weight); if (SIGMAP_CUT(ppi0.M())){ SigmaPlusPiMinusMass->Fill(ppim2g.M(),weight); } } } void JEventProcessor_hyperon::KpPim3GammaAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf,kp_kf,pim_kf; vectorgammas_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -211: // pi- pim_kf=(*locParticleIterator)->Get_P4(); break; case 1: // Gamma gammas_kf.push_back((*locParticleIterator)->Get_P4()); break; default: break; } } } DLorentzVector ppim_kf=proton_kf+pim_kf; DLorentzVector ppimgam1=ppim_kf+gammas_kf[0]; DLorentzVector ppimgam2=ppim_kf+gammas_kf[1]; DLorentzVector ppimgam3=ppim_kf+gammas_kf[2]; DLorentzVector twogam1=gammas_kf[2]+gammas_kf[1]; DLorentzVector twogam2=gammas_kf[2]+gammas_kf[0]; DLorentzVector twogam3=gammas_kf[0]+gammas_kf[1]; TwoGammaMass_kppim3g->Fill(twogam1.M(),weight); TwoGammaMass_kppim3g->Fill(twogam2.M(),weight); TwoGammaMass_kppim3g->Fill(twogam3.M(),weight); PPimMass_kppim3g->Fill(ppim_kf.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ LambdaGammaMass_kppim3g->Fill(ppimgam1.M(),weight); if (PI0_CUT(twogam1.M())){ DLorentzVector ppimpi0=twogam1+ppim_kf; LambdaPi0Mass_kppim3g->Fill(ppimpi0.M(),weight); if (SIGMA0_CUT(ppimgam1.M())){ DLorentzVector ppim3g=ppimgam1+twogam1; Sigma0Pi0Mass_kppim3g->Fill(ppim3g.M(),weight); } } LambdaGammaMass_kppim3g->Fill(ppimgam2.M(),weight); if (PI0_CUT(twogam2.M())){ DLorentzVector ppimpi0=twogam2+ppim_kf; LambdaPi0Mass_kppim3g->Fill(ppimpi0.M(),weight); if (SIGMA0_CUT(ppimgam2.M())){ DLorentzVector ppim3g=ppimgam2+twogam2; Sigma0Pi0Mass_kppim3g->Fill(ppim3g.M(),weight); } } LambdaGammaMass_kppim3g->Fill(ppimgam3.M(),weight); if (PI0_CUT(twogam3.M())){ DLorentzVector ppimpi0=twogam3+ppim_kf; LambdaPi0Mass_kppim3g->Fill(ppimpi0.M(),weight); if (SIGMA0_CUT(ppimgam3.M())){ DLorentzVector ppim3g=ppimgam3+twogam3; Sigma0Pi0Mass_kppim3g->Fill(ppim3g.M(),weight); } } } } void JEventProcessor_hyperon::KpPim4GammaAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf,kp_kf,pim_kf; vectorgammas_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -211: // pi- pim_kf=(*locParticleIterator)->Get_P4(); break; case 1: // Gamma gammas_kf.push_back((*locParticleIterator)->Get_P4()); break; default: break; } } } DLorentzVector ppim_kf=proton_kf+pim_kf; double ppim_mass=ppim_kf.M(); DLorentzVector pair1=gammas_kf[0]+gammas_kf[1]; DLorentzVector pair2=gammas_kf[3]+gammas_kf[2]; double fourgam_mass=(pair1+pair2).M(); DLorentzVector ppim4g_kf=ppim_kf+pair1+pair2; double ppim4g_mass=ppim4g_kf.M(); PPim4GammaMass_kppim4g->Fill(ppim4g_mass,weight); FourGammaMass_kppim4g->Fill(fourgam_mass,weight); PPimMass_kppim4g->Fill(ppim_mass,weight); vector> >doublets; MakeDoublets(doublets); map>>particle_pairs; map>bg_weights; map>>bg_particle_pairs; SetupMesonPairMaps(particle_pairs,bg_particle_pairs,bg_weights); for (unsigned int i=0;iFill(pair1.M(),pair2.M(),weight); MakeMesonPairs(weight,pair1,pair2,particle_pairs,bg_particle_pairs,bg_weights); } if (particle_pairs[Pi0Pi0_].size()>0){ if (LAMBDA_CUT(ppim_mass)){ LambdaPi0Pi0Mass->Fill(ppim4g_mass,weight); } for (unsigned int i=0;iFill(lambdapi01.M2(),lambdapi02.M2(),weight); } } } if (particle_pairs[Pi0Eta_].size()>0){ for (unsigned int i=0;iFill(lambdapi0.M2(),lambdaeta.M2(),weight); } } } for (unsigned int i=0;iFill(ppim4g_mass,bgweight); } DLorentzVector pi0_1=bg_particle_pairs[Pi0Pi0_][i].first; DLorentzVector pi0_2=bg_particle_pairs[Pi0Pi0_][i].second; DLorentzVector lambdapi01=ppim_kf+pi0_1; DLorentzVector lambdapi02=ppim_kf+pi0_2; if (LAMBDA_CUT(ppim_mass)){ LambdaPi0Pi0Dalitz_bg->Fill(lambdapi01.M2(),lambdapi02.M2(),bgweight); } } for (unsigned int i=0;iFill(lambdapi0.M2(),lambdaeta.M2(),bgweight); } } } void JEventProcessor_hyperon::TwoKpTwoPimAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf; vectorkps_kf; vectorpims_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kps_kf.push_back((*locParticleIterator)->Get_P4()); break; case -211: // pi- pims_kf.push_back((*locParticleIterator)->Get_P4()); break; default: break; } } } DLorentzVector kp1pim1=kps_kf[0]+pims_kf[0]; DLorentzVector kp2pim2=kps_kf[1]+pims_kf[1]; DLorentzVector kp1pim2=kps_kf[0]+pims_kf[1]; DLorentzVector kp2pim1=kps_kf[1]+pims_kf[0]; DLorentzVector twokp_twopim=kp1pim1+kp2pim2; DLorentzVector ppim1=proton_kf+pims_kf[0]; DLorentzVector ppim2=proton_kf+pims_kf[1]; DLorentzVector ppim1pim2=ppim1+pims_kf[1]; P2PimMass_2kp->Fill(ppim1pim2.M(),weight); PPimMass_2kp->Fill(ppim1.M(),weight); PPimMass_2kp->Fill(ppim2.M(),weight); if (LAMBDA_CUT(ppim1.M()) || LAMBDA_CUT(ppim2.M())){ LambdaPimMass->Fill(ppim1pim2.M(),weight); } TwoKpTwoPimMass->Fill(twokp_twopim.M(),weight); KpPim_vs_KpPim->Fill(kp1pim1.M(),kp2pim2.M(),weight); KpPim_vs_KpPim->Fill(kp2pim1.M(),kp1pim2.M(),weight); } void JEventProcessor_hyperon::TwoKpTwoPim2GammaAnalysis(DKinFitter *dKinFitter,double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf; vectorkps_kf; vectorpims_kf; vectorgammas_kf; 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_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 321: // K+ kps_kf.push_back((*locParticleIterator)->Get_P4()); break; case -211: // pi- pims_kf.push_back((*locParticleIterator)->Get_P4()); break; case 1: // Gamma gammas_kf.push_back((*locParticleIterator)->Get_P4()); default: break; } } } DLorentzVector kp1pim1=kps_kf[0]+pims_kf[0]; DLorentzVector kp2pim2=kps_kf[1]+pims_kf[1]; DLorentzVector kp1pim2=kps_kf[0]+pims_kf[1]; DLorentzVector kp2pim1=kps_kf[1]+pims_kf[0]; DLorentzVector twokp_twopim=kp1pim1+kp2pim2; DLorentzVector ppim1=proton_kf+pims_kf[0]; DLorentzVector ppim2=proton_kf+pims_kf[1]; DLorentzVector ppim1pim2=ppim1+pims_kf[1]; DLorentzVector twogam_kf=gammas_kf[0]+gammas_kf[1]; double twogam_mass_kf=twogam_kf.M(); TwoGammaMass_2kp2g->Fill(twogam_mass_kf,weight); P2PimMass_2kp->Fill(ppim1pim2.M(),weight); PPimMass_2kp2g->Fill(ppim1.M(),weight); PPimMass_2kp2g->Fill(ppim2.M(),weight); if (LAMBDA_CUT(ppim1.M()) || LAMBDA_CUT(ppim2.M())){ LambdaPimMass_2g->Fill(ppim1pim2.M(),weight); if (PI0_CUT(twogam_mass_kf)){ DLorentzVector ppim1pim2pi0=ppim1pim2+twogam_kf; LambdaPimPi0Mass->Fill(ppim1pim2pi0.M(),weight); if (XIMINUS_CUT(ppim1pim2.M())){ XiMinusPi0Mass->Fill(ppim1pim2pi0.M(),weight); } } } } // Use the PIDFOM scheme to sort tracks according to particle type void JEventProcessor_hyperon::FillParticleVectors(vector&tracks, map >&particles ) const { for (unsigned int j=0;j >probabilities; Particle_t particle_list[8]={Proton,AntiProton,KPlus,KMinus,PiPlus,PiMinus, Positron,Electron}; 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); particles[myParticle].push_back(hyp->Get_TrackTimeBased()); } } } // Create lists of photon and neutron candidates after applying some cuts on // the showers. void JEventProcessor_hyperon::FillNeutralVectors(double t0_rf,const DVector3 &vertex, vector&neutrals, vector&gammaHyps, vector&neutron_v4s, vector&neutron_covs, vector&neutron_x4s ) const{ for (unsigned int i=0;iGet_Hypothesis(Gamma); double tdiff=gamHyp->time()-t0_rf; RFGammaTimeDiff->Fill(tdiff); // Get shower info correspnding to this hypothesis const DNeutralShower *shower=gamHyp->Get_NeutralShower(); // Neutrons? double mn=ParticleMass(Neutron); double dt=gamHyp->t1()-t0_rf; double ct=dt*29.98; double beta=0.; DVector3 diff,pos; if (shower->dDetectorSystem==SYS_FCAL){ pos=(dynamic_cast(shower->dBCALFCALShower))->getPosition(); diff=pos-vertex; beta=diff.Mag()/ct; } else if (shower->dDetectorSystem==SYS_BCAL){ const DBCALShower *bcalshower=dynamic_cast(shower->dBCALFCALShower); pos.SetXYZ(bcalshower->x,bcalshower->y,bcalshower->z); diff=pos-vertex; beta=diff.Mag()/ct; } NeutralBeta->Fill(beta); if (beta0.){ double Eneut=1/sqrt(1.-beta*beta)*ParticleMass(Neutron); double pneut=beta*Eneut; DVector3 dir=diff; dir.SetMag(pneut); DLorentzVector neutron; neutron.SetVect(dir); neutron.SetT(Eneut); neutron_v4s.push_back(neutron); DLorentzVector x4; x4.SetVect(pos); x4.SetT(gamHyp->t1()); neutron_x4s.push_back(x4); double x=diff.x(); double y=diff.y(); double z=diff.z(); double s=diff.Mag(); double fac=1./sqrt(ct*ct-s*s); double fac2=fac*fac; double fac3=fac2*fac; double dpx_dx=mn*fac*(1.+x*x*fac2); double dpx_dy=mn*x*y*fac3; double dpx_dz=mn*x*z*fac3; double dpx_dt=-mn*fac3*x*ct*29.98; double dpy_dx=dpx_dy; double dpy_dz=mn*z*y*fac3; double dpy_dy=mn*fac*(1.+y*y*fac2); double dpy_dt=-mn*fac3*y*ct*29.98; double dpz_dx=dpx_dz; double dpz_dy=dpy_dz; double dpz_dz=mn*fac*(1.+z*z*fac2); double dpz_dt=-mn*fac3*z*ct*29.98; double varx=25.,vary=25.,varz=25.,vart=0.09; TMatrixFSym myMatrix(7); myMatrix(0,0)=dpx_dx*dpx_dx*varx+dpx_dy*dpx_dy*vary+dpx_dz*dpx_dz*varz +dpx_dt*dpx_dt*vart; myMatrix(1,1)=dpy_dx*dpy_dx*varx+dpy_dy*dpy_dy*vary+dpy_dz*dpy_dz*varz +dpy_dt*dpy_dt*vart; myMatrix(2,2)=dpz_dx*dpz_dx*varx+dpz_dy*dpz_dy*vary+dpz_dz*dpz_dz*varz +dpz_dt*dpz_dt*vart; myMatrix(3,3)=varx; myMatrix(4,4)=vary; myMatrix(5,5)=varz; myMatrix(6,6)=vart; neutron_covs.push_back(myMatrix); } else{ RFGammaTimeDiffNoNeutron->Fill(tdiff); // Look for good photons bool got_photon=(fabs(tdiff)lorentzMomentum(); if (shower->dDetectorSystem==SYS_FCAL){ if (gamma_v4.E()dQuality(shower->dBCALFCALShower))->getPosition(); if (pos.Perp()>FCAL_RADIAL_CUT) got_photon=false; if (fabs(pos.X())dDetectorSystem==SYS_BCAL){ if (gamma_v4.E()(shower->dBCALFCALShower); if (bcalshower->z>BCAL_Z_CUT) got_photon=false; } } if (got_photon){ gammaHyps.push_back(gamHyp); } } } } // Run the kinematic fitter requiring energy and momentum conservation void JEventProcessor_hyperon::DoKinematicFit(const DBeamPhoton *beamphoton, map >&particles, vector&gammaHyps, vector&neutron_v4s, vector&neutron_covs, DKinFitUtils_GlueX *dKinFitUtils, DKinFitter *dKinFitter) const { set> InitialParticles, FinalParticles; shared_ptrmyBeam=dKinFitUtils->Make_BeamParticle(beamphoton); shared_ptrmyTarget=dKinFitUtils->Make_TargetParticle(Proton); InitialParticles.insert(myBeam); InitialParticles.insert(myTarget); // Charged particles Particle_t particle_list[7]={KPlus,KMinus,PiPlus,PiMinus,Positron,Electron,Proton}; TLorentzVector vert; bool got_vertex_track=false; for (unsigned int k=0;k<7;k++){ vectormyParticles=particles[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; } } } // Neutral particles for (unsigned int k=0;kmyGamma=dKinFitUtils->Make_DetectedParticle(1,0,0.,vert,gammaHyps[k]->momentum(),0.,gammaHyps[k]->errorMatrix()); FinalParticles.insert(myGamma); } for (unsigned int k=0;kmyNeutron=dKinFitUtils->Make_DetectedParticle(2112,0,ParticleMass(Neutron),vert,neutron_v4s[k].Vect(),0.,make_shared(neutron_covs[k])); FinalParticles.insert(myNeutron); } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); } void JEventProcessor_hyperon::FillThrownHistos(const DMCReaction *reaction, vector&throwns, vector&tagm_geom, vector&tagh_geom) const { double mc_beam_energy=reaction->beam.energy(); double tagged_mc_beam_energy=0.; unsigned int counter=0,column=0; if (tagm_geom[0]->E_to_column(mc_beam_energy,column)){ tagged_mc_beam_energy=0.5*(tagm_geom[0]->getEhigh(column) +tagm_geom[0]->getElow(column)); } else if (tagh_geom[0]->E_to_counter(mc_beam_energy,counter)){ tagged_mc_beam_energy=0.5*(tagh_geom[0]->getEhigh(counter) +tagh_geom[0]->getElow(counter)); } if (tagged_mc_beam_energy>MIN_BEAM_E && tagged_mc_beam_energyFill(tagged_mc_beam_energy); vectorprotons,kps,etas,pi0s,gammas; for (unsigned int i=0;imech==0){ switch(throwns[i]->type){ case KPlus: kps.push_back(throwns[i]->lorentzMomentum()); break; case Proton: protons.push_back(throwns[i]->lorentzMomentum()); break; case Eta: etas.push_back(throwns[i]->lorentzMomentum()); break; case Pi0: pi0s.push_back(throwns[i]->lorentzMomentum()); break; case Gamma: gammas.push_back(throwns[i]->lorentzMomentum()); break; default: break; } } } if (kps.size()==1){ DLorentzVector beam_v4(0.,0.,tagged_mc_beam_energy,tagged_mc_beam_energy); double t_mc=(beam_v4-kps[0]).M2(); MCt->Fill(-t_mc); } } } void JEventProcessor_hyperon::Kp3Pim2PipAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight) const{ const DLorentzVector proton_kf=final_kf[Proton][0]; for (unsigned int i=0;i<3;i++){ DLorentzVector ppim=proton_kf+final_kf[PiMinus][i]; unsigned int ind1=0,ind2=0; switch(i){ case 0: ind1=1; ind2=2; break; case 1: ind1=0; ind2=2; break; case 2: ind1=0; ind2=1; break; default: break; } vector>>pippim_groups; vector>pippim_group; pippim_group.push_back(make_pair(final_kf[PiPlus][0],final_kf[PiMinus][ind1])); pippim_group.push_back(make_pair(final_kf[PiPlus][1],final_kf[PiMinus][ind2])); pippim_groups.push_back(pippim_group); pippim_group.clear(); pippim_group.push_back(make_pair(final_kf[PiPlus][0],final_kf[PiMinus][ind2])); pippim_group.push_back(make_pair(final_kf[PiPlus][1],final_kf[PiMinus][ind1])); pippim_groups.push_back(pippim_group); PPim_k2pip3pim->Fill(ppim.M(),weight); for (unsigned int k=0;kFill(pair1.M(),pair2.M(),weight); if (LAMBDA_CUT(ppim.M())){ PipPim2D_k2pip3pim_LambdaCut->Fill(pair1.M(),pair2.M(),weight); if (R>1.){ PipPim2D_k2pip3pim_Rcut->Fill(pair1.M(),pair2.M(),weight); DLorentzVector lambdapip1=ppim+pippim_groups[k][0].first; DLorentzVector lambdapim1=ppim+pippim_groups[k][0].second; LambdaPipPimDalitz_Rcut->Fill(lambdapim1.M2(),lambdapip1.M2(),weight); DLorentzVector lambdapip2=ppim+pippim_groups[k][1].first; DLorentzVector lambdapim2=ppim+pippim_groups[k][1].second; LambdaPipPimDalitz_Rcut->Fill(lambdapim2.M2(),lambdapip2.M2(),weight); } } } } } void JEventProcessor_hyperon::Kp2PimPipAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight) const{ DLorentzVector proton_kf=final_kf[Proton][0]; DLorentzVector kp_kf=final_kf[KPlus][0]; DLorentzVector pip_kf=final_kf[PiPlus][0]; DLorentzVector ppim1=proton_kf+final_kf[PiMinus][0]; DLorentzVector ppim2=proton_kf+final_kf[PiMinus][1]; DLorentzVector kpim1=kp_kf+final_kf[PiMinus][0]; DLorentzVector kpim2=kp_kf+final_kf[PiMinus][1]; DLorentzVector pippim1=pip_kf+final_kf[PiMinus][0]; DLorentzVector pippim2=pip_kf+final_kf[PiMinus][1]; PPim_vs_Kpim_with_pip->Fill(kpim1.M(),ppim2.M(),weight); PPim_vs_Kpim_with_pip->Fill(kpim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp->Fill(pippim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp->Fill(pippim1.M(),ppim2.M(),weight); if (LAMBDA_CUT(ppim1.M())){ DLorentzVector lambdapip=ppim1+pip_kf; LambdaPip_vs_Kpim->Fill(kpim2.M(),lambdapip.M(),weight); } if (LAMBDA_CUT(ppim2.M())){ DLorentzVector lambdapip=ppim2+pip_kf; LambdaPip_vs_Kpim->Fill(kpim1.M(),lambdapip.M(),weight); } } void JEventProcessor_hyperon::Kp2PimPip2GammaAnalysis(const DLorentzVector &beam_kf,map>&final_kf,double weight) const{ DLorentzVector proton_kf=final_kf[Proton][0]; DLorentzVector kp_kf=final_kf[KPlus][0]; DLoretnzVector pip_kf=final_kf[PiPlus][0]; DLorentzVector ppim1=proton_kf+final_kf[PiMinus][0]; DLorentzVector ppim2=proton_kf+final_kf[PiMinus][1]; DLorentzVector kpim1=kp_kf+final_kf[PiMinus][0]; DLorentzVector kpim2=kp_kf+final_kf[PiMinus][1]; DLorentzVector pippim1=pip_kf+final_kf[PiMinus][0]; DLorentzVector pippim2=pip_kf+final_kf[PiMinus][1]; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; TwoGamma_k2pimpip->Fill(twogam.M(),weight); PPim_vs_Kpim_with_pip2g->Fill(kpim1.M(),ppim2.M(),weight); PPim_vs_Kpim_with_pip2g->Fill(kpim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp2g->Fill(pippim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp2g->Fill(pippim1.M(),ppim2.M(),weight); if (LAMBDA_CUT(ppim1.M())){ DLorentzVector lambdapip=ppim1+pip_kf; LambdaPip_vs_Kpim_with_2g->Fill(kpim2.M(),lambdapip.M(),weight); DLorentzVector lambdapippi0=lambdapip+twogam; DLorentzVector kpimpi0=kpim2+twogam; if (PI0_CUT(twogam.M())){ LambdaPip_vs_KPimPi0->Fill(kpimpi0.M(),lambdapip.M(),weight); LambdaPipPi0_vs_KPim->Fill(kpim2.M(),lambdapippi0.M(),weight); } } if (LAMBDA_CUT(ppim2.M())){ DLorentzVector lambdapip=ppim2+pip_kf; LambdaPip_vs_Kpim_with_2g->Fill(kpim1.M(),lambdapip.M(),weight); DLorentzVector lambdapippi0=lambdapip+twogam; DLorentzVector kpimpi0=kpim1+twogam; if (PI0_CUT(twogam.M())){ LambdaPip_vs_KPimPi0->Fill(kpimpi0.M(),lambdapip.M(),weight); LambdaPipPi0_vs_KPim->Fill(kpim1.M(),lambdapippi0.M(),weight); } } } void JEventProcessor_hyperon::Kp2PimPip4GammaAnalysis(const DLorentzVector &beam_kf,map>&final_kf,double weight) const{ DLorentzVector proton_kf=final_kf[Proton][0]; DLorentzVector kp_kf=final_kf[KPlus][0]; DLoretnzVector pip_kf=final_kf[PiPlus][0]; DLorentzVector ppim1=proton_kf+final_kf[PiMinus][0]; DLorentzVector ppim2=proton_kf+final_kf[PiMinus][1]; DLorentzVector kpim1=kp_kf+final_kf[PiMinus][0]; DLorentzVector kpim2=kp_kf+final_kf[PiMinus][1]; DLorentzVector pippim1=pip_kf+final_kf[PiMinus][0]; DLorentzVector pippim2=pip_kf+final_kf[PiMinus][1]; DLorentzVector fourgam=final_kf[Gamma][0]+final_kf[Gamma][1]+final_kf[Gamma][2]+final_kf[Gamma][3]; FourGamma_k2pimpip->Fill(fourgam.M(),weight); PPim_vs_Kpim_with_pip4g->Fill(kpim1.M(),ppim2.M(),weight); PPim_vs_Kpim_with_pip4g->Fill(kpim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp4g->Fill(pippim2.M(),ppim1.M(),weight); PPim_vs_pippim_with_Kp4g->Fill(pippim1.M(),ppim2.M(),weight); MakeMesonPairs(final_kf[Gamma],weight,FourGamma2d_kppim4g); if (LAMBDA_CUT(ppim1.M())){ DLorentzVector lambdapip=ppim1+pip_kf; LambdaPip_vs_Kpim_with_4g->Fill(kpim2.M(),lambdapip.M(),weight); for (unsigned int i=0;iFill(kpimpi0_1.M(),lambdapippi0_1.M(),weight); DLorentzVector kpimpi0_2=kpim2+pi0_2; DLorentzVector lambdapippi0_2=lambdapip+pi0_1; LambdaPipPi0_vs_KPimPi0->Fill(kpimpi0_2.M(),lambdapippi0_2.M(),weight); } } if (LAMBDA_CUT(ppim2.M())){ DLorentzVector lambdapip=ppim2+pip_kf; LambdaPip_vs_Kpim_with_4g->Fill(kpim1.M(),lambdapip.M(),weight); for (unsigned int i=0;iFill(kpimpi0_1.M(),lambdapippi0_1.M(),weight); DLorentzVector kpimpi0_2=kpim1+pi0_2; DLorentzVector lambdapippi0_2=lambdapip+pi0_1; LambdaPipPi0_vs_KPimPi0->Fill(kpimpi0_2.M(),lambdapippi0_2.M(),weight); } } if (LAMBDA_CUT(ppim1.M())){ DLorentzVector lambdapip=ppim1+pip_kf; for (unsigned int i=0;iFill(kpimpi0_1.M(),lambdapippi0_1.M(),bgweight); DLorentzVector kpimpi0_2=kpim2+pi0_2; DLorentzVector lambdapippi0_2=lambdapip+pi0_1; LambdaPipPi0_vs_KPimPi0_bg->Fill(kpimpi0_2.M(),lambdapippi0_2.M(),bgweight); } } if (LAMBDA_CUT(ppim2.M())){ DLorentzVector lambdapip=ppim2+pip_kf; for (unsigned int i=0;iFill(kpimpi0_1.M(),lambdapippi0_1.M(),bgweight); DLorentzVector kpimpi0_2=kpim1+pi0_2; DLorentzVector lambdapippi0_2=lambdapip+pi0_1; LambdaPipPi0_vs_KPimPi0_bg->Fill(kpimpi0_2.M(),lambdapippi0_2.M(),bgweight); } } } void JEventProcessor_hyperon::MakeDoublets(){ 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); } void JEventProcessor_hyperon::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)); } } }