// $Id$
//
//    File: JEventProcessor_neutron.h
// Created: Wed Jun  1 13:54:40 EDT 2022
// Creator: staylor (on Linux ifarm1802.jlab.org 3.10.0-1160.11.1.el7.x86_64 x86_64)
//

#ifndef _JEventProcessor_neutron_
#define _JEventProcessor_neutron_

#include <JANA/JEventProcessor.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TDirectory.h>
#include <TRACKING/DTrackTimeBased.h>
#include <PID/DNeutralParticleHypothesis.h>
#include <PID/DBeamPhoton.h>
#include <PID/DChargedTrack.h>
#include <PID/DNeutralParticle.h>
#include <KINFITTER/DKinFitter.h>
#include <ANALYSIS/DKinFitUtils_GlueX.h>
#include <TAGGER/DTAGMGeometry.h>
#include <TAGGER/DTAGHGeometry.h>
#include <TRACKING/DMCThrown.h>
#include "../include/mass_cuts.h"


class JEventProcessor_neutron:public jana::JEventProcessor{
public:
JEventProcessor_neutron();
~JEventProcessor_neutron();
const char* className(void){return "JEventProcessor_neutron";}

private:
  jerror_t init(void);						///< Called once at program start.
  jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber);	///< Called everytime a new run number is detected.
  jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber);	///< Called every event.
  jerror_t erun(void);						///< Called everytime run number changes, provided brun has been called.
  jerror_t fini(void);						///< Called after last event of last event source has been processed.

