// $Id$
//
//    File: JEventProcessor_InsertStudies.cc
// Created: Mon Dec  9 14:08:13 EST 2019
// Creator: staylor (on Linux gluon117.jlab.org 3.10.0-957.21.3.el7.x86_64 x86_64)
//

#include "JEventProcessor_InsertStudies.h"
using namespace jana;

#define PI0_CUT(m) (fabs(m-ParticleMass(Pi0))<0.026) // "3.5 sigma" cut
#define TIGHT_PI0_CUT(m) (fabs(m-ParticleMass(Pi0))<0.015) // ~2 sigma cut
#define PI0_VETO_CUT(m) (fabs(m-ParticleMass(Pi0))>0.026) // 3 sigma anticut
#define ETA_CUT(m) (fabs(m-ParticleMass(Eta))<0.04) // "3 sigma" cut?

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


//------------------
// JEventProcessor_InsertStudies (Constructor)
//------------------
JEventProcessor_InsertStudies::JEventProcessor_InsertStudies()
{

}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_InsertStudies::init(void)
{
  MAX_BEAM_E=12.;
  gPARMS->SetDefaultParameter("InsertStudies:MAX_BEAM_E",MAX_BEAM_E);   
  MIN_BEAM_E=3.; 
  gPARMS->SetDefaultParameter("InsertStudies:MIN_BEAM_E",MIN_BEAM_E);
  VETO_BCAL=false; 
  gPARMS->SetDefaultParameter("InsertStudies:VETO_BCAL",VETO_BCAL);
  FCAL_THRESHOLD=0.1;
  gPARMS->SetDefaultParameter("InsertStudies:FCAL_THRESHOLD",FCAL_THRESHOLD);
  BCAL_THRESHOLD=0.05; 
  gPARMS->SetDefaultParameter("InsertStudies:BCAL_THRESHOLD",BCAL_THRESHOLD);
  PROTON_GAMMA_DT_CUT=1.;
  CHISQ_CUT=3.;
  
  gDirectory->mkdir("InsertStudies")->cd();
  gDirectory->mkdir("MC")->cd();
  MCBeam=new TH1F("MCBeam","Thrown beam energy",440,2.8,11.6);
  MCMissingMass=new TH1F("MCMissingMass","Missing mass off proton",3000,0.,3.);
  gDirectory->cd("../");
  gDirectory->mkdir("timing")->cd();
  ProtonGammaTimeDiff=new TH1F("ProtonGammaTimeDiff",";#Deltat [ns]",200,-20,20); 
  BeamProtonTimeDiff=new TH1F("BeamProtonTimeDiff",";#Deltat [ns]",200,-20,20);
  gDirectory->cd("../");
  gDirectory->mkdir("NoProton")->cd();
  BeamTimeDiff=new TH1F("BeamTimeDiff",";#Deltat [ns]",200,-20,20);
  gDirectory->mkdir("TwoGamma")->cd();
  TwoGammaChiSqMissingProton=new TH1F("TwoGammaChiSqMissingProton",";#chi^{2}/ndf",1000,0,100);
  Measured2GammaMissingProton=new TH1F("Measured2GammaMissingProton",";M(2#gamma) [GeV]",
			      3000,0.,3.); 
  TwoGammaMissingProton=new TH1F("TwoGammaMissingProton",";M(2#gamma) [GeV]",
			     3000,0.,3.);
  MissingMass2g=new TH1F("MissingMass2g",";Missing mass [GeV]",
			     3000,0.,3.);
  
  gDirectory->cd("../");

  gDirectory->mkdir("FourGamma")->cd();
  FourGammaChiSqMissingProton=new TH1F("FourGammaChiSqMissingProton",";#chi^{2}/ndf",1000,0,100);
  Measured4GammaMissingProton=new TH1F("Measured4GammaMissingProton",";M(4#gamma) [GeV]",
			      3000,0.,3.); 
  FourGammaMissingProton=new TH1F("FourGammaMissingProton",";M(4#gamma) [GeV]",
			     3000,0.,3.);
  MissingMass=new TH1F("MissingMass",";Missing mass [GeV]",
			     3000,0.,3.);

  FourGamma2DMissingProton=new TH2F("FourGamma2DMissingProton","M(2#gamma pair 1) vs M(2#gamma pair 2)",
				    1000,0.,1.,1000,0,1);
  
  gDirectory->cd("../");
  gDirectory->cd("../");
  gDirectory->mkdir("FourGamma")->cd();
  FourGammaChiSq=new TH2F("FourGammaChiSq",";#chi^{2}/ndf",1000,0,100,2,-0.5,1.5);
  Measured4GammaMass=new TH2F("Measured4GammaMass",";M(4#gamma) [GeV]",
			      3000,0.,3.,5,-0.5,4.5); 
  FourGammaMass=new TH2F("FourGammaMass",";M(4#gamma) [GeV]",
			 3000,0.,3.,5,-0.5,4.5); 
  TwoPi0NoOverlap=new TH1F("TwoPi0NoOverlap",";M(2#pi^{0}) [GeV]",
			   3000,0.,3.);
  TwoEtaNoOverlap=new TH1F("TwoEtaNoOverlap",";M(2#eta) [GeV]",
			   3000,0.,3.); 
  EtaPi0NoOverlap=new TH1F("EtaPi0NoOverlap",";M(#eta#pi^{0}) [GeV]",
			   3000,0.,3.);
  
  Pi02gNoVeto=new TH2F("Pi02gNoVeto","With no veto on 2nd #pi^{0};M(#pi^{0}2#gamma) [GeV]",
		       3000,0.,3.,5,-0.5,4.5); 
  Pi02gWithVeto=new TH2F("Pi02gWithVeto","With veto on 2nd #pi^{0};M(#pi^{0}2#gamma) [GeV]",
			 3000,0.,3.,5,-0.5,4.5);  
  TwoGammaOffPi0=new TH2F("TwoGammaOffPi0","2#gamma mass recoiling off #pi^{0}};M(2#gamma) [GeV]",
			  1000,0.,1.,5,-0.5,4.5);
  
  FourGamma2D=new TH2F("FourGamma2D","M(2#gamma pair 1) vs M(2#gamma pair 2)",
		       1000,0.,1.,1000,0,1);
  FourGamma2DNoOverlap=new TH2F("FourGamma2DNoOverlap","M(2#gamma pair 1) vs M(2#gamma pair 2)",
		       1000,0.,1.,1000,0,1);
  Pi02gVsPi0g=new TH2F("Pi02gVsPi0g","M(#pi^{0}2#gamma) vs M(#pi^{0}#gamma)",
		       1000,0.,1.,2000,0,2);
  Pi02gVsPi0gWithVeto=new TH2F("Pi02gVsPi0gWithVeto",
			       "M(#pi^{0}2#gamma) vs M(#pi^{0}#gamma)",
			       1000,0.,1.,2000,0,2);
  
  Pi02gNoVetoNoOverlap=new TH2F("Pi02gNoVetoNoOverlap","With no veto on 2nd #pi^{0};M(#pi^{0}2#gamma) [GeV]",
		       3000,0.,3.,5,-0.5,4.5); 
  Pi02gWithVetoNoOverlap=new TH2F("Pi02gWithVetoNoOverlap","With veto on 2nd #pi^{0};M(#pi^{0}2#gamma) [GeV]",
			 3000,0.,3.,5,-0.5,4.5);  
  TwoGammaOffPi0NoOverlap=new TH2F("TwoGammaOffPi0NoOverlap","2#gamma mass recoiling off #pi^{0}};M(2#gamma) [GeV]",
			  1000,0.,1.,5,-0.5,4.5);
  Pi02gVsPi0gNoOverlap=new TH2F("Pi02gVsPi0gNoOverlap","M(#pi^{0}2#gamma) vs M(#pi^{0}#gamma)",
		       1000,0.,1.,2000,0,2);
  Pi02gVsPi0gWithVetoNoOverlap=new TH2F("Pi02gVsPi0gWithVetoNoOverlap",
			       "M(#pi^{0}2#gamma) vs M(#pi^{0}#gamma)",
			       1000,0.,1.,2000,0,2);
  

  
  gDirectory->cd("../"); 
  
  gDirectory->mkdir("SixGamma")->cd();
  SixGammaChiSq=new TH2F("SixGammaChiSq",";#chi^{2}/ndf",1000,0,100,2,-0.5,1.5);
  Measured6GammaMass=new TH1F("Measured6GammaMass",";M(6#gamma) [GeV]",
			      3000,0.,3.);  
  SixGammaMass=new TH1F("SixGammaMass",";M(6#gamma) [GeV]",3000,0.,3.); 
  ThreePi0Mass=new TH1F("ThreePi0Mass",";M(3#pi^{0}) [GeV]",3000,0.,3.);
  Eta2Pi0Mass=new TH1F("Eta2Pi0Mass",";M(#eta2#pi^{0}) [GeV]",3000,0.,3.);
  ThreePi0MassNoOverlap=new TH1F("ThreePi0MassNoOverlap",";M(3#pi^{0}) [GeV]",3000,0.,3.);
  Eta2Pi0MassNoOverlap=new TH1F("Eta2Pi0MassNoOverlap",";M(#eta2#pi^{0}) [GeV]",3000,0.,3.);
  FourGamma2D_6g=new TH2F("FourGamma2D_6g",";M(2#gamma pair 1) vs  M(2#gamma pair 2)",1000,0.,1.,1000,0,1);
  Pi0gVsPi0g=new TH2F("Pi0gVsPi0g","M#(#pi^{0}#gamma, pair 1) vs M(#pi^{0}#gamma, pair 2)",1000,0,2,1000,0,2);
  Pi0gVsPi0gNoOverlap=new TH2F("Pi0gVsPi0gNoOverlap","M#(#pi^{0}#gamma, pair 1) vs M(#pi^{0}#gamma, pair 2)",1000,0,2,1000,0,2);
  FourGamma2D_6g_NoOverlap=new TH2F("FourGamma2D_6g_NoOverlap",";M(2#gamma pair 1) vs  M(2#gamma pair 2)",1000,0.,1.,1000,0,1);

  TwoEtaPi0Mass=new TH1F("TwoEtaPi0Mass",";M(2#eta#pi^{0}) [GeV]",3000,0.,3.);
  TwoEtaPi0MassNoOverlap=new TH1F("TwoEtaPi0MassNoOverlap",";M(2#eta#pi^{0}) [GeV]",3000,0.,3.);
  TwoEtaPi0Dalitz=new TH2F("TwoEtaPi0Dalitz","M^2(#eta_2#pi0) vs M^2(#eta_1#pi0)",100,0,4,100,0,4);

  

  gDirectory->cd("../");


  gDirectory->cd("../");

  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_InsertStudies::brun(JEventLoop *eventLoop, int32_t runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_InsertStudies::evnt(JEventLoop *loop, uint64_t eventnumber)
{ vector<const DMCThrown*>throwns;
  loop->Get(throwns); 

  vector<const DChargedTrack*>tracks; 
  loop->Get(tracks);

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

  vector<const DBeamPhoton*>beamphotons;
  loop->Get(beamphotons);

  vector<const DTAGMGeometry *>tagm_geom;
  loop->Get(tagm_geom); 
  vector<const DTAGHGeometry *>tagh_geom;
  loop->Get(tagh_geom);

  // Check for thrown events
  const DMCReaction* reaction=NULL;
  if (throwns.size()>0){
    loop->GetSingle(reaction);
  }
  
  japp->WriteLock("myInsertStudiesLock");

  // Handle thrown events
  if (reaction!=NULL&& tagm_geom.size()>0 && tagh_geom.size()>0){
    FillThrownHistos(reaction,throwns,tagm_geom,tagh_geom);
  }

  // Assign charge track by PID
  map<Particle_t,vector<const DTrackTimeBased *> >particles;
  FillParticleVectors(tracks,particles);

  // Find t0 (at "vertex") for event
  double t0_rf=-1000.;
  if (tracks.size()>0){
    t0_rf=tracks[0]->Get_BestTrackingFOM()->t0();
  }
  else if (neutrals.size()>0){
    t0_rf=neutrals[0]->Get_Hypothesis(Gamma)->t0();
  }

  // Get gamma candidates:  check that they are in time with the charged 
  // particles at the vertex and count how many in-time photons are in the 
  // insert.
  vector<const DNeutralParticleHypothesis*>gammaHyps;
  unsigned int insert_count=FillGammaVector(t0_rf,neutrals,gammaHyps);

  if (particles.size()==0 && gammaHyps.size()>0){ // no protons, at least one photon candidate
    DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop);
    DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils);   
    dKinFitUtils->Reset_NewEvent();
    dKinFitter->Reset_NewEvent();
    
    for (unsigned int i=0;i<beamphotons.size();i++){
      if (beamphotons[i]->lorentzMomentum().E()<MIN_BEAM_E
	  || beamphotons[i]->lorentzMomentum().E()>MAX_BEAM_E) continue;
       double dt_st=beamphotons[i]->time()-t0_rf;

      BeamTimeDiff->Fill(dt_st);
	
      double weight=1.;
      bool got_beam_photon=false;
      if (fabs(dt_st)>6.012&&fabs(dt_st)<18.036){
	weight=-1./6.;
	got_beam_photon=true;
      }
      if (fabs(dt_st)<2.004){
	got_beam_photon=true;
      }
      if (got_beam_photon==false) continue;
      
      //--------------------------------
      // Kinematic fit
      //--------------------------------
      dKinFitter->Reset_NewFit();
      DoKinematicFitMissingProton(beamphotons[i],particles,gammaHyps,dKinFitUtils,dKinFitter);
      double ChiSq=dKinFitter->Get_ChiSq();
      unsigned int NDoF=dKinFitter->Get_NDF();
      double ReducedChiSq=ChiSq/double(NDoF);
      switch(gammaHyps.size()){
      case 2:
	{
	  TwoGammaChiSqMissingProton->Fill(ReducedChiSq);
	  DLorentzVector twogam;
	  for (unsigned int j=0;j<2;j++){
	    twogam+=gammaHyps[j]->lorentzMomentum();
	  }
	  DLorentzVector target_v4(0,0,0,ParticleMass(Proton));
	  DLorentzVector beam_v4=beamphotons[i]->lorentzMomentum();
	  DLorentzVector missing_v4=beam_v4+target_v4-twogam;
	  Measured2GammaMissingProton->Fill(twogam.M(),weight);
	  MissingMass2g->Fill(missing_v4.M(),weight);
	  if (ReducedChiSq<CHISQ_CUT){	   
	    TwoGammaAnalysis(insert_count,dKinFitter,weight,true);
	  }
	}
	break;
      case 4:
	{
	  FourGammaChiSqMissingProton->Fill(ReducedChiSq);
	  DLorentzVector fourgam;
	  for (unsigned int j=0;j<4;j++){
	    fourgam+=gammaHyps[j]->lorentzMomentum();
	  }
	  DLorentzVector target_v4(0,0,0,ParticleMass(Proton));
	  DLorentzVector beam_v4=beamphotons[i]->lorentzMomentum();
	  DLorentzVector missing_v4=beam_v4+target_v4-fourgam;
	  Measured4GammaMissingProton->Fill(fourgam.M(),weight);
	  MissingMass->Fill(missing_v4.M(),weight);
	  if (ReducedChiSq<CHISQ_CUT){	   
	    FourGammaAnalysis(insert_count,dKinFitter,weight,true);
	  }
	}
	break;
      default:
	break;
      }
    } 

    delete dKinFitter;
    delete dKinFitUtils;
  }


  if (particles[Proton].size()==1){
    DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop);
    DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils);   
    dKinFitUtils->Reset_NewEvent();
    dKinFitter->Reset_NewEvent();

    for (unsigned int i=0;i<beamphotons.size();i++){
      if (beamphotons[i]->lorentzMomentum().E()<MIN_BEAM_E
	  || beamphotons[i]->lorentzMomentum().E()>MAX_BEAM_E) continue;
      
      double dt_st=beamphotons[i]->time()-t0_rf;

      BeamProtonTimeDiff->Fill(dt_st);
	
      double weight=1.;
      bool got_beam_photon=false;
      if (fabs(dt_st)>6.012&&fabs(dt_st)<18.036){
	weight=-1./6.;
	got_beam_photon=true;
      }
      if (fabs(dt_st)<2.004){
	got_beam_photon=true;
      }
      if (got_beam_photon==false) continue;
      
      //--------------------------------
      // Kinematic fit
      //--------------------------------
      dKinFitter->Reset_NewFit();
      DoKinematicFit(beamphotons[i],particles,gammaHyps,dKinFitUtils,dKinFitter);
      double ChiSq=dKinFitter->Get_ChiSq();
      unsigned int NDoF=dKinFitter->Get_NDF();
      double ReducedChiSq=ChiSq/double(NDoF);
      if (particles[PiPlus].size()==0 && particles[PiMinus].size()==0
	  && particles[Electron].size()==0 && particles[Positron].size()==0){
	switch(gammaHyps.size()){
	case 4:
	  FourGammaChiSq->Fill(ReducedChiSq,(weight>0)?0.:1.);
	  if (ReducedChiSq<CHISQ_CUT){
	    DLorentzVector fourgam;
	    for (unsigned int j=0;j<4;j++){
	      fourgam+=gammaHyps[j]->lorentzMomentum();
	    }
	    Measured4GammaMass->Fill(fourgam.M(),insert_count,weight);
	    FourGammaAnalysis(insert_count,dKinFitter,weight);

	  }
	  break;
	case 6: 
	  SixGammaChiSq->Fill(ReducedChiSq,(weight>0)?0.:1.);
	  if (ReducedChiSq<CHISQ_CUT){
	    DLorentzVector sixgam;
	    for (unsigned int j=0;j<6;j++){
	      sixgam+=gammaHyps[j]->lorentzMomentum();
	    }
	    Measured6GammaMass->Fill(sixgam.M(),weight);
	    SixGammaAnalysis(dKinFitter,weight);

	  }
	  break;
	default:
	  break;
	}
      } 
    }

    delete dKinFitter;
    delete dKinFitUtils;
  }

  japp->Unlock("myInsertStudiesLock");
  
  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_InsertStudies::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_InsertStudies::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}



