// $Id$
//
//    File: DEventProcessor_strangeness.h
// Created: Thu Jul  7 08:30:20 EDT 2016
// Creator: staylor (on Linux gluon111.jlab.org 2.6.32-431.20.3.el6.x86_64 x86_64)
//

#ifndef _JEventProcessor_strangeness_
#define _JEventProcessor_strangeness_
//#include <RooDataSet.h>
//#include <RooRealVar.h>
#include <JANA/JEventProcessor.h>
#include <JANA/JApplication.h>

#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>

#include <TRIGGER/DL1MCTrigger.h>
#include <TRIGGER/DL1Trigger.h>
#include <TRIGGER/DTrigger.h>
#include <PID/DChargedTrack.h>
#include <PID/DNeutralParticleHypothesis.h>
#include <PID/DChargedTrackHypothesis.h>
#include <PID/DNeutralParticle.h>
#include <PID/DNeutralShower.h>
#include <PID/DBeamPhoton.h>
#include <PID/DMCReaction.h>
#include <TRACKING/DMCThrown.h>
#include <TRACKING/DTrackTimeBased.h>
#include <FCAL/DFCALShower.h>
#include <BCAL/DBCALShower.h>
#include <TOF/DTOFPoint.h>
#include <PID/DDetectorMatches.h>
#include <PID/DVertex.h>
#include <DVector3.h>
#include <DLorentzVector.h>

#include <KINFITTER/DKinFitter.h>
#include <ANALYSIS/DKinFitUtils_GlueX.h>

using namespace jana;
using namespace std;

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

 private:
 enum particle_pair_t{
    Pi0Pi0_,
    Pi0Eta_,
    EtaEta_,
    Pi0EtaPrime_,
  };

