// $Id$
//
//    File: JEventProcessor_survey.h
// Created: Mon Dec  6 14:28:07 EST 2021
// Creator: staylor (on Linux ifarm1901.jlab.org 3.10.0-1160.11.1.el7.x86_64 x86_64)
//

#ifndef _JEventProcessor_survey_
#define _JEventProcessor_survey_

#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 <TOF/DTOFPoint.h>
#include <KINFITTER/DKinFitter.h>
#include <ANALYSIS/DKinFitUtils_GlueX.h>
#include <TAGGER/DTAGMGeometry.h>
#include <TAGGER/DTAGHGeometry.h>
#include <HDDM/DEventHitStatistics.h>
#include "../include/mass_cuts.h"

class JEventProcessor_survey:public jana::JEventProcessor{
 public:
  JEventProcessor_survey();
  ~JEventProcessor_survey();
  const char* className(void){return "JEventProcessor_survey";}
  
 private:
  enum particle_pair_t{
    Pi0Pi0_,
    Pi0Eta_,
    EtaEta_,
    Pi0EtaPrime_,
    EtaEtaPrime_,
  };

  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 FillChargedParticleVectors(vector<const DChargedTrack *>&tracks,
				  map<Particle_t,vector<const DTrackTimeBased *>>&chargedParticles
				  ) const;
  void FillNeutralParticleVectors(double t0_rf,const DVector3 &vertex,
				  vector<const DTOFPoint*>&tofs,
				  vector<const DNeutralParticle*>&neutrals,
				  map<Particle_t,vector<const DNeutralParticleHypothesis*>>&neutralParticles,
				  int &num_in_insert,
				  int &num_tof_match
				  ) const;
  void DoKinematicFit(const DBeamPhoton *beamphoton,
		      map<Particle_t,vector<const DTrackTimeBased *> >&chargedParticles,
		      map<Particle_t,vector<const DNeutralParticleHypothesis*>>&neutralParticles,
		      DKinFitUtils_GlueX *dKinFitUtils,
		      DKinFitter *dKinFitter,
		      bool isMissing=false) const;
  void GetKF4vectors(DKinFitter *dKinFitter,DLorentzVector &beam_kf,
		     DLorentzVector &missing_kf,
		     map<Particle_t,vector<DLorentzVector>>&final_kf) const;
  static bool SortParticleProbability(const pair<double,Particle_t>a,
				      const pair<double,Particle_t>b){
    return a.first>b.first;
  }
  bool MakeMesonPairs(vector<unsigned int>&paired_doublets,
		      vector<DLorentzVector>&gammas_kf,
		      double weight,TH2F *histo);
  void FillPairedMesonMassHistos(double weight,double mass,
				 TH1F *myPi0Pi0,TH1F *myPi0Pi0bg,
				 TH1F *myPi0Eta,TH1F *myPi0Etabg,
				 TH1F *myEtaEta,TH1F *myEtaEtabg
				 );
  void FillGFHisto(double weight,const DLorentzVector &beam,
		   const DLorentzVector &meson,const DLorentzVector &analyzer,
		   TH2F *histo) const;
  
