// $Id$
//
//    File: JEventProcessor_hyperon.h
// Created: Wed Aug 30 09:03:49 EDT 2017
// Creator: staylor (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64)
//

#ifndef _JEventProcessor_hyperon_
#define _JEventProcessor_hyperon_

#include <JANA/JEventProcessor.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <RooDataSet.h>
#include <RooRealVar.h>

#include <TRIGGER/DL1MCTrigger.h>
#include <TRIGGER/DL1Trigger.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 <ANALYSIS/DAnalysisUtilities.h>
#include <KINFITTER/DKinFitter.h>
#include <ANALYSIS/DKinFitUtils_GlueX.h>
#include <TAGGER/DTAGMGeometry.h>
#include <TAGGER/DTAGHGeometry.h>

using namespace jana;
using namespace std;

class JEventProcessor_hyperon:public jana::JEventProcessor{
 public:
  JEventProcessor_hyperon();
  ~JEventProcessor_hyperon();
  const char* className(void){return "JEventProcessor_hyperon";}
  
 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.
  
  bool MakeMesonPairs(double weight,DLorentzVector pair1,DLorentzVector pair2,
		     map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>&particle_pairs,
		     map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>&bg_particle_pairs,
		     map<particle_pair_t,vector<double>>&bg_weights);
 void SetupMesonPairMaps(map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>&particle_pairs,
			 map<particle_pair_t,vector<pair<DLorentzVector,DLorentzVector>>>&bg_particle_pairs,
			 map<particle_pair_t,vector<double>>&bg_weights);

 
  void MakeDoublets(vector<vector<pair<unsigned int, unsigned int > > >&doublets) 
    const;

  void FillParticleVectors(vector<const DChargedTrack *>&tracks,
			   map<Particle_t, vector<const DTrackTimeBased *> >&particles
			   ) const;
  void FillNeutralVectors(double t0_rf,const DVector3 &vertex,
			  vector<const DNeutralParticle*>&neutrals,
			  vector<const DNeutralParticleHypothesis*>&gammaHyps,
			  vector<DLorentzVector>&neutron_v4s,
			  vector<TMatrixFSym>&neutron_covs,
			  vector<DLorentzVector>&neutron_x4s
			  ) const;
  void DoKinematicFit(const DBeamPhoton *beamphoton,
		      map<Particle_t,vector<const DTrackTimeBased *> >&particles,
		      vector<const DNeutralParticleHypothesis*>&gammaHyps,
		      vector<DLorentzVector>&neutron_v4s,
		      vector<TMatrixFSym>&neutron_covs,
		      DKinFitUtils_GlueX *dKinFitUtils,
		      DKinFitter *dKinFitter) const;
  void FillThrownHistos(const DMCReaction *reaction,
			vector<const DMCThrown*>&throwns,
			vector<const DTAGMGeometry *>&tagm_geom,
			vector<const DTAGHGeometry *>&tagh_geom) const;
  
  void KpPimAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPimEpEmAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPimEpEm1GAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPimEpEm2GAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPimGammaAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPim2GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPim3GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPim4GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void KpPim6GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void  TwoKpTwoPimAnalysis(DKinFitter *dKinFitter,double weight);
  void TwoKpTwoPim2GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void FillKpPim4GHistos(const DLorentzVector &pair1, const DLorentzVector &pair2,const DLorentzVector &kp_kf,const DLorentzVector &pim_kf, const DLorentzVector &proton_kf, double weight) const;
  double BackgroundFit(const DBeamPhoton *beam,
		       const DTrackTimeBased *proton,
		       const DTrackTimeBased *pip,
		       const DTrackTimeBased *pim);
  void KpPipPimNAnalysis(double t0_rf,DKinFitter *dKinFitter,
			 const DLorentzVector &neutron_x4,double weight);
  void KpPipPimN1GAnalysis(DKinFitter *dKinFitter,
			   const DLorentzVector &neutron_x4,double weight);
  void KpPipPimN2GAnalysis(DKinFitter *dKinFitter,
			   const DLorentzVector &neutron_x4,double weight);
  void TwoKpKmN2GAnalysis(DKinFitter *dKinFitter,
			   const DLorentzVector &neutron_x4,double weight);
  void KpN2GammaAnalysis(DKinFitter *dKinFitter,double weight);
  void Kp3Pim2PipAnalysis(DKinFitter *dKinFitter,DAnalysisUtilities *dAnalysisUtitilities,double weight);
  void Kp2PimPipAnalysis(DKinFitter *dKinFitter,DAnalysisUtilities *dAnalysisUtitilities,double weight);
  void Kp2PimPip2GammaAnalysis(DKinFitter *dKinFitter,DAnalysisUtilities *dAnalysisUtitilities,double weight);
  void Kp2PimPip4GammaAnalysis(DKinFitter *dKinFitter,DAnalysisUtilities *dAnalysisUtitilities,double weight);
  
