#ifndef fiTQun_h_seen
#define fiTQun_h_seen

#include "fiTQun_shared.h"

#include "hemi/table.h"

#include <TFile.h>
#include <TH3F.h>
#include <TTree.h>

class fiTQun {
 private:

  static int fNPMT; 

  static fiTQun* staticthis;

  static int fParConst;
  static double t0bound[2];

  double Phi0_local[nPID]; // Get rid of these?
  double dEdx_local[nPID]; // Get rid of these?
  
  bool flgPMTdraw;
  bool fignrbkPMT;
  int ndraws;

  int nHitTDraws;
  
  int cring;
  
  int PIDring[nRngMX];//Particle type of each ring
  double mom_fortime[nRngMX];
  
  TRotation RingDirTfmMtx[nRngMX];
  
  static int flgTQcorr;
  
  bool flgRingDirTfmFit;
  
  bool flgdirgrdsd;
  
  int flgBndry;
 
  double QEEffCorr;
  double QEEffCorrVal;
  
  double WAttenLrecp[2];
  
  static double ftdrk;   // Set these as static. They are used at pre-fitting. 
  static double twnd[2]; // This way pre-fitter and time window definitions need to 
  static double mu_dark; // be called only once and not for every fiTQun derived fitter.
    
  double buff1Rpars[fiTQun_shared::nSnglTrkParams];//Single ring fit parameters
  
  double tOffset;
  
  int cpdfnpar[nPMT_max];
  double cpdfthr[nPMT_max][2];
  double cpdfcoeff[nPMT_max][6];//coefficients for gaussian upper/lower tail
  double cpdfpar[nPMT_max][6];//charge PDF parameters for each PMT
  
  double tcpdfpar[nRngMX][2*ntpdfmupar];//corrected time likelihood parameters
  
  bool UsePi0MassConstraint;

  bool mMaskForPrefit;


  static double fpunhit(double mu, double *par);//unhit probability
  static double fphit(double mu, double *par);//hit probability
  static double flogchrgPDF(double mu, int nparam, double *par, double *thr, double *coeff);//charge likelihood
  
  double Omega(double R);//solid angle factor of PMT
  double T(double R, double Zsrc, double Zpmt);//transmission function
//  double epsilon(double coseta);//PMT angular acceptance
  double J(double R, double coseta, double Zsrc, double Zpmt);
  
  void SetRingDirTfmMtx(int iring, double &theta, double &phi);
  void TfmRingDirToGlobal(int iring, double &theta, double &phi);
  
  static void OneRingWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void TwoRingWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void TwoColinearRingWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void OnepipWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void VtxPrefitWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void MultiRingWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void MRVtxEWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  
  static void FitQEEffWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  static void FitMuWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  
  void clamp(TH3F *h, float &x, float &y, float &z);

  virtual double ftsct(int icab, int iring);
  virtual double ftsct(int icab, int iring, bool LUT);
  virtual double ftdir(int icab, int iring);
  virtual double ftdir(int icab, int iring, bool LUT);

  double flogtimePDF(int icab);
  double totqpnlty(int nHittmp);
  void SetCPDFParams();
  double VtxEstmtr(double* X, int nHittmp=-1, int* icabtmp=NULL, double* tHittmp=NULL, double tMintmp=0., double tMaxtmp=-1., double* chrgtmp=NULL, bool mask = false);
  double VtxEstmtr(double* X, bool mask) {return VtxEstmtr(X,-1,NULL,NULL,0.,-1,NULL,mask);};
  
  double GetScatRatio(double* t2s, double zs, double rs, double phipos, int icab, double ct, double srcphi);
  
  int GetHitInfo();
  
  int GetnHitqTot(int& tmpnhit, double& tmptotq);
  
  int SnglTrk(int iPID, double *X, int iring, int flgUpst=0);
//  int UpstTrk(int iPID, int iring, double *pippar, bool fpipdown=false);

  int NSeg;
  int SegIndex[nRngMX];
  int MuMomentumCorr;

  
  //#########################################################
  
  static const int nBinQT_cos = 50;
  static const int nBinQT_Tres = 40;
  
  static const double Rang_Tres_coarse[2];// = {-50,350};
  static const double Rang_Tres_fine[2];// = {-15,65};
  
  float qObsTot;
  float qObsArr[2][nBinQT_Tres][nBinQT_cos];
  
  int nHitTot;
  int nHitArr[2][nBinQT_Tres][nBinQT_cos];
  
  float mudirTot;
  float mudirArr[nBinQT_cos];
  float musctTot;
  float musctArr[nBinQT_cos];
  
  float muArr[2][nBinQT_cos];
  float tPDFArr[2][nBinQT_Tres][nBinQT_cos];
  
  float DlnLArrQ[2][nBinQT_Tres][nBinQT_cos];
  float DlnLArrT[2][nBinQT_Tres][nBinQT_cos];
  
//  double TQdist[80][180];
  
  //#########################################################
  