jerror_t 
JEventProcessor_InsertStudies::FillParticleVectors(vector<const DChargedTrack *>&tracks,
						   map<Particle_t, vector<const DTrackTimeBased *> >&particles
						   ) const {
  for (unsigned int j=0;j<tracks.size();j++){  
    const DChargedTrackHypothesis *hyp=NULL;
    vector<pair<double,Particle_t> >probabilities;
    Particle_t particle_list[8]={Proton,AntiProton,KPlus,KMinus,PiPlus,PiMinus,
				 Positron,Electron};
    for (unsigned int i=0;i<8;i++){
      if ((hyp=tracks[j]->Get_Hypothesis(particle_list[i]))!=NULL){
	probabilities.push_back(make_pair(hyp->Get_FOM(),particle_list[i]));
      } 
    }
    if (probabilities.size()>0){
      sort(probabilities.begin(),probabilities.end(),SortParticleProbability);
      Particle_t myParticle=probabilities[0].second;
      hyp=tracks[j]->Get_Hypothesis(myParticle);
      particles[myParticle].push_back(hyp->Get_TrackTimeBased());
    }
  }
    
  return NOERROR;
}

unsigned int 
JEventProcessor_InsertStudies::FillGammaVector(double t0_rf,
					       vector<const DNeutralParticle*>&neutrals,
					       vector<const DNeutralParticleHypothesis*>&gammaHyps) const{
  unsigned int insert_count=0;
  bool vetoed_bcal=false;
  for (unsigned int i=0;i<neutrals.size();i++){
    const DNeutralParticleHypothesis *gamHyp=neutrals[i]->Get_Hypothesis(Gamma);
    double tdiff=gamHyp->time()-t0_rf;

    ProtonGammaTimeDiff->Fill(tdiff);

    if (fabs(tdiff)<PROTON_GAMMA_DT_CUT){  
      DLorentzVector gamma_v4=gamHyp->lorentzMomentum();
      const DNeutralShower *shower=NULL;
      shower=gamHyp->Get_NeutralShower();
      if (shower->dDetectorSystem==SYS_FCAL){
	if (gamma_v4.E()<FCAL_THRESHOLD){
	  continue;
	}
	DVector3 pos=(dynamic_cast<const DFCALShower*>(shower->dBCALFCALShower))->getPosition();
	if (fabs(pos.X())<50.16 && fabs(pos.Y())<50.16){
	  insert_count++;
	}	 
      }
      else if (shower->dDetectorSystem==SYS_BCAL){
	if (gamma_v4.E()<BCAL_THRESHOLD){
	  continue;
	}
	else if (VETO_BCAL){
	  vetoed_bcal=true;
	  break;
	}
      }
      gammaHyps.push_back(gamHyp);
    }
  }
  if (vetoed_bcal) gammaHyps.clear();
  return insert_count;
}