  static bool SortParticleProbability(const pair<double,Particle_t>a,
				      const pair<double,Particle_t>b){
    return a.first>b.first;
  }
  
  void MakeDoublets();
  void MakeMesonPairs(vector<DLorentzVector>&gammas_kf,
		      double weight,TH2F *histo);

  TH2F *BeamPhotonTimeDiff,*EOverPPlus,*EOverPPlus_BCAL,*PositronTOF;
  TH2F *KpK0sqVsNK0sq,*LambdaPipPimDalitz_Rcut;
  TH1F *KpPipPimNCL;
  TH1F *PPim_k2pip3pim;
  TH2F *PipPim2D_k2pip3pim,*PipPim2D_k2pip3pim_LambdaCut;
  TH1F *LambdaPip_k2pip3pim_Rcut,*PPimDOCA_k2pip3pim;
  TH2F *PipPim2D_k2pip3pim_Rcut,*PPimMassVsR_k2pip3pim;

  TH1F *RFGammaTimeDiff,*PPi0_with_K,*SigPlusPiMinusMass;
  TH1F *RFGammaTimeDiffNoNeutron;
  TH2F *NeutrondTvsZ,*NeutrondXvsZ,*NeutrondYvsZ;

  TH1F *MissingMassOffKaon,*MissingPChiSq,*MissingMassOffKaonCut;
  TH1F *MissingMassOffKPi,*MissingMassOffKaonKF;
  TH2F *MissingProton,*MissingProtonRecon,*MissingProtonReconCut;
  TH1F *ReconMissingPDiff,*ReconMissingThetaDiff,*ReconMissingPhiDiff;
  TH2F *LambdaGammaThetaVsMass,*Sigma0GammaThetaVsMass;

  TH1F *PPimMass_with_Kepem2g,*EpEmMass_with_Kppim2g;
  TH1F *LambdaEpEm_2g,*TwoGMass_with_Kppimepem;
  TH1F *PPimMass_with_Kepem1g,*EpEmMass_with_Kppim1g;
  TH1F *LambdaEpEm_1g,*EpEm1GMass_with_Kppim;
  TH1F *PPimMass_with_K,*PPimMass_with_Kepem,*EpEmMass_with_Kppim;
  TH1F *LambdaEpEm;
  TH1F *PPimGammaMass_with_K,*PPimMass_with_Kg;
  TH1F *LambdaGammaMass,*PPimGammaMass_with_Kg,*PPimMass_with_K2g;
  TH1F *LambdaGammaMass_K2g,*TwoGammaMass_Kppim,*LambdaPi0Mass;
  TH1F *LambdaEtaMass,*Sigma0GammaMass,*SigmaPlusPiMinusMass;
  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 *PPim4GammaMass_kppim4g,*FourGammaMass_kppim4g;
    TH2F *FourGamma2d_kppim4g,*LambdaPi0Pi0Dalitz,*LambdaEtaPi0Dalitz;
    TH1F *PPimMass_kppim4g,*LambdaPi0Pi0Mass,*KpPipPimMissingMass;
    TH2F *KpBeamPhotonTimeDiff,*NPipSqVsNPimSq,*NPipSqVsNPimSq_kf;
    TH1F *NeutralBeta,*KpPipPimNMissingMassSq,*KpPipPimNMissingMassSqCut;
    TH1F *KinFitChiSq,*KinFitCL,*NPipPimMass;
    TH2F *NPip_vs_Kpim,*NKp_vs_pippim,*ProtonTOFdtWithFDCdEdx;
    TH1F *SigmaMinusPipMass,*SigmaPlusPimMass;
    TH1F *TwoGammaMassNoProton_kf,*UnusedEnergy;
    TH2F *LambdaGammaMassVsE,*LambdaGammaMassVst;
    TH2F *K_mom_vs_theta_1g;
    TH1F *XiMinusPi0Mass;
    TH1F *MCBeam,*MCt;
    