 public:

  fiTQun(int anpmt = nPMT_max);
  ~fiTQun(){;}
  
  void InitEvent(int ievttmp);
  
  void SetTimeWindow(int itwndtmp, double tStart=0., double tEnd=-1.);
  
  void SetTQCorrflg(int tmp) {flgTQcorr=tmp;}
  
  void SetRingDirTfmFit(bool tmp) {flgRingDirTfmFit=tmp;}
  
  void SetDirsdflg(bool flg) {flgdirgrdsd=flg;}

  void SetQEEffCorr(double tmp) {QEEffCorrVal=tmp;}
  
  double GetElossmax(int iPID, double mom);
  
  double DirGridSrch(int iPID, double *X);

  bool SetMaskFlgForPrefit(bool mask){ mMaskForPrefit = mask;}
  void DoVtxPrefit(double *X);
  void DoPeakSearch(double *Xtmp, bool flgInWindow=false,bool mask = false, bool clear = true);
  void DoIngateFit(int flgInGateDcyFit, bool flgRePrefit, int npeaks, float *tpeaksf, double *X = 0);
//  void DoDecayESearch(int icluster_good, int ipeak, int npeak, float* tpeaksf, double * Xtmp = 0); // //the three paremeter will be updated. the information will contain previous peak (not only electron).
  void DoDecayESearch(); // //the three paremeter will be updated. the information will contain previous peak (not only electron).
  void MaskIngates();
  void MaskPMT(double *par=NULL, bool fMaskPeak=false, bool thresH = false);
  
  int IsPeakInCluster(double *X, double tStarttmp, double tEndtmp);
  void SetOneRingSeed(int iPID, double *X, bool fSeedVtx=false, bool fSeedDir=false, double *Xnoshft=NULL);
  double SetPi0Seed(double* fitpars, int& PCflg, bool seedpi0mass, bool translatevertex);
  double SetTwoColinearRingSeed(double* fitpars);
//  void SetPipSeed(double *pippar1, double *pippar2);
  void SetElossSeed(int iPID, double *pippar);
  
  double OneRingFit(int iPID, double *X, int& PCflg, int *flgparfix=NULL, double *arStps=NULL, int iring=-1);
  double TwoRingFit(int iPID1, int iPID2, bool usepi0mass, double *X, int& PCflg, int *flgparfix=NULL, double *arStps=NULL);
  
