// $Id$ // // File: JEventProcessor_antiproton.cc // Created: Tue Jun 14 14:01:44 EDT 2022 // Creator: staylor (on Linux ifarm1802.jlab.org 3.10.0-1160.11.1.el7.x86_64 x86_64) // #include "JEventProcessor_antiproton.h" using namespace jana; // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_antiproton()); } } // "C" //------------------ // JEventProcessor_antiproton (Constructor) //------------------ JEventProcessor_antiproton::JEventProcessor_antiproton() { } //------------------ // ~JEventProcessor_antiproton (Destructor) //------------------ JEventProcessor_antiproton::~JEventProcessor_antiproton() { } //------------------ // init //------------------ jerror_t JEventProcessor_antiproton::init(void) { japp->CreateLock("mykinlock"); CL_CUT=0.01; gPARMS->SetDefaultParameter("ANTIPROTON:CL_CUT",CL_CUT); PROTON_GAMMA_DT_CUT=2.; gPARMS->SetDefaultParameter("ANTIPROTON:PROTON_GAMMA_DT_CUT", PROTON_GAMMA_DT_CUT); PI0_CUT_VALUE=0.025; ETA_CUT_VALUE=0.05; ETAPRIME_CUT_VALUE=0.05; NUM_SIGMA_BG=2; gPARMS->SetDefaultParameter("ANTIPROTON:NUM_SIGMA_BG",NUM_SIGMA_BG); SPLIT_CUT=0.5; gPARMS->SetDefaultParameter("ANTIPROTON:SPLIT_CUT",SPLIT_CUT); FCAL_POS_CUT=10.0; gPARMS->SetDefaultParameter("ANTIPROTON:FCAL_POS_CUT",FCAL_POS_CUT); FCAL_RADIAL_CUT=104.0; gPARMS->SetDefaultParameter("ANTIPROTON:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT); BCAL_Z_CUT=384.0; gPARMS->SetDefaultParameter("ANTIPROTON:BCAL_Z_CUT",BCAL_Z_CUT); FCAL_THRESHOLD=0.1; gPARMS->SetDefaultParameter("ANTIPROTON:FCAL_THRESHOLD",FCAL_THRESHOLD); BCAL_THRESHOLD=0.05; gPARMS->SetDefaultParameter("ANTIPROTON:BCAL_THRESHOLD",BCAL_THRESHOLD); gDirectory->mkdir("antiproton")->cd(); gDirectory->mkdir("AntiP2P")->cd(); { PantiP= new TH1F("PantiP",";M(p#bar{p}) [GeV]",1600,1.8,5.8); } gDirectory->cd("../"); gDirectory->mkdir("AntiP2P1G")->cd(); { PantiP_with_1g= new TH1F("PantiP_with_1g",";M(p#bar{p}) [GeV]",1600,1.8,5.8); PantiPGamma= new TH1F("PantiPGamma",";M(p#bar{p}#gamma) [GeV]",1600,1.8,5.8); } gDirectory->cd("../"); gDirectory->mkdir("AntiP2P2G")->cd(); { PantiP_with_2g= new TH1F("PantiP_with_2g",";M(p#bar{p}) [GeV]",1600,1.8,5.8); TwoGamma_with_ppap=new TH1F("TwoGamma_with_ppap",";M(2#gamma)[GeV]",1600,0,4); PantiPPi0= new TH1F("PantiPPi0",";M(p#bar{p}#pi^{0}) [GeV]",1600,1.8,5.8); PantiPEta= new TH1F("PantiPEta",";M(p#bar{p}#eta) [GeV]",1600,1.8,5.8); } gDirectory->cd("../"); gDirectory->mkdir("PipPimAntiP2P")->cd(); { PPiMinusVsAntiPPiPlus_kf=new TH2F("PPiMinusVsAntiPPiPlus_kf", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,2.9,400,0.9,2.9); LambdaLambdaBarMass= new TH1F("LambdaLambdaBarMass", ";M(#Lambda#bar{Lambda}) [GeV]",1600,2.,6.); } gDirectory->cd("../"); gDirectory->mkdir("PipPimAntiP2P1G")->cd(); { PPiMinusVsAntiPPiPlus_with_1g=new TH2F("PPiMinusVsAntiPPiPlus_with_1g", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,2.9,400,0.9,2.9); Sigma0AntiLambdaMass=new TH1F("Sigma0AntiLambdaMass", ";M(#Sigma^{0}#bar{#Lambda}) [GeV]", 1600,1.,5.); AntiSigma0LambdaMass=new TH1F("AntiSigma0LambdaMass", ";M(#bar{#Sigma^{0}}#Lambda) [GeV]", 1600,1.,5.); AntiLambdaGammaMass=new TH1F("AntiLambdaGammaMass", ";M(#bar{#Lambda}#gamma) [GeV]",1600,1.,5.); LambdaGammaMass_with_AntiLambda=new TH1F("LambdaGammaMass_with_AntiLambda", ";M(#Lambda#gamma) [GeV]", 1600,1.,5.); CosThetaAntiLambdaRest_1g=new TH1F("CosThetaAntiLambdaRest_1g", ";cos(#theta(#bar{p} #gamma))",20,-1,1); CosThetaLambdaRest_1g=new TH1F("CosThetaLambdaRest_1g", ";cos(#theta(p #gamma))",20,-1,1); PPiMinusVsAntiPPiPlus_with_1g_NU=new TH2F("PPiMinusVsAntiPPiPlus_with_1g_NU", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,2.9,400,0.9,2.9); Sigma0AntiLambdaMass_NU=new TH1F("Sigma0AntiLambdaMass_NU", ";M(#Sigma^{0}#bar{#Lambda}) [GeV]", 1600,1.,5.); AntiSigma0LambdaMass_NU=new TH1F("AntiSigma0LambdaMass_NU", ";M(#bar{#Sigma^{0}}#Lambda) [GeV]", 1600,1.,5.); AntiLambdaGammaMass_NU=new TH1F("AntiLambdaGammaMass_NU", ";M(#bar{#Lambda}#gamma) [GeV]",1600,1.,5.); LambdaGammaMass_with_AntiLambda_NU=new TH1F("LambdaGammaMass_with_AntiLambda_NU", ";M(#Lambda#gamma) [GeV]", 1600,1.,5.); CosThetaAntiLambdaRest_1g_NU=new TH1F("CosThetaAntiLambdaRest_1g_NU", ";cos(#theta(#bar{p} #gamma))",20,-1,1); CosThetaLambdaRest_1g_NU=new TH1F("CosThetaLambdaRest_1g_NU", ";cos(#theta(p #gamma))",20,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("PipPimAntiP2P2G")->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(#bar{p} #gamma))",20,-1,1); CosThetaLambdaRest=new TH1F("CosThetaLambdaRest", ";cos(#theta(p #gamma))",20,-1,1); PPim_vs_aPPip_NU=new TH2F("PPim_vs_aPPip_NU", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,1.9,400,0.9,1.9); AntiLambdaGamma_vs_LambdaGamma_NU=new TH2F("AntiLambdaGamma_vs_LambdaGamma_NU", ";M(#Lambda#gamma) [GeV];M(#bar{#Lambda}#gamma) [GeV]", 400,0.9,1.9,400,0.9,1.9); Sigma0AntiSigma0_NU=new TH1F("Sigma0AntiSigma0_NU", ";M(#Sigma^{0}#bar{#Sigma}^{0}) [GeV]", 800,2.,4.); CosThetaAntiLambdaRest_NU=new TH1F("CosThetaAntiLambdaRest_NU", ";cos(#theta(#bar{p} #gamma))",20,-1,1); CosThetaLambdaRest_NU=new TH1F("CosThetaLambdaRest_NU", ";cos(#theta(p #gamma))",20,-1,1); } gDirectory->cd("../"); gDirectory->mkdir("TwoPip2PimAntiP2P")->cd(); { PPiMinusVsAntiPPiPlus_2pip2pim=new TH2F("PPiMinusVsAntiPPiPlus_2pip2pim", ";M(#bar{p}#pi^{+}) [GeV];M(p#pi^{-}) [GeV]", 400,0.9,2.9,400,0.9,2.9); P2PiMinusVsAntiP2PiPlus=new TH2F("P2PiMinusVsAntiP2PiPlus_2", ";M(#bar{#Lambda}#pi^{+}) [GeV];M(#Lambda#pi^{-}) [GeV]",400,0.9,2.9,400,0.9,2.9); XiAntiXi=new TH1F("XiAntiXi",";M(#Xi#bar{#Xi}) [GeV]",400,2.5,4.5); } gDirectory->cd("../"); gDirectory->cd("../"); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_antiproton::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_antiproton::evnt(JEventLoop *loop, uint64_t eventnumber) { vectortracks; loop->Get(tracks); if (tracks.size()==0) return RESOURCE_UNAVAILABLE; // Assign charge track by PID map>chargedParticles; FillChargedParticleVectors(tracks,chargedParticles); if (chargedParticles[AntiProton].size()==0) return RESOURCE_UNAVAILABLE; vectorbeamphotons; loop->Get(beamphotons); if (beamphotons.size()==0) return RESOURCE_UNAVAILABLE; vectorneutrals; loop->Get(neutrals); japp->WriteLock("myantilock"); // Find t0 (at "vertex") for event and assign list of neutral particles double t0_rf=tracks[0]->Get_BestTrackingFOM()->t0(); DVector3 vertex=tracks[0]->Get_BestTrackingFOM()->position(); // Assign neutrals by ID double unused_energy=0.; vectorgammaParticles; FillGammaParticleVector(unused_energy,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;itime()-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(); if (CL>CL_CUT){ DLorentzVector beam_kf; map>final_kf; GetKF4vectors(dKinFitter,beam_kf,final_kf); // Particle counts 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 (numProton==2&&numAntiProton==1){ if (numPiPlus==0&&numPiMinus==0){ switch(numGamma){ case 0: AntiProton2ProtonAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 1: AntiProton2Proton1GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 2: AntiProton2Proton2GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } if (numPiPlus==1&&numPiMinus==1){ switch(numGamma){ case 0: PipPimAntiProton2ProtonAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 1: PipPimAntiProton2Proton1GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; case 2: PipPimAntiProton2Proton2GammaAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } else if (numPiPlus==2&&numPiMinus==2){ switch(numGamma){ case 0: TwoPip2PimAntiProton2ProtonAnalysis(unused_energy,beam_kf,final_kf,weight); break; default: break; } } } } // CL cut } // loop over beam photons //delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; japp->Unlock("myantilock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_antiproton::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_antiproton::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_antiproton::FillChargedParticleVectors(vector&tracks,map>&chargedParticles) const { Particle_t particle_list[4]={PiPlus,PiMinus,Proton,AntiProton}; for (unsigned int j=0;j >probabilities; for (unsigned int i=0;i<4;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()); } } } // Get list of final state photons // Create lists of photon candidates after applying some cuts on the showers. void JEventProcessor_antiproton::FillGammaParticleVector(double &unused_energy, double t0_rf, const DVector3 &vertex, vector&neutrals, vector&gammaParticles ) const{ for (unsigned int i=0;iGet_Hypothesis(Gamma); double tdiff=gamHyp->time()-t0_rf; DLorentzVector gamma_v4=gamHyp->lorentzMomentum(); // Get shower info correspnding to this hypothesis const DNeutralShower *shower=gamHyp->Get_NeutralShower(); // Look for good photons bool got_photon=(fabs(tdiff)dDetectorSystem==SYS_FCAL){ if (gamma_v4.E()dQuality(shower->dBCALFCALShower))->getPosition(); pos-=mFCALCenter; 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){ gammaParticles.push_back(gamHyp); } else if (fabs(tdiff) >&chargedParticles, vector&gammaParticles, 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[4]={PiPlus,PiMinus,Proton,AntiProton}; for (unsigned int k=0;k<4;k++){ vectormyParticles=chargedParticles[particle_list[k]]; for (unsigned int m=0;mmyParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]); FinalParticles.insert(myParticle); } } // Final state photons for (unsigned int k=0;kmyGamma=dKinFitUtils->Make_DetectedParticle(22,0,0.,vert4,gammaHyp->momentum(),0.,gammaHyp->errorMatrix()); FinalParticles.insert(myGamma); } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); } void JEventProcessor_antiproton::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()); } } } //------------------------------------------------------------------------------ // AntiProton 2Proton Ngamma //------------------------------------------------------------------------------ void JEventProcessor_antiproton::AntiProton2ProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector antip_p1=final_kf[AntiProton][0]+final_kf[Proton][0]; DLorentzVector antip_p2=final_kf[AntiProton][0]+final_kf[Proton][1]; PantiP->Fill(antip_p1.M(),weight); PantiP->Fill(antip_p2.M(),weight); } void JEventProcessor_antiproton::AntiProton2Proton1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector antip_p1=final_kf[AntiProton][0]+final_kf[Proton][0]; DLorentzVector antip_p2=final_kf[AntiProton][0]+final_kf[Proton][1]; DLorentzVector antip_p1_g=antip_p1+final_kf[Gamma][0]; DLorentzVector antip_p2_g=antip_p2+final_kf[Gamma][0]; PantiP_with_1g->Fill(antip_p1.M(),weight); PantiP_with_1g->Fill(antip_p2.M(),weight); PantiPGamma->Fill(antip_p1_g.M(),weight); PantiPGamma->Fill(antip_p2_g.M(),weight); } void JEventProcessor_antiproton::AntiProton2Proton2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector antip_p1=final_kf[AntiProton][0]+final_kf[Proton][0]; DLorentzVector antip_p2=final_kf[AntiProton][0]+final_kf[Proton][1]; DLorentzVector twogam=final_kf[Gamma][0]+final_kf[Gamma][1]; double twogam_mass_kf=twogam.M(); DLorentzVector antip_p1_2g=antip_p1+twogam; DLorentzVector antip_p2_2g=antip_p2+twogam; TwoGamma_with_ppap->Fill(twogam_mass_kf,weight); PantiP_with_2g->Fill(antip_p1.M(),weight); PantiP_with_2g->Fill(antip_p2.M(),weight); if (PI0_CUT(twogam_mass_kf,PI0_CUT_VALUE)){ PantiPPi0->Fill(antip_p2_2g.M(),weight); PantiPPi0->Fill(antip_p2_2g.M(),weight); } if (ETA_CUT(twogam_mass_kf,ETA_CUT_VALUE)){ PantiPEta->Fill(antip_p2_2g.M(),weight); PantiPEta->Fill(antip_p2_2g.M(),weight); } } //------------------------------------------------------------------------------ // AntiProton 2Proton pi+ pi- Ngamma //------------------------------------------------------------------------------ void JEventProcessor_antiproton::PipPimAntiProton2ProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector proton1_pim=final_kf[PiMinus][0]+final_kf[Proton][0]; DLorentzVector proton2_pim=final_kf[PiMinus][0]+final_kf[Proton][1]; DLorentzVector antiproton_pip=final_kf[PiPlus][0]+final_kf[AntiProton][0]; PPiMinusVsAntiPPiPlus_kf->Fill(antiproton_pip.M(),proton1_pim.M(),weight); PPiMinusVsAntiPPiPlus_kf->Fill(antiproton_pip.M(),proton2_pim.M(),weight); if (LAMBDA_CUT(antiproton_pip.M())){ if (LAMBDA_CUT(proton1_pim.M())){ DLorentzVector lamlambar=antiproton_pip+proton1_pim; LambdaLambdaBarMass->Fill(lamlambar.M(),weight); } if (LAMBDA_CUT(proton2_pim.M())){ DLorentzVector lamlambar=antiproton_pip+proton2_pim; LambdaLambdaBarMass->Fill(lamlambar.M(),weight); } } } void JEventProcessor_antiproton::PipPimAntiProton2Proton1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector gamma_kf=final_kf[Gamma][0]; DLorentzVector proton1_pim=final_kf[PiMinus][0]+final_kf[Proton][0]; DLorentzVector proton2_pim=final_kf[PiMinus][0]+final_kf[Proton][1]; DLorentzVector antiproton_pip=final_kf[PiPlus][0]+final_kf[AntiProton][0]; PPiMinusVsAntiPPiPlus_with_1g->Fill(antiproton_pip.M(),proton1_pim.M(),weight); PPiMinusVsAntiPPiPlus_with_1g->Fill(antiproton_pip.M(),proton2_pim.M(),weight); bool got_lambda1=LAMBDA_CUT(proton1_pim.M()); bool got_lambda2=LAMBDA_CUT(proton2_pim.M()); if (LAMBDA_CUT(antiproton_pip.M()) && (got_lambda1 || got_lambda2)){ DLorentzVector AntiLambdaGamma=antiproton_pip+gamma_kf; DVector3 beta_antilambda=(1./(antiproton_pip.E()))*antiproton_pip.Vect(); DLorentzVector antiproton_in_antilambda_rest=final_kf[AntiProton][0]; antiproton_in_antilambda_rest.Boost(-beta_antilambda); AntiLambdaGammaMass->Fill(AntiLambdaGamma.M(),weight); if (SIGMA0_CUT(AntiLambdaGamma.M())){ DLorentzVector antisiglam1=AntiLambdaGamma+proton1_pim; DLorentzVector antisiglam2=AntiLambdaGamma+proton2_pim; if (LAMBDA_CUT(proton1_pim.M())){ AntiSigma0LambdaMass->Fill(antisiglam1.M(),weight); } if (LAMBDA_CUT(proton2_pim.M())){ AntiSigma0LambdaMass->Fill(antisiglam2.M(),weight); } DLorentzVector gamma_in_antilambda_rest=gamma_kf; gamma_in_antilambda_rest.Boost(-beta_antilambda); DVector3 pvec_antiproton=antiproton_in_antilambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_antilambda_rest.Vect(); double costheta_gamma_antip=pvec_gamma1.Dot(pvec_antiproton) /(pvec_gamma1.Mag()*pvec_antiproton.Mag()); CosThetaAntiLambdaRest_1g->Fill(costheta_gamma_antip,weight); } if (got_lambda1){ DLorentzVector LambdaGamma=proton1_pim+gamma_kf; LambdaGammaMass_with_AntiLambda->Fill(LambdaGamma.M(),weight); DVector3 beta_lambda=(1./(proton1_pim.E()))*proton1_pim.Vect(); DLorentzVector proton_in_lambda_rest=final_kf[Proton][0]; proton_in_lambda_rest.Boost(-beta_lambda); if (SIGMA0_CUT(LambdaGamma.M())){ DLorentzVector sigantilam=LambdaGamma+antiproton_pip; Sigma0AntiLambdaMass->Fill(sigantilam.M(),weight); DLorentzVector gamma_in_lambda_rest=gamma_kf; gamma_in_lambda_rest.Boost(-beta_lambda); DVector3 pvec_proton=proton_in_lambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_lambda_rest.Vect(); double costheta_gamma_p=pvec_gamma1.Dot(pvec_proton) /(pvec_gamma1.Mag()*pvec_proton.Mag()); CosThetaLambdaRest_1g->Fill(costheta_gamma_p,weight); } } if (got_lambda2){ DLorentzVector LambdaGamma=proton2_pim+gamma_kf; LambdaGammaMass_with_AntiLambda->Fill(LambdaGamma.M(),weight); DVector3 beta_lambda=(1./(proton2_pim.E()))*proton2_pim.Vect(); DLorentzVector proton_in_lambda_rest=final_kf[Proton][0]; proton_in_lambda_rest.Boost(-beta_lambda); if (SIGMA0_CUT(LambdaGamma.M())){ DLorentzVector sigantilam=LambdaGamma+antiproton_pip; Sigma0AntiLambdaMass->Fill(sigantilam.M(),weight); DLorentzVector gamma_in_lambda_rest=gamma_kf; gamma_in_lambda_rest.Boost(-beta_lambda); DVector3 pvec_proton=proton_in_lambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_lambda_rest.Vect(); double costheta_gamma_p=pvec_gamma1.Dot(pvec_proton) /(pvec_gamma1.Mag()*pvec_proton.Mag()); CosThetaLambdaRest_1g->Fill(costheta_gamma_p,weight); } } } if (unused_energy<0.01){ PPiMinusVsAntiPPiPlus_with_1g_NU->Fill(antiproton_pip.M(),proton1_pim.M(),weight); PPiMinusVsAntiPPiPlus_with_1g_NU->Fill(antiproton_pip.M(),proton2_pim.M(),weight); if (LAMBDA_CUT(antiproton_pip.M()) && (got_lambda1 || got_lambda2)){ DLorentzVector AntiLambdaGamma=antiproton_pip+gamma_kf; DVector3 beta_antilambda=(1./(antiproton_pip.E()))*antiproton_pip.Vect(); DLorentzVector antiproton_in_antilambda_rest=final_kf[AntiProton][0]; antiproton_in_antilambda_rest.Boost(-beta_antilambda); AntiLambdaGammaMass_NU->Fill(AntiLambdaGamma.M(),weight); if (SIGMA0_CUT(AntiLambdaGamma.M())){ DLorentzVector antisiglam1=AntiLambdaGamma+proton1_pim; DLorentzVector antisiglam2=AntiLambdaGamma+proton2_pim; if (LAMBDA_CUT(proton1_pim.M())){ AntiSigma0LambdaMass_NU->Fill(antisiglam1.M(),weight); } if (LAMBDA_CUT(proton2_pim.M())){ AntiSigma0LambdaMass_NU->Fill(antisiglam2.M(),weight); } DLorentzVector gamma_in_antilambda_rest=gamma_kf; gamma_in_antilambda_rest.Boost(-beta_antilambda); DVector3 pvec_antiproton=antiproton_in_antilambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_antilambda_rest.Vect(); double costheta_gamma_antip=pvec_gamma1.Dot(pvec_antiproton) /(pvec_gamma1.Mag()*pvec_antiproton.Mag()); CosThetaAntiLambdaRest_1g_NU->Fill(costheta_gamma_antip,weight); } if (got_lambda1){ DLorentzVector LambdaGamma=proton1_pim+gamma_kf; LambdaGammaMass_with_AntiLambda_NU->Fill(LambdaGamma.M(),weight); DVector3 beta_lambda=(1./(proton1_pim.E()))*proton1_pim.Vect(); DLorentzVector proton_in_lambda_rest=final_kf[Proton][0]; proton_in_lambda_rest.Boost(-beta_lambda); if (SIGMA0_CUT(LambdaGamma.M())){ DLorentzVector sigantilam=LambdaGamma+antiproton_pip; Sigma0AntiLambdaMass_NU->Fill(sigantilam.M(),weight); DLorentzVector gamma_in_lambda_rest=gamma_kf; gamma_in_lambda_rest.Boost(-beta_lambda); DVector3 pvec_proton=proton_in_lambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_lambda_rest.Vect(); double costheta_gamma_p=pvec_gamma1.Dot(pvec_proton) /(pvec_gamma1.Mag()*pvec_proton.Mag()); CosThetaLambdaRest_1g_NU->Fill(costheta_gamma_p,weight); } } if (got_lambda2){ DLorentzVector LambdaGamma=proton2_pim+gamma_kf; LambdaGammaMass_with_AntiLambda_NU->Fill(LambdaGamma.M(),weight); DVector3 beta_lambda=(1./(proton2_pim.E()))*proton2_pim.Vect(); DLorentzVector proton_in_lambda_rest=final_kf[Proton][0]; proton_in_lambda_rest.Boost(-beta_lambda); if (SIGMA0_CUT(LambdaGamma.M())){ DLorentzVector sigantilam=LambdaGamma+antiproton_pip; Sigma0AntiLambdaMass_NU->Fill(sigantilam.M(),weight); DLorentzVector gamma_in_lambda_rest=gamma_kf; gamma_in_lambda_rest.Boost(-beta_lambda); DVector3 pvec_proton=proton_in_lambda_rest.Vect(); DVector3 pvec_gamma1=gamma_in_lambda_rest.Vect(); double costheta_gamma_p=pvec_gamma1.Dot(pvec_proton) /(pvec_gamma1.Mag()*pvec_proton.Mag()); CosThetaLambdaRest_1g_NU->Fill(costheta_gamma_p,weight); } } } } } void JEventProcessor_antiproton::PipPimAntiProton2Proton2GammaAnalysis(double unused_energy,const 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 (unused_energy<0.01){ PPim_vs_aPPip_NU->Fill(appip_mass,ppim_mass,weight); } if (LAMBDA_CUT(appip_mass)&&LAMBDA_CUT(ppim_mass)){ FillSigma0Histos(unused_energy,ppim,appip,gamma1,gamma2,proton,antiproton,weight); FillSigma0Histos(unused_energy,ppim,appip,gamma2,gamma1,proton,antiproton,weight); } } } void JEventProcessor_antiproton::TwoPip2PimAntiProton2ProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf, map>&final_kf, double weight) const{ DLorentzVector proton1_pim1=final_kf[PiMinus][0]+final_kf[Proton][0]; DLorentzVector proton2_pim1=final_kf[PiMinus][0]+final_kf[Proton][1]; DLorentzVector proton1_pim2=final_kf[PiMinus][1]+final_kf[Proton][0]; DLorentzVector proton2_pim2=final_kf[PiMinus][1]+final_kf[Proton][1]; DLorentzVector antiproton_pip1=final_kf[PiPlus][0]+final_kf[AntiProton][0]; DLorentzVector antiproton_pip2=final_kf[PiPlus][1]+final_kf[AntiProton][0]; DLorentzVector antiproton_2pip=antiproton_pip1+final_kf[PiPlus][1]; DLorentzVector proton1_2pim=proton1_pim1+final_kf[PiMinus][1]; DLorentzVector proton2_2pim=proton2_pim1+final_kf[PiMinus][1]; PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip1.M(),proton1_pim1.M(), weight); if (LAMBDA_CUT(antiproton_pip1.M()) && LAMBDA_CUT(proton1_pim1.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton1_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton1_2pim.M()>1.3 && proton1_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton1_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip1.M(),proton2_pim1.M(), weight); if (LAMBDA_CUT(antiproton_pip1.M()) && LAMBDA_CUT(proton2_pim1.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton2_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton2_2pim.M()>1.3 && proton2_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton2_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip2.M(),proton1_pim1.M(), weight); if (LAMBDA_CUT(antiproton_pip2.M()) && LAMBDA_CUT(proton1_pim1.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton1_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton1_2pim.M()>1.3 && proton1_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton1_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip2.M(),proton2_pim1.M(), weight); if (LAMBDA_CUT(antiproton_pip2.M()) && LAMBDA_CUT(proton2_pim1.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton2_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton2_2pim.M()>1.3 && proton2_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton2_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip1.M(),proton1_pim2.M(), weight); if (LAMBDA_CUT(antiproton_pip1.M()) && LAMBDA_CUT(proton1_pim2.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton1_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton1_2pim.M()>1.3 && proton1_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton1_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip1.M(),proton2_pim2.M(), weight); if (LAMBDA_CUT(antiproton_pip1.M()) && LAMBDA_CUT(proton2_pim2.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton2_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton2_2pim.M()>1.3 && proton2_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton2_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip2.M(),proton1_pim2.M(), weight); if (LAMBDA_CUT(antiproton_pip2.M()) && LAMBDA_CUT(proton1_pim2.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton1_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton1_2pim.M()>1.3 && proton1_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton1_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } PPiMinusVsAntiPPiPlus_2pip2pim->Fill(antiproton_pip2.M(),proton2_pim2.M(), weight); if (LAMBDA_CUT(antiproton_pip2.M()) && LAMBDA_CUT(proton2_pim2.M())){ P2PiMinusVsAntiP2PiPlus->Fill(antiproton_2pip.M(),proton2_2pim.M(),weight); if (antiproton_2pip.M()>1.3 && antiproton_2pip.M()<1.345 && proton2_2pim.M()>1.3 && proton2_2pim.M()<1.345){ DLorentzVector xiantixi=antiproton_2pip+proton2_2pim; XiAntiXi->Fill(xiantixi.M(),weight); } } } void JEventProcessor_antiproton::FillSigma0Histos(double unused_energy,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 (unused_energy<0.01){ AntiLambdaGamma_vs_LambdaGamma_NU->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); } if (unused_energy<0.01){ if (SIGMA0_CUT(alg_mass) && SIGMA0_CUT(lg_mass)){ DLorentzVector lg_alg=lg+alg; double lg_alg_mass=lg_alg.M(); Sigma0AntiSigma0_NU->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_NU->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_NU->Fill(cosTheta_pg,weight); } } }