// $Id$
//
//    File: JEventProcessor_EpEm.h
// Created: Tue Oct 21 16:50:20 EDT 2014
// Creator: staylor (on Linux gluon05.jlab.org 2.6.32-358.18.1.el6.x86_64 x86_64)
//

#ifndef _JEventProcessor_EpEm_
#define _JEventProcessor_EpEm_
//#include <RooDataSet.h>
//#include <RooRealVar.h>

#include <JANA/JEventProcessor.h>
#include "TRACKING/DTrackFinder.h"
#include <PID/DParticleID.h>
#include <TRACKING/DTrackTimeBased.h>
#include <PID/DChargedTrack.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <KINFITTER/DKinFitter.h>
#include <ANALYSIS/DKinFitUtils_GlueX.h>
#include <TRACKING/DMagneticFieldStepper.h>
#include <HDGEOMETRY/DGeometry.h>
#include <TRACKING/DReferenceTrajectory.h>
#include <ANALYSIS/DAnalysisUtilities.h>

class JEventProcessor_EpEm:public jana::JEventProcessor{
 public:
  JEventProcessor_EpEm();
  ~JEventProcessor_EpEm();
  const char* className(void){return "JEventProcessor_EpEm";}
  
 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 GetKF4vectors(DKinFitter *dKinFitter,DLorentzVector &beam_kf,
		     map<Particle_t,vector<DLorentzVector>>&final_kf) const;

  void MakeDoublets(vector<vector<pair<unsigned int, unsigned int > > >&doublets) const;
  void MakeTriplets(vector<vector<pair<unsigned int, unsigned int > > >&triplets) const;
  void FillNeutralParticleVectors(double t0_rf,const DVector3 &vertex,
				  vector<const DNeutralParticle*>&neutrals,
				  map<Particle_t,vector<const DNeutralParticleHypothesis*>>&neutralParticles) const;
  void FillChargedParticleVectors(vector<const DChargedTrack *>&tracks,map<Particle_t,vector<const DTrackTimeBased *>>&chargedParticles) const;
  void DoKinematicFit(double t0_rf,const DVector3 &vertex,
		      const DBeamPhoton *beamphoton,
		      map<Particle_t,vector<const DTrackTimeBased *> >&chargedParticles,
		      map<Particle_t,vector<const DNeutralParticleHypothesis*>>&neutralParticles,
		      DKinFitUtils_GlueX *dKinFitUtils,
		      DKinFitter *dKinFitter) const;
  
