// $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"

// Routine used to create our JEventProcessor
#include <JANA/JApplication.h>
#include <JANA/JFactory.h>
extern "C"{
  void InitPlugin(JApplication *app){
    InitJANAPlugin(app);
    app->Add(new JEventProcessor_antiproton());
  }
} // "C"

//------------------
// ~JEventProcessor_antiproton (Destructor)
//------------------
JEventProcessor_antiproton::~JEventProcessor_antiproton()
{

}

//------------------
// Init
//------------------
void JEventProcessor_antiproton::Init(void)
{
  auto app = GetApplication();
  
  CL_CUT=0.01;
  app->SetDefaultParameter("ANTIPROTON:CL_CUT",CL_CUT);   

  PROTON_GAMMA_DT_CUT=2.; 
  app->SetDefaultParameter("ANTIPROTON:PROTON_GAMMA_DT_CUT",
			      PROTON_GAMMA_DT_CUT);
  MIN_MASS_CUT=0.0;
  app->SetDefaultParameter("ANTIPROTON:MIN_MASS_CUT",MIN_MASS_CUT); 

  PI0_CUT_VALUE=0.025;
  ETA_CUT_VALUE=0.05;
  ETAPRIME_CUT_VALUE=0.05;
  NUM_SIGMA_BG=2; 
  app->SetDefaultParameter("ANTIPROTON:NUM_SIGMA_BG",NUM_SIGMA_BG); 

  SPLIT_CUT=0.;
  app->SetDefaultParameter("ANTIPROTON:SPLIT_CUT",SPLIT_CUT); 
  FCAL_POS_CUT=10.0;
  app->SetDefaultParameter("ANTIPROTON:FCAL_POS_CUT",FCAL_POS_CUT); 
  FCAL_RADIAL_CUT=104.0;
  app->SetDefaultParameter("ANTIPROTON:FCAL_RADIAL_CUT",FCAL_RADIAL_CUT);
  BCAL_Z_CUT=384.0;
  app->SetDefaultParameter("ANTIPROTON:BCAL_Z_CUT",BCAL_Z_CUT); 
  FCAL_THRESHOLD=0.1;
  app->SetDefaultParameter("ANTIPROTON:FCAL_THRESHOLD",FCAL_THRESHOLD);
  BCAL_THRESHOLD=0.05;
  app->SetDefaultParameter("ANTIPROTON:BCAL_THRESHOLD",BCAL_THRESHOLD);

  lockService = app->GetService<JLockService>();
  lockService->RootWriteLock();
  
  gDirectory->mkdir("antiproton")->cd();
  
  RFGammaTimeDiff=new TH1F("RFGammaTimeDiff","#deltat between RF and #gamma at vertex",200,-10,10);

  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("AntiP2P4G")->cd();
  {
    FourGamma2d_kf=new TH2F("FourGamma2d_kf",";M(2#gamma,pair1)[GeV];M(2#gamma,pair2)[GeV]",500,0,1,500,0,1); 
    AntiPPi0_vs_PPi0=new TH2F("AntiPPi0_vs_PPi0",
			      ";M(p#pi^{0}) [GeV];M(#bar{p}#pi^{0}) [GeV]",
			      400,0.9,2.9,400,0.9,2.9);
    SigpAntiSigp= new TH1F("SigpAntiSigp",";M(#Sigma^{+}#bar{#Sigma}^{-}) [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("../");

  lockService->RootUnLock();
  
  // Make a container to keep track of doublets of photons
  vector<pair<unsigned int,unsigned int>>pairs;
  pairs.push_back(make_pair(0,1));
  pairs.push_back(make_pair(2,3));
  doublets.push_back(pairs); 
  pairs.clear();
  pairs.push_back(make_pair(0,2));
  pairs.push_back(make_pair(1,3));
  doublets.push_back(pairs); 
  pairs.clear(); 
  pairs.push_back(make_pair(0,3));
  pairs.push_back(make_pair(2,1));
  doublets.push_back(pairs);

  // set up maps for keeping track of pairs of mesons
  vector<pair<DLorentzVector,DLorentzVector>>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);
  particle_pairs.emplace(EtaX_,my_pair_vector);
  particle_pairs.emplace(Pi0X_,my_pair_vector);

  bachelor_gammas.emplace(EtaX_,my_pair_vector);
  bachelor_gammas.emplace(Pi0X_,my_pair_vector);

  lockService->CreateLock("myantilock");
}

//------------------
// BeginRun
//------------------
void JEventProcessor_antiproton::BeginRun(const std::shared_ptr<const JEvent>& event)
{
  const DGeometry *geom=DEvent::GetDGeometry(event);
  double x0,y0,z0;
  geom->GetFCALPosition(x0,y0,z0);
  mFCALCenter.SetXYZ(x0,y0,z0);
}

//------------------
// Process
//------------------
void JEventProcessor_antiproton::Process(const std::shared_ptr<const JEvent>& event)
{
  vector<const DChargedTrack*>tracks; 
  event->Get(tracks);
  if (tracks.size()==0) return;

  // Assign charge track by PID
  map<Particle_t,vector<const DTrackTimeBased *>>chargedParticles;
  FillChargedParticleVectors(tracks,chargedParticles);
  if (chargedParticles[AntiProton].size()==0) return;
  
  vector<const DBeamPhoton*>beamphotons;
  event->Get(beamphotons);
  if (beamphotons.size()==0) return;

  vector<const DNeutralParticle*>neutrals;
  event->Get(neutrals);

  lockService->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.;
  vector<DNeutralParticleHypothesis>gammaParticles;
  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<pair<const DBeamPhoton*,double>>beamPhotons;
  for (unsigned int i=0;i<beamphotons.size();i++){
    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(event);
  DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(event);
  DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils);   
  dKinFitUtils->Reset_NewEvent();
  dKinFitter->Reset_NewEvent();
    
  for (unsigned int i=0;i<beamPhotons.size();i++){
    double weight=beamPhotons[i].second;

    //--------------------------------
    // Kinematic fit w/ nothing missing
    //--------------------------------
    if (DoKinematicFit(t0_rf,vertex,beamPhotons[i].first,chargedParticles,
		       gammaParticles,dKinFitUtils,dKinFitter)){
      double CL = dKinFitter->Get_ConfidenceLevel();
      if (CL>CL_CUT){
	DLorentzVector beam_kf;
	map<Particle_t,vector<DLorentzVector>>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;
	    case 4:
	      AntiProton2Proton4GammaAnalysis(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;
  
  lockService->Unlock("myantilock");
}

//------------------
// EndRun
//------------------
void JEventProcessor_antiproton::EndRun(void)
{
}

//------------------
// Finish
//------------------
void JEventProcessor_antiproton::Finish(void)
{
}

// Use the PIDFOM scheme to sort tracks according to particle type
void
JEventProcessor_antiproton::FillChargedParticleVectors(vector<const DChargedTrack *>&tracks,map<Particle_t,vector<const DTrackTimeBased *>>&chargedParticles) const {
  Particle_t particle_list[4]={PiPlus,PiMinus,Proton,AntiProton};
  for (unsigned int j=0;j<tracks.size();j++){  
    const DChargedTrackHypothesis *hyp=NULL;
    vector<pair<double,Particle_t> >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<const DNeutralParticle*>&neutrals,
						    vector<DNeutralParticleHypothesis>&gammaParticles
						    ) const{
  vector<const DNeutralParticleHypothesis*>tempFCALGammas;
  for (unsigned int i=0;i<neutrals.size();i++){
    const DNeutralParticleHypothesis *gamHyp=neutrals[i]->Get_Hypothesis(Gamma);
    double tdiff=gamHyp->time()-t0_rf;
    DLorentzVector gamma_v4=gamHyp->lorentzMomentum(); 

    RFGammaTimeDiff->Fill(tdiff);

    // Get shower info correspnding to this hypothesis
    const DNeutralShower *shower=gamHyp->Get_NeutralShower();

    // Look for good photons
    bool got_photon=(fabs(tdiff)<PROTON_GAMMA_DT_CUT);
    if (got_photon){      
      if (shower->dDetectorSystem==SYS_FCAL){
	if (shower->dQuality<SPLIT_CUT) got_photon=false;
	DVector3 pos=(dynamic_cast<const DFCALShower*>(shower->dBCALFCALShower))->getPosition();
	DVector3 pos_rel=pos-mFCALCenter;
	
	if (pos_rel.Perp()>FCAL_RADIAL_CUT) got_photon=false;
	if (fabs(pos_rel.X())<FCAL_POS_CUT && fabs(pos_rel.Y())<FCAL_POS_CUT){
	  got_photon=false;
	}
	if (got_photon){
	  tempFCALGammas.push_back(gamHyp);
	}
      }
      else if (shower->dDetectorSystem==SYS_BCAL){
	const DBCALShower *bcalshower=dynamic_cast<const DBCALShower*>(shower->dBCALFCALShower);
	if (gamma_v4.E()<BCAL_THRESHOLD){
	  got_photon=false;
	}
	if (bcalshower->z>BCAL_Z_CUT) got_photon=false;
	if (got_photon){
	  //got_bcal=true; 
	  gammaParticles.push_back(*gamHyp);
	}
      }
    }
  }
  double minMass=MIN_MASS_CUT;
  unsigned int numfcals=tempFCALGammas.size();
  unsigned int min_i=numfcals,min_j=numfcals;
  for (unsigned int i=0;i<numfcals;i++){
    for (unsigned int j=i+1;j<numfcals;j++){
      DLorentzVector twogam=tempFCALGammas[i]->lorentzMomentum()
	+tempFCALGammas[j]->lorentzMomentum();
      //      TwoGammaMassFCAL->Fill(twogam.M());
      if (twogam.M()<minMass){
	minMass=twogam.M();
	min_i=i;
	min_j=j;
      }
    }
  }
  for (unsigned int i=0;i<numfcals;i++){
    if (i!=min_i && i!=min_j && tempFCALGammas[i]->lorentzMomentum().E()>FCAL_THRESHOLD){
      gammaParticles.push_back(*tempFCALGammas[i]);
    }
  }
  if (min_i<numfcals && min_j<numfcals){
    DNeutralParticleHypothesis extraGamma;
    DVector3 momentum=tempFCALGammas[min_i]->momentum()
	    +tempFCALGammas[min_j]->momentum();
    if (momentum.Mag()>FCAL_THRESHOLD){
      TMatrixFSym tempMatrix(7);
      for (unsigned i=0;i<7;i++){
	for (unsigned int j=0;j<7;j++){
	  tempMatrix(i,j)=(*tempFCALGammas[min_i]->errorMatrix())(i,j)
	    +(*tempFCALGammas[min_j]->errorMatrix())(i,j);
	}
      }
      extraGamma.setMomentum(momentum);
      extraGamma.setErrorMatrix(const_pointer_cast<const TMatrixFSym>(make_shared<TMatrixFSym>(tempMatrix)));
      gammaParticles.push_back(extraGamma);
    }
  }
}


// Run the kinematic fitter requiring energy and momentum conservation
bool JEventProcessor_antiproton::DoKinematicFit(double t0_rf,const DVector3 &vertex,

						 const DBeamPhoton *beamphoton,
						 map<Particle_t,vector<const DTrackTimeBased *> >&chargedParticles,
						 vector<DNeutralParticleHypothesis>&gammaParticles,
						 DKinFitUtils_GlueX *dKinFitUtils,
						 DKinFitter *dKinFitter) const {
  set<shared_ptr<DKinFitParticle>> InitialParticles, FinalParticles;
  
  dKinFitter->Reset_NewFit();
  shared_ptr<DKinFitParticle>myBeam=dKinFitUtils->Make_BeamParticle(beamphoton);
  shared_ptr<DKinFitParticle>myTarget=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++){ 
    vector<const DTrackTimeBased *>myParticles=chargedParticles[particle_list[k]];
    for (unsigned int m=0;m<myParticles.size();m++){
      shared_ptr<DKinFitParticle>myParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]);
      FinalParticles.insert(myParticle);
    } 
  }

  // Final state photons
  for (unsigned int k=0;k<gammaParticles.size();k++){
    DNeutralParticleHypothesis gammaHyp=gammaParticles[k];
    shared_ptr<DKinFitParticle>myGamma=dKinFitUtils->Make_DetectedParticle(22,0,0.,vert4,gammaHyp.momentum(),0.,gammaHyp.errorMatrix());
    FinalParticles.insert(myGamma);
  }

  // make constraint
  shared_ptr<DKinFitConstraint_P4>locP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles);
  // set constraint
  dKinFitter->Add_Constraint(locP4Constraint);
  
  // PERFORM THE KINEMATIC FIT
  return dKinFitter->Fit_Reaction();

}

void JEventProcessor_antiproton::GetKF4vectors(DKinFitter *dKinFitter,
					       DLorentzVector &beam_kf,
					       map<Particle_t,vector<DLorentzVector>>&final_kf) const{
  
  set<shared_ptr<DKinFitParticle>>myParticles=dKinFitter->Get_KinFitParticles();
  set<shared_ptr<DKinFitParticle>>::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());
    }
  }
}