  void TwoGammaAnalysis(DLorentzVector beam_kf,
			map<Particle_t,vector<DLorentzVector>>&final_kf,
			double weight) const;
  void FourGammaAnalysis(int num_fdc,int num_in_insert,int num_tof_match,
			 DLorentzVector beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&final_kf,
			 double weight);
  void SixGammaAnalysis(DLorentzVector beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&final_kf,
			 double weight) const;
  void NeutronTwoGammaAnalysis(DLorentzVector beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
  void NeutronFourGammaAnalysis(DLorentzVector beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight);
  void EpEmAnalysis(DLorentzVector beam_kf,
		    map<Particle_t,vector<DLorentzVector>>&final_kf,
		    double weight) const;
  void EpEmTwoGammaAnalysis(DLorentzVector beam_kf,
			    map<Particle_t,vector<DLorentzVector>>&final_kf,
			    double weight) const;
  void PipPimTwoGammaAnalysis(DLorentzVector beam_kf,
			      map<Particle_t,vector<DLorentzVector>>&final_kf,
			      double weight) const;
  void PipPimFourGammaAnalysis(DLorentzVector beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			double weight);  
  void KpKmTwoGammaAnalysis(DLorentzVector beam_kf,
			    map<Particle_t,vector<DLorentzVector>>&final_kf,
			    double weight) const;
  void KpKmFourGammaAnalysis(DLorentzVector beam_kf,
			     map<Particle_t,vector<DLorentzVector>>&final_kf,
			     double weight);
  void KpKmPipPimTwoGammaAnalysis(DLorentzVector beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,
				 double weight) const;
  void KpKmPipPimFourGammaAnalysis(DLorentzVector beam_kf,
				   map<Particle_t,vector<DLorentzVector>>&final_kf,
				   double weight);
  void FillSigma0Histos(const DLorentzVector &lambda,
			const DLorentzVector &lambda_bar,
			DLorentzVector gamma1,
			DLorentzVector gamma2,
			DLorentzVector proton,
			DLorentzVector antiproton,
			double weight) const ;
  void AntiProtonPipPim2GammaAnalysis(DLorentzVector beam_kf,
				      map<Particle_t,vector<DLorentzVector>>&final_kf,
				      double weight) const;
  void NeutronKpPipPimAnalysis(DLorentzVector beam_kf,
			       map<Particle_t,vector<DLorentzVector>>&final_kf,
			       double weight) const;
    
  map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>particle_pairs;
  map<particle_pair_t,vector<double>>bg_weights; 
  map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>bg_particle_pairs;
  vector<vector<pair<unsigned int,unsigned int>> >doublets;
  
  double PI0_CUT_VALUE,ETA_CUT_VALUE,ETAPRIME_CUT_VALUE,NUM_SIGMA_BG;
  double FCAL_THRESHOLD,SPLIT_CUT,FCAL_RADIAL_CUT,MIN_BEAM_E,MAX_BEAM_E;
  double BCAL_THRESHOLD,FCAL_POS_CUT,BCAL_Z_CUT,BETA_CUT,CL_CUT,TOF_FCAL_CUT;
  DVector3 mFCALCenter;

  TH1F *MC_meson_mass,*MC_beam_E,*MC_beam_E_all,*MC_t;

  TH1F *NeutralBeta,*GammaRFTimeDiff,*BeamPhotonTimeDiff;
  TH1F *ExclusiveCL,*MissingProtonCL;
  
  TH2F *PipPi0_NPi0_Dalitz,*PipEta_NEta_Dalitz,*FourGamma2d_NoMesonPairs;
  TH1F *TwoGammaMass,*TwoGammaMass_NPip,*TwoGammaMass_KpKm;
  TH2F *Pi02g_vs_Pi0g,*Pi02g_vs_2g;
  TH1F *FourGammaNoPairs_FDCcut,*FourGammaNoPairs_TOFcut;
  TH1F *FourGammaMass,*FourGammaNoPairs,*FourGammaMass_NPip,*FourGammaMass_KpKm;
  TH2F *FourGamma2d,*FourGamma2d_NPip,*FourGamma2d_KpKm;
  TH1F *KpKmMass,*KpKmMass_4g,*TwoGammaMass_PipPim,*PipPimPi0,*PipPimEta;
  TH2F *TwoGamma_vs_KpKm,*FourGamma2d_PipPim;
  TH1F *PipPimMass_4g,*FourGammaMass_PipPim;
  TH1F *Pi0Pi0,*Pi0Pi0bg,*Pi0Eta,*Pi0Etabg,*EtaEta,*EtaEtabg;
  TH1F *Pi0Pi0_NPip,*Pi0Pi0bg_NPip,*Pi0Eta_NPip,*Pi0Etabg_NPip;
  TH1F *EtaEta_NPip,*EtaEtabg_NPip;
  TH1F *Pi0Pi0_PipPim,*Pi0Pi0bg_PipPim,*Pi0Eta_PipPim,*Pi0Etabg_PipPim;
  TH1F *EtaEta_PipPim,*EtaEtabg_PipPim;
  TH1F *Pi0Pi0_KpKm,*Pi0Pi0bg_KpKm,*Pi0Eta_KpKm,*Pi0Etabg_KpKm;
  TH1F *EtaEta_KpKm,*EtaEtabg_KpKm;
  TH1F *Pi0Pi0_Phi_bg,*Pi0Pi0bg_Phi_bg,*Pi0Eta_Phi_bg,*Pi0Etabg_Phi_bg;
  TH1F *EtaEta_Phi_bg,*EtaEtabg_Phi_bg;
  TH1F *Pi0Pi0_Phi,*Pi0Pi0bg_Phi,*Pi0Eta_Phi,*Pi0Etabg_Phi;
  TH1F *EtaEta_Phi,*EtaEtabg_Phi;
  TH1F *EpEmMass,*PipPimPi0_KpKm,*PipPimEta_KpKm;
  TH2F *TwoGamma_vs_EpEm,*PPim_vs_aPPip,*AntiLambdaGamma_vs_LambdaGamma;
  TH2F *Pi0CosThetaGF,*Pi0CosThetaGF_bg,*EtaCosThetaGF,*EtaCosThetaGF_bg;
  TH2F *Pi0CosThetaGF_n,*Pi0CosThetaGF_n_bg,*EtaCosThetaGF_n,*EtaCosThetaGF_n_bg;
  TH2F *PipPim_vs_nKp,*NPim_vs_KpPip,*NPip_vs_KpPim;
  TH1F *NKs,*SigmapPim,*SigmamPip;
  TH1F *Sigma0AntiSigma0;
  TH2F *EtaPhiDalitz,*FourGamma2d_KpKmPipPim,*TwoGamma_vs_KpKm_with_PipPim;
  TH1F *SixGammaMass,*KpKmMass_4gPipPim,*FourGammaMass_KpKmPipPim;
  TH1F *Pi0Pi0_KpKmPipPim,*Pi0Pi0bg_KpKmPipPim,*Pi0Eta_KpKmPipPim,*Pi0Etabg_KpKmPipPim;
  TH1F *EtaEta_KpKmPipPim,*EtaEtabg_KpKmPipPim;
  TH2F *EtaPrimeCosThetaGF_with_Phi,*EtaCosThetaGF_with_Phi,*OmegaCosThetaGF_with_Phi;
  TH2F *EtaPrimeCosThetaGF_with_Phi_bg,*EtaCosThetaGF_with_Phi_bg,*OmegaCosThetaGF_with_Phi_bg;
  TH1F *PipPimPi0_KpKm_bg,*PipPimEta_KpKm_bg;
  TH1F *CosThetaLambdaRest,*CosThetaAntiLambdaRest;
  TH1F *NFDCPseudos,*TofFcaldt;
  TH2F *TOFdE_with_FCAL;
  TH2F *TofFcaldxdy;
  TH1F * FourGamma_Insertcut,*FourGammaNoPairs_Insertcut;
};

#endif // _JEventProcessor_survey_