jerror_t init(void);						///< Called once at program start.
  jerror_t brun(jana::JEventLoop* locEventLoop, int32_t locRunNumber);	///< Called every time a new run number is detected.
  jerror_t evnt(jana::JEventLoop* locEventLoop, uint64_t locEventNumber);	///< Called every event.
  jerror_t erun(void);						///< Called every time run number changes, provided brun has been called.
  jerror_t fini(void);						///< Called after last event of last event source has been processed.
  
  //For non-custom reaction-independent histograms, it is recommended that you simultaneously run the monitoring_hists plugin instead of defining them here
  
  //For making EventStore skims
  int Get_FileNumber(JEventLoop* locEventLoop) const;
  ofstream dEventStoreSkimStream;

  static bool SortParticleProbability(const pair<double,Particle_t>a,
				      const pair<double,Particle_t>b){
    return a.first>b.first;
  }
 
  bool DoKinematicFit(double t0_rf,const DVector3 &vertex,
		      const DBeamPhoton *beamphoton,
		      map<Particle_t,vector<const DTrackTimeBased *> >&chargedParticles,
		      vector<const DNeutralParticleHypothesis*>&gammaParticles,
		      DKinFitUtils_GlueX *dKinFitUtils,
		      DKinFitter *dKinFitter) const;
  void FillChargedParticleVectors(vector<const DChargedTrack *>&tracks,map<Particle_t,vector<const DTrackTimeBased *>>&chargedParticles) const;
  void FillGammaParticleVector(double &unused_energy,double t0_rf,
			       const DVector3 &vertex,
			       vector<const DNeutralParticle*>&neutrals,
			       vector<const DNeutralParticleHypothesis*>&gammaParticles
			       ) const;
  void GetKF4vectors(DKinFitter *dKinFitter,DLorentzVector &beam_kf,
		     map<Particle_t,vector<DLorentzVector>>&final_kf) const;
  bool MakeMesonPairs(vector<DLorentzVector>&gammas_kf,double weight,
		      TH2F *histo);

  void KpKmProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const; 
  void KpKm1GammaProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void KpKm2GammaProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void KpKm3GammaProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void KpKm4GammaProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight);
  void KpKm6GammaProtonAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void FillKpKm4GHistos(double unused_energy,const DLorentzVector &pair1, 
			const DLorentzVector &pair2,
			const DLorentzVector &kp_kf, 
			const DLorentzVector &km_kf, double weight) const;
  void KpKmPipPim4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight);
  void TwoKpTwoKmAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const;
  void TwoKpKmPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const; 
  void TwoKpKmPim1GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const; 
  void TwoKpKm2PimPipAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			      map<Particle_t,vector<DLorentzVector>>&final_kf,
			      double weight) const;
  void KpKmPipPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,
				double weight) const;
  void KpKmPipPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,
			  double weight) const;
  void TwoKminusTwoPiPlusAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				  map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
		     map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPimGammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPim3GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPim4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void TwoKpTwoPimAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void TwoKpTwoPim2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void Kp2PimPipAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void Kp2PimPip2GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void Kp2PimPip4GammaAnalysis(double unused_energy,const DLorentzVector &beam_kf,map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void Kp3Pim2PipAnalysis(double unused_energy,const DLorentzVector &beam_kf,
			  map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;

  void MakeDoublets();
  void MakeTriplets();

  TH2F *FourGamma2d_KpKmPipPim,*PKmPip_vs_KpPim2Pi0,*PKmPip_vs_KpPim2Pi0_Phicut;
  TH1F *PipPimPi0_with_PhiPi0,*PipPimPi0_with_PhiPi0_bg;
  TH1F *PipPimPi0_with_PhiPi0_2pi0bg,*KpKmMass_with_2pi0pippim_bg;
  TH2F *PKmPip_vs_KpPim2Pi0_bg,*PKmPip_vs_KpPim2Pi0_Phicut_bg;
  TH1F *KpKmMass_with_2pi0pippim;
  TH1F *NeutralBeta,*KinFitCL,*KinFitChiSq;
  TH2F *BeamPhotonTimeDiff,*EOverPPlus,*EOverPPlus_BCAL,*PositronTOF;
  TH2F *ProtonThetaVsP,*ProtonRvsZ,*ProtondEdxSC;
  TH2F *ProtondEdxFDC,*ProtondEdx,*ProtondEdxCutFDC,*ProtondEdxCut,*ProtonTOFdt;
  TH2F *ProtonTOFdtNodEdx;
  TH2F *ProtonBCALdtWithdEdxCut,*ProtonBCALdtWithFDCdEdxCut;
  TH2F *ProtonTOFdtWithdEdxCut,*ProtonFCALdtWithdEdxCut;
  TH2F *ProtonFCALdt;
  TH2F *EOverPMinus,*EOverPMinus_BCAL,*ElectronTOF;
  TH2F *KPlusdEdx,*KPlusTOF,*KPlusBCALTOF;
  TH2F *KMinusdEdx,*KMinusTOF,*KMinusBCALTOF;
  TH1F *ProtonGammaTimeDiff,*PhiGammaMass;
  TH1F *KpKm2GammaMass,*TwoGammaMass_KpKm,*KpKmMass_2g,*PhiPi0Mass,*PhiEtaMass;
  TH1F *PhiEtaPrimeMass;
  TH1F *KpKmPi0Mass,*KpKmPi0Mass_baryon_cut,*KstarKm,*KstarKp,*KpKmEtaMass;
  TH2F *KpKmPi0Dalitz,*KpKmMass_vs_2GammaMass,*PPi0_vs_KpKm;
  TH1F *KpKm_with_Pi0,*KpKm3GammaMass,*TwoGammaMass_kpkm3g,*ThreeGammaMass_KpKm;
  TH1F *PhiOmegaMass,*PhiOmegaMass_with_pi0;
  TH2F *KpKmMass_vs_3GammaMass,*CLCorrelation,*KpKmVertex;
  TH1F *KpKmReducedChiSq;
  TH1F *KpKmMass,*KmPMass;
  TH2F *KpKmPDalitz,*KpKmMass_vs_Ebeam,*KpKm_vs_t[48],*KpCosThetaH_vs_t[48];
  TH2F *KpCosThetaVsPhi[48];
  TH1F *KpKm4GammaMass,*FourGammaMass_KpKm,*TwoPi0Mass_KpKm,*EtaPi0Mass_KpKm;
  TH1F *KpKmMass_4g,*KpKmMass_with_2pi0,*KpKmMass_with_etapi0,*KstarKstar;
  TH2F *FourGamma2d_KpKm,*KpPi0_vs_KmPi0,*PhiPi0Pi0Dalitz;
  TH2F *FourGamma2d_bg;
  TH1F *KpKm6GammaMass,*SixGammaMass_KpKm,*KpKmMass_6g,*KpKmEtaMass_6g_kf;
  TH2F *KpKm_vs_6g;
  TH1F *PhiEta_6g,*F1pi0_KpKm;
  TH1F *TwoGammaMass_KpKm6g_kf,*TwoPi0EtaMass_KpKm6g_kf,*TwoEtaPi0Mass_KpKm6g_kf;
  TH1F *ThreeEtaMass_KpKm6g_kf,*Pi0Gamma_kpkm,*ThreePi0Mass_KpKm6g_kf;
  TH1F *TwoKpTwoKmMass,*P2KmMass_2kp,*PKmMass_2kp,*TwoPhiMass;
  TH1F *PKmMass_2kp_NoPhis;
  TH2F *KpKm_vs_KpKm,*PKmPipMass_vs_KpPimPi0Mass;
  TH2F *PKmPipMass_vs_KpPimMass,*KmPipMass_vs_KpPimMass,*KmKpMass_vs_PipPimMass;
  TH1F *PhiPipPimMass;
  TH1F *TwoGamma_PKpKmPipPim;
  TH2F *KpKmMass_vs_PipPimPi0Mass,*KpKmMass_vs_PipPimEtaMass;
  TH1F *PhiOmegaMass_pippim,*PhiEtaMass_pippim,*PhiEtaPrimeMass_pippim;
  TH2F *ProtonBCALdt;
  
  TH1F *PhiKpKmMass;
  TH1F *KpKmGammaMass,*KmPMass_1g,*KpKmMass_1g;
  TH2F *KpKmGammaMass_vs_Ebeam,*KpThetaVsPhi_lab;
  TH2F *MCProtonThetaVsP,*MCMissingMassVsEbeam,*KpCosThetaVsPhi_mc[48];
  TH1F *MCProtonP,*MCProtonAcceptedP,*MCBeam,*MCBeamTagged,*MCMissingMass;
  TH2F *KpCosThetaH_vs_t_mc[48];
  TH1F *TwoKpKmPimMass_1g,*PPimKmMass_2kp_1g,*PPimMass_2kpkm_1g,*LambdaKmMass_2kpkm_1g;
  TH1F *LambdaGamma_2kpkm_1g,*Sigma0KmMass_2kpkm_1g;
  TH1F *TwoKpKmPimMass,*PPimKmMass_2kp,*PPimMass_2kpkm,*LambdaKmMass_2kpkm;
  TH1F *KpKmPi0_2pi0,*KpKmPi0_etapi0,*KpKmPi0_2pi0_no_phi;
  TH2F *KmPipVsKmPip,*MCBeamCorrelation;
  TH1F *TrigBits,*UnusedEnergy,*TwoKpMissingMass,*TwoKpMissingMassCut;
  TH2F *EtaPi0VsKpKm,*TwoPi0VsKpKm;
  TH1F *KPlusdEdxCL,*KPlusBCALdtCL,*KPlusTOFdtCL,*KPlusPIDCL; 
  TH1F *ProtondEdxCL,*ProtonBCALdtCL,*ProtonTOFdtCL,*ProtonPIDCL;
  TH1F *PiMinusdEdxCL,*PiMinusBCALdtCL,*PiMinusTOFdtCL,*PiMinusPIDCL;
  TH2F *DPiMinusdEdx,*PiMinusdEdxSC,*PiMinusTOFdt; 
  TH1F *PiPlusdEdxCL,*PiPlusBCALdtCL,*PiPlusTOFdtCL,*PiPlusPIDCL;
  TH2F *DPiPlusdEdx,*PiPlusdEdxSC,*PiPlusTOFdt; 
  TH2F *DKPlusdEdx,*KPlusdEdxSC,*KPlusBCALdt;
  TH2F *DProtondEdx,*PiPlusBCALdt,*PiMinusBCALdt,*PiMinusdEdx,*KPlusTOFdt;
  TH2F *ProtondEdxAmp,*PiMinusdEdxFDC,*KPlusTOFdtCut,*KPlusBCALdtCut;
  TH2F *KPlusdEdxFDCCut,*KPlusdEdxFDC,*KMinusdEdxFDC,*PiPlusdEdx;
  TH2F *PiPlusdEdxFDC,*KPlusdEdxCut,*PiPlusdEdxCut;
  TH2F *DKMinusdEdx,*KMinusdEdxSC,*KMinusBCALdt,*KMinusTOFdt;
  TH1F *KMinusdEdxCL,*KMinusBCALdtCL,*KMinusTOFdtCL,*KMinusPIDCL; 
  TH2F *DProtondEdxSC;
  TH1F *ProtondEdxSCCL;
  TH1F *KMinusdEdxFDCCL,*KPlusdEdxFDCCL;
  TH2F *DKMinusdEdxFDC,*DKPlusdEdxFDC; 
  TH1F *ProtondEdxFDCCL;
  TH2F *DProtondEdxFDC;  
  TH1F *PiPlusdEdxFDCCL;
  TH2F *DPiPlusdEdxFDC;  
  TH1F *PiMinusdEdxFDCCL;
  TH2F *PKm_vs_t;
  TH2F *DPiMinusdEdxFDC,*PiPlusdEdxFDCcut;
  TH2F *PipPimMass_vs_PPimMass;
  TH1F *LambdaKmMass_2kpkm2pimpip,*LambdaKmMass_with_Ks;
  
  TH1F *F1pi0_KpKm_bg,*TwoPi0Mass_KpKm_bg,*EtaPi0Mass_KpKm_bg;
  TH2F *EtaPi0VsKpKm_bg;
  TH1F *KpKmMass_with_2pi0_bg,*KpKmMass_with_etapi0_bg,*KpKmPi0_2pi0_bg;
  TH1F *KpKmPi0_2pi0_no_phi_bg,*KpKmPi0_etapi0_bg,*KstarKstar_bg;
  TH2F *KpPi0_vs_KmPi0_bg,*PhiPi0Pi0Dalitz_bg;
    
  TH1F *PhiA0Mass,*PhiA0Mass_bg,*PhiA2Mass,*PhiA2Mass_bg;
  
  TH1F *PPim_k2pip3pim;
  TH2F *PipPim2D_k2pip3pim,*PipPim2D_k2pip3pim_LambdaCut;

  TH1F *PPimMass_with_K,*KPimMass_with_p;
  TH1F *PPimGammaMass_with_K,*PPimMass_with_Kg;
  TH1F *LambdaGammaMass,*PPimGammaMass_with_Kg,*PPimMass_with_K2g;
  TH1F *LambdaGammaMass_K2g,*TwoGammaMass_Kppim,*LambdaPi0Mass;
  TH1F *LambdaGammaMass_K2g_NU,*LambdaGammaMass_NU;
  TH1F *LambdaEtaMass,*Sigma0GammaMass,*SigmaPlusPiMinusMass;
  TH1F *LambdaEtaMass_NU,*Sigma0GammaMass_NU,*LambdaPi0Mass_NU;
  TH1F *N2gMass,*TwoGammaMass_NK,*NPi0Mass,*NEtaMass; 
  TH1F *TwoGammaMass_2kp2g,*PPimMass_2kp2g,*LambdaPimMass_2g;
  TH1F *LambdaPimPi0Mass,*TwoKpTwoPimMass,*P2PimMass_2kp;
  TH1F *PPimMass_2kp,*LambdaPimMass;
  TH2F *KpPim_vs_KpPim,*CLCorrelation2;
    
  TH1F *TwoGammaMass_kppim3g,*PPimMass_kppim3g,*LambdaGammaMass_kppim3g;
  TH1F *LambdaPi0Mass_kppim3g,*Sigma0Pi0Mass_kppim3g;
  TH1F *TwoGammaMass_kppim3g_NU,*PPimMass_kppim3g_NU,*LambdaGammaMass_kppim3g_NU;
  TH1F *LambdaPi0Mass_kppim3g_NU,*Sigma0Pi0Mass_kppim3g_NU;
  
  TH1F *PPim4GammaMass_kppim4g,*FourGammaMass_kppim4g;
  TH2F *FourGamma2d_kppim4g,*LambdaPi0Pi0Dalitz,*LambdaEtaPi0Dalitz;
  TH1F *PPimMass_kppim4g,*LambdaPi0Pi0Mass;
  TH1F *TwoGammaMassNoProton_kf;
  TH2F *LambdaGammaMassVsE,*LambdaGammaMassVst;
  TH2F *K_mom_vs_theta_1g,*LambdaGammaThetaVsMass,*Sigma0GammaThetaVsMass;
  TH1F *XiMinusPi0Mass,*PPi0_with_K;

  TH1F *TwoGammaMass_2kp2g_NU;
  TH1F *PPimMass_2kp2g_NU;
  TH1F *LambdaPimMass_2g_NU;
  TH1F *LambdaPimPi0Mass_NU;
  TH1F *XiMinusPi0Mass_NU;

  TH1F *TwoKpKmPimMass_1g_NU;
  TH1F *PPimKmMass_2kp_1g_NU;
  TH1F *PPimMass_2kpkm_1g_NU;
  TH1F *LambdaKmMass_2kpkm_1g_NU;
  TH1F *LambdaGamma_2kpkm_1g_NU;
  TH1F *Sigma0KmMass_2kpkm_1g_NU;
  
  TH2F *PPim_vs_Kpim_with_pip,*PPim_vs_pippim_with_Kp;
  TH2F *LambdaPip_vs_Kpim;

  TH1F *TwoGamma_k2pimpip;
  TH2F *PPim_vs_Kpim_with_pip2g,*PPim_vs_pippim_with_Kp2g;
  TH2F *LambdaPip_vs_Kpim_with_2g;
    
  TH1F *FourGamma_k2pimpip;
  TH2F *PPim_vs_Kpim_with_pip4g,*PPim_vs_pippim_with_Kp4g;
  TH2F *LambdaPip_vs_Kpim_with_4g;

  TH2F *LambdaPipPi0_vs_KPim,*LambdaPip_vs_KPimPi0;
  TH2F *LambdaPipPi0_vs_KPimPi0,*FourGamma2d_k2pimpip;
  

  // RooRealVar *RooMass,*RooE,*RooT,*RooWeight;
  //RooDataSet *RooData;
    double BETA_CUT,NUM_SIGMA_BG;
    bool VETO_GAP_EVENTS;
    bool EMULATE_TRIGGER,FILL_ROO_DATASET;
  double FCAL_POS_CUT;
  double BCAL_THRESHOLD,FCAL_THRESHOLD,FCAL_THETA_CUT;
  double PI0_VETO_CUT,TRACK_CL_CUT;
  double DEDX_CUT[3],FDC_DEDX_CUT[3];
  double CL_CUT,BG_CL_CUT,MIN_BEAM_E,MAX_BEAM_E;
  int DEBUG_LEVEL; 
  double PROTON_GAMMA_DT_CUT;
  bool DROP_ONE_CELL_BCAL_SHOWERS;
  double UNUSED_ENERGY_CUT,SPLIT_CUT,BCAL_Z_CUT,FCAL_RADIAL_CUT;
  double PROTON_CL_CUT,PIM_CL_CUT,PIP_CL_CUT,KP_CL_CUT,KM_CL_CUT,ANTIPROTON_CL_CUT;
  double PI0_CUT_VALUE,ETA_CUT_VALUE,ETAPRIME_CUT_VALUE;
  DVector3 mFCALCenter;

  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_strangeness_