// Run the kinematic fitter requiring energy and momentum conservation
void JEventProcessor_InsertStudies::DoKinematicFit(const DBeamPhoton *beamphoton,
						   map<Particle_t,vector<const DTrackTimeBased *> >&particles,
						   vector<const DNeutralParticleHypothesis*>&gammaHyps,
						   DKinFitUtils_GlueX *dKinFitUtils,
						   DKinFitter *dKinFitter) const {
  set<shared_ptr<DKinFitParticle>> InitialParticles, FinalParticles;
  
  shared_ptr<DKinFitParticle>myBeam=dKinFitUtils->Make_BeamParticle(beamphoton);
  shared_ptr<DKinFitParticle>myTarget=dKinFitUtils->Make_TargetParticle(Proton);
  
  InitialParticles.insert(myBeam);  
  InitialParticles.insert(myTarget);
  
  // Charged particles
  Particle_t particle_list[5]={Proton,PiPlus,PiMinus,Positron,Electron};
  TLorentzVector vert;
  bool got_vertex_track=false;
  for (unsigned int k=0;k<5;k++){ 
    vector<const DTrackTimeBased *>myParticles=particles[particle_list[k]];
    for (unsigned int m=0;m<myParticles.size();m++){
      shared_ptr<DKinFitParticle>myParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]);
      FinalParticles.insert(myParticle);
      if (got_vertex_track==false){
	vert.SetVect(myParticles[m]->position());
	vert.SetT(myParticles[m]->t0());
	got_vertex_track=true;
      }
    } 
  }
  // Neutral particles
  for (unsigned int k=0;k<gammaHyps.size();k++){
    shared_ptr<DKinFitParticle>myGamma=dKinFitUtils->Make_DetectedParticle(1,0,0.,vert,gammaHyps[k]->momentum(),0.,gammaHyps[k]->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
  dKinFitter->Fit_Reaction();
}

void JEventProcessor_InsertStudies::TwoGammaAnalysis(unsigned int insert_count,
						     DKinFitter *dKinFitter,
						     double weight,bool missing_proton) const {
  DLorentzVector beam_kf;
  DLorentzVector proton_kf;
  vector<DLorentzVector>gammas_kf;
  
  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){
      switch ((*locParticleIterator)->Get_PID()){
      case 2212:  // Proton
	proton_kf=(*locParticleIterator)->Get_P4();
	break;
      case 1:  // Gamma
	gammas_kf.push_back((*locParticleIterator)->Get_P4());
	break;
      default:
	break;
      }
    }
  }

  //  double Ebeam=beam_kf.E();
  DLorentzVector twogam_kf=gammas_kf[0]+gammas_kf[1];
  double twogam_mass_kf=twogam_kf.M();
  
  if (missing_proton){
    TwoGammaMissingProton->Fill(twogam_mass_kf,weight);
  }

}