void DoKinematicFit(double t0_rf,const DVector3 &vertex,
		    const DBeamPhoton *beamphoton,
		    vector<const DNeutralParticleHypothesis *>&neutronHyps,
			map<Particle_t,vector<const DTrackTimeBased *> >&chargedParticles,
		      vector<const DNeutralParticleHypothesis*>&gammas,
		      DKinFitUtils_GlueX *dKinFitUtils,
		      DKinFitter *dKinFitter,
		      bool isMissing) const;
  void FillChargedParticleVectors(vector<const DChargedTrack *>&tracks,
				    map<Particle_t, vector<const DTrackTimeBased *> >&particles
				    ) const;

  void FillNeutralParticleVectors(double &unused_energy,double t0_rf,
				  const DVector3 &vertex,
				  vector<const DNeutralParticle*>&neutrals,
				  map<Particle_t,vector<const DNeutralParticleHypothesis*>>&neutralParticles) const; 
  static bool SortParticleProbability(const pair<double,Particle_t>a,
					const pair<double,Particle_t>b){
return a.first>b.first;
  }
  void GetKF4vectors(DKinFitter *dKinFitter,DLorentzVector &beam_kf,
		     DLorentzVector &missing_kf,
		     map<Particle_t,vector<DLorentzVector>>&final_kf) const;
  void FillGJHisto(double weight,const DLorentzVector &beam,
		   const DLorentzVector &meson,const DLorentzVector &analyzer,
		   TH2F *histo) const;
  void MakeMesonPairs(vector<DLorentzVector>&gammas_kf,double weight,
		      TH2F *histo);
  void MakeTriplets();

  void NeutronTwoGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
  void NeutronThreeGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
  void NeutronFourGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight);
  void NeutronSixGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
  void MissingPionAnalysis(const DTrackTimeBased *pion_track,
			   const DLorentzVector &missing_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,
			   double weight) const;
  void NeutronKpPipPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
  void NeutronKpPipPim1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				     map<Particle_t,vector<DLorentzVector>>&final_kf,
				     double weight) const;
  void NeutronKpPipPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				     map<Particle_t,vector<DLorentzVector>>&final_kf,
				     double weight) const;

  void Neutron2PipPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			      map<Particle_t,vector<DLorentzVector>>&final_kf,
			      double weight) const;
  void Neutron2PipPim1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				    map<Particle_t,vector<DLorentzVector>>&final_kf,
				    double weight) const;
  void Neutron2PipPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				    map<Particle_t,vector<DLorentzVector>>&final_kf,
				    double weight) const;
  void Neutron2PipPim4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				    map<Particle_t,vector<DLorentzVector>>&final_kf,
				    double weight);
  void TwoKpKmN2GAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const;
  void KpN2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void TwoNeutronAnalysis(const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const;
  void TwoNeutronPipPimAnalysis(const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void PAntiPNPipAnalysis(const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const;
  void PAntiPNPip2GAnalysis(const DLorentzVector &beam_kf,
			    map<Particle_t,vector<DLorentzVector>>&final_kf,
			    double weight) const;
  void TwoProtonPim2GAnalysis(const DLorentzVector &beam_kf,
			      map<Particle_t,vector<DLorentzVector>>&final_kf,
			      double weight) const;

  
  double PI0_CUT_VALUE,ETA_CUT_VALUE,ETAPRIME_CUT_VALUE,NUM_SIGMA_BG;
  double FCAL_THRESHOLD,SPLIT_CUT,FCAL_RADIAL_CUT;
  double BCAL_THRESHOLD,FCAL_POS_CUT,BCAL_Z_CUT,BETA_CUT,CL_CUT;
  DVector3 mFCALCenter;

  TH2F *MCM_vs_E,*MCM_vs_t,*MC_EtaCosThetaGJ_n;

  TH1F *UnusedEnergy;
  TH1F *NeutralBeta,*GammaRFTimeDiff,*BeamPhotonTimeDiff,*MissingPionChiSq;
  TH1F *MissingPionCL,*ExclusiveCL,*ExclusiveChiSq,*TwoGammaMass_NPip;
  TH2F *PipPi0_NPi0_Dalitz,*Pi0CosThetaGJ_n,*PipEta_NEta_Dalitz,*EtaCosThetaGJ_n;
  TH2F *Pi0CosThetaGJ_n_bg,*EtaCosThetaGJ_n_bg,*DeltaPOverP_vs_P;
  TH1F *TwoGammaMass_missingPip,*NKs,*SigmamPip,*SigmapPim;
  TH2F *PipPim_vs_nKp,*NPim_vs_KpPip,*NPip_vs_KpPim;
  TH2F *DeltaPhi_vs_P,*DeltaTheta_vs_P,*DeltaPOverP_vs_theta;
  TH1F *FourGammaMass_NPip,*Pi0Pi0_NPip,*EtaEta_NPip,*Pi0Eta_NPip;
  TH2F *FourGamma2d_NPip,*PipPim_vs_NPip,*PipPim_vs_NPip_1G;
  TH1F *NeutronPipMass_2pip,*NeutronPipMass_2pip1G,*PipPim1G_N;
  TH1F *SixGammaMass_NPi,*Pip6G_N,*ThreePi0Mass_NPi,*PipEta_N6g,*TwoPi0EtaMass_NPi;
  TH1F *TwoEtaPi0Mass_NPi,*ThreeEtaMass_NPi;
  TH2F *TwoEtaPi0Dalitz_NPi,*Pi0gVsPi0g_NPi;
  TH3F *SixGamma3D_NPi;
  TH1F *TwoPipPimPi0_N,*PipPimPi0_n,*EtaPip_N,*OmegaPip_N2g;
  TH2F *NPi0_vs_2PipPim,*PipPim_vs_NPip_with_Pi0;
  TH2F *SigmaPlusPi0_SigmaPlusK0_Dalitz,*PipPim_vs_NPip_with_Eta;
  TH2F *SigmaPlusEta_SigmaPlusK0_Dalitz;
  TH1F *PipPimEta_n,*TwoPipPimEta_N,*EtaPrimePip_N,*EtaPrimePip_N_NU;
  TH1F *PipPimPi0_N_with_Pi0,*EtaPi0_N_eta_3pi,*PipPimPi0_N_with_Eta;
  TH1F *PipPimEta_N_with_Pi0;
  TH2F *Omega2PI_N_Dalitz,*PipPim_vs_NPip_with_2Pi0,*PipPim_vs_NPip_with_EtaPi0;
  TH1F *FourGammaMass_N3Pi,*TwoPi0_N3Pi,*EtaPi0_N3Pi;
  TH2F *FourGamma2d_N3Pi,*NEta_vs_2PipPim;
  TH1F *TwoGammaMass_n2pip,*NeutronPipMass_2pip2g;
  TH1F *TwoGammaMass_n2pip_NU;
  TH1F *TwoGamma_2kpkmn,*NPi0_2kpkm,*SigmaPlusGamma,*SigmaMinusGamma;
  TH2F *NPipSqVsNPimSq_1G,*NPipSqVsNPimSq_2G;
  TH1F *TwoGamma_with_KNPipPim;
  TH2F *SigmaPlusPi0_vs_KPim,*SigmaPlusPi0_SigmaPlusPim_Dalitz;
  TH2F *SigmaPlusPim_vs_KPi0,*SigmaMinusPip_vs_KPi0,*SigmaMinusPi0_vs_KPip;
  TH2F *SigmaMinusPi0_SigmaMinusPip_Dalitz,*SigmaPlusEta_vs_KPim;
  TH1F *SigmamPipPi0,*SigmapPimPi0;
  TH2F *SigmaPlusEta_SigmaPlusPim_Dalitz,*SigmaPlusPim_vs_KEta;
  TH2F *SigmaMinusPip_vs_KEta;
  TH2F *SigmaMinusEta_vs_KPip,*SigmaMinusEta_SigmaMinusPip_Dalitz;
  TH1F *SigmamPipEta,*SigmapPimEta;
  TH1F *TwoGammaMass_NK,*N2gMass,*NPi0Mass,*NEtaMass;
  TH2F *ThreeGam_vs_NPip;
  TH1F *OmegaPip_n3g;
  TH2F *PipEtaMass_vs_E,*PipEtaMass_vs_t,*PipPi0Mass_vs_E,*PipPi0Mass_vs_t;
  TH2F *PipEtaPrimeMass_vs_E,*PipEtaPrimeMass_vs_t,*EtaPrimeCosThetaGJ_n;
  TH1F *TwoNeutronMass,*AntiPNPip;
  TH2F *NantiN_vs_t,*NantiN_vs_E;
  TH2F *NeutronPip_vs_NeutronPim,*AntiPNPip_vs_2g,*NPip_vs_AntiPPi0;
  TH2F *NPi0_vs_AntiPPiP;

  TH2F *TwoNeutronMomentum2D,*NeutronPVsTheta,*ProtonPVsThetaWithN;
  TH2F *TwoNeutronVertex,*PNDalitz;
  TH1F *TwoGammaMass_with_2pnpim;
  TH2F *PPim_vs_AntiNPi0,*PPi0_vs_AntiNPim;
 
  enum particle_pair_t{
    Pi0Pi0_,
    Pi0Eta_,
    EtaEta_,
    Pi0EtaPrime_,
  };
  map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>particle_pairs;
  vector<vector<pair<unsigned int,unsigned int>> >doublets;
  vector<vector<pair<unsigned int,unsigned int>> > triplets;
 

};

#endif // _JEventProcessor_neutron_