  void EplusEminusAnalysis(const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void EplusEminusGammaAnalysis(const DLorentzVector &beam_kf,
				map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void EplusEminus2GammaAnalysis(const DLorentzVector &beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight); 
  void EplusEminus3GammaAnalysis(const DLorentzVector &beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void EplusEminus4GammaAnalysis(const DLorentzVector &beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void EplusEminusPiplusPiminusAnalysis(const DLorentzVector &beam_kf,
					map<Particle_t,vector<DLorentzVector>>&final_kf,double weight); 
  void TwoEplusTwoEminusAnalysis(const DLorentzVector &beam_kf,
				 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void NEpEmPipAnalysis(const DLorentzVector &beam_kf,
			map<Particle_t,vector<DLorentzVector>>&final_kf,double weight);
  void KpPimEpEmAnalysis(const DLorentzVector &beam_kf,
			 map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPimEpEm1GAnalysis(const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  void KpPimEpEm2GAnalysis(const DLorentzVector &beam_kf,
			   map<Particle_t,vector<DLorentzVector>>&final_kf,double weight) const;
  
  const DParticleID* pid_algorithm;
  double GetProtonProbability(const DChargedTrackHypothesis *hyp) const;
  double GetPiPlusProbability(const DChargedTrackHypothesis *hyp) const; 
  double GetPiMinusProbability(const DChargedTrackHypothesis *hyp) const;
  double GetElectronProbability(const DChargedTrackHypothesis *hyp) const;
  double GetPositronProbability(const DChargedTrackHypothesis *hyp) const;
  static bool SortParticleProbability(const pair<double,Particle_t>a,
				     const pair<double,Particle_t>b){
    return a.first>b.first;
  }

  TH1F *MCMissingMass,*UnusedEnergy;
  TH2F *MCMissingMassVsEbeam,*PimRvsZ,*PipRvsZ,*ProtonRvsZ;

  TH1F *TwoEplusTwoEminusCL,*TwoEplusTwoEminusMass_kf;
  TH1F *EplusEminusMass_1g;
  TH2F *MCProtonThetaVsP,*ProtonThetaVsP;
  TH1F *MCProtonAcceptedP;

  TH1F *EplusEminusMissingGammaCL;
  TH1F *EplusEminusMass_with_2g,*EplusEminusMass_PipPim;
  TH1F *EplusEminusMass_with_pi0,*EplusEminusMass_with_eta;
  TH1F *EplusEminusMass_PipPim_eta_cut,*EplusEminusMass_PipPim_etaprime_cut;

  TH1F *EpEmDOCA,*EplusEminusMass_Rcut;
  TH2F *EpEmMassVsR; 

  TH1F *EpEmDOCA_n,*EplusEminusMass_npip_Rcut,*Nepem_pip_Rcut;
  TH1F *EpEmPip_n_Rcut,*Npip_epem_Rcut;
  TH2F *EpEmMassVsR_n; 

  TH1F *EpEmDOCA_1g,*EplusEminus1GammaMass_Rcut,*EplusEminusMass_1g_Rcut;
  TH2F *EpEmMassVsR_1g;

  TH1F *EpEmDOCA_2g,*EplusEminusMass_with_2g_Rcut,*TwoGammaMass_eplus_eminus_Rcut;
  TH1F *EplusEminusEtaMass_Rcut,*EplusEminusPi0Mass_Rcut;
  TH1F *EplusEminusMass_with_pi0_Rcut,*EplusEminusMass_with_eta_Rcut;
  TH2F *EpEmMassVsR_2g;

  TH1F *EpEmDOCA_3g,*EplusEminusMass_with_3g_Rcut,*EplusEminus3GammaMass_Rcut;

  TH1F *EpEmDOCA_4g,*EplusEminusMass_with_4g_Rcut,*EplusEminus4GammaMass_Rcut;
  TH1F *EplusEminusMass_with_4g,*EplusEminus4GammaMass_kf;
  TH1F *FourGammaMass_epem,*TwoPi0Mass_epem;
  TH2F *FourGamma2d_epem;

  TH1F *EpEmDOCA_pippim,*EplusEminusPiplusPiminusMass_Rcut;
  TH2F *EpEmMassVsR_pippim;

  TH2F *BeamPhotonTimeDiff,*PiMinusdEdxFDCCut;
  TH1F *EpEmMissingMassSq,*FCAL2ShowerDistance;
  TH2F *ProtondEdxFDC,*SCdIddT,*SCdIddTLowLevel;
  TH2F *SCdIddTLowLevelHasFADC;
   TH2F *ProtondEdx,*NormProtondEdx;
   TH1F *SCt;
  TH2F *dProtondEdx,*ProtondEdxCut,*ProtondEdxCutFDC;
  TH2F *ProbProtondEdx,*ProbKPlusdEdx,*ProbKMinusdEdx;

  TH2F *ProtonTOFdtNodEdx,*ProtonSCtCorr,*ProtonSCt;

  TH1F *TimeDiffTaggerTracks,*TDiffGammaProton,*NumGammaPerEvent,*NumGammaPerEventMM2Cut; 
  TH2F *ProtonTOFdt,*ProtonBCALdt,*ProtonBCALdtForward,*ProtonBCALdtWithdEdxCut;
  TH2F *ProtonBCALdtWithFDCdEdxCut,*ProtonFCALdt;
  TH2F *ProtonTOFdtWithdEdxCut,*ProtonFCALdtWithdEdxCut;

  TH2F *PiPlusdEdxFDC,*PiPlusTOFwithFDCdEdxCut,*PiPlusBCALTOFwithFDCdEdxCut;

  TH1F *TrigMask;

  TH2F *ProbProtonTOF,*NormProtonTOF;
  TH2F *PiPlusTOFwithdEdxCut,*ProtondEdxSC,*ProtondEdxSCNoCDCdEdx;

  TH2F *ProbPiMinusTOF,*NormPiMinusTOF;
  TH2F *ProbPiPlusTOF,*NormPiPlusTOF,*PiPlusBCALTOFwithdEdxCut;
  TH2F *ProtondEdxFDCWithSC;

  TH2F *NormPiPlusdEdx,*NormPiMinusdEdx,*PiPlusdEdxSC;
  TH2F *ProbPiPlusdEdx,*ProbPiMinusdEdx;

  TH2F *dPiPlusdEdx,*dPiMinusdEdx;

  TH2F *PiPlusdEdx,*PiMinusdEdx,*PiPlusdEdxCut,*PiMinusdEdxCut;
  TH2F *PiPlusdEdxFDCCut;

  TH2F *MCt_vs_E,*MCt_vs_Ecoarse;
  TH1F *MCBeam,*MCBeamTagged;

  TH1F *EplusEminusPiplusPiminusCL,*EplusEminusPiplusPiminusMass;
  TH1F *EplusEminusPiplusPiminusMass_cut;
  TH1F *TwoGammaMass_epem3g,* EplusEminus3GammaMass_kf,* EplusEminus3GammaCL;
  TH1F *EplusEminusMass_with_3g;
  TH1F *EplusEminusPi0GammaMass_kf;
  TH2F *DProtondEdxSC;
  TH1F *ProtondEdxSCCL;

  TH2F *EOverPMinus,*EOverPPlus,*EOverPMinus_BCAL,*EOverPPlus_BCAL;
  TH2F *EOverPMinus_cut,*EOverPPlus_cut,*EOverPMinus_BCAL_cut,*EOverPPlus_BCAL_cut;
  TH2F *Pim_EOverP_with_timecut,*Pip_EOverP_with_timecut;
  TH1F *TwoGammaMass_eplus_eminus,*EplusEminusPi0Mass_kf;
  TH1F *EplusEminusPi0Mass_cut,*EplusEminusEtaMass_cut;
  TH2F *DEOverPMinus_BCAL,*DEOverPPlus_BCAL;
  TH2F *DEOverPMinus,*DEOverPPlus;

  TH1F *Pi0Gamma_epem2g,*EtaGamma_epem2g;
  TH1F *Pi0Gamma_epem2g_cut,*EtaGamma_epem2g_cut;
  TH1F *EplusEminus1GammaMass_with_2g;
  TH1F *EpEmGammaMass;
  TH1F *ProtonGammaTimeDiff;
 TH1F *NumSCHits_for_3gam,*DTimeSC;

 TH2F *ReconPhiTheta;

 TH1F *TDiff2Gamma,*TDiff2GammaBeam,*EDiff2GammaBeam;
 
 double SPLIT_CUT;
 
  TH1F *MCProtonAcceptance,*MCProtonP,*BeamEnergy,*EtaCosTheta_phi_mc;
 
 TH1F *MCt;
 TH2F *MCdeltat;
 
 TH2F *EplusEminusMomenta;
 TH2F *EplusEminusNoHadrons_dEcut;
 TH1F *EplusEminusMass,*EplusEminusNoHadrons,*EplusEminusBeamDeltaE;
 TH1F *EplusEminusGamma,*TwoGamma_eplus_eminus,*TwoGammaEplusEminus;
 TH1F *EplusEminusGammaWithProton,*EplusEminusMassWithProton;
 TH1F *TwoGamma_eplus_eminus_with_proton,*TwoGammaEplusEminusWithProton;
 
 TH1F *EplusEminusGammaCL,*EplusEminus2GammaCL,*EplusEminusCL;
 TH1F *EplusEminus1GammaMass_kf,*EplusEminus2GammaMass_kf,*EplusEminusMass_kf;
 TH1F *EplusEminusEtaMass_kf;

 TH1F *Npip_epem,*EplusEminusMass_npip,*Nepem_pip,*EpEmPip_n;
 TH1F *TwoEplusTwoEminusMass_Rcut;
 TH2F *EpEmDOCA2,*EpEm_vs_EpEm,*EpEm_vs_EpEm_Rcut;

 TH2F *ElectronTOF,*PositronTOF;

 TH1F *Eta_sinphicosphi,*EtaPrime_sinphicosphi,*Pi0Gamma_with_epem;
 TH1F *OmegaPi0Mass_epem,*PiPlusPiMinusMass_epem;
 TH1F *NeutralBeta,*KinFitCL,*KinFitChiSq;

 TH1F *PPimMass_with_Kepem,*LambdaEpEm,*EpEmMass_with_Kppim;
 TH1F *PPimMass_with_Kepem1g,*LambdaEpEm_1g,*EpEmMass_with_Kppim1g;
 TH1F *EpEm1GMass_with_Kppim;
 TH1F * PPimMass_with_Kepem2g,*LambdaEpEm_2g,*EpEmMass_with_Kppim2g;
 TH1F *TwoGMass_with_Kppimepem;

 // RooRealVar *RooMass,*RooE,*RooT,*RooWeight;
 //RooDataSet *RooData;

 bool EMULATE_TRIGGER,FILL_ROO_DATASET;
 double FCAL_POS_CUT,BETA_CUT;
 double BCAL_THRESHOLD,FCAL_THRESHOLD;
 double PI0_VETO_CUT;
 int MIN_CLUSTERS_IN_INSERT;
 double BCAL_DTRACK_CUT;
 double DEDX_CUT[3],DEDX_CUT_FDC[3];
 double TRACK_CL_CUT,CL_CUT,MIN_BEAM_E,MAX_BEAM_E;
 bool INSERT_ONLY,FCAL_ONLY;
 int DEBUG_LEVEL;
 double UNUSED_ENERGY_CUT,BCAL_Z_CUT,FCAL_RADIAL_CUT;
 double PROTON_GAMMA_DT_CUT;
 bool DROP_ONE_CELL_BCAL_SHOWERS;
 double PI0_CUT_VALUE,ETA_CUT_VALUE,ETAPRIME_CUT_VALUE,NUM_SIGMA_BG;

 const DMagneticFieldMap *bfield; 
 DGeometry *dgeom;
 double PROTON_CL_CUT,PIM_CL_CUT,PIP_CL_CUT,ELECTRON_CL_CUT,POSITRON_CL_CUT;
 bool VETO_GAP_EVENTS;
 DVector3 mFCALCenter;

};

#endif // _JEventProcessor_EpEm_