void JEventProcessor_InsertStudies::FourGammaAnalysis(unsigned int insert_count,
						      DKinFitter *dKinFitter,
						      double weight,bool missing_proton) const {
  DLorentzVector beam_kf;
  DLorentzVector proton_kf;
  vector<DLorentzVector>gammas_kf;
  
  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){
      switch ((*locParticleIterator)->Get_PID()){
      case 2212:  // Proton
	proton_kf=(*locParticleIterator)->Get_P4();
	break;
      case 1:  // Gamma
	gammas_kf.push_back((*locParticleIterator)->Get_P4());
	break;
      default:
	break;
      }
    }
  }

  //  double Ebeam=beam_kf.E();
  DLorentzVector fourgam_kf=gammas_kf[0]+gammas_kf[1]+gammas_kf[2]+gammas_kf[3];
  double fourgam_mass_kf=fourgam_kf.M();
 
  vector<vector<pair<unsigned int,unsigned int>> >doublets;
  MakeDoublets(doublets);

  if (missing_proton){
    FourGammaMissingProton->Fill(fourgam_mass_kf,weight);
  }
  else{
    FourGammaMass->Fill(fourgam_mass_kf,insert_count,weight);
  }
  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];

    if (missing_proton){
      FourGamma2DMissingProton->Fill(pair1.M(),pair2.M(),weight);
    }
    
    vector<DLorentzVector>wrong_pairs;
    wrong_pairs.push_back(gammas_kf[doublets[i][0].first]
			  +gammas_kf[doublets[i][1].first]); 
    wrong_pairs.push_back(gammas_kf[doublets[i][0].first]
			  +gammas_kf[doublets[i][1].second]);  
    wrong_pairs.push_back(gammas_kf[doublets[i][0].second]
			  +gammas_kf[doublets[i][1].first]); 
    wrong_pairs.push_back(gammas_kf[doublets[i][0].second]
			  +gammas_kf[doublets[i][1].second]);
    bool pion_overlap=false;
    for (unsigned int j=0;j<wrong_pairs.size();j++){
      if (TIGHT_PI0_CUT(wrong_pairs[j].M())){
	pion_overlap=true;
	break;
      }
    }
    if (pion_overlap==false){
      if (missing_proton){

      }
      else{
	FourGamma2DNoOverlap->Fill(pair1.M(),pair2.M(),weight);
      }
    }
    
    FourGamma2D->Fill(pair1.M(),pair2.M(),weight);
    if (pion_overlap==false){
      if (PI0_CUT(pair1.M())&&PI0_CUT(pair2.M())){
	if (missing_proton){

	}
	else{
	  TwoPi0NoOverlap->Fill(fourgam_mass_kf,weight);
	}
      } 
      if (ETA_CUT(pair1.M())&&ETA_CUT(pair2.M())){
	if (missing_proton){

	}
	else{
	  TwoEtaNoOverlap->Fill(fourgam_mass_kf,weight);
	}
      } 
      if ((PI0_CUT(pair1.M())&&ETA_CUT(pair2.M()))
	  || (PI0_CUT(pair2.M())&&ETA_CUT(pair1.M()))){
	if (missing_proton){

	}
	else{
	  EtaPi0NoOverlap->Fill(fourgam_mass_kf,weight);
	}
      }
    }      
    if (PI0_CUT(pair1.M())){
      if (missing_proton){

      }
      else {
	Pi02gNoVeto->Fill(fourgam_mass_kf,insert_count,weight);
      }

      DLorentzVector gamma1=gammas_kf[doublets[i][1].first]; 
      DLorentzVector gamma2=gammas_kf[doublets[i][1].second];
      DLorentzVector pi0g_1=pair1+gamma1;
      DLorentzVector pi0g_2=pair1+gamma2;
      if (missing_proton){

      }
      else {
	Pi02gVsPi0g->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	Pi02gVsPi0g->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	
	if (PI0_VETO_CUT(pair2.M())){
	  Pi02gWithVeto->Fill(fourgam_mass_kf,insert_count,weight);
	  TwoGammaOffPi0->Fill(pair2.M(),insert_count,weight);
	  Pi02gVsPi0gWithVeto->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	  Pi02gVsPi0gWithVeto->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	}
	if (pion_overlap==false){ 
	  Pi02gNoVetoNoOverlap->Fill(fourgam_mass_kf,insert_count,weight);
	  Pi02gVsPi0gNoOverlap->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	  Pi02gVsPi0gNoOverlap->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	
	  if (PI0_VETO_CUT(pair2.M())){
	    Pi02gWithVetoNoOverlap->Fill(fourgam_mass_kf,insert_count,weight);
	    TwoGammaOffPi0NoOverlap->Fill(pair2.M(),insert_count,weight);
	    Pi02gVsPi0gWithVetoNoOverlap->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	    Pi02gVsPi0gWithVetoNoOverlap->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	  }
	}
      }
    } 
    if (PI0_CUT(pair2.M())){
      if (missing_proton){
      }
      else{
	Pi02gNoVeto->Fill(fourgam_mass_kf,insert_count,weight); 
      }

      DLorentzVector gamma1=gammas_kf[doublets[i][0].first]; 
      DLorentzVector gamma2=gammas_kf[doublets[i][0].second];
      DLorentzVector pi0g_1=pair1+gamma1;
      DLorentzVector pi0g_2=pair1+gamma2;
      if (missing_proton){

      }
      else {
	Pi02gVsPi0g->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	Pi02gVsPi0g->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	
	if (PI0_VETO_CUT(pair1.M())){
	  Pi02gWithVeto->Fill(fourgam_mass_kf,insert_count,weight);
	  TwoGammaOffPi0->Fill(pair1.M(),insert_count,weight); 
	  Pi02gVsPi0gWithVeto->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	  Pi02gVsPi0gWithVeto->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	}
	if (pion_overlap==false){
	  Pi02gNoVetoNoOverlap->Fill(fourgam_mass_kf,insert_count,weight);
	  Pi02gVsPi0gNoOverlap->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	  Pi02gVsPi0gNoOverlap->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	  
	  if (PI0_VETO_CUT(pair1.M())){
	    Pi02gWithVetoNoOverlap->Fill(fourgam_mass_kf,insert_count,weight);
	    TwoGammaOffPi0NoOverlap->Fill(pair1.M(),insert_count,weight);
	    Pi02gVsPi0gWithVetoNoOverlap->Fill(pi0g_1.M(),fourgam_mass_kf,weight);
	    Pi02gVsPi0gWithVetoNoOverlap->Fill(pi0g_2.M(),fourgam_mass_kf,weight);
	  }
	}
      }

    }
    
  }

}

