// $Id$ // // File: JEventProcessor_EpEm.cc // Created: Tue Oct 21 16:50:20 EDT 2014 // Creator: staylor (on Linux gluon05.jlab.org 2.6.32-358.18.1.el6.x86_64 x86_64) // #include "JEventProcessor_EpEm.h" using namespace jana; #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #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_EpEm()); } } // "C" //------------------ // JEventProcessor_EpEm (Constructor) //------------------ JEventProcessor_EpEm::JEventProcessor_EpEm() { } //------------------ // ~JEventProcessor_EpEm (Destructor) //------------------ JEventProcessor_EpEm::~JEventProcessor_EpEm() { } //------------------ // init //------------------ jerror_t JEventProcessor_EpEm::init(void) { japp->CreateLock("myEpEmXlock"); BETA_CUT=0.9; PI0_CUT_VALUE=0.025; ETA_CUT_VALUE=0.05; ETAPRIME_CUT_VALUE=0.05; NUM_SIGMA_BG=2; EMULATE_TRIGGER=true; gPARMS->SetDefaultParameter("EPEM:EMULATE_TRIGGER",EMULATE_TRIGGER); VETO_GAP_EVENTS=true; gPARMS->SetDefaultParameter("EPEM:VETO_GAP_EVENTS",VETO_GAP_EVENTS); INSERT_ONLY=false; gPARMS->SetDefaultParameter("EPEM:INSERT_ONLY",INSERT_ONLY); FCAL_ONLY=false; gPARMS->SetDefaultParameter("EPEM:FCAL_ONLY",FCAL_ONLY); DROP_ONE_CELL_BCAL_SHOWERS=true; gPARMS->SetDefaultParameter("EPEM:DROP_ONE_CELL_BCAL_SHOWERS",DROP_ONE_CELL_BCAL_SHOWERS); PROTON_GAMMA_DT_CUT=2.; gPARMS->SetDefaultParameter("EPEM:PROTON_GAMMA_DT_CUT", PROTON_GAMMA_DT_CUT); PIM_CL_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:PIM_CL_CUT",PIM_CL_CUT); PIP_CL_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:PIP_CL_CUT",PIM_CL_CUT); PROTON_CL_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:PROTON_CL_CUT",PROTON_CL_CUT); ELECTRON_CL_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:ELECTRON_CL_CUT",ELECTRON_CL_CUT); POSITRON_CL_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:POSITRON_CL_CUT",POSITRON_CL_CUT); BCAL_DTRACK_CUT=0.0; gPARMS->SetDefaultParameter("EPEM:BCAL_DTRACK_CUT",BCAL_DTRACK_CUT); SPLIT_CUT=0.5; gPARMS->SetDefaultParameter("EPEM:SPLIT_CUT",SPLIT_CUT); PI0_VETO_CUT=0.042; // 3 sigma? gPARMS->SetDefaultParameter("EPEM:PI0_VETO_CUT",PI0_VETO_CUT); DEBUG_LEVEL=0; gPARMS->SetDefaultParameter("EPEM:DEBUG_LEVEL",DEBUG_LEVEL); CL_CUT=0.01; gPARMS->SetDefaultParameter("EPEM:CL_CUT",CL_CUT); TRACK_CL_CUT=0.; gPARMS->SetDefaultParameter("EPEM:TRACK_CL_CUT",TRACK_CL_CUT); FCAL_POS_CUT=10.0; gPARMS->SetDefaultParameter("EPEM:FCAL_POS_CUT",FCAL_POS_CUT); FCAL_RADIAL_CUT=104.0; gPARMS->SetDefaultParameter("EPEM:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT); BCAL_Z_CUT=384.0; gPARMS->SetDefaultParameter("EPEM:BCAL_Z_CUT",BCAL_Z_CUT); FCAL_THRESHOLD=0.1; gPARMS->SetDefaultParameter("EPEM:FCAL_THRESHOLD",FCAL_THRESHOLD); BCAL_THRESHOLD=0.05; gPARMS->SetDefaultParameter("EPEM:BCAL_THRESHOLD",BCAL_THRESHOLD); MIN_BEAM_E=0.0; gPARMS->SetDefaultParameter("EPEM:MIN_BEAM_E",MIN_BEAM_E); MAX_BEAM_E=12.0; gPARMS->SetDefaultParameter("EPEM:MAX_BEAM_E",MAX_BEAM_E); MIN_CLUSTERS_IN_INSERT=0; gPARMS->SetDefaultParameter("EPEM:MIN_CLUSTERS_IN_INSERT", MIN_CLUSTERS_IN_INSERT); UNUSED_ENERGY_CUT=0.01; gPARMS->SetDefaultParameter("EPEM:UNUSED_ENERGY_CUT",UNUSED_ENERGY_CUT); FILL_ROO_DATASET=false; gPARMS->SetDefaultParameter("EPEM:FILL_ROO_DATASET", FILL_ROO_DATASET); // japp->RootWriteLock(); /* if (FILL_ROO_DATASET){ RooMass=new RooRealVar("m","",0.,3.,"GeV"); RooT=new RooRealVar("t","",0.,3.,"GeV^{-2}"); RooE=new RooRealVar("E","",3.,12.,"GeV"); RooWeight=new RooRealVar("weight","",-1.,1.,""); RooData=new RooDataSet("MassData","MassData",RooArgSet(*RooMass,*RooT,*RooE,*RooWeight)); } */ // ... fill historgrams or trees ... gDirectory->mkdir("EpEm")->cd(); TrigMask=new TH1F("TrigMask","Trigger bits",15,-0.5,14.5); gDirectory->mkdir("MC")->cd(); { MCProtonThetaVsP=new TH2F("MCProtonThetaVsP","Thrown #theta vs p for proton",200,0,2,90,0,90); MCProtonThetaVsP->Sumw2(); MCProtonThetaVsP->SetXTitle("p [GeV/c]"); MCProtonThetaVsP->SetYTitle("#theta [degrees]"); MCProtonAcceptedP=new TH1F("MCProtonAcceptedP","Accepted proton",200,0,2); MCProtonAcceptedP->Sumw2(); MCBeam = new TH1F("MCBeam","Photon energy",48,2.8,12.4); MCBeamTagged = new TH1F("MCBeamTagged","Photon energy",48,2.8,12.4); // MCProtonAcceptance=new TH1F("MCProtonAcceptance","proton acceptance", //200,0,2); MCProtonP=new TH1F("MCProtonP","thrown proton distribution", 200,0,2); MCProtonP->Sumw2(); MCt=new TH1F("MCt","t distribution",1000,0,2.); MCt_vs_E=new TH2F("MCt_vs_E","t dist",27,6.4,11.8,250,0,5); MCt_vs_Ecoarse=new TH2F("MCt_vs_Ecoarse","t dist",8,7.,11.0,50,0,2); MCdeltat=new TH2F("MCdeltat","t residual",100,-2,0,100,-1,1); MCMissingMass =new TH1F("MCMissingMass","mX",1600,0,4.); MCMissingMassVsEbeam=new TH2F("MCMissingMassVsEbeam", "MX vs E#gamma", 48,2.8,12.4,1600,0,4.); //48,2.8,12.4,300,0,3); } gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); { FCAL2ShowerDistance=new TH1F("FCAL2ShowerDistance","distance between showers", 500,0,50); ProtonGammaTimeDiff=new TH1F("ProtonGammaTimeDiff","#deltat between proton and #gamma at vertex",200,-10,10); BeamPhotonTimeDiff=new TH2F("BeamPhotonTimeDiff","t(tag)-t(rf) at vertex", 48,2.8,12.4,1000,-50,50); UnusedEnergy=new TH1F("UnusedEnergy","unused shower energy",100,0,0.5); 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); } gDirectory->cd("../"); gDirectory->mkdir("e+e-")->cd(); { EplusEminusMass_kf=new TH1F("EplusEminusMass_kf",";M(e^{+}e^{-}) [GeV]",1600,0,4.); EpEmDOCA=new TH1F("EpEmDOCA",";DOCA [cm]",1000,0,10); EpEmMassVsR=new TH2F("EpEmMassVsR",";R [cm];M(e^{+}e^{-}) [GeV]",100,0,20, 1000,0.9,3.9); EplusEminusMass_Rcut=new TH1F("EplusEminusMass_Rcut",";M(e^{+}e^{-}) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("e+e-gamma")->cd(); { EplusEminus1GammaMass_kf=new TH1F("EplusEminus1GammaMass_kf",";M(e^{+}e^{-}#gamma) [GeV]",1600,0,4.); EplusEminusMass_1g=new TH1F("EplusEminusMass_1g","M(e^{+}e^{-}) [GeV]",1600,0,4.); EpEmDOCA_1g=new TH1F("EpEmDOCA_1g",";DOCA [cm]",1000,0,10); EpEmMassVsR_1g=new TH2F("EpEmMassVsR_1g",";R [cm];M(e^{+}e^{-}) [GeV]",100,0,20, 1000,0.9,3.9); EplusEminus1GammaMass_Rcut=new TH1F("EplusEminus1GammaMass_Rcut",";M(e^{+}e^{-}#gamma) [GeV]",1600,0,4.); EplusEminusMass_1g_Rcut=new TH1F("EplusEminusMass_1g_Rcut","M(e^{+}e^{-}) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("e+e-2gamma")->cd(); { EplusEminus2GammaMass_kf=new TH1F("EplusEminus2GammaMass_kf","e+e-2#gamma",1600,0,4.); TwoGammaMass_eplus_eminus=new TH1F("TwoGammaMass_eplus_eminus","M(2#gamma)",1600,0,4.); EplusEminusPi0Mass_kf=new TH1F("EplusEminusPi0Mass_kf","e+e-#pi^0",1600,0,4.); EplusEminus1GammaMass_with_2g=new TH1F("EplusEminus1GammaMass_with_2g","e+e-#gamma",1600,0,4.); EplusEminusEtaMass_kf=new TH1F("EplusEminusEtaMass_kf","e+e-#eta",1600,0,4.); EplusEminusPi0Mass_cut=new TH1F("EplusEminusPi0Mass_cut","e+e-#pi^0",1600,0,4.); EplusEminusEtaMass_cut=new TH1F("EplusEminusEtaMass_cut","e+e-#eta",1600,0,4.); EplusEminusMass_with_2g=new TH1F("EplusEminusMass_with_2g","e+e-",1600,0,4.); EplusEminusMass_with_pi0=new TH1F("EplusEminusMass_with_pi0","e+e-",1600,0,4.); EplusEminusMass_with_eta=new TH1F("EplusEminusMass_with_eta","e+e-",1600,0,4.); EtaGamma_epem2g=new TH1F("EtaGamma_epem2g","e+e-2#gamma",1600,0,4.); Pi0Gamma_epem2g=new TH1F("Pi0Gamma_epem2g","e+e-2#gamma",1600,0,4.); EtaGamma_epem2g_cut=new TH1F("EtaGamma_epem2g_cut","e+e-2#gamma",1600,0,4.); Pi0Gamma_epem2g_cut=new TH1F("Pi0Gamma_epem2g_cut","e+e-2#gamma",1600,0,4.); OmegaPi0Mass_epem=new TH1F("OmegaPi0Mass_epem","e+e-2#gamma",1600,0,4.); EplusEminusMass_with_pi0_Rcut=new TH1F("EplusEminusMass_with_pi0_Rcut","e+e-",1600,0,4.); EplusEminusMass_with_eta_Rcut=new TH1F("EplusEminusMass_with_eta_Rcut","e+e-",1600,0,4.); EplusEminusEtaMass_Rcut=new TH1F("EplusEminusEtaMass_Rcut","e+e-#eta",1600,0,4.); EplusEminusPi0Mass_Rcut=new TH1F("EplusEminusPi0Mass_Rcut","e+e-#pi^0",1600,0,4.); TwoGammaMass_eplus_eminus_Rcut=new TH1F("TwoGammaMass_eplus_eminus_Rcut","M(2#gamma)",1600,0,4.); EplusEminusMass_with_2g_Rcut=new TH1F("EplusEminusMass_with_2g_Rcut","e+e-",1600,0,4.); EpEmDOCA_2g=new TH1F("EpEmDOCA_2g",";DOCA [cm]",1000,0,10); EpEmMassVsR_2g=new TH2F("EpEmMassVsR_2g",";R [cm];M(e^{+}e^{-}) [GeV]",100,0,20, 1000,0.9,3.9); } gDirectory->cd("../"); gDirectory->mkdir("e+e-3gamma")->cd(); { EplusEminusPi0GammaMass_kf=new TH1F("EplusEminusPi0GammaMass_kf","e+e-3#gamma",1600,0,4.); TwoGammaMass_epem3g=new TH1F("TwoGammaMass_epem3g","M(2#gamma)",1600,0,4.); Pi0Gamma_with_epem=new TH1F("Pi0Gamma_with_epem","M(3#gamma)",1600,0,4.); EplusEminusMass_with_3g=new TH1F("EplusEminusMass_with_3g","e+e-",1600,0,4.); EplusEminus3GammaMass_kf=new TH1F("EplusEminus3GammaMass_kf","e+e-3#gamma",1600,0,4.); EplusEminus3GammaMass_Rcut=new TH1F("EplusEminus3GammaMass_Rcut","e+e-3#gamma",1600,0,4.); EplusEminusMass_with_3g_Rcut=new TH1F("EplusEminusMass_with_3g_Rcut","e+e-",1600,0,4.); EpEmDOCA_3g=new TH1F("EpEmDOCA_3g",";DOCA [cm]",1000,0,10); } gDirectory->cd("../"); gDirectory->mkdir("e+e-4gamma")->cd(); { EplusEminusMass_with_4g=new TH1F("EplusEminusMass_with_4g","e+e-",1600,0,4.); EplusEminus4GammaMass_kf=new TH1F("EplusEminus4GammaMass_kf","e+e-4#gamma",1600,0,4.); FourGamma2d_epem=new TH2F("FourGamma2d_epem",";M(2#gamma,pair2) [GeV];M(2#gamma,pair1) [GeV]",500,0,1,500,0,1); FourGammaMass_epem=new TH1F("FourGammaMass_epem","; M(4#gamma) [GeV]",1600,0,4.); TwoPi0Mass_epem=new TH1F("TwoPi0Mass_epem","; M(2#pi^{0}) [GeV]",1600,0,4.); EplusEminus4GammaMass_Rcut=new TH1F("EplusEminus4GammaMass_Rcut","e+e-4#gamma",1600,0,4.); EplusEminusMass_with_4g_Rcut=new TH1F("EplusEminusMass_with_4g_Rcut","e+e-",1600,0,4.); EpEmDOCA_4g=new TH1F("EpEmDOCA_4g",";DOCA [cm]",1000,0,10); } gDirectory->cd("../"); gDirectory->mkdir("e+e-pi+pi-")->cd(); { EplusEminusPiplusPiminusMass=new TH1F("EplusEminusPiplusPiminusMass","e+e-pi+ppi-",1600,0,4.); EplusEminusPiplusPiminusMass_cut=new TH1F("EplusEminusPiplusPiminusMass_cut","e+e-pi+ppi-",1600,0,4.); EplusEminusMass_PipPim=new TH1F("EplusEminusMass_PipPim","e+e-",1600,0,4.); PiPlusPiMinusMass_epem=new TH1F("PiPlusPiMinusMass_epem","#pi+#pi-",1600,0,4.); EplusEminusMass_PipPim_etaprime_cut=new TH1F("EplusEminusMass_PipPim_etaprime_cut","e+e-, #eta' cut",1600,0,4.); EplusEminusMass_PipPim_eta_cut=new TH1F("EplusEminusMass_PipPim_eta_cut","e+e-, #eta cut",1600,0,4.); EtaPrime_sinphicosphi=new TH1F("EtaPrime_sinphicosphi","sin(#phi)cos(#phi)",20,-0.5,0.5); Eta_sinphicosphi=new TH1F("Eta_sinphicosphi","sin(#phi)cos(#phi)",20,-0.5,0.5); EpEmDOCA_pippim=new TH1F("EpEmDOCA_pippim",";DOCA [cm]",1000,0,10); EpEmMassVsR_pippim=new TH2F("EpEmMassVsR_pippim",";R [cm];M(e^{+}e^{-}) [GeV]",100,0,20, 1000,0.9,3.9); EplusEminusPiplusPiminusMass_Rcut=new TH1F("EplusEminusPiplusPiminusMass_Rcut","M(e^{+}e^{-}#pi^{+}#pi^{-}) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("2e+2e-")->cd(); { TwoEplusTwoEminusMass_kf=new TH1F("TwoEplusTwoEminusMass_kf","2e+2e-",1600,0,4.); EpEm_vs_EpEm=new TH2F("EpEm_vs_EpEm",";M(e^{+}e^{-}, pair 1) [GeV];M(e^+}e^{-}, pair 2) [GEV]",400,0,1,400,0,1); EpEm_vs_EpEm_Rcut=new TH2F("EpEm_vs_EpEm_Rcut",";M(e^{+}e^{-}, pair 1) [GeV];M(e^+}e^{-}, pair 2) [GEV]",400,0,1,400,0,1); EpEmDOCA2=new TH2F("EpEmDOCA2",";DOCA, pair1 [cm]; DOCA, pair 2 [cm]", 1000,0,10,1000,0,10); TwoEplusTwoEminusMass_Rcut=new TH1F("TwoEplusTwoEminusMass_Rcut", ";M(2e^{+}2e^{-}) [GeV]",1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("ne+e-pi+")->cd(); { EplusEminusMass_npip=new TH1F("EplusEminusMass_npip",";M(e^{+}e^{-}) [GeV]",1600,0,4.); Nepem_pip=new TH1F("Nepem_pip",";M(ne^{+}e^{-}) [GeV]",1600,0.9,4.9); EpEmPip_n=new TH1F("EpEmPip_n",";M(#pi^{+}e^{+}e^{-}) [GeV]",1600,0.,4.); Npip_epem=new TH1F("Npip_epem",";M(n#pi^{+}) [GeV]",1600,0.9,4.9); EpEmDOCA_n=new TH1F("EpEmDOCA_n",";DOCA [cm]",1000,0,10); EpEmMassVsR_n=new TH2F("EpEmMassVsR_n",";R [cm];M(e^{+}e^{-}) [GeV]",100,0,20, 1000,0.9,3.9); EplusEminusMass_npip_Rcut=new TH1F("EplusEminusMass_npip_Rcut",";M(e^{+}e^{-}) [GeV]",1600,0,4.); Nepem_pip_Rcut=new TH1F("Nepem_pip_Rcut",";M(ne^{+}e^{-}) [GeV]",1600,0.9,4.9); EpEmPip_n_Rcut=new TH1F("EpEmPip_n_Rcut",";M(#pi^{+}e^{+}e^{-}) [GeV]",1600,0.,4.); Npip_epem_Rcut=new TH1F("Npip_epem_Rcut",";M(n#pi^{+}) [GeV]",1600,0.9,4.9); } gDirectory->cd("../"); gDirectory->mkdir("KpPimEpEm")->cd(); { PPimMass_with_Kepem= new TH1F("PPimMass_with_Kepem",";M(p#pi-) [GeV]",1600,0,4.); LambdaEpEm= new TH1F("LambdaEpEm",";M#Lambdae+e-) [GeV]",1600,0,4.); EpEmMass_with_Kppim=new TH1F("EpEmMass_with_Kppim",";M(e+e-) [GeV]", 1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpPimEpEm1G")->cd(); { PPimMass_with_Kepem1g= new TH1F("PPimMass_with_Kepem1g",";M(p#pi-) [GeV]",1600,0,4.); LambdaEpEm_1g= new TH1F("LambdaEpEm_1g",";M#Lambdae+e-) [GeV]",1600,0,4.); EpEmMass_with_Kppim1g=new TH1F("EpEmMass_with_Kppim1g",";M(e+e-) [GeV]", 1600,0,4.); EpEm1GMass_with_Kppim=new TH1F("EpEm1GMass_with_Kppim",";M(e+e-#gamma) [GeV]", 1600,0,4.); } gDirectory->cd("../"); gDirectory->mkdir("KpPimEpEm2G")->cd(); { PPimMass_with_Kepem2g= new TH1F("PPimMass_with_Kepem2g",";M(p#pi-) [GeV]",1600,0,4.); LambdaEpEm_2g= new TH1F("LambdaEpEm_2g",";M#Lambdae+e-) [GeV]",1600,0,4.); EpEmMass_with_Kppim2g=new TH1F("EpEmMass_with_Kppim2g",";M(e+e-) [GeV]", 1600,0,4.); TwoGMass_with_Kppimepem=new TH1F("TwoGMass_with_Kppimepem",";M(2#gamma) [GeV]", 1600,0,4.); } gDirectory->cd("../"); BeamEnergy=new TH1F("BeamEnergy","E(beam)",48,2.8,12.4); gDirectory->cd("../"); japp->RootUnLock(); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_EpEm::brun(JEventLoop *eventLoop, int32_t runnumber) { DApplication* dapp = dynamic_cast(eventLoop->GetJApplication()); dgeom = dapp->GetDGeometry(runnumber); bfield = dapp->GetBfield(runnumber); double x0,y0,z0; dgeom->GetFCALPosition(x0,y0,z0); mFCALCenter.SetXYZ(x0,y0,z0); //const DMagneticFieldMap* locMagneticFieldMap = dapp->GetBfield(runnumber); //return NOERROR; vector pid_algorithms; eventLoop->Get(pid_algorithms); if(pid_algorithms.size()<1){ _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<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->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); //return NOERROR; static int thrown_count=0; // static int recon_count=0; /* japp->WriteLock("mylock"); dKinFitUtils->Reset_NewEvent(loop->GetJEvent().GetEventNumber()); dKinFitter->Reset_NewEvent(); japp->Unlock("mylock"); */ vectortagm_geom; loop->Get(tagm_geom); vectortagh_geom; loop->Get(tagh_geom); vectorneutrals; loop->Get(neutrals); vectorbeamphotons; loop->Get(beamphotons); vectortracks; loop->Get(tracks); vectorthrowns; loop->Get(throwns); vectorschits; loop->Get(schits); const DMCReaction* reaction=NULL; vector triggers; if (throwns.size()>0){ loop->GetSingle(reaction); //loop->Get(triggers); } japp->RootWriteLock(); double mc_beam_energy=0; double tagged_mc_beam_energy=0.; if (reaction!=NULL){ mc_beam_energy=reaction->beam.energy(); MCBeam->Fill(mc_beam_energy); if (tagm_geom.size()&&tagh_geom.size()){ unsigned int counter=0,column=0; if ( tagm_geom[0]->E_to_column(mc_beam_energy,column)){ double Egam=0.5*(tagm_geom[0]->getEhigh(column) +tagm_geom[0]->getElow(column)); MCBeamTagged->Fill(Egam); tagged_mc_beam_energy=Egam; // cout << Egam-mc_beam_energy << endl; } else if (tagh_geom[0]->E_to_counter(mc_beam_energy,counter)){ double Egam=0.5*(tagh_geom[0]->getEhigh(counter) +tagh_geom[0]->getElow(counter)); //cout << Egam << " " << mc_beam_energy << endl; MCBeamTagged->Fill(Egam); tagged_mc_beam_energy=Egam; //cout << Egam-mc_beam_energy << endl; } } } double diff_min=1e6; int itrue=-1; for (unsigned int i=0;ilorentzMomentum().E(); BeamEnergy->Fill(beam_E); //cout << beamphotons[i]->lorentzMomentum().E() << endl; if (mc_beam_energy>0. && fabs(mc_beam_energy-beam_E)-1){ if (diff_min<0.015){ //cout << " .... " << beamphotons[itrue]->lorentzMomentum().E()- tagged_mc_beam_energy << endl; //MCBeamTagged->Fill(beamphotons[itrue]->lorentzMomentum().E()); } } double mc_mom=0.; double t_mc=0.; vectormcprotons,mcetas,mcgammas,mcpi0s,mcKps,mcKms; if (tagged_mc_beam_energy>0.){ for (unsigned int i=0;imech==0){ switch(throwns[i]->type){ case Proton: { mcprotons.push_back(throwns[i]->lorentzMomentum()); DVector3 mom=throwns[i]->momentum(); MCProtonThetaVsP->Fill(mom.Mag(),180./M_PI*mom.Theta()); MCProtonP->Fill(mom.Mag()); thrown_count++; mc_mom=mom.Mag(); DLorentzVector target(0,0,0,ParticleMass(Proton)); t_mc=(throwns[i]->lorentzMomentum()-target).M2(); MCt->Fill(-t_mc); double myBeam=reaction->beam.energy(); MCt_vs_E->Fill(myBeam,-t_mc); if (myBeam>7. && myBeam<11.){ MCt_vs_Ecoarse->Fill(myBeam,-t_mc); } DLorentzVector beam(0,0,myBeam,myBeam); DLorentzVector missing=target+beam-throwns[i]->lorentzMomentum(); MCMissingMass->Fill(missing.M()); MCMissingMassVsEbeam->Fill(myBeam,missing.M()); // printf("Thrown: p %f count %d\n",mom.Mag(),thrown_count); } 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; case KPlus: mcKps.push_back(throwns[i]->lorentzMomentum()); break; case KMinus: mcKms.push_back(throwns[i]->lorentzMomentum()); break; default: break; } } } } if (EMULATE_TRIGGER) { loop->Get(triggers); if (triggers.size()>0){ int trig_mask=triggers[0]->Get_L1TriggerBits(); // cout << "trigger: " << trig_mask <Fill(trig_mask); if (!((trig_mask&0x1)||(trig_mask&0x4))){ japp->RootUnLock(); //cout << "Not a valid trigger: " << trig_mask <RootUnLock(); // Assign charge track by PID map>chargedParticles; FillChargedParticleVectors(tracks,chargedParticles); if (chargedParticles[Electron].size()==0) return RESOURCE_UNAVAILABLE; if (chargedParticles[Positron].size()==0) return RESOURCE_UNAVAILABLE; japp->WriteLock("myEpEmXlock"); if (chargedParticles[Proton].size()==1&&throwns.size()>0&&mc_mom>0.){ MCProtonAcceptedP->Fill(mc_mom); } // 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 map>neutralParticles; FillNeutralParticleVectors(t0_rf,vertex,neutrals,neutralParticles); // 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; 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(); KinFitCL->Fill(CL); if (CL>CL_CUT){ DLorentzVector beam_kf; map>final_kf; GetKF4vectors(dKinFitter,beam_kf,final_kf); // Particle counts unsigned int numElectron=final_kf[Electron].size(); unsigned int numPositron=final_kf[Positron].size(); unsigned int numPiPlus=final_kf[PiPlus].size(); unsigned int numPiMinus=final_kf[PiMinus].size(); unsigned int numKPlus=final_kf[KPlus].size(); unsigned int numKMinus=final_kf[KMinus].size(); unsigned int numGamma=final_kf[Gamma].size(); unsigned int numNeutron=final_kf[Neutron].size(); unsigned int numProton=final_kf[Proton].size(); unsigned int numAntiProton=final_kf[AntiProton].size(); if (numProton==1 && numAntiProton==0 && numNeutron==0){ if (numPositron==1 && numElectron==1){ // No kaons if (numKPlus==0 && numKMinus==0){ // no extra charge pions if (numPiPlus==0 && numPiMinus==0){ switch(numGamma){ case 0: EplusEminusAnalysis(beam_kf,final_kf,weight); break; case 1: EplusEminusGammaAnalysis(beam_kf,final_kf,weight); break; case 2: EplusEminus2GammaAnalysis(beam_kf,final_kf,weight); break; case 3: EplusEminus3GammaAnalysis(beam_kf,final_kf,weight); break; case 4: EplusEminus4GammaAnalysis(beam_kf,final_kf,weight); break; default: break; } } // with pi+pi- if (numPiMinus==1 && numPiPlus==1){ switch(numGamma){ case 0: EplusEminusPiplusPiminusAnalysis(beam_kf,final_kf,weight); break; default: break; } } } // no kaons // With kaon if (numKPlus==1 && numPiMinus==1 && numKMinus==0 && numProton==1 && numPiPlus==0 && numNeutron==0){ switch(numGamma){ case 0: KpPimEpEmAnalysis(beam_kf,final_kf,weight); break; case 1: KpPimEpEm1GAnalysis(beam_kf,final_kf,weight); break; case 2: KpPimEpEm2GAnalysis(beam_kf,final_kf,weight); break; default: break; } } } // 2e+ 2e- if (numPositron==2 && numElectron==2 && numPiPlus==0 && numPiMinus==0 && numKPlus==0 && numKMinus==0){ if (numGamma==0){ TwoEplusTwoEminusAnalysis(beam_kf,final_kf,weight); } } } if (numProton==0 && numNeutron==1 && numAntiProton==0){ if (numPositron==1 && numElectron==1){ if (numKPlus==0 && numKMinus==0 && numPiPlus==1 && numPiMinus==0){ NEpEmPipAnalysis(beam_kf,final_kf,weight); } } } } } //delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; japp->Unlock("myEpEmXlock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_EpEm::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_EpEm::fini(void) { /* if (FILL_ROO_DATASET){ RooData->Write(); } */ // Called before program exit after event processing is finished. return NOERROR; } void JEventProcessor_EpEm::EplusEminus2GammaAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector epem2g_kf=epem_kf+twogam; double epem2g_mass_kf=epem2g_kf.M(); EplusEminus2GammaMass_kf->Fill(epem2g_mass_kf,weight); TwoGammaMass_eplus_eminus->Fill(twogam.M(),weight); EplusEminusMass_with_2g->Fill(epem_kf.M(),weight); if (OMEGA_CUT(epem_kf.M())){ if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ OmegaPi0Mass_epem->Fill(epem2g_mass_kf,weight); } } if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ EplusEminusMass_with_pi0->Fill(epem_kf.M(),weight); } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ EplusEminusMass_with_eta->Fill(epem_kf.M(),weight); } DLorentzVector epemg_1=epem_kf+final_kf[Gamma][0]; DLorentzVector epemg_2=epem_kf+final_kf[Gamma][1]; EplusEminus1GammaMass_with_2g->Fill(epemg_1.M(),weight); if (PI0_CUT(epemg_1.M(),PI0_CUT_VALUE)){ Pi0Gamma_epem2g->Fill(epem2g_mass_kf,weight); } if (ETA_CUT(epemg_1.M(),ETA_CUT_VALUE)){ EtaGamma_epem2g->Fill(epem2g_mass_kf,weight); } EplusEminus1GammaMass_with_2g->Fill(epemg_2.M(),weight); if (PI0_CUT(epemg_2.M(),PI0_CUT_VALUE)){ Pi0Gamma_epem2g->Fill(epem2g_mass_kf,weight); } if (ETA_CUT(epemg_2.M(),ETA_CUT_VALUE)){ EtaGamma_epem2g->Fill(epem2g_mass_kf,weight); } if (PI0_ANTI_CUT(twogam.M(),PI0_VETO_CUT)){ if (PI0_CUT(epemg_2.M(),PI0_CUT_VALUE)){ Pi0Gamma_epem2g_cut->Fill(epem2g_mass_kf,weight); } if (ETA_CUT(epemg_2.M(),ETA_CUT_VALUE)){ EtaGamma_epem2g_cut->Fill(epem2g_mass_kf,weight); } if (PI0_CUT(epemg_1.M(),PI0_CUT_VALUE)){ Pi0Gamma_epem2g_cut->Fill(epem2g_mass_kf,weight); } if (ETA_CUT(epemg_1.M(),ETA_CUT_VALUE)){ EtaGamma_epem2g_cut->Fill(epem2g_mass_kf,weight); } } if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ EplusEminusPi0Mass_kf->Fill(epem2g_kf.M(),weight); } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ EplusEminusEtaMass_kf->Fill(epem2g_kf.M(),weight); } if (epem_kf.M()>0.07){ if (PI0_CUT(twogam.M(),PI0_CUT_VALUE)){ EplusEminusPi0Mass_cut->Fill(epem2g_kf.M(),weight); } if (ETA_CUT(twogam.M(),ETA_CUT_VALUE)){ EplusEminusEtaMass_cut->Fill(epem2g_kf.M(),weight); } } } void JEventProcessor_EpEm::EplusEminus3GammaAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector threegam=final_kf[Gamma][0]+final_kf[Gamma][1]+final_kf[Gamma][2]; DLorentzVector epem3g_kf=epem_kf+threegam; double epem3g_mass_kf=epem3g_kf.M(); EplusEminus3GammaMass_kf->Fill(epem3g_mass_kf,weight); EplusEminusMass_with_3g->Fill(epem_kf.M(),weight); DLorentzVector pair1=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector pair2=final_kf[Gamma][1]+final_kf[Gamma][2]; DLorentzVector pair3=final_kf[Gamma][0]+final_kf[Gamma][2]; TwoGammaMass_epem3g->Fill(pair1.M(),weight); if (PI0_CUT(pair1.M(),PI0_CUT_VALUE)&&PI0_ANTI_CUT(pair2.M(),PI0_VETO_CUT) && PI0_ANTI_CUT(pair3.M(),PI0_VETO_CUT)){ Pi0Gamma_with_epem->Fill(threegam.M(),weight); EplusEminusPi0GammaMass_kf->Fill(epem3g_mass_kf,weight); } TwoGammaMass_epem3g->Fill(pair2.M(),weight); if (PI0_CUT(pair2.M(),PI0_CUT_VALUE)&&PI0_ANTI_CUT(pair1.M(),PI0_VETO_CUT) && PI0_ANTI_CUT(pair3.M(),PI0_VETO_CUT)){ Pi0Gamma_with_epem->Fill(threegam.M(),weight); EplusEminusPi0GammaMass_kf->Fill(epem3g_mass_kf,weight); } TwoGammaMass_epem3g->Fill(pair3.M(),weight); if (PI0_CUT(pair3.M(),PI0_CUT_VALUE)&&PI0_ANTI_CUT(pair2.M(),PI0_VETO_CUT) && PI0_ANTI_CUT(pair1.M(),PI0_VETO_CUT)){ Pi0Gamma_with_epem->Fill(threegam.M(),weight); EplusEminusPi0GammaMass_kf->Fill(epem3g_mass_kf,weight); } } void JEventProcessor_EpEm::EplusEminusGammaAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector epem1g_kf=epem_kf+final_kf[Gamma][0]; double epem1g_mass_kf=epem1g_kf.M(); EplusEminus1GammaMass_kf->Fill(epem1g_mass_kf,weight); EplusEminusMass_1g->Fill(epem_kf.M(),weight); } void JEventProcessor_EpEm::EplusEminusAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; EplusEminusMass_kf->Fill(epem_kf.M(),weight); } void JEventProcessor_EpEm::EplusEminus4GammaAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector proton_kf=final_kf[Proton][0]; DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector fourgam=final_kf[Gamma][0]+final_kf[Gamma][1]+final_kf[Gamma][2]+final_kf[Gamma][3]; double fourgam_mass=fourgam.M(); FourGammaMass_epem->Fill(fourgam_mass,weight); DLorentzVector epem4g_kf=epem_kf+fourgam; double epem4g_mass_kf=epem4g_kf.M(); EplusEminus4GammaMass_kf->Fill(epem4g_mass_kf,weight); EplusEminusMass_with_4g->Fill(epem_kf.M(),weight); vector> >doublets; MakeDoublets(doublets); bool got_meson_pairs=false; map>>particle_pairs; 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); vectorpaired_doublets(doublets.size()); for (unsigned int i=0;iFill(pair1_mass,pair2_mass,weight); if (PI0_CUT(pair1_mass,PI0_CUT_VALUE)&&PI0_CUT(pair2_mass,PI0_CUT_VALUE)){ particle_pairs[Pi0Pi0_].push_back(make_pair(pair1,pair2)); paired_doublets[i]=1; got_meson_pairs=true; } if (ETA_CUT(pair1_mass,ETA_CUT_VALUE) && ETA_CUT(pair2_mass,ETA_CUT_VALUE)){ particle_pairs[EtaEta_].push_back(make_pair(pair1,pair2)); paired_doublets[i]=1; got_meson_pairs=true; } if (ETA_CUT(pair1_mass,ETA_CUT_VALUE) && PI0_CUT(pair2_mass,PI0_CUT_VALUE)){ particle_pairs[Pi0Eta_].push_back(make_pair(pair2,pair1)); paired_doublets[i]=1; got_meson_pairs=true; } if (ETA_CUT(pair2_mass,ETA_CUT_VALUE) && PI0_CUT(pair1_mass,PI0_CUT_VALUE)){ particle_pairs[Pi0Eta_].push_back(make_pair(pair1,pair2)); paired_doublets[i]=1; got_meson_pairs=true; } if (PI0_CUT(pair1_mass,PI0_CUT_VALUE) && ETAPRIME_CUT(pair2_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair1,pair2)); paired_doublets[i]=1; got_meson_pairs=true; } if (PI0_CUT(pair2_mass,PI0_CUT_VALUE) && ETAPRIME_CUT(pair1_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[Pi0EtaPrime_].push_back(make_pair(pair2,pair1)); paired_doublets[i]=1; got_meson_pairs=true; } if (ETA_CUT(pair1_mass,ETA_CUT_VALUE) && ETAPRIME_CUT(pair2_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[EtaEtaPrime_].push_back(make_pair(pair1,pair2)); paired_doublets[i]=1; got_meson_pairs=true; } if (ETA_CUT(pair2_mass,ETA_CUT_VALUE) && ETAPRIME_CUT(pair1_mass,ETAPRIME_CUT_VALUE)){ particle_pairs[EtaEtaPrime_].push_back(make_pair(pair2,pair1)); paired_doublets[i]=1; got_meson_pairs=true; } } if (particle_pairs[Pi0Pi0_].size()>0){ TwoPi0Mass_epem->Fill(fourgam_mass,weight); } } void JEventProcessor_EpEm::TwoEplusTwoEminusAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector twoeptwoem_kf=final_kf[Electron][0]+final_kf[Positron][0]+final_kf[Electron][1]+final_kf[Positron][1]; TwoEplusTwoEminusMass_kf->Fill(twoeptwoem_kf.M(),weight); DLorentzVector pair11=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector pair22=final_kf[Electron][1]+final_kf[Positron][1]; DLorentzVector pair12=final_kf[Electron][0]+final_kf[Positron][1]; DLorentzVector pair21=final_kf[Electron][1]+final_kf[Positron][0]; EpEm_vs_EpEm->Fill(pair11.M(),pair22.M(),weight); EpEm_vs_EpEm->Fill(pair12.M(),pair21.M(),weight); } void JEventProcessor_EpEm::EplusEminusPiplusPiminusAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector pippim=final_kf[PiMinus][0]+final_kf[PiPlus][0]; DLorentzVector pippim_epem=epem_kf+pippim; EplusEminusPiplusPiminusMass->Fill(pippim_epem.M(),weight); EplusEminusMass_PipPim->Fill(epem_kf.M(),weight); PiPlusPiMinusMass_epem->Fill(pippim.M(),weight); DVector3 beta=(1./pippim_epem.E())*pippim_epem.Vect(); DLorentzVector pip_kf=final_kf[PiPlus][0]; pip_kf.Boost(-beta); DLorentzVector ep_kf=final_kf[Positron][0]; ep_kf.Boost(-beta); DLorentzVector pim_kf=final_kf[PiMinus][0]; pim_kf.Boost(-beta); DLorentzVector em_kf=final_kf[Electron][0]; em_kf.Boost(-beta); DVector3 epemhat=(em_kf.Vect()).Cross(ep_kf.Vect()); epemhat.SetMag(1.); DVector3 pippimhat=(pim_kf.Vect()).Cross(pip_kf.Vect()); pippimhat.SetMag(1.); epem_kf.Boost(-beta); DVector3 zhat=epem_kf.Vect(); zhat.SetMag(1.); double sinphicosphi=(epemhat.Cross(pippimhat)).Dot(zhat)*epemhat.Dot(pippimhat); // (epem_kf+pippim).Print(); if (pippim.M()>0.65 && pippim.M()<0.85){ EplusEminusPiplusPiminusMass_cut->Fill(pippim_epem.M(),weight); } if (ETAPRIME_CUT(pippim_epem.M(),ETAPRIME_CUT_VALUE)){ EtaPrime_sinphicosphi->Fill(sinphicosphi,weight); EplusEminusMass_PipPim_etaprime_cut->Fill(epem_kf.M(),weight); } if (ETA_CUT(pippim_epem.M(),ETA_CUT_VALUE)){ Eta_sinphicosphi->Fill(sinphicosphi,weight); EplusEminusMass_PipPim_eta_cut->Fill(epem_kf.M(),weight); } } // Use the PIDFOM scheme to sort tracks according to particle type void JEventProcessor_EpEm::FillChargedParticleVectors(vector&tracks,map>&chargedParticles) const { Particle_t particle_list[8]={Electron,Positron,Proton,AntiProton,KPlus,KMinus,PiPlus,PiMinus}; 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()); } } } // Sort neutral particles into lists of photons and neutrons void JEventProcessor_EpEm::FillNeutralParticleVectors(double t0_rf, const DVector3 &vertex, vector&neutrals, map>&neutralParticles) const { for (unsigned int i=0;iGet_Hypothesis(Gamma); // Get shower info correspnding to this hypothesis const DNeutralShower *shower=gamHyp->Get_NeutralShower(); // Compute beta double ct=29.98*(gamHyp->t1()-t0_rf); DVector3 diff; if (shower->dDetectorSystem==SYS_FCAL){ const DFCALShower *fcalshower=dynamic_cast(shower->dBCALFCALShower); DVector3 pos=fcalshower->getPosition(); diff=pos-vertex; } 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; 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())dDetectorSystem==SYS_BCAL){ if (E(shower->dBCALFCALShower); if (bcalshower->z>BCAL_Z_CUT) continue; } neutralParticles[Gamma].push_back(gamHyp); } } } } void JEventProcessor_EpEm::GetKF4vectors(DKinFitter *dKinFitter, DLorentzVector &beam_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_DetectedParticle){ Particle_t myPID=PDGtoPType((*locParticleIterator)->Get_PID()); final_kf[myPID].push_back((*locParticleIterator)->Get_P4()); } } } // Run the kinematic fitter requiring energy and momentum conservation void JEventProcessor_EpEm::DoKinematicFit(double t0_rf,const DVector3 &vertex, const DBeamPhoton *beamphoton, map >&chargedParticles, map>&neutralParticles, DKinFitUtils_GlueX *dKinFitUtils, DKinFitter *dKinFitter) const { set> InitialParticles, FinalParticles; dKinFitter->Reset_NewFit(); shared_ptrmyBeam=dKinFitUtils->Make_BeamParticle(beamphoton); shared_ptrmyTarget=dKinFitUtils->Make_TargetParticle(Proton); InitialParticles.insert(myBeam); InitialParticles.insert(myTarget); // Vertex info DLorentzVector vert4; vert4.SetVect(vertex); vert4.SetT(t0_rf); // Charged particles Particle_t particle_list[8]={Electron,Positron,Proton,AntiProton,PiPlus,PiMinus,KPlus,KMinus}; 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); } } // Final state photons for (unsigned int k=0;kmyGamma=dKinFitUtils->Make_DetectedParticle(22,0,0.,vert4,gammaHyp->momentum(),0.,gammaHyp->errorMatrix()); FinalParticles.insert(myGamma); } // final state neutrons for (unsigned int k=0;kmyNeutron=dKinFitUtils->Make_DetectedParticle(2112,0,ParticleMass(Neutron),vert4,neutronHyp->momentum(),0.,neutronHyp->errorMatrix()); 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_EpEm::NEpEmPipAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight){ DLorentzVector epem_kf=final_kf[Electron][0]+final_kf[Positron][0]; DLorentzVector neutron_pip=final_kf[Neutron][0]+final_kf[PiPlus][0]; DLorentzVector epemn=final_kf[Neutron][0]+epem_kf; DLorentzVector epempip=epem_kf+final_kf[PiPlus][0]; Npip_epem->Fill(neutron_pip.M(),weight); Nepem_pip->Fill(epemn.M(),weight); EpEmPip_n->Fill(epempip.M(),weight); EplusEminusMass_npip->Fill(epem_kf.M(),weight); } //---------------------------------------------------------------------------- // K+ p pi- e+ e- Ngamma //---------------------------------------------------------------------------- void JEventProcessor_EpEm::KpPimEpEmAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight) const { DLorentzVector epem_kf=final_kf[Positron][0]+final_kf[Electron][0]; DLorentzVector ppim_kf=final_kf[Proton][0]+final_kf[PiMinus][0]; DLorentzVector ppimepem=ppim_kf+epem_kf; PPimMass_with_Kepem->Fill(ppim_kf.M(),weight); EpEmMass_with_Kppim->Fill(epem_kf.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ LambdaEpEm->Fill(ppimepem.M(),weight); } } void JEventProcessor_EpEm::KpPimEpEm1GAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight) const { DLorentzVector epem_kf=final_kf[Positron][0]+final_kf[Electron][0]; DLorentzVector epemg=epem_kf+final_kf[Gamma][0]; DLorentzVector ppim_kf=final_kf[Proton][0]+final_kf[PiMinus][0]; DLorentzVector ppimepem=ppim_kf+epem_kf; PPimMass_with_Kepem1g->Fill(ppim_kf.M(),weight); EpEmMass_with_Kppim1g->Fill(epem_kf.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ LambdaEpEm_1g->Fill(ppimepem.M(),weight); } EpEm1GMass_with_Kppim->Fill(epemg.M(),weight); } void JEventProcessor_EpEm::KpPimEpEm2GAnalysis(const DLorentzVector &beam_kf, map>&final_kf,double weight) const { DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; DLorentzVector epem_kf=final_kf[Positron][0]+final_kf[Electron][0]; DLorentzVector ppim_kf=final_kf[Proton][0]+final_kf[PiMinus][0]; DLorentzVector ppimepem=ppim_kf+epem_kf; PPimMass_with_Kepem2g->Fill(ppim_kf.M(),weight); EpEmMass_with_Kppim2g->Fill(epem_kf.M(),weight); if (LAMBDA_CUT(ppim_kf.M())){ LambdaEpEm_2g->Fill(ppimepem.M(),weight); } TwoGMass_with_Kppimepem->Fill(twogam.M(),weight); } void JEventProcessor_EpEm::MakeDoublets(vector > >&doublets) const { 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); }