bool JEventProcessor_antiproton::MakeMesonPairs(vector<DLorentzVector>&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();
  particle_pairs[Pi0X_].clear();
  particle_pairs[EtaX_].clear();
  bachelor_gammas[Pi0X_].clear();
  bachelor_gammas[EtaX_].clear();

  // Find pairs of mesons
  for (unsigned int i=0;i<doublets.size();i++){
    DLorentzVector pair1=gammas_kf[doublets[i][0].first]
      +gammas_kf[doublets[i][0].second];
    DLorentzVector pair2=gammas_kf[doublets[i][1].first]
      +gammas_kf[doublets[i][1].second];
    double pair1_mass=pair1.M();
    double pair2_mass=pair2.M();
    histo->Fill(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));
    }
    else if (pi0_1_cut
	     && PI0_ANTI_CUT(pair2_mass,PI0_VETO_CUT)
	     ){
      particle_pairs[Pi0X_].push_back(make_pair(pair1,pair2));
      bachelor_gammas[Pi0X_].push_back(make_pair(gammas_kf[doublets[i][1].first],
						 gammas_kf[doublets[i][1].second]));
    }
    else if (pi0_2_cut
	     && PI0_ANTI_CUT(pair1_mass,PI0_VETO_CUT)
	     ){
      particle_pairs[Pi0X_].push_back(make_pair(pair2,pair1));
      bachelor_gammas[Pi0X_].push_back(make_pair(gammas_kf[doublets[i][0].first],
						 gammas_kf[doublets[i][0].second]));
    }
    else if (eta_1_cut){
      particle_pairs[EtaX_].push_back(make_pair(pair1,pair2));
      bachelor_gammas[EtaX_].push_back(make_pair(gammas_kf[doublets[i][1].first],
						 gammas_kf[doublets[i][1].second]));
    }
    else if (eta_2_cut){
      particle_pairs[EtaX_].push_back(make_pair(pair2,pair1));
      bachelor_gammas[EtaX_].push_back(make_pair(gammas_kf[doublets[i][0].first],
						 gammas_kf[doublets[i][0].second]));
    }
  }

  return (particle_pairs[Pi0Pi0_].size()>0 
	  || particle_pairs[Pi0Eta_].size()>0
	  || particle_pairs[EtaEta_].size()>0
	  || particle_pairs[Pi0EtaPrime_].size()>0
	  );
}