void JEventProcessor_InsertStudies::SixGammaAnalysis(DKinFitter *dKinFitter,double weight) const {
  DLorentzVector beam_kf;
  DLorentzVector proton_kf;
  vector<DLorentzVector>gammas_kf;
  
  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){
      switch ((*locParticleIterator)->Get_PID()){
      case 2212:  // Proton
	proton_kf=(*locParticleIterator)->Get_P4();
	break;
      case 1:  // Gamma
	gammas_kf.push_back((*locParticleIterator)->Get_P4());
	break;
      default:
	break;
      }
    }
  }
  vector<vector<pair<unsigned int,unsigned int>> > triplets;
  MakeTriplets(triplets);
  DLorentzVector sixgam_kf;
  for (unsigned int i=0;i<6;i++){
    sixgam_kf+=gammas_kf[i];
  }
  double sixgammamass=sixgam_kf.M();
  SixGammaMass->Fill(sixgammamass,weight);
  for (unsigned int i=0;i<triplets.size();i++){
    DLorentzVector pair1=gammas_kf[triplets[i][0].first]
      +gammas_kf[triplets[i][0].second];
    DLorentzVector pair2=gammas_kf[triplets[i][1].first]
      +gammas_kf[triplets[i][1].second];
    DLorentzVector pair3=gammas_kf[triplets[i][2].first]
      +gammas_kf[triplets[i][2].second];
    
    FourGamma2D_6g->Fill(pair1.M(),pair2.M(),weight); 
    FourGamma2D_6g->Fill(pair1.M(),pair3.M(),weight);
    FourGamma2D_6g->Fill(pair2.M(),pair3.M(),weight);

    vector<DLorentzVector>wrong_pairs;
    wrong_pairs.push_back(gammas_kf[triplets[i][0].first]
			  +gammas_kf[triplets[i][1].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][0].first]
			  +gammas_kf[triplets[i][1].second]);  
    wrong_pairs.push_back(gammas_kf[triplets[i][0].second]
			  +gammas_kf[triplets[i][1].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][0].second]
			  +gammas_kf[triplets[i][1].second]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][0].first]
			  +gammas_kf[triplets[i][2].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][0].first]
			  +gammas_kf[triplets[i][2].second]);  
    wrong_pairs.push_back(gammas_kf[triplets[i][0].second]
			  +gammas_kf[triplets[i][2].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][0].second]
			  +gammas_kf[triplets[i][2].second]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][1].first]
			  +gammas_kf[triplets[i][2].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][1].first]
			  +gammas_kf[triplets[i][2].second]);  
    wrong_pairs.push_back(gammas_kf[triplets[i][1].second]
			  +gammas_kf[triplets[i][2].first]); 
    wrong_pairs.push_back(gammas_kf[triplets[i][1].second]
			  +gammas_kf[triplets[i][2].second]);
    bool pion_overlap=false;
    for (unsigned int j=0;j<wrong_pairs.size();j++){
      if (TIGHT_PI0_CUT(wrong_pairs[j].M())){
	pion_overlap=true;
	break;
      }
    }
    if (pion_overlap==false){
      FourGamma2D_6g_NoOverlap->Fill(pair1.M(),pair2.M(),weight); 
      FourGamma2D_6g_NoOverlap->Fill(pair1.M(),pair3.M(),weight);
      FourGamma2D_6g_NoOverlap->Fill(pair2.M(),pair3.M(),weight);
    }

    vector<DLorentzVector>pi0s;
    vector<DLorentzVector>etas;
    vector<unsigned int>unused_pair;
    for (unsigned int j=0;j<triplets[i].size();j++){
      pair1=gammas_kf[triplets[i][j].first]+gammas_kf[triplets[i][j].second];
      if (PI0_CUT(pair1.M())){
	pi0s.push_back(pair1);
      }
      else if (ETA_CUT(pair1.M())){
	etas.push_back(pair1);
      }
      else {
	unused_pair.push_back(j);
      }
    }
    if (pi0s.size()==3){
      ThreePi0Mass->Fill(sixgammamass,weight);
      if (pion_overlap==false){
	ThreePi0MassNoOverlap->Fill(sixgammamass,weight);
      }
    }
    if (pi0s.size()==1 && etas.size()==2){
      TwoEtaPi0Mass->Fill(sixgammamass,weight);
      if (pion_overlap==false){
	TwoEtaPi0MassNoOverlap->Fill(sixgammamass,weight);

	DLorentzVector eta1pi0=pi0s[0]+etas[0];
	DLorentzVector eta2pi0=pi0s[0]+etas[1];
	TwoEtaPi0Dalitz->Fill(eta1pi0.M2(),eta2pi0.M2(),weight);

      }

    }
    if (pi0s.size()==2){
      if (etas.size()==1){
	Eta2Pi0Mass->Fill(sixgammamass,weight);
	if (pion_overlap==false){
	  Eta2Pi0MassNoOverlap->Fill(sixgammamass,weight);
	}
      }
      else{
	DLorentzVector gamma1=gammas_kf[triplets[i][unused_pair[0]].first];
	DLorentzVector gamma2=gammas_kf[triplets[i][unused_pair[0]].second];
	DLorentzVector pi0_1_g_1=pi0s[0]+gamma1;
	DLorentzVector pi0_2_g_2=pi0s[1]+gamma2;
	Pi0gVsPi0g->Fill(pi0_1_g_1.M(),pi0_2_g_2.M(),weight);	
	DLorentzVector pi0_1_g_2=pi0s[0]+gamma2;
	DLorentzVector pi0_2_g_1=pi0s[1]+gamma1;
	Pi0gVsPi0g->Fill(pi0_1_g_2.M(),pi0_2_g_1.M(),weight);
	if (pion_overlap==false){
	  Pi0gVsPi0gNoOverlap->Fill(pi0_1_g_1.M(),pi0_2_g_2.M(),weight);	
	  Pi0gVsPi0gNoOverlap->Fill(pi0_1_g_2.M(),pi0_2_g_1.M(),weight);
	}
      }
    }
  }

}