  double MultiRingFit(int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], int& PCflg, bool flgFixVtx=true, bool flgFixDir=true, int flgConvfit=0);
  
  double MRVtxEFit(int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], int& PCflg, bool flgFixTime=false, bool flgFixVtx=false, bool flgFixE=false);
  
  double Do1RFit(int iPID, double *X, int& PCflg, int flgSeed=1, double *X0=NULL, int ipeak=0);
  double DoPi0Fit(int iPID1, int iPID2, bool usepi0mass, int seedtype, double *X, int& PCflg, int defaultIndex, double *X0=NULL);
  double TwoColinearRingFit(int iPID1, int iPID2, double *X, int& PCflg, int *flgparfix=NULL, double *arStps=NULL);
  double DoUpst1RFit(int iPID, double *X, int& PCflg, double *X0=NULL, int ipeak=0);
  double OnepipFit(double *X, int& PCflg, int *flgparfix=NULL, double *arStps=NULL, bool fpipdown=false);
  
  void SortRingArr(int& nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams]);
  int SortEvis(int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double *RingEvis, int *idxSrtEvis);
  double CheckFakeRing(int& nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams]);
  
  void DoMoreMRFit(int flgMoreMRfit, int flgModeSeqFit=1, int flgMSfit=1, int *timeArr=NULL);
  void DoSeqMRFit(int idSeedMRFit, int flgModeSeqFit=1);
  void DoMRVtxEFit(int idSeedMRFit);
  void DoMRFit(int flgConvfit=0, int flgMRfit=1);
  
  int GetMRfitID(int nringtmp, int *TypeArr);
  int OneMoreRing(int nringtmp, int *TypeArr, int nTypeNew, int *TypeArrnew, double *lnLArrNew, int flgFitlevel, int flgConvfit=0);
  void SearchNewRing(int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], int nTypeNew, int *TypeArrnew, double ParArrnew[][fiTQun_shared::nSnglTrkParams], double *arlogL);
  
  int IsThereNewRing(double DlnL, double Enew, int nringnew=-1, bool flgMuLike=false);
  
  void Resetmu();
  void Resetmuring(int iring);
  void ChkNonNegmu(int iring);
  void DrawPreddist(char *chrHistnm=NULL, bool fSmpl=false);
  void DrawPMTs(char *cSufx=NULL);

  void DrawHitTDis(double* X, int nHittmp, int* icabtmp, double* tHittmp, double tMintmp, double tMaxtmp, double* chrgtmp, bool mask, char* hittname);
  void SetActivePMTs(int iflgMode=1);

  void TimeLEval(int niter);
  void TestTbl(int iPID, double R0, double costh0, double mom);
  double GetOneRngnglogL(int iPID, double *X, int& PCflg, int iring=-1);
  double GetTwoRngnglogL(int iPID1, int iPID2, double *X, int& PCflg);
  double GetTwoColinearRngnglogL(int iPID1, int iPID2, double *X, int& PCflg);
  double GetMRnglogL(int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams]);
  int Get1Rmudist(int iPID, double *X, double *muarr=NULL);
  int Get1Rmudist(int iPID, double *X, double *muarr, bool direct, bool indirect);
  double GetTotmu() {return fiTQun_shared::mutot;}
  
  void FitQEEff(int iPDG, double *X);
  double FitMu();
  
  void ReadSharedParams(int fScat=1, int *fRngtmp=NULL, bool fFitCProf=true, bool fFittProf=true, int iDetGeom=1) {
    fiTQun_shared::Get()->ReadParameters(fScat,fRngtmp,fFitCProf,fFittProf,iDetGeom);
    for (int iPID=0; iPID<nPID; iPID++) fiTQun_shared::fRngTyp[iPID]=fiTQun_shared::Get()->GetRingType(iPID);
  }
  
  int ClrInfo(int iMod) {return fiTQun_shared::Get()->ClrInfo(iMod);}
  
  void GetN50(double *par, int &N50, double &Q50);
  
  double GetDwall(double *pos);
  double GetToHKArc(double *pos, double *dir, int iquadrant);
  double GetToWall(double *pos, double *dir);
  
  void SetTreeQTdist(TTree *tmpTree);
  void GetQTdist_1Re(double d_shift=0., TH2D *hTQdist=NULL);
  void GetQTdist_1R(int iPID, double *X, TH2D *hTQdist=NULL);
  void Get1R2Rdist(int *ibincosArr, int ibinTresArr[][2], double *tauArr);
  void GetQTbin(double qObs, double Tres, int &qBin, int &TresBin);
  
  int GetEPipPID(double lnLpip, double lnLe, double momE) {
    if (lnLpip-lnLe>momE*fiTQun_shared::Get()->GetPIDCutEPip(1)+fiTQun_shared::Get()->GetPIDCutEPip(0)) return fiTQun_parameters::ie;// e
    else return fiTQun_parameters::ipip;// pi+
  }
  
  int GetMuPipPID(double lnLmu, double lnLpip, double Eloss) {
    if (lnLmu-lnLpip>Eloss*fiTQun_shared::Get()->GetPIDCutMuPip(1)+fiTQun_shared::Get()->GetPIDCutMuPip(0)) return fiTQun_parameters::ipip;//pi+
    else return fiTQun_parameters::imu;
  }
  
  // GetPhi0PhotMom returns the momentum of the second photon in a pi0 decay given
  // the direction and momentum of the first photon and the direction of the second photon
  double GetPi0PhotMom(double mom1, double dir1x, double dir1y, double dir1z, double dir2x, double dir2y, double dir2z);
  double GetPi0PhotMom(double mom1, double theta1, double phi1, double theta2, double phi2);
  double GetPi0PhotMom(double mom1, double cos_openingangle);
  double GetPhotAngle(double theta1, double phi1, double theta2, double phi2);
  double GetTruePi0Params(double* tPi0Params);


  //Multi-segment muon fit/////////////////////////////////////////////////////////////////////
  void SearchNewSegment(int    nringtmp,
                        int    iseg,
                        int*   segIndx,
                        int*   TypeArr,
                        int    segPID,
                        double setParArr[][fiTQun_shared::nSnglTrkParams],
                        double segEloss,
                        double arlogL);
  static void MultiSegWrapper(int& nDim, double* gout, double& result, double par[], int flg);
  double MultiSegFit(int nringtmp,
                     int nseg,
		     int* segIndx,
                     int* TypeArr,
                     int segPID, 
                     double ParArr[][fiTQun_shared::nSnglTrkParams],
                     int* flgFixArr,
                     int& PCflg);
  void DoMSFit(int flgFixVtx, int flg1R, double ParArrOrig[][fiTQun_shared::nSnglTrkParams]=NULL,
               int* TypeArrMR=NULL, int nringMR=1, int iMER=0, int MRfitID=1);
  double SegDirGridSrch(int iPID, double *X, double dSeg);
  //set internal parameters
  void SetNSeg(int nseg) {NSeg=nseg;}
  void SetMuMomCorr(int mumomcorr) {MuMomentumCorr=mumomcorr;} //toggle empirical momentum correction
  double GetSegLength(double pseg, double eloss);
  double GetSegTime(double pseg, double eloss);
//  void DoPG(double mom);
  ///////////////////////////////////////////////////////////////////////////////////////////////

 protected:

  fiTQun_shared * fqshared;
  int ievent,itwnd;
  double nglogL();//negative log likelihood
  int OneRing(int iPID, double *X, int iring);
  
};

#endif