//------------------------------------------------------------------------------
// AntiProton 2Proton Ngamma
//------------------------------------------------------------------------------
void JEventProcessor_antiproton::AntiProton2ProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&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<Particle_t,vector<DLorentzVector>>&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<Particle_t,vector<DLorentzVector>>&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);
  }

}

void JEventProcessor_antiproton::AntiProton2Proton4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&final_kf,
						      double weight){
  bool got_meson_pairs=MakeMesonPairs(final_kf[Gamma],weight,FourGamma2d_kf);
  DLorentzVector antiproton=final_kf[AntiProton][0];
  DLorentzVector proton1=final_kf[Proton][0];
  DLorentzVector proton2=final_kf[Proton][1];
  if (got_meson_pairs){
    if (particle_pairs[Pi0Pi0_].size()==1){
      DLorentzVector pi0_1=particle_pairs[Pi0Pi0_][0].first; 
      DLorentzVector pi0_2=particle_pairs[Pi0Pi0_][0].second;

      DLorentzVector ap_pi0_1=antiproton+pi0_1;
      double m_ap_pi1=ap_pi0_1.M();
      DLorentzVector p1_pi0_2=proton1+pi0_2;
      double m_p1_pi2=p1_pi0_2.M();
      DLorentzVector p2_pi0_2=proton2+pi0_2;
      double m_p2_pi2=p2_pi0_2.M();
      DLorentzVector ap_p1_4pi=ap_pi0_1+p1_pi0_2;
      double m_ap_p1_4pi=ap_p1_4pi.M();
      DLorentzVector ap_p2_4pi=ap_pi0_1+p2_pi0_2;
      double m_ap_p2_4pi=ap_p2_4pi.M();
      if (m_ap_pi1>1.16 && m_ap_pi1<1.21){
	if (m_p1_pi2>1.16 && m_p1_pi2<1.21){
	  SigpAntiSigp->Fill(m_ap_p1_4pi,weight);
	}
	if (m_p2_pi2>1.16 && m_p2_pi2<1.21){
	  SigpAntiSigp->Fill(m_ap_p2_4pi,weight);
	}
      }
      AntiPPi0_vs_PPi0->Fill(m_p1_pi2,m_ap_pi1,weight); 
      AntiPPi0_vs_PPi0->Fill(m_p2_pi2,m_ap_pi1,weight);
      
      DLorentzVector ap_pi0_2=antiproton+pi0_2;
      double m_ap_pi2=ap_pi0_2.M();
      DLorentzVector p1_pi0_1=proton1+pi0_1;
      double m_p1_pi1=p1_pi0_1.M();
      DLorentzVector p2_pi0_1=proton2+pi0_1;
      double m_p2_pi1=p2_pi0_1.M();
      AntiPPi0_vs_PPi0->Fill(m_p1_pi1,m_ap_pi2,weight); 
      AntiPPi0_vs_PPi0->Fill(m_p2_pi1,m_ap_pi2,weight);
      if (m_ap_pi2>1.16 && m_ap_pi2<1.21){
	if (m_p1_pi1>1.16 && m_p1_pi1<1.21){
	  SigpAntiSigp->Fill(m_ap_p1_4pi,weight);
	}
	if (m_p2_pi1>1.16 && m_p2_pi1<1.21){
	  SigpAntiSigp->Fill(m_ap_p2_4pi,weight);
	}
      }
    }
  }

}



//------------------------------------------------------------------------------
// AntiProton 2Proton pi+ pi- Ngamma
//------------------------------------------------------------------------------
void JEventProcessor_antiproton::PipPimAntiProton2ProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&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<Particle_t,vector<DLorentzVector>>&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<Particle_t,vector<DLorentzVector>>&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<Particle_t,vector<DLorentzVector>>&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);
    }
    
  }
}