void JEventProcessor_InsertStudies::MakeTriplets(vector<vector<pair<unsigned int, unsigned int > > >&triplets) const { 
  vector<pair<unsigned int,unsigned int>>pairs;
  pairs.push_back(make_pair(0,1));
  pairs.push_back(make_pair(2,3));
  pairs.push_back(make_pair(4,5));
  triplets.push_back(pairs); 
  pairs.clear();
  pairs.push_back(make_pair(0,1));
  pairs.push_back(make_pair(2,4));
  pairs.push_back(make_pair(3,5));
  triplets.push_back(pairs);
  pairs.clear(); 
  pairs.push_back(make_pair(0,1));
  pairs.push_back(make_pair(2,5));
  pairs.push_back(make_pair(3,4));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,2));
  pairs.push_back(make_pair(1,3));
  pairs.push_back(make_pair(4,5));
  triplets.push_back(pairs);
  pairs.clear(); 
  pairs.push_back(make_pair(0,2));
  pairs.push_back(make_pair(1,4));
  pairs.push_back(make_pair(5,3));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,2));
  pairs.push_back(make_pair(1,5));
  pairs.push_back(make_pair(3,4));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,3));
  pairs.push_back(make_pair(2,1));
  pairs.push_back(make_pair(4,5));
  triplets.push_back(pairs);
  pairs.clear();  
  pairs.push_back(make_pair(0,3));
  pairs.push_back(make_pair(2,4));
  pairs.push_back(make_pair(1,5));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,3));
  pairs.push_back(make_pair(2,5));
  pairs.push_back(make_pair(4,1));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,4));
  pairs.push_back(make_pair(1,3));
  pairs.push_back(make_pair(2,5));
  triplets.push_back(pairs);
  pairs.clear(); 
  pairs.push_back(make_pair(0,4));
  pairs.push_back(make_pair(1,5));
  pairs.push_back(make_pair(2,3));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,4));
  pairs.push_back(make_pair(1,2));
  pairs.push_back(make_pair(3,5));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,5));
  pairs.push_back(make_pair(1,3));
  pairs.push_back(make_pair(4,2));
  triplets.push_back(pairs);
  pairs.clear();  
  pairs.push_back(make_pair(0,5));
  pairs.push_back(make_pair(1,4));
  pairs.push_back(make_pair(3,2));
  triplets.push_back(pairs);
  pairs.clear();
  pairs.push_back(make_pair(0,5));
  pairs.push_back(make_pair(1,2));
  pairs.push_back(make_pair(4,3));
  triplets.push_back(pairs);
}