    TH2F *NPipSqVsNPimSq_1G,*NPipSqVsNPimSq_2G;
    TH1F *TwoGamma_with_KNPipPim;
    TH2F *SigmaPlusPi0_vs_KPim,*SigmaPlusPim_vs_KPi0;
    TH2F *SigmaMinusPi0_vs_KPip,*SigmaMinusPip_vs_KPi0;
    TH2F *SigmaPlusPi0_SigmaPlusPim_Dalitz,*SigmaMinusPi0_SigmaMinusPip_Dalitz; 
    TH1F *SigmamPipPi0,*SigmapPimPi0;
    TH2F *SigmaPlusEta_vs_KPim,*SigmaPlusPim_vs_KEta;
    TH2F *SigmaMinusEta_vs_KPip,*SigmaMinusPip_vs_KEta;
    TH2F *SigmaPlusEta_SigmaPlusPim_Dalitz,*SigmaMinusEta_SigmaMinusPip_Dalitz; 
    TH1F *SigmamPipEta,*SigmapPimEta;
	
    TH1F *TwoGamma_2kpkmn,*NPi0_2kpkm;
    TH1F *SigmaPlusGamma,*SigmaMinusGamma;

    TH1F *PPimDOCA_k2pimpip;
    TH2F *PPim_vs_Kpim_with_pip,*PPim_vs_pippim_with_Kp;
    TH2F *PPimMassVsR_k2pimpip,*LambdaPip_vs_Kpim;

    TH1F *PPimDOCA_k2pimpip2g,*TwoGamma_k2pimpip;
    TH2F *PPim_vs_Kpim_with_pip2g,*PPim_vs_pippim_with_Kp2g;
    TH2F *PPimMassVsR_k2pimpip2g,*LambdaPip_vs_Kpim_with_2g;
    
    TH1F *PPimDOCA_k2pimpip4g,*FourGamma_k2pimpip;
    TH2F *PPim_vs_Kpim_with_pip4g,*PPim_vs_pippim_with_Kp4g;
    TH2F *PPimMassVsR_k2pimpip4g,*LambdaPip_vs_Kpim_with_4g;

    TH2F *LambdaPipPi0_vs_KPim,*LambdaPip_vs_KPimPi0;
    TH2F *LambdaPipPi0_vs_KPimPi0,*FourGamma2d_k2pimpip;
    
    TH2F *LambdaPi0Pi0Dalitz_bg,*LambdaEtaPi0Dalitz_bg,*FourGamma2d_bg;
    TH1F *LambdaPi0Pi0Mass_bg;
    TH2F *LambdaPipPi0_vs_KPimPi0_bg;
    
    RooRealVar *RooMass;
    RooDataSet *RooData,*RooAccData;

    const DParticleID* pid_algorithm;
    DVector3 mFCALCenter;

    bool EMULATE_TRIGGER;
    double FCAL_POS_CUT,BETA_CUT,NUM_SIGMA_BG;
    double BCAL_THRESHOLD,FCAL_THRESHOLD,FCAL_THETA_CUT;
    double PI0_VETO_CUT,TRACK_CL_CUT;
    double FCAL_RADIAL_CUT,BCAL_Z_CUT;    
    double CL_CUT,BG_CL_CUT,MIN_BEAM_E,MAX_BEAM_E;
    int DEBUG_LEVEL;
    double KAON_GAMMA_DT_CUT;
    bool DROP_ONE_CELL_BCAL_SHOWERS;
    double UNUSED_ENERGY_CUT,SPLIT_CUT;
    double PROTON_CL_CUT,PIM_CL_CUT,PIP_CL_CUT,KP_CL_CUT;
    double PI0_CUT_VALUE,ETA_CUT_VALUE,ETAPRIME_CUT_VALUE;

    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;
};

#endif // _JEventProcessor_hyperon_