void JEventProcessor_InsertStudies::MakeDoublets(vector<vector<pair<unsigned int, unsigned int > > >&doublets) const { 
  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);
}

void JEventProcessor_InsertStudies::FillThrownHistos(const DMCReaction *reaction,
						     vector<const DMCThrown*>&throwns,
						     vector<const DTAGMGeometry *>&tagm_geom,
						     vector<const DTAGHGeometry *>&tagh_geom) const {
  double mc_beam_energy=reaction->beam.energy();
  double tagged_mc_beam_energy=0.;  
  unsigned int counter=0,column=0;
  if (tagm_geom[0]->E_to_column(mc_beam_energy,column)){
    tagged_mc_beam_energy=0.5*(tagm_geom[0]->getEhigh(column)
			       +tagm_geom[0]->getElow(column));
  }
  else if (tagh_geom[0]->E_to_counter(mc_beam_energy,counter)){
    tagged_mc_beam_energy=0.5*(tagh_geom[0]->getEhigh(counter)
			       +tagh_geom[0]->getElow(counter));
  }
  if (tagged_mc_beam_energy>MIN_BEAM_E && tagged_mc_beam_energy<MAX_BEAM_E){
    MCBeam->Fill(tagged_mc_beam_energy);

    vector<DLorentzVector>protons,etas,pi0s,gammas;
    for (unsigned int i=0;i<throwns.size();i++){
      if (throwns[i]->mech==0){     
	switch(throwns[i]->type){
	case Proton:
	  protons.push_back(throwns[i]->lorentzMomentum());
	  break;
	case Eta:
	  etas.push_back(throwns[i]->lorentzMomentum());
	  break;
	case Pi0:
	  pi0s.push_back(throwns[i]->lorentzMomentum());
	  break;
	case Gamma:
	  gammas.push_back(throwns[i]->lorentzMomentum());
	  break;
	default:
	  break;
	}
      }
    }
    if (protons.size()==1){
      DLorentzVector beam_v4(0.,0.,tagged_mc_beam_energy,tagged_mc_beam_energy);
      DLorentzVector target_v4(0.,0.,0.,ParticleMass(Proton));
      DLorentzVector proton_v4=protons[0];
      DLorentzVector missing_v4=beam_v4+target_v4-proton_v4;
      MCMissingMass->Fill(missing_v4.M());

    }
  }
}

// Run the kinematic fitter requiring energy and momentum conservation
void JEventProcessor_InsertStudies::DoKinematicFitMissingProton(const DBeamPhoton *beamphoton,
						   map<Particle_t,vector<const DTrackTimeBased *> >&particles,
						   vector<const DNeutralParticleHypothesis*>&gammaHyps,
						   DKinFitUtils_GlueX *dKinFitUtils,
						   DKinFitter *dKinFitter) const {
  set<shared_ptr<DKinFitParticle>> InitialParticles, FinalParticles;
  
  shared_ptr<DKinFitParticle>myBeam=dKinFitUtils->Make_BeamParticle(beamphoton);
  shared_ptr<DKinFitParticle>myTarget=dKinFitUtils->Make_TargetParticle(Proton);
  
  InitialParticles.insert(myBeam);  
  InitialParticles.insert(myTarget);

  // Missing proton
  shared_ptr<DKinFitParticle>myParticle=dKinFitUtils->Make_MissingParticle(Proton);
  
  // Charged particles
  Particle_t particle_list[4]={PiPlus,PiMinus,Positron,Electron};
  TLorentzVector vert(0,0,65.,0.);
  bool got_vertex_track=false;
  for (unsigned int k=0;k<4;k++){ 
    vector<const DTrackTimeBased *>myParticles=particles[particle_list[k]];
    for (unsigned int m=0;m<myParticles.size();m++){
      shared_ptr<DKinFitParticle>myParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]);
      FinalParticles.insert(myParticle);
      if (got_vertex_track==false){
	vert.SetVect(myParticles[m]->position());
	vert.SetT(myParticles[m]->t0());
	got_vertex_track=true;
      }
    } 
  }
  // Neutral particles
  for (unsigned int k=0;k<gammaHyps.size();k++){
    shared_ptr<DKinFitParticle>myGamma=dKinFitUtils->Make_DetectedParticle(1,0,0.,vert,gammaHyps[k]->momentum(),0.,gammaHyps[k]->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
  dKinFitter->Fit_Reaction();
}