#include "fiTQun_shared.h"
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <stdio.h>
#include <time.h>

#include <TF1.h>

#include <TH1D.h>
#include <TH2D.h>
#include <TH3F.h>
#include <TH3D.h>
#include <TVector3.h>

#include "fQChrgPDF.h"
#include "TRuntimeParameters.hxx"
#include "fQEventManager.h"

// Put compiler flag around sk info
#ifndef NOSKLIBRARIES
#include "skheadC.h"
#include "skparmC.h"
#include "geopmtC.h"
#include "sktqC.h"
#include "skbadcC.h"
#else
#include "WCSimWrap.h"
#endif
/// Some of variables from following are 
/// used in the code, so keep it.
#include "spliTChanOutC.h"
#include "fitqunoutC.h"

using namespace fiTQun_parameters;

extern"C" {
#ifndef NOSKLIBRARIES
  void kzbloc_(char *,int *,int);
  
  void fortconsts_();
  
  extern struct{
    int isk23[MAXPM];
    float qefactor[MAXPM];
  } pmtinfcmn_;
#endif
  
  void clrfqcmns_(int *);
}

fiTQun_shared* fiTQun_shared::staticthis = NULL;
int fiTQun_shared::fMaxNPMT = 11146;

const double pi=3.1415926535898;
const double costhmax=1.-1e-12;

const int fiTQun_shared::PIDarr[nPID]={22,11,13,211,321,2212,48};//Supported particle types(PDG code)
const TString fiTQun_shared::PIDnames[nPID]={"g","e","mu","pip","kp","p","CG"};//Short Particle Names
const double fiTQun_shared::m0[nPID]={0.,0.5110,105.6584,139.5702,493.677,938.2720,0.};//Particle rest mass
double fiTQun_shared::dEdx[nPID]={1.,1.,2.14786,2.13069,2.06859,2.04043,1.};//dE/dx for upstream track calculation

const int fiTQun_shared::nTrkPar[nTrkTyp]={8,7,7,7,7,7,7,12,11};//Number of fit parameters for each track type
const int fiTQun_shared::npar[3]={4,6,4};//number of charge pdf parameters(Lo,Mid,Hi)
double fiTQun_shared::Phi0[nPID]={0.983,0.973,0.898,0.921,0.921,0.921,1.004};

double fiTQun_shared::expbinwrecp;
double fiTQun_shared::expbins[nexpbinm1+1];
double fiTQun_shared::expTbl[nexpbinm1+1];

double fiTQun_shared::epsilon_binwrecp;
double fiTQun_shared::epsilon_bins[nEpsilonBins+1];
double fiTQun_shared::epsilon_Tbl [nEpsilonBins+1];

int fiTQun_shared::flgHit[nPMT_max];
int fiTQun_shared::flgMask[nPMT_max];
int fiTQun_shared::ihitPMT[nPMT_max];
int fiTQun_shared::nhitPMT;
double fiTQun_shared::tsigmvtx;

double fiTQun_shared::chrg[nPMT_max];
double fiTQun_shared::mutot;
double fiTQun_shared::mu[nRngMX][nPMT_max];
double fiTQun_shared::muscat[nRngMX][nPMT_max];
double fiTQun_shared::tHit[nPMT_max];
double fiTQun_shared::tau[nRngMX][nPMT_max];
double fiTQun_shared::nglnLarr[nPMT_max];
double fiTQun_shared::nglnLtarr[nPMT_max];
double fiTQun_shared::costheta0[nPMT_max];
double fiTQun_shared::qtot_fwd;
double fiTQun_shared::mutot_fwd;

int fiTQun_shared::fRngTyp
[nPID];
int fiTQun_shared::nring;
int fiTQun_shared::PIDring[nRngMX];//Particle type of each ring

double fiTQun_shared::inDirtPDFpars[3];

// LUT time PDFs
hemi::Table3D<float> * fiTQun_shared::tpdf_direct[nPID];
hemi::Table3D<float> * fiTQun_shared::tpdf_indirect[nPID];

int    fiTQun_shared::npfitr;       // Number of pre-fit iterations. Maximum is 10.
double fiTQun_shared::sgmar[10]; // Width (2*sigma^2) of the gaussian for the vertex estimator
double fiTQun_shared::grdszVtx;  // Spatial step size for vertex grid search
double fiTQun_shared::grdsztime; // Temporal step size for vertex grid search

int    fiTQun_shared::nRingMaxTrials; // Number of rings to fit with maximum number of ring hypothesis (further rings fit only with one ring hypothesis -- e-like)

fiTQun_shared::fiTQun_shared( int aMaxNPmt ){
  
  flgLoaded=false;
  
  fMaxNPMT = aMaxNPmt;
  
  SK_Ver=1;

  WCSim_Config = "NULL";
  WCSim_PMTType = "NULL";
  
  // Set default indirect time likelihood parameters to SK values
  // gamma:
  inDirtPDFpars[0] = 25;
  // sigma:
  inDirtPDFpars[1] = 8;
  // t0:
  inDirtPDFpars[2] = 25;
  
  // Default values for vertex pre-fit:
  npfitr = 4;
  sgmar[0] = 500.;
  sgmar[1] = 200.;
  sgmar[2] = 80.;
  sgmar[3] = 32.;
  grdszVtx = 400.;
  grdsztime = 15.;

  // Default value for nRingMaxTrials
  nRingMaxTrials = 4;

  WAttenL[0]=7200.;
  WAttenL[1]=7900.;
  
  const double xmax=-10.;
  for (int i=0; i<=nexpbinm1; i++) {
    expbins[i]=i*xmax/nexpbinm1;
    expTbl[i]=exp(expbins[i]);
  }
  expbinwrecp=1./(expbins[1]-expbins[0]);
  
  const int fRngTyp_Def[nPID]={0,1,1,2,0,2,0};//Default ring type for each particle type (1:normal, 2;upstream track)
  for (int iPID=0; iPID<nPID; iPID++) {
    fRngTyp[iPID]=fRngTyp_Def[iPID];
  }
  
  fQuiet = 1;
  MSFitMethod=0;

  // Read shared parameters on construction
#ifdef NOSKLIBRARIES
  TString strDetName = "WCSim";
#else
  SK_Ver = skheadg_.sk_geometry;
  TString strDetName = Form("SK%d",SK_Ver);
#endif

  fQuiet = TRuntimeParameters::Get().GetParameterI("fiTQun.fQuiet");

  MSFitMethod =  TRuntimeParameters::Get().GetParameterI("fiTQun.MSFitMethod");
  std::cout<< "Multi-segment fit method: "<<MSFitMethod<<std::endl;

  MSElossMin = TRuntimeParameters::Get().GetParameterD("fiTQun.MSElossMin");

  MSThetRes = TRuntimeParameters::Get().GetParameterD("fiTQun.MSThetRes");

  twnd_def[0] = TRuntimeParameters::Get().GetParameterD("fiTQun.DefaultTimeWindowStart");
  twnd_def[1] = TRuntimeParameters::Get().GetParameterD("fiTQun.DefaultTimeWindowEnd");

  flgscat = TRuntimeParameters::Get().GetParameterI("fiTQun.UseScatteredLight");

  PrefitVtxMinDwall = TRuntimeParameters::Get().GetParameterD("fiTQun.PrefitVtxMinDwall");

  PkFlatThr = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.PeakThr%s",strDetName.Data()));
  std::cout << "Peak finder threshold=" << PkFlatThr << std::endl;

  chrgfact = 1;
	
  thrTres[0] = TRuntimeParameters::Get().GetParameterD("fiTQun.ThrResTimeL");
  thrTres[1] = TRuntimeParameters::Get().GetParameterD("fiTQun.ThrResTimeR");
  std::cout << "ThrTres=[" << thrTres[0] << "," << thrTres[1] << "]" << std::endl;

  flgMaskDeadTimePMT = (bool) TRuntimeParameters::Get().GetParameterI("fiTQun.MaskDeadTimePMT");

  flgTVarCorr = TRuntimeParameters::Get().GetParameterI("fiTQun.CorrDetTVar");

  GainCorrFact = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.GainCorrFact%s",strDetName.Data()));
  std::cout << "before this" << std:: endl;
  //  std::cout << "MaxAng[deg]=" << TotqMaxAng << std::endl;

  ctotqpnl = TRuntimeParameters::Get().GetParameterD("fiTQun.TotalChargeConstraintWeight");
  std::cout << "Total charge constraint: " << ctotqpnl << std::endl;
  double TotqMaxAng = TRuntimeParameters::Get().GetParameterD("fiTQun.TotalChargeConstraintAngle");
  std::cout << "MaxAng[deg]=" << TotqMaxAng << std::endl;
  this->SetTotqMaxAng(TotqMaxAng);
    
  flgBiasCorr = TRuntimeParameters::Get().GetParameterI("fiTQun.DoBiasCorrection");

  flgtime = TRuntimeParameters::Get().GetParameterI("fiTQun.UseTimeLikelihood");

  // Read in indirect time pdf parameters
  inDirtPDFpars[0] = TRuntimeParameters::Get().GetParameterD("fiTQun.inDirtPDFGamma");
  inDirtPDFpars[1] = TRuntimeParameters::Get().GetParameterD("fiTQun.inDirtPDFSigma");
  inDirtPDFpars[2] = TRuntimeParameters::Get().GetParameterD("fiTQun.inDirtPDFt0");
  std::cout << "Indirect time PDF parameters: gamma " << fiTQun_shared::inDirtPDFpars[0] << "; sigma " << fiTQun_shared::inDirtPDFpars[1] << "; t0 " << fiTQun_shared::inDirtPDFpars[2] << std::endl;
  
  // what PDF mode should we use for the time PDF? TRUE is LUT, false is fitted. FIXME: generalise
  flgTimePDFmode = TRuntimeParameters::Get().GetParameterI("fiTQun.UseLookupTableTimePDF");
  
  this->SetWAttL(TRuntimeParameters::Get().GetParameterD(Form("fiTQun.WaterAttenuationLength%s",strDetName.Data())));

  nwtr = TRuntimeParameters::Get().GetParameterD("fiTQun.WaterRefractiveIndex");

  darkrate = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.DarkRate%s",strDetName.Data()));
  
  flgdraw = TRuntimeParameters::Get().GetParameterI("fiTQun.DrawFitStages");
  std::cout << "Draw = " << flgdraw << std::endl;

  int flgDeactPMT = TRuntimeParameters::Get().GetParameterI("fiTQun.PhotoCoverageFracDenom");
  flgDeactPMT += 10*TRuntimeParameters::Get().GetParameterI("fiTQun.DeactCapPMTs");
  if (flgDeactPMT>1) {
    this->SetActivePMTs(flgDeactPMT);
  }

  std::string outZBSFileName = TRuntimeParameters::Get().GetOutZBSFileName();
  char* coutZBSFileName = (char*)outZBSFileName.c_str();
  int outZBSFileLength;
  for (int i=0; i<256; i++) {
    if (coutZBSFileName[i]=='\0') {
      outZBSFileLength=i;
      break;
    }
  }

  this->Setofilenm(coutZBSFileName);

  QEEff = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.QEEff%s",strDetName.Data()));

  this->SetPhi0(ig,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrGamma1RFit"));
  this->SetPhi0(ie,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrElectron1RFit"));
  this->SetPhi0(imu,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrMuon1RFit"));
  this->SetPhi0(ipip,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrPionPlus1RFit"));
  this->SetPhi0(ikp,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrKaonPlus1RFit"));
  this->SetPhi0(ip,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrProton1RFit"));
  this->SetPhi0(icg,TRuntimeParameters::Get().GetParameterD("fiTQun.QECorrConeGenerator1RFit"));

  this->SetdEdx(ig,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxGamma"));
  this->SetdEdx(ie,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxElectron"));
  this->SetdEdx(imu,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxMuon"));
  this->SetdEdx(ipip,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxPionPlus"));
  this->SetdEdx(ikp,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxKaonPlus"));
  this->SetdEdx(ip,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxProton"));
  this->SetdEdx(icg,TRuntimeParameters::Get().GetParameterD("fiTQun.dEdxConeGenerator"));


  QmuConvFact = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.QmuConvFact%s",strDetName.Data()));

  PIDCut_epip[0] = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.PIDCutEPipa0%s",strDetName.Data()));
  PIDCut_epip[1] = TRuntimeParameters::Get().GetParameterD(Form("fiTQun.PIDCutEPipa1%s",strDetName.Data()));
  
  PIDCut_mupip[0] =TRuntimeParameters::Get().GetParameterD(Form("fiTQun.PIDCutMuPip0%s",strDetName.Data()));
  PIDCut_mupip[1] =TRuntimeParameters::Get().GetParameterD(Form("fiTQun.PIDCutMuPip1%s",strDetName.Data()));

  nPi0SeedScanPoints = TRuntimeParameters::Get().GetParameterI("fiTQun.nPi0SeedScanPoints");
  pi0PreFitVtxShiftThr = TRuntimeParameters::Get().GetParameterD("fiTQun.pi0PreFitVtxShiftThr");
  pi0PreFitVtxShiftStepSize = TRuntimeParameters::Get().GetParameterD("fiTQun.pi0PreFitVtxShiftStepSize");
  pi0PreFitVtxShiftNsteps = TRuntimeParameters::Get().GetParameterI("fiTQun.pi0PreFitVtxShiftNsteps");

  RCCut_a[0]=TRuntimeParameters::Get().GetParameterD(Form("fiTQun.RCCuta0%s",strDetName.Data()));
  RCCut_a[1]=TRuntimeParameters::Get().GetParameterD(Form("fiTQun.RCCuta1%s",strDetName.Data()));
  RCCut_1Ea0=TRuntimeParameters::Get().GetParameterD(Form("fiTQun.RCCut1Ea0%s",strDetName.Data()));
  RCCut_1Mua0=TRuntimeParameters::Get().GetParameterD(Form("fiTQun.RCCut1Mua0%s",strDetName.Data()));
  
  nSegMax = TRuntimeParameters::Get().GetParameterI("fiTQun.nSegMax");

  MSScattSig = TRuntimeParameters::Get().GetParameterD("fiTQun.MSScattSig");

  // Number of rings where maximum number of hypothesis are tried. Further rings are fitted with e-like assumption only. Default is 4.
  if( TRuntimeParameters::Get().HasParameter("fiTQun.nRingMaxTrials") )
    nRingMaxTrials = TRuntimeParameters::Get().GetParameterI("fiTQun.nRingMaxTrials");
  
}

void fiTQun_shared::SetTimePDFMode (bool flg){
  flgTimePDFmode = flg;
  if (flgTimePDFmode) {
    // if we are using the LUT mode, then we need to load the tables
    fiTQun_shared::Get()->LoadLookupTableTimePDF();
  }
}



void fiTQun_shared::ReadParameters(int fScat, int *fRngtmp, bool fFitCProf, bool fFittProf, int iDetGeom){
  
  // All fiTQun constant root files must be placed in fQconstdir
  
  char *ctmp = getenv("FITQUN_ROOT");
  if (ctmp!=NULL) {
    fQconstdir=ctmp;
    fQconstdir+="/";
  }
  else {
    ctmp = getenv("ATMPD_ROOT");
    if (ctmp!=NULL) {
      fQconstdir=ctmp;
      fQconstdir+="/src/recon/fitqun/";
    }
    else {
      std::cout << "ATMPD_ROOT not set!" << std::endl;
      exit(-1);
    }
  }
  
  fQconstdir+="const/";
  
  if (fRngtmp!=NULL) {
    for (int iPID=0; iPID<nPID; iPID++) {
      fRngTyp[iPID]=fRngtmp[iPID];
    }
  }
  
  std::cout << "fiTQun_shared::ReadParameters ..." << std::endl;
  for (int iPID=0; iPID<nPID; iPID++) {
    if (fRngTyp[iPID]) {
      std::cout << " Enabling " << PIDnames[iPID] << " ring hypothesis";
      if (fRngTyp[iPID]==2) std::cout << " - upstream track";
      std::cout << std::endl;
    }
  }
  
#ifndef NOSKLIBRARIES
  SK_Ver = skheadg_.sk_geometry;
#endif
  
  std::cout << "SK version: " << SK_Ver << std::endl;
  
  strSKVer=Form("_sk%d",SK_Ver);

  TFile *TestFile = new TFile(Form("%s/%s_sk%d.root",fQconstdir.Data(),"cPDFpar",SK_Ver));
  if (TestFile->IsZombie()) {// new constants exist!
    std::cout << "Constants not found in " << fQconstdir << std::endl;
    exit(-1);
  }
  delete TestFile;
  
  std::cout << "Reading fiTQun constants from: " << fQconstdir << std::endl;
  
  TFile* tf;
  
  TH2D *hscattbl2d;
  TH3D *hscattbl3d;
  SctTbl2d[0]=NULL;
  SctTbl3d=NULL;
  
  if (fScat==3) {
    std::cout << "Scattering table mode 3 is no longer supported!" << std::endl;
    exit(-1);
  }
  
  
  // If we're using WCSim AND the WCSim detector configuration has been specified
  // in the parameters file load the appropriate PMT charge response. Otherwise
  // stick to default behaviour which is to determine the SK version automatically
  // from SKDETSIM input or to use SK1 for WCSim.
  TString ScatTable_Config = strSKVer;
  
#ifdef NOSKLIBRARIES
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ){
    WCSim_Config = TRuntimeParameters::Get().GetParameterS("fiTQun.WCSimConfig");
    ScatTable_Config = Form("_%s", WCSim_Config.Data());
  }
#endif
  
  if (fScat==1 || fScat==4) {//scattered light enabled
    
    if (fRngTyp[icg]) {//CG Reconstruction
      tf = new TFile(Form("%s/%s%s.root",fQconstdir.Data(),"fiTQun_scattablesF_CG",ScatTable_Config.Data()));//CG scattering table
      std::cout << "Opening CG scattering table file..." << std::endl;
    }
    else {
      tf = new TFile(Form("%s/%s%s.root",fQconstdir.Data(),"fiTQun_scattablesF",ScatTable_Config.Data()));
      std::cout << "Opening scattering table file..." << std::endl;
    }
    if (tf->IsZombie()) {
      std::cout << "Scattering table file does not exist!" << std::endl;
      DeleteInstances();
      exit(-1);
    }
    
    tf->GetObject("topscattable",top_scattable);
    tf->GetObject("botscattable",bot_scattable);
    tf->GetObject("sidescattable",side_scattable);
    
    tf->Close();
    delete tf;
    
    top_scattable->PrintMinElement();
    bot_scattable->PrintMinElement();
    side_scattable->PrintMinElement();

    top_scattable->fillbininfo();
    bot_scattable->fillbininfo();
    side_scattable->fillbininfo();
    
  }
  if (fScat==2 || fScat==3 || fScat==4) {//3d scattering table
    
    std::cout << "Opening 3D scattering table file..." << std::endl;
    tf = new TFile(Form("%s/%s%s.root",fQconstdir.Data(),"fiTQun_scattable3d",ScatTable_Config.Data()));
    if (tf->IsZombie()) {
      std::cout << "3D scattering table file does not exist!" << std::endl;
      DeleteInstances();
      exit(-1);
    }
    
    tf->GetObject("hscattable3D",hscattbl3d);
    
    nsct3dRbinm1=hscattbl3d->GetNbinsX()-1;
    nsct3dwlbinm1=hscattbl3d->GetNbinsY()-1;
    nsct3dthbinm1=hscattbl3d->GetNbinsZ()-1;
    
    sct3dRbins = new double[nsct3dRbinm1+1];
    sct3dwlbins = new double[nsct3dwlbinm1+1];
    sct3dthbins = new double[nsct3dthbinm1+1];
    
    SctTbl3d = new float**[nsct3dRbinm1+1];
    for (int iR=0; iR<=nsct3dRbinm1; iR++) {
      sct3dRbins[iR]=hscattbl3d->GetXaxis()->GetBinCenter(iR+1);
      SctTbl3d[iR] = new float*[nsct3dwlbinm1+1];
      for (int iwl=0; iwl<=nsct3dwlbinm1; iwl++) {
        sct3dwlbins[iwl]=hscattbl3d->GetYaxis()->GetBinCenter(iwl+1);
        SctTbl3d[iR][iwl] = new float[nsct3dthbinm1+1];
        for (int ith=0; ith<=nsct3dthbinm1; ith++) {
          sct3dthbins[ith]=hscattbl3d->GetZaxis()->GetBinCenter(ith+1);
          SctTbl3d[iR][iwl][ith]=hscattbl3d->GetBinContent(iR+1,iwl+1,ith+1);
        }
      }
    }
    sct3dRbinw=sct3dRbins[1]-sct3dRbins[0];
    sct3dwlbinw=sct3dwlbins[1]-sct3dwlbins[0];
    sct3dthbinw=sct3dthbins[1]-sct3dthbins[0];
    delete hscattbl3d;
    tf->Close();
    delete tf;
    std::cout << "3D scattering table successfully loaded!" << std::endl;
    
    if (fScat!=3) {
      nsct2dRbin=nsct3dRbinm1+1;
      nsct2dthbin=nsct3dthbinm1+1;    
      sct2dRbins = new double[nsct2dRbin];
      sct2dthbins = new double[nsct2dthbin];
      for (int iR=0; iR<nsct2dRbin; iR++) sct2dRbins[iR]=sct3dRbins[iR];
      for (int ith=0; ith<nsct2dthbin; ith++) sct2dthbins[ith]=sct3dthbins[ith];
      
      for (int istp=0; istp<3; istp++) {
        SctTbl2d[istp] = new float[nsct2dRbin*nsct2dthbin];
      }
      
      sct2dRbinwrecp=1./(sct2dRbins[1]-sct2dRbins[0]);
      sct2dthbinwrecp=1./(sct2dthbins[1]-sct2dthbins[0]);
    }
    
  }
  if (!(fScat>0 && fScat<=4)) {
    std::cout << "fiTQun_shared: Scattered light disabled!" << std::endl;
  }

  TransDetToGlbl.SetXYZ(0.,0.,0.);// Detector centre in global coordinate system
  
  TVector3 DetZinGlbl(0.,0.,1.);// Detector local Z-axis in global coordinate system
  TVector3 DetXinGlbl(1.,0.,0.);// Detector local X-axis in global coordinate system

  /// put this into compiler flag
#ifndef NOSKLIBRARIES
  
  fortconsts_();//fill commom blocks for SK constants
  
  PMTRadius = 25.4;
  for (int icab=0; icab<fMaxNPMT; icab++) {
    for (int j=0; j<3; j++) {
      PMTpos[icab][j]=geopmt_.xyzpm[icab][j];
      PMTdir[icab][j]=geopmt_.dxyzpm[icab][j];
    }
    PMTtype[icab]=pmtinfcmn_.isk23[icab]-2;//set PMT type(0:Old, 1:New)
    if (!(PMTtype[icab]==0 || PMTtype[icab]==1)) PMTtype[icab]=-1;//Missing PMTs
//    if (combad_.ibad[icab]!=0) PMTtype[icab]=-1;//ignore bad channels - do this event-by-event!
    
    PMTQE[icab]=pmtinfcmn_.qefactor[icab];
    
  }
  /// else compiler flage
#else
  /// now fill from wcsim root file

  WCSimWrap * wc  = WCSimWrap::Get();
  float ** geopmtxyz = wc->PMTposs();
  float ** geopmtdxyz = wc->PMTdirs();
  int * geopmttype = wc->PMTtypes();
  float * geopmtqe = wc->PMTQEs();
//  float * geopmtr = wc->PMTRs();
//  float * geopmtphi = wc->PMTphis();
  PMTRadius = wc->PMTRad();
  std::cout << "PMTRadius was found to be " << PMTRadius << std::endl;
  for (int icab=0; icab<fMaxNPMT; icab++) {
    //std::cout<<"<fiTQun_shared::ReadParameters> icab="<<icab<<" / "<<fMaxNPMT<<std::endl;
    //std::cout<<"\tx,y,z="
    //	     <<geopmtxyz[icab][0]<<", "
    //	     <<geopmtxyz[icab][1]<<", "
    //	     <<geopmtxyz[icab][2]<<std::endl;

    for (int j=0; j<3; j++) {
      PMTpos[icab][j]=geopmtxyz[icab][j];
      PMTdir[icab][j]=geopmtdxyz[icab][j];
    }
    PMTtype[icab]=geopmttype[icab];//set PMT type(0:Old, 1:New)
    if (!(PMTtype[icab]==0 || PMTtype[icab]==1)) PMTtype[icab]=-1;//Dead PMT
    
    PMTQE[icab]=geopmtqe[icab];

  }

  TransDetToGlbl.SetXYZ(wc->LocCoordOffset(0),wc->LocCoordOffset(1),wc->LocCoordOffset(2));
  
  DetZinGlbl.SetXYZ(wc->LocCoordZdir(0),wc->LocCoordZdir(1),wc->LocCoordZdir(2));
  DetXinGlbl.SetXYZ(wc->LocCoordXdir(0),wc->LocCoordXdir(1),wc->LocCoordXdir(2));
  
#endif
  ///
    
  RotDetToGlbl.SetZAxis(DetZinGlbl,DetXinGlbl);// Rotation matrix which transforms detector local frame to global frame
  
  double PMTmaxR=0.;
  for (int j=0; j<3; j++) PMTmaxPos[j]=0.;
  for (int icab=0; icab<fMaxNPMT; icab++) {

    // Transform PMT position&direction to detector local frame
    TVector3 tubepos = TfmDetToGlblCoord(TVector3(PMTpos[icab]),true,true);
    tubepos.GetXYZ(PMTpos[icab]);
    
    TVector3 tubedir = TfmDetToGlblCoord(TVector3(PMTdir[icab]),false,true);
    tubedir.SetMag(1.);
    tubedir.GetXYZ(PMTdir[icab]);
    
    PMT_R[icab] = tubepos.Perp();
    PMTphi[icab] = tubepos.Phi();
    
    for (int j=0; j<3; j++) if (PMTmaxPos[j]<fabs(PMTpos[icab][j])) PMTmaxPos[j]=fabs(PMTpos[icab][j]);
    if (PMTmaxR<PMT_R[icab]) PMTmaxR=PMT_R[icab];
  }
  
  DetGeomType=iDetGeom;
  DetGeomR=PMTmaxR;
  DetGeomL=PMTmaxPos[2];

  std::cout << "Setting detector geometry: nPMT=" << fMaxNPMT << ", Type=" << DetGeomType << ", R="<< DetGeomR << ", L=" << DetGeomL << std::endl;

  std::cout << Form("Local coord. origin in global coord.: [%.1f,%.1f,%.1f]",TransDetToGlbl.X(),TransDetToGlbl.Y(),TransDetToGlbl.Z()) << std::endl;
  std::cout << Form("Local coord. X-axis in global coord.: [%.3f,%.3f,%.3f]",RotDetToGlbl.XX(),RotDetToGlbl.YX(),RotDetToGlbl.ZX()) << std::endl;
  std::cout << Form("Local coord. Y-axis in global coord.: [%.3f,%.3f,%.3f]",RotDetToGlbl.XY(),RotDetToGlbl.YY(),RotDetToGlbl.ZY()) << std::endl;
  std::cout << Form("Local coord. Z-axis in global coord.: [%.3f,%.3f,%.3f]",RotDetToGlbl.XZ(),RotDetToGlbl.YZ(),RotDetToGlbl.ZZ()) << std::endl;
  
  // Set Additional Hyper-K variables (indexed by quadrant)
  DetGeomHKR[0] = 2909.5; // first quadrant
  DetGeomHKR[1] = 2909.5; // second quadrant
  DetGeomHKR[2] = 2709.5; // third quadrant
  DetGeomHKR[3] = 2709.5; // fourth quadrant
  DetGeomHKArcXCenter[0] = -800.; // first quadrant
  DetGeomHKArcXCenter[1] =  800.; // second quadrant
  DetGeomHKArcXCenter[2] =  600.; // third quadrant
  DetGeomHKArcXCenter[3] = -600.; // fourth quadrant

  // The code assumes that the maximum x-(y-)value and the minimum x-(y-) value
  // are equidistant from the origin.
  DetGeomHKXMax = DetGeomHKR[0] + DetGeomHKArcXCenter[0]; // maximum x-value of ID wall
  DetGeomHKYMax = 2109.5; // maximum y-value of ID wall
  DetGeomHKXMaxAtYMax = DetGeomHKArcXCenter[0] + sqrt(DetGeomHKR[0]*DetGeomHKR[0] - DetGeomHKYMax*DetGeomHKYMax);
  DetGeomHKXMaxAtYMin = DetGeomHKArcXCenter[3] + sqrt(DetGeomHKR[3]*DetGeomHKR[3] - DetGeomHKYMax*DetGeomHKYMax);

  LoadProfiles(fFitCProf, fFittProf);
  
  flgLoaded=true;

  // Read in angular response function:
  // Check if angular response file path was given in parameters funciton
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ){
    // If yes, open it and populate array
    WCSim_PMTType = TRuntimeParameters::Get().GetParameterS("fiTQun.WCSimConfig");
    
    TString angResp_FileName = Form("%s/%s_%s.root",fQconstdir.Data(),"angResp", WCSim_PMTType.Data());
    
    TFile * angRespIn = new TFile(angResp_FileName.Data());

    if (! angRespIn->IsOpen() ){
      std::cerr << "Error: PMT angular response file not found. Exiting..." << std::endl;
      std::cerr << angResp_FileName.Data() << std::endl;
      exit(-1);
    }
    
    std::cout << "Reading PMT angular response curve from " << angResp_FileName.Data() << std::endl;

    TF1 * angRespFunc = (TF1*) angRespIn->Get("angResp");
    
    double cosetamax=1.;
    for (int i = 0; i <= nEpsilonBins; i++) {
      epsilon_bins[i] = i*cosetamax/nEpsilonBins;
      epsilon_Tbl[i] = angRespFunc->Eval(epsilon_bins[i]);
    }
    epsilon_binwrecp = 1./(epsilon_bins[1]-epsilon_bins[0]);

    angRespIn->Close();
  } else {
    // If not populate sample "old" function and populate array
    std::cout << "Warning: PMT type not specified. Using default (SK) angular response function." << std::endl;
    
    double cosetamax=1.;
    for (int i = 0; i <= nEpsilonBins; i++) {
      epsilon_bins[i] = i*cosetamax/nEpsilonBins;
      
      // "Old" epsilon function
      if (epsilon_bins[i] < 0.196688){
	double val = -0.183124+epsilon_bins[i]*(3.4469-epsilon_bins[i]*4.6833);
	if (val > 0.) {
	  epsilon_Tbl[i] = val;
	} else {
	  epsilon_Tbl[i] = 0.;
	}
      } else if (epsilon_bins[i] < 0.986394) {
	epsilon_Tbl[i] = -0.207593+epsilon_bins[i]*(4.03995+epsilon_bins[i]*(-8.99093+epsilon_bins[i]*(10.6753-epsilon_bins[i]*4.5175)));
      } else {
	epsilon_Tbl[i] = 0.879163+epsilon_bins[i]*0.120837;
      }
    }
    epsilon_binwrecp = 1./(epsilon_bins[1]-epsilon_bins[0]);
  }
}

TVector3 fiTQun_shared::TfmDetToGlblCoord(TVector3 tmpVct, bool flgPos, bool flgInverse){
  
  if (flgInverse) {
    if (flgPos) tmpVct-=TransDetToGlbl;
    tmpVct*=RotDetToGlbl.Inverse();
  }
  else {
    tmpVct*=RotDetToGlbl;
    if (flgPos) tmpVct+=TransDetToGlbl;
  }
  
  return tmpVct;
  
}

void fiTQun_shared::Boundlogmu(double& logmu) {
  if (!(logmu>=tpdfmurang[0] && logmu<tpdfmurang[1])) {
    if (logmu<tpdfmurang[0]) logmu=tpdfmurang[0];//ensure logmu is in the fitted range
    else logmu=tpdfmurang[1];
  }
}


void fiTQun_shared::IntrpScat3Dto2D(double *dwall) {
  for (int istp=0; istp<3; istp++) {
    double dwalltmp=dwall[istp];
    
    double diwl=nsct3dwlbinm1*(dwalltmp-sct3dwlbins[0])/(sct3dwlbins[nsct3dwlbinm1]-sct3dwlbins[0]);
    int iwl=int(diwl);
    int ipwl=iwl+1;
    if (!(diwl < nsct3dwlbinm1 && diwl > 0.)) {
      if (diwl >= nsct3dwlbinm1) {
        ipwl = nsct3dwlbinm1;
        iwl = ipwl;
      }
      else {
        iwl = 0;
        ipwl = 0;
      }
    }
    
    double yd=(dwalltmp-sct3dwlbins[iwl])/sct3dwlbinw;
    
    for (int iR=0; iR<nsct2dRbin; iR++) {
      for (int ith=0; ith<nsct2dthbin; ith++) {
        SctTbl2d[istp][iR*nsct2dthbin+ith]=SctTbl3d[iR][iwl][ith]*(1.-yd)+SctTbl3d[iR][ipwl][ith]*yd;
      }
    }
  }
}

double fiTQun_shared::GetScatRatio2D(int istp, double R, double costh) {
  
  double dR=(R-sct2dRbins[0])*sct2dRbinwrecp;
  int iR=(int)dR;
  double dth=(costh-sct2dthbins[0])*sct2dthbinwrecp;
  int ith=(int)dth;
  
  if (iR > nsct2dRbin-2) {
    iR = nsct2dRbin - 2;
  }
  else if (iR < 0) {
    iR = 0;
  }
  
  if (ith > nsct2dthbin-2) {
    ith = nsct2dthbin - 2;
  }
  else if (ith < 0) {
    ith = 0;
  }
  
  double xd=dR-(double)iR;
  double yd=dth-(double)ith;
  
  int ibase=iR*nsct2dthbin+ith;
  double w1=SctTbl2d[istp][ibase]*(1.-yd)+SctTbl2d[istp][ibase+1]*yd;
  ibase+=nsct2dthbin;
  double w2=SctTbl2d[istp][ibase]*(1.-yd)+SctTbl2d[istp][ibase+1]*yd;
  
  return w1*(1.-xd)+w2*xd;
}

double fiTQun_shared::explintrp(double x) {
  double dix=x*expbinwrecp;
  
  if (!(dix > 0. && dix < nexpbinm1)) {
    return exp(x);
  }
  
  int ix=(int)dix;
  double xd=(x-expbins[ix])*expbinwrecp;
  return expTbl[ix]*(1.-xd)+expTbl[ix+1]*xd;
}

double fiTQun_shared::epsilon(double coseta) {
  double dix=coseta*epsilon_binwrecp;
  
  if (!(dix >= 0. && dix < nEpsilonBins)) {
    return 0;
  }
  
  int ix=(int)dix;
  double xd=(coseta-epsilon_bins[ix])*epsilon_binwrecp;
  return epsilon_Tbl[ix]*(1.-xd)+epsilon_Tbl[ix+1]*xd;
}

double fiTQun_shared::pCherenkovThr(int iPID) {
  const double nphase=1.334;
  return m0[iPID]/sqrt(nphase*nphase-1.);
}

//double fiTQun_shared::WAttenL(double Zsrc, double Zpmt) {
//  double Zbound = -1100.;
//  if (Zsrc > Zbound) {
//    if (Zpmt > Zbound) {
//    } else {
//    }
//  } else {
//    if (Zpmt > Zbound) {
//    } else {
//    }
//  }
//}


fiTQun_shared::~fiTQun_shared(){
  DeleteInstances();
  std::cout << "fiTQun_shared destructed!!" << std::endl;
}

void fiTQun_shared::DeleteInstances(){
  
  DeleteProfiles();
  
  if (SctTbl2d[0]) {
    delete[] sct2dRbins;
    delete[] sct2dthbins;
    for (int istp=0; istp<3; istp++) {
      delete[] SctTbl2d[istp];
    }
  }
  
  if (SctTbl3d) {
    delete[] sct3dRbins;
    delete[] sct3dwlbins;
    delete[] sct3dthbins;
    for (int iR=0; iR<=nsct3dRbinm1; iR++) {
      for (int iwl=0; iwl<=nsct3dwlbinm1; iwl++) {
        delete[] SctTbl3d[iR][iwl];
      }
      delete[] SctTbl3d[iR];
    }
    delete[] SctTbl3d;
  }
  
}


void fiTQun_shared::SaveTimeWindow(int itwnd, double *twindow){
  fitquntwnd_.fqtwnd[itwnd][0] = twindow[0];
  fitquntwnd_.fqtwnd[itwnd][1] = twindow[1];
}

void fiTQun_shared::GetPreFit(int itwnd, double* fitparams, double *twindow){
  //  fitparams[0] = fitquntwnd_.fqtwnd_prftpos[itwnd][0];
  //  fitparams[1] = fitquntwnd_.fqtwnd_prftpos[itwnd][1];
  //  fitparams[2] = fitquntwnd_.fqtwnd_prftpos[itwnd][2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitquntwnd_.fqtwnd_prftpos[itwnd]),true,true);
  tmppos.GetXYZ(fitparams);
  
  fitparams[3] = fitquntwnd_.fqtwnd_prftt0 [itwnd];
  if (twindow) {
    twindow[0] = fitquntwnd_.fqtwnd[itwnd][0];
    twindow[1] = fitquntwnd_.fqtwnd[itwnd][1];
  }
}

void fiTQun_shared::ClearPeaks(int itwnd){
  fitquntwnd_.fqtwnd_npeak[itwnd] = 0;
  for (int ipeak=0; ipeak<MAXNPEAK; ipeak++) {
    fitquntwnd_.fqtwnd_peakt0   [itwnd][ipeak] = 0.;
    fitquntwnd_.fqtwnd_peakiness[itwnd][ipeak] = 0.;
  }
}

void fiTQun_shared::SavePeak(int itwnd, double peakt0, double peakiness){
  if (fitquntwnd_.fqtwnd_npeak[itwnd] >= MAXNPEAK) {
    std::cout << "In fiTQun_shared::SavePeak: cannot save more than " << MAXNPEAK << " peak information." << std::endl;
    std::cout << "     ---> The following will not be saved (itwnd, peakt0, peakiness) = (" << itwnd << "," << peakt0 << "," << peakiness << ")." << std::endl;
    return;
  }
  
  int i=fitquntwnd_.fqtwnd_npeak[itwnd];
  fitquntwnd_.fqtwnd_peakt0   [itwnd][i] = peakt0;
  fitquntwnd_.fqtwnd_peakiness[itwnd][i] = peakiness;
  fitquntwnd_.fqtwnd_npeak[itwnd]++;
}

void fiTQun_shared::GetPeaks(int itwnd, int& npeaks, double *peakt0, double *peakiness){
  npeaks=fitquntwnd_.fqtwnd_npeak[itwnd];
  
  for (int ipeak=0; ipeak<npeaks; ipeak++) {
    peakt0[ipeak]=fitquntwnd_.fqtwnd_peakt0[itwnd][ipeak];
    peakiness[ipeak]=fitquntwnd_.fqtwnd_peakiness[itwnd][ipeak];
  }
}

void fiTQun_shared::SaveSubEvtInfo(int itwnd, int ipeak, int nHPMT, double Totalq){
  int ise=ipeak;
  for (int idx=0; idx<itwnd; idx++) {
    ise+=fitquntwnd_.fqtwnd_npeak[idx];
  }
  if (ise >= MAXNPEAK) {
    std::cout << "In fiTQun_shared::SaveSubEvtInfo: cannot save more than " << MAXNPEAK << " sub-event information." << std::endl;
    std::cout << "     ---> The following will not be saved (ise, nHPMT, Totalq) = (" << ise << "," << nHPMT << "," << Totalq << ")." << std::endl;
    return;
  }
  
  fitqun1r_.fqitwnd  [ise] = itwnd;
  fitqun1r_.fqipeak  [ise] = ipeak;
  fitqun1r_.fqnhitpmt[ise] = nHPMT;
  fitqun1r_.fqtotq   [ise] = Totalq;
  
  if (fitqun1r_.fqnse<ise+1) fitqun1r_.fqnse=ise+1;
}

void fiTQun_shared::SaveDefaultSnglTrkFit(double* fitparams, int itwnd, int ipeak, int iPID, double totmu, double nll, int PCflg, int *n50, double *q50){
  
  int ise=ipeak;
  for (int idx=0; idx<itwnd; idx++) {
    ise+=fitquntwnd_.fqtwnd_npeak[idx];
  }
  if (ise >= MAXNPEAK) {
    std::cout << "In fiTQun_shared::SaveDefaultSnglTrkFit: cannot save more than " << MAXNPEAK << " one-ring information." << std::endl;
    return;
  }
  
  //  fitqun1r_.fq1rpos  [ise][iPID][0] = fitparams[0];
  //  fitqun1r_.fq1rpos  [ise][iPID][1] = fitparams[1];
  //  fitqun1r_.fq1rpos  [ise][iPID][2] = fitparams[2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitqun1r_.fq1rpos[ise][iPID]);
  
  fitqun1r_.fq1rt0   [ise][iPID]    = fitparams[3];
  
  //fitqun1r_.fq1rdir  [ise][iPID][0] = sin(fitparams[4])*cos(fitparams[5]);
  //  fitqun1r_.fq1rdir  [ise][iPID][1] = sin(fitparams[4])*sin(fitparams[5]);
  //  fitqun1r_.fq1rdir  [ise][iPID][2] = cos(fitparams[4]);
  TVector3 orgdir;
  orgdir.SetMagThetaPhi(1.,fitparams[4],fitparams[5]);
  TVector3 tmpdir = TfmDetToGlblCoord(orgdir,false,false);
  tmpdir.GetXYZ(fitqun1r_.fq1rdir[ise][iPID]);
  
  fitqun1r_.fq1rmom  [ise][iPID]    = fitparams[6];
  
  if (fRngTyp[iPID]==2) {// upstream-track
    fitqun1r_.fq1rdconv[ise][iPID]  = 0.;
    fitqun1r_.fq1reloss[ise][iPID]  = fitparams[7];
  }
  else {
    fitqun1r_.fq1rdconv[ise][iPID]  = fitparams[7];
    fitqun1r_.fq1reloss[ise][iPID]  = 0.;
  }
  
  fitqun1r_.fq1rtotmu[ise][iPID]    = totmu;
  fitqun1r_.fq1rnll  [ise][iPID]    = nll;
  fitqun1r_.fq1rpcflg[ise][iPID]    = PCflg;
  
  if (iPID==ie && n50 && q50) {//save n50
    fitqun1r_.fqn50[ise]=*n50;
    fitqun1r_.fqq50[ise]=*q50;
  }
  
  //std::cout << "fq1rpos["<< ise << "][" << iPID << "][0] = " << fitqun1r_.fq1rpos[ise][iPID][0] << std::endl;
  //std::cout << "-lnL=" << fitqun1r_.fq1rnll[ise][iPID]  << "; totq, totmu = " << fitqun1r_.fq1rtotmu[ise][iPID] << ", " << fitqun1r_.fq1rtotmu[ise][iPID] << ";  Time=" << fitqun1r_.fq1rt0   [ise][iPID] << std::endl;
}

void fiTQun_shared::GetDefaultSnglTrkFit(double* fitparams, int itwnd, int ipeak, int iPID, double& totmu, double& nll, int& PCflg){
  int ise=ipeak;
  for (int idx=0; idx<itwnd; idx++) {
    ise+=fitquntwnd_.fqtwnd_npeak[idx];
  }
  if (ise >= MAXNPEAK) {
    std::cout << "In fiTQun_shared::GetDefaultSnglTrkFit: ise >= MAXNPEAK!" << std::endl;
    return;
  }
  
  //  fitparams[0] = fitqun1r_.fq1rpos  [ise][iPID][0];
  //  fitparams[1] = fitqun1r_.fq1rpos  [ise][iPID][1];
  //  fitparams[2] = fitqun1r_.fq1rpos  [ise][iPID][2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitqun1r_.fq1rpos[ise][iPID]),true,true);
  tmppos.GetXYZ(fitparams);
  
  fitparams[3] = fitqun1r_.fq1rt0   [ise][iPID];
  
  //  fitparams[4] = acos (fitqun1r_.fq1rdir[ise][iPID][2]);
  //  fitparams[5] = atan2(fitqun1r_.fq1rdir[ise][iPID][1],fitqun1r_.fq1rdir[ise][iPID][0]);
  TVector3 tmpdir = TfmDetToGlblCoord(TVector3(fitqun1r_.fq1rdir[ise][iPID]),false,true);
  fitparams[4] = tmpdir.Theta();
  fitparams[5] = tmpdir.Phi();

  fitparams[6] = fitqun1r_.fq1rmom  [ise][iPID];
  
  if (fRngTyp[iPID]==2) {// upstream-track
    fitparams[7] = fitqun1r_.fq1reloss[ise][iPID];
  }
  else {
    fitparams[7] = fitqun1r_.fq1rdconv[ise][iPID];
  }
  
  totmu        = fitqun1r_.fq1rtotmu[ise][iPID];
  nll          = fitqun1r_.fq1rnll  [ise][iPID];
  PCflg        = fitqun1r_.fq1rpcflg[ise][iPID];
}

void fiTQun_shared::SaveTestSnglTrkFit(double* fitparams, int istage, int itwnd, int iPID, double totmu, double nll, int PCflg){
  if (fitqun1rtest_.fqtestn1r >= FQTESTMAX1R) {
    std::cout << "In fiTQun_shared::SaveTestSnglTrkFit: cannot save more than " << FQTESTMAX1R << " test 1 ring fits." << std::endl;
    std::cout << "     ---> This test fit will not be saved (itwnd, iPID, totmu, nll) = (" << itwnd << "," << iPID << "," << totmu << "," << nll << ")." << std::endl;
    return;
  }

  fitqun1rtest_.fqtest1rstage[fitqun1rtest_.fqtestn1r]    = istage;
  fitqun1rtest_.fqtest1rse   [fitqun1rtest_.fqtestn1r]    = itwnd;
  fitqun1rtest_.fqtest1rpid  [fitqun1rtest_.fqtestn1r]    = iPID;

  //  fitqun1rtest_.fqtest1rpos  [fitqun1rtest_.fqtestn1r][0] = fitparams[0];
  //  fitqun1rtest_.fqtest1rpos  [fitqun1rtest_.fqtestn1r][1] = fitparams[1];
  //  fitqun1rtest_.fqtest1rpos  [fitqun1rtest_.fqtestn1r][2] = fitparams[2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitqun1rtest_.fqtest1rpos[fitqun1rtest_.fqtestn1r]);

  fitqun1rtest_.fqtest1rt0   [fitqun1rtest_.fqtestn1r]    = fitparams[3];

  //  fitqun1rtest_.fqtest1rdir  [fitqun1rtest_.fqtestn1r][0] = sin(fitparams[4])*cos(fitparams[5]);
  //  fitqun1rtest_.fqtest1rdir  [fitqun1rtest_.fqtestn1r][1] = sin(fitparams[4])*sin(fitparams[5]);
  //  fitqun1rtest_.fqtest1rdir  [fitqun1rtest_.fqtestn1r][2] = cos(fitparams[4]);
  TVector3 orgdir;
  orgdir.SetMagThetaPhi(1.,fitparams[4],fitparams[5]);
  TVector3 tmpdir = TfmDetToGlblCoord(orgdir,false,false);
  tmpdir.GetXYZ(fitqun1rtest_.fqtest1rdir[fitqun1rtest_.fqtestn1r]);

  fitqun1rtest_.fqtest1rmom  [fitqun1rtest_.fqtestn1r]    = fitparams[6];
  
  if (fRngTyp[iPID]==2) {// upstream-track
    fitqun1rtest_.fqtest1rdconv[fitqun1rtest_.fqtestn1r]  = 0.;
    fitqun1rtest_.fqtest1reloss[fitqun1rtest_.fqtestn1r]  = fitparams[7];
  }
  else {
    fitqun1rtest_.fqtest1rdconv[fitqun1rtest_.fqtestn1r]  = fitparams[7];
    fitqun1rtest_.fqtest1reloss[fitqun1rtest_.fqtestn1r]  = 0.;
  }
  
  fitqun1rtest_.fqtest1rtotmu[fitqun1rtest_.fqtestn1r]    = totmu;
  fitqun1rtest_.fqtest1rnll  [fitqun1rtest_.fqtestn1r]    = nll;
  fitqun1rtest_.fqtest1rpcflg[fitqun1rtest_.fqtestn1r]    = PCflg;

  fitqun1rtest_.fqtestn1r++;
}

void fiTQun_shared::SaveDefaultPi0Fit(double* fitparams, int index, double totmu, double nll, int PCflg){
  //  fitqunpi0_.fqpi0pos[index][0] = fitparams[0];
  //  fitqunpi0_.fqpi0pos[index][1] = fitparams[1];
  //  fitqunpi0_.fqpi0pos[index][2] = fitparams[2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitqunpi0_.fqpi0pos[index]);

  fitqunpi0_.fqpi0t0[index] = fitparams[3];

  //  fitqunpi0_.fqpi0dir1[index][0] = sin(fitparams[4])*cos(fitparams[5]);
  //  fitqunpi0_.fqpi0dir1[index][1] = sin(fitparams[4])*sin(fitparams[5]);
  //  fitqunpi0_.fqpi0dir1[index][2] = cos(fitparams[4]);
  TVector3 orgdir1;
  orgdir1.SetMagThetaPhi(1.,fitparams[4],fitparams[5]);
  TVector3 tmpdir1 = TfmDetToGlblCoord(orgdir1,false,false);
  tmpdir1.GetXYZ(fitqunpi0_.fqpi0dir1[index]);

  fitqunpi0_.fqpi0mom1[index] = fitparams[6];
  fitqunpi0_.fqpi0dconv1[index] = fitparams[7];

  //  fitqunpi0_.fqpi0dir2[index][0] = sin(fitparams[8])*cos(fitparams[9]);
  //  fitqunpi0_.fqpi0dir2[index][1] = sin(fitparams[8])*sin(fitparams[9]);
  //  fitqunpi0_.fqpi0dir2[index][2] = cos(fitparams[8]);
  TVector3 orgdir2;
  orgdir2.SetMagThetaPhi(1.,fitparams[8],fitparams[9]);
  TVector3 tmpdir2 = TfmDetToGlblCoord(orgdir2,false,false);
  tmpdir2.GetXYZ(fitqunpi0_.fqpi0dir2[index]);

  fitqunpi0_.fqpi0mom2[index] = fitparams[10];
  fitqunpi0_.fqpi0dconv2[index] = fitparams[11];
  
  fitqunpi0_.fqpi0mass[index] = sqrt( 2.*fitqunpi0_.fqpi0mom1[index]*fitqunpi0_.fqpi0mom2[index]
				      *(1 - tmpdir1.Dot(tmpdir2)) );
  fitqunpi0_.fqpi0photangle[index] = acos(tmpdir1.Dot(tmpdir2));
  TVector3 momtot = fitqunpi0_.fqpi0mom1[index]*tmpdir1 + fitqunpi0_.fqpi0mom2[index]*tmpdir2;
  fitqunpi0_.fqpi0momtot[index] = momtot.Mag();
  TVector3 dirtot = momtot.Unit();
  fitqunpi0_.fqpi0dirtot[index][0] = dirtot.X();
  fitqunpi0_.fqpi0dirtot[index][1] = dirtot.Y();
  fitqunpi0_.fqpi0dirtot[index][2] = dirtot.Z();

  fitqunpi0_.fqpi0totmu[index] = totmu;
  fitqunpi0_.fqpi0nll[index] = nll;
  fitqunpi0_.fqpi0pcflg[index] = PCflg;
}

void fiTQun_shared::SaveTestPi0Fit(double* fitparams, int istage, double totmu, double nll, int PCflg){
  if (fitqunpi0test_.fqtestnpi0 >= FQTESTMAXPI0) {
    std::cout << "In fiTQun_shared::SaveTestPi0Fit: cannot save more than " << FQTESTMAXPI0 << " test 1 ring fits." << std::endl;
    std::cout << "     ---> This test fit will not be saved (istage = " << istage << ")" << std::endl;
    return;
  }

  fitqunpi0test_.fqtestpi0stage[fitqunpi0test_.fqtestnpi0]  = istage;

  //  fitqunpi0test_.fqtestpi0pos[fitqunpi0test_.fqtestnpi0][0] = fitparams[0];
  //  fitqunpi0test_.fqtestpi0pos[fitqunpi0test_.fqtestnpi0][1] = fitparams[1];
  //  fitqunpi0test_.fqtestpi0pos[fitqunpi0test_.fqtestnpi0][2] = fitparams[2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitqunpi0test_.fqtestpi0pos[fitqunpi0test_.fqtestnpi0]);
  
  fitqunpi0test_.fqtestpi0t0[fitqunpi0test_.fqtestnpi0] = fitparams[3];

  //  fitqunpi0test_.fqtestpi0dir1[fitqunpi0test_.fqtestnpi0][0] = sin(fitparams[4])*cos(fitparams[5]);
  //  fitqunpi0test_.fqtestpi0dir1[fitqunpi0test_.fqtestnpi0][1] = sin(fitparams[4])*sin(fitparams[5]);
  //  fitqunpi0test_.fqtestpi0dir1[fitqunpi0test_.fqtestnpi0][2] = cos(fitparams[4]);
  TVector3 orgdir1;
  orgdir1.SetMagThetaPhi(1.,fitparams[4],fitparams[5]);
  TVector3 tmpdir1 = TfmDetToGlblCoord(orgdir1,false,false);
  tmpdir1.GetXYZ(fitqunpi0test_.fqtestpi0dir1[fitqunpi0test_.fqtestnpi0]);
  
  fitqunpi0test_.fqtestpi0mom1[fitqunpi0test_.fqtestnpi0] = fitparams[6];
  fitqunpi0test_.fqtestpi0dconv1[fitqunpi0test_.fqtestnpi0] = fitparams[7];

  //  fitqunpi0test_.fqtestpi0dir2[fitqunpi0test_.fqtestnpi0][0] = sin(fitparams[8])*cos(fitparams[9]);
  //  fitqunpi0test_.fqtestpi0dir2[fitqunpi0test_.fqtestnpi0][1] = sin(fitparams[8])*sin(fitparams[9]);
  //  fitqunpi0test_.fqtestpi0dir2[fitqunpi0test_.fqtestnpi0][2] = cos(fitparams[8]);
  TVector3 orgdir2;
  orgdir2.SetMagThetaPhi(1.,fitparams[8],fitparams[9]);
  TVector3 tmpdir2 = TfmDetToGlblCoord(orgdir2,false,false);
  tmpdir2.GetXYZ(fitqunpi0test_.fqtestpi0dir2[fitqunpi0test_.fqtestnpi0]);
  
  fitqunpi0test_.fqtestpi0mom2[fitqunpi0test_.fqtestnpi0] = fitparams[10];
  fitqunpi0test_.fqtestpi0dconv2[fitqunpi0test_.fqtestnpi0] = fitparams[11];

  fitqunpi0test_.fqtestpi0mass[fitqunpi0test_.fqtestnpi0] = sqrt( 2.*fitparams[6]*fitparams[10]
			       *(1 - tmpdir1.Dot(tmpdir2)) );
  fitqunpi0test_.fqtestpi0photangle[fitqunpi0test_.fqtestnpi0] = acos(tmpdir1.Dot(tmpdir2));
  TVector3 momtot = fitparams[6]*tmpdir1 + fitparams[10]*tmpdir2;
  fitqunpi0test_.fqtestpi0momtot[fitqunpi0test_.fqtestnpi0] = momtot.Mag();
  TVector3 dirtot = momtot.Unit();
  fitqunpi0test_.fqtestpi0dirtot[fitqunpi0test_.fqtestnpi0][0] = dirtot.X();
  fitqunpi0test_.fqtestpi0dirtot[fitqunpi0test_.fqtestnpi0][1] = dirtot.Y();
  fitqunpi0test_.fqtestpi0dirtot[fitqunpi0test_.fqtestnpi0][2] = dirtot.Z();

  fitqunpi0test_.fqtestpi0totmu[fitqunpi0test_.fqtestnpi0] = totmu;
  fitqunpi0test_.fqtestpi0nll[fitqunpi0test_.fqtestnpi0] = nll;
  fitqunpi0test_.fqtestpi0pcflg[fitqunpi0test_.fqtestnpi0] = PCflg;

  fitqunpi0test_.fqtestnpi0++;
}

void fiTQun_shared::ClearMRFit() {
  for (int iFit=0; iFit<FQMAXNMRFIT; iFit++) {
    SaveMRFit(0,0,NULL,NULL,0.,0.,0,iFit);
  }
  fitqunmr_.fqnmrfit=0;
}

void fiTQun_shared::SaveMRFit(int idFit, int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double totmu, double nll, int PCflg, int iFit){
  
  if (iFit<0) {
    iFit=0;
    while (iFit<fitqunmr_.fqnmrfit) {
      if (fitqunmr_.fqmrifit[iFit] == idFit) {
        break;
      }
      iFit++;
    }
  }
  
  if (iFit >= FQMAXNMRFIT) {
    std::cout << "In fiTQun_shared::SaveMRFit: cannot save more than " << FQMAXNMRFIT << " multi-ring fits." << std::endl;
    std::cout << "     ---> This MR fit will not be saved (idFit = " << idFit << ")" << std::endl;
    return;
  }
  if (!(nringtmp>=0 && nringtmp<=FQMAXNRING)) {
    std::cout << "In fiTQun_shared::SaveMRFit: nring must be [0-" << FQMAXNRING << "]!" << std::endl;
    exit(-1);
  }
  
  if (fitqunmr_.fqnmrfit<iFit+1) fitqunmr_.fqnmrfit=iFit+1;
  
  fitqunmr_.fqmrifit[iFit] = idFit;
  fitqunmr_.fqmrnring[iFit] = nringtmp;
  fitqunmr_.fqmrpcflg[iFit] = PCflg;
  fitqunmr_.fqmrnll[iFit] = nll;
  fitqunmr_.fqmrtotmu[iFit] = totmu;
  for (int iring=0; iring<nringtmp; iring++) {
    fitqunmr_.fqmrpid[iFit][iring] = TypeArr[iring];
    fitqunmr_.fqmrmom[iFit][iring] = ParArr[iring][6];
    if (fRngTyp[TypeArr[iring]]==2) {// upstream-track
      fitqunmr_.fqmrdconv[iFit][iring] = 0.;
      fitqunmr_.fqmreloss[iFit][iring] = ParArr[iring][7];
    }
    else {
      fitqunmr_.fqmrdconv[iFit][iring] = ParArr[iring][7];
      fitqunmr_.fqmreloss[iFit][iring] = 0.;
    }
    fitqunmr_.fqmrt0[iFit][iring] = ParArr[iring][3];

    //    for (int ix=0; ix<3; ix++) {
    //      fitqunmr_.fqmrpos[iFit][iring][ix] = ParArr[iring][ix];
    //    }
    TVector3 tmppos = TfmDetToGlblCoord(TVector3(ParArr[iring]),true,false);
    tmppos.GetXYZ(fitqunmr_.fqmrpos[iFit][iring]);

    //    fitqunmr_.fqmrdir[iFit][iring][0] = sin(ParArr[iring][4])*cos(ParArr[iring][5]);
    //    fitqunmr_.fqmrdir[iFit][iring][1] = sin(ParArr[iring][4])*sin(ParArr[iring][5]);
    //    fitqunmr_.fqmrdir[iFit][iring][2] = cos(ParArr[iring][4]);
    TVector3 orgdir;
    orgdir.SetMagThetaPhi(1.,ParArr[iring][4],ParArr[iring][5]);
    TVector3 tmpdir = TfmDetToGlblCoord(orgdir,false,false);
    tmpdir.GetXYZ(fitqunmr_.fqmrdir[iFit][iring]);
  }
  
  for (int iring=nringtmp; iring<FQMAXNRING; iring++) {// clear the rest
    fitqunmr_.fqmrpid[iFit][iring] = 0;
    fitqunmr_.fqmrmom[iFit][iring] = 0.;
    fitqunmr_.fqmrdconv[iFit][iring] = 0.;
    fitqunmr_.fqmreloss[iFit][iring] = 0.;
    fitqunmr_.fqmrt0[iFit][iring] = 0.;
    for (int ix=0; ix<3; ix++) {
      fitqunmr_.fqmrpos[iFit][iring][ix] = 0.;
      fitqunmr_.fqmrdir[iFit][iring][ix] = 0.;
    }
  }
  
}

int fiTQun_shared::GetMRFit(int idFit, int &nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double &totmu, double &nll, int &PCflg, int iFit){
  
  if (iFit<0) {
    iFit=0;
    while (1) {
      if (iFit>=fitqunmr_.fqnmrfit) {
        //      std::cout << "In fiTQun_shared::GetMRFit: cannot find result of (idFit = " << idFit << ")" << std::endl;
        return -1;
      }
      
      if (fitqunmr_.fqmrifit[iFit] == idFit) {
        break;
      }
      iFit++;
    }
  }
  
  nringtmp = fitqunmr_.fqmrnring[iFit];
  PCflg = fitqunmr_.fqmrpcflg[iFit];
  nll = fitqunmr_.fqmrnll[iFit];
  totmu = fitqunmr_.fqmrtotmu[iFit];
  for (int iring=0; iring<nringtmp; iring++) {
    TypeArr[iring] = fitqunmr_.fqmrpid[iFit][iring];
    ParArr[iring][6] = fitqunmr_.fqmrmom[iFit][iring];
    if (fRngTyp[TypeArr[iring]]==2) {// upstream-track
      ParArr[iring][7] = fitqunmr_.fqmreloss[iFit][iring];
    }
    else {
      ParArr[iring][7] = fitqunmr_.fqmrdconv[iFit][iring];
    }
    ParArr[iring][3] = fitqunmr_.fqmrt0[iFit][iring];

    //    for (int ix=0; ix<3; ix++) {
    //      ParArr[iring][ix] = fitqunmr_.fqmrpos[iFit][iring][ix];
    //    }
    TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitqunmr_.fqmrpos[iFit][iring]),true,true);
    tmppos.GetXYZ(ParArr[iring]);
    
    //    ParArr[iring][4] = acos(fitqunmr_.fqmrdir[iFit][iring][2]);
    //    ParArr[iring][5] = atan2(fitqunmr_.fqmrdir[iFit][iring][1],fitqunmr_.fqmrdir[iFit][iring][0]);
    TVector3 tmpdir = TfmDetToGlblCoord(TVector3(fitqunmr_.fqmrdir[iFit][iring]),false,true);
    ParArr[iring][4] = tmpdir.Theta();
    ParArr[iring][5] = tmpdir.Phi();
  }
  
  return iFit;
  
}

int fiTQun_shared::MoveMRFit(int idFit, int iDest){// iDest<0 will remove the result
  
  if (iDest>=fitqunmr_.fqnmrfit) {
    std::cout << "In fiTQun_shared::MoveMRFit: Invalid destination iDest = " << iDest << std::endl;
    return -1;
  }
  
  bool flgRemove=false;
  if (iDest<0) {
    iDest=fitqunmr_.fqnmrfit-1;
    flgRemove=true;
  }
  
  int bufnring, bufTypeArr[FQMAXNRING], bufPCflg;
  double bufParArr[FQMAXNRING][fiTQun_shared::nSnglTrkParams], buftotmu, bufnlL;
  
  int iOrig = GetMRFit(idFit,bufnring,bufTypeArr,bufParArr,buftotmu,bufnlL,bufPCflg);
  if (iOrig<0) {
    std::cout << "In fiTQun_shared::MoveMRFit: Cannot find result. idFit = " << idFit << std::endl;
    return -1;
  }
  
  int iShift=0;
  if (iDest<iOrig) {
    iShift=1;
  }
  else if (iDest>iOrig) {
    iShift=-1;
  }
  
  if (iShift) {
    for (int iFit=iOrig-iShift; iShift*(iFit-iDest)>=0; iFit-=iShift) {
      int tmpnring, tmpTypeArr[FQMAXNRING], tmpPCflg;
      double tmpParArr[FQMAXNRING][fiTQun_shared::nSnglTrkParams], tmptotmu, tmpnlL;
      
      int idtmp=fitqunmr_.fqmrifit[iFit];
      GetMRFit(idtmp,tmpnring,tmpTypeArr,tmpParArr,tmptotmu,tmpnlL,tmpPCflg,iFit);
      SaveMRFit(idtmp,tmpnring,tmpTypeArr,tmpParArr,tmptotmu,tmpnlL,tmpPCflg,iFit+iShift);
    }
    
    SaveMRFit(idFit,bufnring,bufTypeArr,bufParArr,buftotmu,bufnlL,bufPCflg,iDest);
  }
  
  if (flgRemove) {
    SaveMRFit(0,0,NULL,NULL,0.,0.,0,iDest);
    fitqunmr_.fqnmrfit--;
  }
  
  return 0;
  
}

void fiTQun_shared::GetListMRFit(int &nMRFit, int *idList, double *nllList, int *iSort){
  nMRFit=fitqunmr_.fqnmrfit;
  for (int iFit=0; iFit<nMRFit; iFit++) {
    idList[iFit]=fitqunmr_.fqmrifit[iFit];
    nllList[iFit]=fitqunmr_.fqmrnll[iFit];
    
    iSort[iFit]=iFit;
    for (int jFit=iFit-1; jFit>=0; jFit--) {
      if (nllList[iSort[jFit]]>nllList[iSort[jFit+1]]) {
        int iBuf=iSort[jFit];
        iSort[jFit]=iSort[jFit+1];
        iSort[jFit+1]=iBuf;
      }
      else break;
    }
  }
}

void fiTQun_shared::SaveDefaultPDK_MuGammaFit(double* fitparams, int index, double totmu, double nll, int PCflg){
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitqunpmg_.fqpmgpos1[index]);

  fitqunpmg_.fqpmgt01[index] = fitparams[3];

  TVector3 tmppos2 = TfmDetToGlblCoord(TVector3(&(fitparams[7])),true,false);
  tmppos2.GetXYZ(fitqunpmg_.fqpmgpos2[index]);

  fitqunpmg_.fqpmgt02[index] = fitparams[10];

  TVector3 orgdir1;
  orgdir1.SetMagThetaPhi(1.,fitparams[4],fitparams[5]);
  TVector3 tmpdir1 = TfmDetToGlblCoord(orgdir1,false,false);
  tmpdir1.GetXYZ(fitqunpmg_.fqpmgdir1[index]);

  fitqunpmg_.fqpmgmom1[index] = fitparams[6];

  TVector3 orgdir2;
  orgdir2.SetMagThetaPhi(1.,fitparams[11],fitparams[12]);
  TVector3 tmpdir2 = TfmDetToGlblCoord(orgdir2,false,false);
  tmpdir2.GetXYZ(fitqunpmg_.fqpmgdir2[index]);

  fitqunpmg_.fqpmgmom2[index] = fitparams[13];

  fitqunpmg_.fqpmgtotmu[index] = totmu;
  fitqunpmg_.fqpmgnll[index] = nll;
  fitqunpmg_.fqpmgpcflg[index] = PCflg;
}
void fiTQun_shared::LoadLookupTableTimePDF() {
  std::cout << "loading Lookup Table Time PDFs...";

  TString LUTtPDF_Config = strSKVer;

  // If we're using WCSim AND the WCSim detector configuration has been specified
  // in the parameters file load the appropriate PMT charge response. Otherwise
  // stick to default behaviour which is to determine the SK version automatically
  // from SKDETSIM input or to use SK1 for WCSim.
#ifdef NOSKLIBRARIES
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ){
    WCSim_Config = TRuntimeParameters::Get().GetParameterS("fiTQun.WCSimConfig");
    LUTtPDF_Config = Form("_%s", WCSim_Config.Data());
  }
#endif

  TString object_name = "hist_tpdf";
  TString fq_root = "";

  // make sure we can find the files
  if (getenv("FITQUN_ROOT")==NULL)
    {
      printf("ERROR! FITQUN_ROOT is not set. See fiTQun::LoadLookupTableTimePDF()\n");
      exit(-1);
    } else {
    fq_root = getenv("FITQUN_ROOT");
  }

  for (int ipid = 0; ipid < nPID; ++ipid) {
    if (!fRngTyp[ipid]) continue;
    // direct time PDF lookup table
    TString LUTtpdf_dir_name = Form("%s/%d_LUTtpdfhist_dir%s.root",fQconstdir.Data(),PIDarr[ipid], LUTtPDF_Config.Data());
    TFile *f = new TFile(LUTtpdf_dir_name.Data(), "OPEN");
    TH3F *h = (TH3F*)f->Get(object_name.Data());
    if (!h) { std::cerr << "could not open object " << object_name << " from " << LUTtpdf_dir_name << std::endl; exit(-1);}
    tpdf_direct[ipid] = new hemi::Table3D<float>(h->GetNbinsX(), h->GetNbinsY(), h->GetNbinsZ(),
					      // low edges
					      h->GetXaxis()->GetBinCenter(1),
					      h->GetYaxis()->GetBinCenter(1),
					      h->GetZaxis()->GetBinCenter(1),
					      // upper edges
					      h->GetXaxis()->GetBinCenter(h->GetNbinsX()),
					      h->GetYaxis()->GetBinCenter(h->GetNbinsY()),
					      h->GetZaxis()->GetBinCenter(h->GetNbinsZ()));

    for (int z = 1; z <= h->GetNbinsZ(); ++z)
      for (int y = 1; y <= h->GetNbinsY(); ++y)
	for (int x = 1; x <= h->GetNbinsX(); ++x) {
	  int idx = hemi::index(x-1, y-1, z-1, h->GetNbinsX(), h->GetNbinsY());
	  tpdf_direct[ipid]->writeOnlyHostPtr()[idx] = h->GetBinContent(x, y, z);
	}
    delete f;

    // indirect time PDF lookup table
    TString LUTtpdf_sct_name = Form("%s/%d_LUTtpdfhist_sct%s.root",fQconstdir.Data(),PIDarr[ipid], LUTtPDF_Config.Data());
    f = new TFile(LUTtpdf_sct_name.Data(), "OPEN");
    h = (TH3F*)f->Get(object_name.Data());
    if (!h) { std::cerr << "could not open object " << object_name << " from " << LUTtpdf_sct_name << std::endl; exit(-1);}
    tpdf_indirect[ipid] = new hemi::Table3D<float>(h->GetNbinsX(), h->GetNbinsY(), h->GetNbinsZ(),
						   // low edges                         
						   h->GetXaxis()->GetBinCenter(1),
						   h->GetYaxis()->GetBinCenter(1),
						   h->GetZaxis()->GetBinCenter(1),
						   // upper edges           
						   h->GetXaxis()->GetBinCenter(h->GetNbinsX()),
						   h->GetYaxis()->GetBinCenter(h->GetNbinsY()),
						   h->GetZaxis()->GetBinCenter(h->GetNbinsZ()));

    for (int z = 1; z <= h->GetNbinsZ(); ++z)
      for (int y = 1; y <= h->GetNbinsY(); ++y)
        for (int x = 1; x <= h->GetNbinsX(); ++x) {
          int idx = hemi::index(x-1, y-1, z-1, h->GetNbinsX(), h->GetNbinsY());
          tpdf_indirect[ipid]->writeOnlyHostPtr()[idx] = h->GetBinContent(x, y, z);
        }
    
    h = NULL; // dont let this pointer dangle
    delete f;
  }

  std::cout << " done" << std::endl;
}


void fiTQun_shared::LoadProfiles(bool fFitCProf, bool fFittProf){

  TString cPDF_FileName = Form("%s/%s%s.root",fQconstdir.Data(),"cPDFpar",strSKVer.Data());

  // If we're using WCSim AND the PMT type has been specified in the parameters file
  // load the appropriate PMT charge response. Otherwise stick to default behaviour
  // which is to determine the SK version automatically from SKDETSIM input or to
  // use SK1 for WCSim.
#ifdef NOSKLIBRARIES
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimPMTType") ){
    WCSim_PMTType = TRuntimeParameters::Get().GetParameterS("fiTQun.WCSimPMTType");
    cPDF_FileName = Form("%s/%s_%s.root",fQconstdir.Data(),"cPDFpar", WCSim_PMTType.Data());
  }
#endif

  char * LoadParamsArg = const_cast<char*> (cPDF_FileName.Data()); // Cast to non-const
  fQChrgPDF::Get()->LoadParams(LoadParamsArg);//read charge PDF parameters
  
  fUseFitCprof=fFitCProf;
  R0bins=NULL;
  th0bins=NULL;
  
  //pdg range and time//
  TFile *fpdg   = new TFile(fQconstdir+"gpdgRange.root");
  TFile* fpdgt = new TFile(fQconstdir+"gpdgTime.root");
  if (fpdg->IsZombie() || fpdgt->IsZombie()) {
    std::cout << "PDG range table missing!" << std::endl;
    exit(-1);
  }
  gPDGRange = new TGraph(*(TGraph*)fpdg->Get("gpdgRange"));
  std::cout<<gPDGRange->Eval(1000.)<<std::endl;
  gPDGTime = new TGraph(*(TGraph*)fpdgt->Get("gpdgTime"));
  delete fpdg;
  delete fpdgt;
  ////////////
  
  //Initialize Cherenkov profiles
  for (int iPID=0; iPID<nPID; iPID++) {//loop over particle types
    if (!fRngTyp[iPID]) continue;
    
    //Initiatize integral table In
    if (fUseFitCprof) {//Load fit Cherenkov profile
      LoadCProf(iPID);
    }
    else {
      std::cout << "Loading Cherenkov profile: "  << PIDnames[iPID] << std::endl;

      TString cProfFileName = Form(fQconstdir+"CProf_%d.root",PIDarr[iPID]);
      
      if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ) {
        cProfFileName = Form(fQconstdir+"CProf_%d_WCSim.root",PIDarr[iPID]);
      }
      
      TFile *fcprof = new TFile(cProfFileName);
      if (fcprof->IsZombie()) {
        std::cout << "Cherenkov profile is missing!" << std::endl;
        delete fcprof;
        DeleteInstances();
        exit(-1);
      }
      
      for (int i=0; i<3; i++) {//I_n
        TH3F *hI3d = (TH3F*)fcprof->Get(Form("hI3d_%d",i));
        if (i==0) {//Set integral table bins
          if (R0bins==NULL) {//R0, costh0 bins same for all particle types
            nR0bin=hI3d->GetNbinsX();
            nth0bin=hI3d->GetNbinsY();
            R0bins = new double[nR0bin];
            th0bins = new double[nth0bin];
            for (int j=0; j<nR0bin; j++) R0bins[j]=hI3d->GetXaxis()->GetBinLowEdge(j+1);
            for (int j=0; j<nth0bin; j++) th0bins[j]=hI3d->GetYaxis()->GetBinLowEdge(j+1);
            R0binwrecp=1./(R0bins[1]-R0bins[0]);
            th0binwrecp=1./(th0bins[1]-th0bins[0]);
          }
          
          nmom[iPID]=hI3d->GetNbinsZ();
          mombins[iPID] = new double[nmom[iPID]];
          for (int j=0; j<nmom[iPID]; j++) mombins[iPID][j]=hI3d->GetZaxis()->GetBinLowEdge(j+1);
        }
        InTbl[iPID][i] = new float**[nR0bin];
        for (int j=0; j<nR0bin; j++) {
          InTbl[iPID][i][j] = new float*[nth0bin];
          for (int k=0; k<nth0bin; k++) {
            InTbl[iPID][i][j][k] = new float[nmom[iPID]];
            for (int l=0; l<nmom[iPID]; l++) {
              InTbl[iPID][i][j][k][l]=hI3d->GetBinContent(j+1,k+1,l+1);
            }
          }
        }
      }
      
      for (int i=0; i<2; i++) {//isotropic source integrale table
        TH1D *hIiso = (TH1D*)fcprof->Get(Form("hI_iso_%d",i+1));
        int nmomtmp=hIiso->GetNbinsX();
        gIniso[iPID][i] = new TGraph(nmomtmp);
        for (int j=0; j<nmomtmp; j++) {
          gIniso[iPID][i]->SetPoint(j,hIiso->GetXaxis()->GetBinLowEdge(j+1),hIiso->GetBinContent(j+1));
        }
      }
      gNphot[iPID] = new TGraph(*(TGraph*)fcprof->Get("gNphot"));
      int ngpts=gNphot[iPID]->GetN();
      gNphtmom[iPID] = new TGraph(ngpts);
      for (int i=0; i<ngpts; i++) {
        Double_t gx,gy;
        gNphot[iPID]->GetPoint(i,gx,gy);
        gNphtmom[iPID]->SetPoint(i,gy,gx);
      }
      
      gSmax[iPID] = new TGraph(*(TGraph*)fcprof->Get("gsthr"));
      
      delete fcprof;
    }
  }
  
  fUseFittpdf=fFittProf;
  //Initialize time likelihood parameters

  TString tPDF_Config = strSKVer;

  // If we're using WCSim AND the WCSim detector configuration has been specified
  // in the parameters file load the appropriate PMT charge response. Otherwise
  // stick to default behaviour which is to determine the SK version automatically
  // from SKDETSIM input or to use SK1 for WCSim.
#ifdef NOSKLIBRARIES
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ){
    WCSim_Config = TRuntimeParameters::Get().GetParameterS("fiTQun.WCSimConfig");
    tPDF_Config = Form("_%s", WCSim_Config.Data());
  }
#endif
  
  for (int iPID=0; iPID<nPID; iPID++) {//loop over particle types
    if (!fRngTyp[iPID]) continue;
    TFile *ftprof = new TFile(Form("%s/%d_tpdfpar%s.root",fQconstdir.Data(),PIDarr[iPID],tPDF_Config.Data()));
    if (ftprof->IsZombie()) {
      std::cout << "Time likelihood parameter file missing!" << std::endl;
      delete ftprof;
      DeleteInstances();
      exit(-1);
    }
    
    TH1D *htpdfinfo = (TH1D*)ftprof->Get("htpdfinfo");
    int ntmp=(int)(htpdfinfo->GetBinContent(1)+0.5);
    tpdfmurang[0]=htpdfinfo->GetBinContent(2);
    tpdfmurang[1]=htpdfinfo->GetBinContent(3);
    int ntmp2=(int)(htpdfinfo->GetBinContent(4)+0.5);
    tpdfprang[iPID][0]=htpdfinfo->GetBinContent(5);
    tpdfprang[iPID][1]=htpdfinfo->GetBinContent(6);
    
    if (ntmp!=ntpdfmupar || ntmp2!=ntpdfppar) {
      std::cout << "Incompatible time likelihood parameter file!" << std::endl;
      delete ftprof;
      DeleteInstances();
      exit(-1);
    }
    
    delete htpdfinfo;
    
    if (fUseFittpdf) {
      std::cout << "Loading fit time profile parameters: "  << PIDnames[iPID] << std::endl;
      
      TH2D *htpdfparmn = (TH2D*)ftprof->Get("htpdfparmn");
      TH2D *htpdfparsg = (TH2D*)ftprof->Get("htpdfparsg");
      
      for (int imupar=0; imupar<ntpdfmupar; imupar++) {
        for (int ippar=0; ippar<ntpdfppar; ippar++) {
          tpdfpar[iPID][0][imupar][ippar]=htpdfparmn->GetBinContent(imupar+1,ippar+1);
          tpdfpar[iPID][1][imupar][ippar]=htpdfparsg->GetBinContent(imupar+1,ippar+1);
        }
      }
    }
    else {
      std::cout << "Loading time profile parameters: "  << PIDnames[iPID] << std::endl;
      for (int i=0; i<ntpdfmupar; i++) {
        gtpar[iPID][0][i] = new TGraphErrors(*(TGraphErrors*)ftprof->Get(Form("gtcmnpar_%d",i)));
        gtpar[iPID][1][i] = new TGraphErrors(*(TGraphErrors*)ftprof->Get(Form("gtcsgpar_%d",i)));
      }
    }
    
    delete ftprof;
    
    std::cout << "Time likelihood parameter " << PIDnames[iPID] << ": pmin=" << tpdfprang[iPID][0] << ", pmax=" << tpdfprang[iPID][1] << std::endl;
  }
  
}

void fiTQun_shared::LoadCProf(int iPID) { //Initialize Cherenkov profile parameters
  
  std::cout << "Loading fit Cherenkov profile parameters: "  << PIDnames[iPID] << std::endl;

  TString cProfFileName = Form(fQconstdir+"CProf_%d_fit.root",PIDarr[iPID]);
  
  if ( TRuntimeParameters::Get().HasParameter("fiTQun.WCSimConfig") ) {
    cProfFileName = Form(fQconstdir+"CProf_%d_fit_WCSim.root",PIDarr[iPID]);
  }
  
  TFile *fcprof = new TFile(cProfFileName);
  if (fcprof->IsZombie()) {
    std::cout << "Cherenkov profile is missing!" << std::endl;
    delete fcprof;
    DeleteInstances();
    exit(-1);
  }
  
  TH1D *hprofinf = (TH1D*)fcprof->Get("hprofinf");
  int npars=(int)(hprofinf->GetBinContent(1)+0.5);
  int nsectmax=(int)(hprofinf->GetBinContent(2)+0.5);
  momofst[iPID]=hprofinf->GetBinContent(3);
  momminstp[iPID]=hprofinf->GetBinContent(4);
  momminmax[iPID][0]=hprofinf->GetBinContent(5);
  momminmax[iPID][1]=hprofinf->GetBinContent(6);
  delete hprofinf;
  
  if (npars!=5) {
    std::cout << "npars!=5" << std::endl;
    delete fcprof;
    DeleteInstances();
    exit(-1);
  }
  
  //Initiatize integral table In
  for (int i=0; i<3; i++) {//I_n
    TH3D *hI3d = (TH3D*)fcprof->Get(Form("hI3d_par_%d",i));
    TH2D *hInsect = (TH2D*)fcprof->Get(Form("hI3d_nsect_%d",i));
    if (i==0) {//Set integral table bins
      if (R0bins==NULL) {//R0, costh0 bins same for all particle types
        nR0bin=hI3d->GetNbinsX();
        nth0bin=hI3d->GetNbinsY();
        R0bins = new double[nR0bin];
        th0bins = new double[nth0bin];
        for (int j=0; j<nR0bin; j++) R0bins[j]=hI3d->GetXaxis()->GetBinLowEdge(j+1);
        for (int j=0; j<nth0bin; j++) th0bins[j]=hI3d->GetYaxis()->GetBinLowEdge(j+1);
        R0binwrecp=1./(R0bins[1]-R0bins[0]);
        th0binwrecp=1./(th0bins[1]-th0bins[0]);
      }
      nmom[iPID]=(int)((momminmax[iPID][1]-momminmax[iPID][0])/momminstp[iPID]+0.5);//number of momentum bins
      mombins[iPID] = new double[nmom[iPID]+1];
      for (int j=0; j<=nmom[iPID]; j++) mombins[iPID][j]=momminmax[iPID][0]+j*momminstp[iPID];
      
      int nInParPtrs=nR0bin*nth0bin*nmom[iPID];
      InParPtr[iPID] = new double*[nInParPtrs*3];
    }
    int nsecttot=(int)(hInsect->Integral()+0.5);//total number of sections
    InPar[iPID][i] = new double[nsecttot*npars];//array which contains polynomial coefficients
    int ictr=0;
    for (int j=0; j<nR0bin; j++) {
      for (int k=0; k<nth0bin; k++) {
        int nsect=(int)(hInsect->GetBinContent(j+1,k+1)+0.5);//number of sections in this R0,th0 bin
        double *params = new double[(npars+1)*nsect+1];
        for (int ipar=0; ipar<=(npars+1)*nsect; ipar++) {
          params[ipar]=hI3d->GetBinContent(j+1,k+1,ipar+1);
        }
        int l=0;
        for (int isect=0; isect<nsect; isect++) {
          for (int ipar=0; ipar<npars; ipar++) {
            InPar[iPID][i][ictr+ipar]=params[isect*npars+ipar];
          }
          while (mombins[iPID][l]+1e-6<params[nsect*npars+isect+1]) {//momentum bin is in isect
            InParPtr[iPID][((l*3+i)*nR0bin+j)*nth0bin+k]=&InPar[iPID][i][ictr];
            l++;
          }
          ictr+=npars;
        }
        delete params;
      }
    }
    delete hInsect;
    delete hI3d;
  }
  
  for (int i=0; i<2; i++) {//isotropic source integrale table
    TH1D *hIiso = (TH1D*)fcprof->Get(Form("hI_iso_%d",i+1));
    int nmomtmp=hIiso->GetNbinsX();
    gIniso[iPID][i] = new TGraph(nmomtmp);
    for (int j=0; j<nmomtmp; j++) {
      gIniso[iPID][i]->SetPoint(j,hIiso->GetXaxis()->GetBinLowEdge(j+1),hIiso->GetBinContent(j+1));
    }
  }
  gNphot[iPID] = new TGraph(*(TGraph*)fcprof->Get("gNphot"));
  int ngpts=gNphot[iPID]->GetN();
  gNphtmom[iPID] = new TGraph(ngpts);
  for (int i=0; i<ngpts; i++) {
    Double_t gx,gy;
    gNphot[iPID]->GetPoint(i,gx,gy);
    gNphtmom[iPID]->SetPoint(i,gy,gx);
  }
  
  gSmax[iPID] = new TGraph(*(TGraph*)fcprof->Get("gsthr"));
  
  Load1dpars(fcprof,"gNphot_pars",(npars1dcprof+1)*nsect1dcprof+1,nphotpar[iPID]);
  Load1dpars(fcprof,"gsthr_pars",(npars1dcprof+1)*nsect1dcprof+1,smaxpar[iPID]);
  Load1dpars(fcprof,"gI_iso_1_pars",(npars1dcprof+1)*nsect1dcprof+1,Inisopar[iPID][0]);
  Load1dpars(fcprof,"gI_iso_2_pars",(npars1dcprof+1)*nsect1dcprof+1,Inisopar[iPID][1]);
  
  delete fcprof;
  
}

void fiTQun_shared::Load1dpars(TFile *fin, const char *ctmp, int npartmp, double *par){
  TH1D *hparam = (TH1D*)fin->Get(ctmp);
  for (int ipar=0; ipar<npartmp; ipar++) {
    par[ipar]=hparam->GetBinContent(ipar+1);
  }
}

double fiTQun_shared::fconpoly(double xval, int npars, int nsect, double *par) {
  int ibase=npars*nsect;//index of the first boundary
  
  int isect;
  for (isect=0; isect<nsect; isect++) {
    if (xval<par[ibase+isect+1] || isect==nsect-1) {
      break;
    }
  }
  
  double xpow=1.;
  double tmpval=0.;
  for (int ipar=0; ipar<npars; ipar++) {
    tmpval+=par[isect*npars+ipar]*xpow;
    xpow*=xval;
  }
  
  return tmpval;
}

void fiTQun_shared::SetTypemom(int iPID, double mom, double *In, double& nphot, double& smax) {
  
  if (!fRngTyp[iPID]) {
    std::cout << "Particle profile not loaded: PDG=" << fiTQun_shared::PIDarr[iPID] << std::endl;
    exit(-1);
  }
  
  if (fUseFitCprof) {
    //bound xtmp!!!
    double xtmp=mom-momofst[iPID];
    if (!(xtmp>1.)) xtmp=1.;
    xtmp=log(xtmp);
    
    if (In) {
      int idx=(int)((xtmp-momminmax[iPID][0])/momminstp[iPID]);
      if (idx<0) idx=0;
      if (idx>nmom[iPID]-1) idx=nmom[iPID]-1;
      InPtrCur = &InParPtr[iPID][idx*3*nR0bin*nth0bin];
      for (int i=0; i<5; i++) {
        lnmom[i]=pow(xtmp,i);
      }
      
      In[0]=exp(fconpoly(xtmp,npars1dcprof,nsect1dcprof,Inisopar[iPID][0]));
      In[1]=exp(fconpoly(xtmp,npars1dcprof,nsect1dcprof,Inisopar[iPID][1]));
    }
    
    nphot=exp(fconpoly(xtmp,npars1dcprof,nsect1dcprof,nphotpar[iPID]));
    smax=exp(fconpoly(xtmp,npars1dcprof,nsect1dcprof,smaxpar[iPID]));
  }
  else {
    for (imomcur=1; imomcur<nmom[iPID]; imomcur++) {
      if (mom<mombins[iPID][imomcur]) break;
    }
    if (imomcur==nmom[iPID]) imomcur--;
    imomcur--;
    
    mombfrac=(mom-mombins[iPID][imomcur])/(mombins[iPID][imomcur+1]-mombins[iPID][imomcur]);
    for (int i=0; i<3; i++) {
      InTblcur[i]=InTbl[iPID][i];
    }
    
    nphot=gNphot[iPID]->Eval(mom);
    smax=gSmax[iPID]->Eval(mom);
    
    In[0] = gIniso[iPID][0]->Eval(mom);
    In[1] = gIniso[iPID][1]->Eval(mom);
  }
  
}

double fiTQun_shared::EvalPdgRange(double mom){
  return gPDGRange->Eval(mom);
}

double fiTQun_shared::EvalPdgTime(double mom){
  return gPDGTime->Eval(mom);
}

void fiTQun_shared::EvalIn(double *In, double R0, double costh0) {//This must not be called before SetTypemom
  
  if (!fUseFitCprof) {
    InterpolateIn(In,R0,costh0);
    return;
  }
  
  //  double i1,i2,j1,j2,w1,w2;
  
  double dR0=(R0-R0bins[0])*R0binwrecp;
  int iR0=(int)dR0;
  double dth0=(costh0-th0bins[0])*th0binwrecp;
  int ith0=(int)dth0;
  
  if (iR0 > nR0bin-2) {
    iR0 = nR0bin-2;
  }
  else if (iR0 < 0) {
    iR0 = 0;
  }
  
  if (ith0 > nth0bin-2) {
    ith0 = nth0bin-2;
  }
  else if (ith0 < 0) {
    ith0 = 0;
  }
  
  double yd=dth0-(double)ith0;
  double xd=dR0-(double)iR0;
  
  int ibase=iR0*nth0bin+ith0;
  for (int i=0; i<3; i++) {
    double i1 = InPtrCur[ibase][0]*lnmom[0]+InPtrCur[ibase][1]*lnmom[1]+InPtrCur[ibase][2]*lnmom[2]+InPtrCur[ibase][3]*lnmom[3]+InPtrCur[ibase][4]*lnmom[4];
    double i2 = InPtrCur[ibase+1][0]*lnmom[0]+InPtrCur[ibase+1][1]*lnmom[1]+InPtrCur[ibase+1][2]*lnmom[2]+InPtrCur[ibase+1][3]*lnmom[3]+InPtrCur[ibase+1][4]*lnmom[4];
    int ibase2=ibase+nth0bin;
    double j1 = InPtrCur[ibase2][0]*lnmom[0]+InPtrCur[ibase2][1]*lnmom[1]+InPtrCur[ibase2][2]*lnmom[2]+InPtrCur[ibase2][3]*lnmom[3]+InPtrCur[ibase2][4]*lnmom[4];
    double j2 = InPtrCur[ibase2+1][0]*lnmom[0]+InPtrCur[ibase2+1][1]*lnmom[1]+InPtrCur[ibase2+1][2]*lnmom[2]+InPtrCur[ibase2+1][3]*lnmom[3]+InPtrCur[ibase2+1][4]*lnmom[4];
    
    double w1=i1*(1.-yd)+i2*yd;
    double w2=j1*(1.-yd)+j2*yd;
    
    In[i]=w1*(1.-xd)+w2*xd;
    
    ibase+=nR0bin*nth0bin;
  }
}

void fiTQun_shared::InterpolateIn(double *In, double R0, double costh0) {//This must not be called before SetTypemom
  double i1,i2,j1,j2,w1,w2;
  
  int iR0=(int)((R0-R0bins[0])*R0binwrecp);
  int ipR0=iR0+1;
  int ith0=(int)((costh0-th0bins[0])*th0binwrecp);
  int ipth0=ith0+1;
  
  int ipmom=imomcur+1;
  
  if (iR0 >= nR0bin-1) {
    ipR0 = nR0bin - 1;
    iR0 = ipR0 - 1;
  }
  else if (iR0 < 0) {
    iR0 = 0;
    ipR0 = 1;
  }
  
  if (ith0 >= nth0bin-1) {
    ipth0 = nth0bin - 1;
    ith0 = ipth0 - 1;
  }
  else if (ith0 < 0) {
    ith0 = 0;
    ipth0 = 1;
  }
  
  double xd=(R0-R0bins[iR0])*R0binwrecp;
  double yd=(costh0-th0bins[ith0])*th0binwrecp;
  double zd=mombfrac;
  
  for (int i=0; i<3; i++) {
    i1 = InTblcur[i][iR0][ith0][imomcur]*(1.-zd)+InTblcur[i][iR0][ith0][ipmom]*zd;
    i2 = InTblcur[i][iR0][ipth0][imomcur]*(1.-zd)+InTblcur[i][iR0][ipth0][ipmom]*zd;
    j1 = InTblcur[i][ipR0][ith0][imomcur]*(1.-zd)+InTblcur[i][ipR0][ith0][ipmom]*zd;
    j2 = InTblcur[i][ipR0][ipth0][imomcur]*(1.-zd)+InTblcur[i][ipR0][ipth0][ipmom]*zd;
    
    w1=i1*(1.-yd)+i2*yd;
    w2=j1*(1.-yd)+j2*yd;
    
    In[i]=w1*(1.-xd)+w2*xd;
  }
}

void fiTQun_shared::GetTimeCoeff(int iPID, double momtmp, double *tcmupar) {
  if (!(momtmp>=tpdfprang[iPID][0] && momtmp<tpdfprang[iPID][1])) {
    if (momtmp<tpdfprang[iPID][0]) momtmp=tpdfprang[iPID][0];
    else momtmp=tpdfprang[iPID][1];
  }
  
  for (int imupar=0; imupar<ntpdfmupar; imupar++) {
    if (fUseFittpdf) {
      tcmupar[2*imupar]=0.;
      tcmupar[2*imupar+1]=0.;
      double ppow=1.;
      for (int ippar=0; ippar<ntpdfppar; ippar++) {
        tcmupar[2*imupar]+=tpdfpar[iPID][0][imupar][ippar]*ppow;
        tcmupar[2*imupar+1]+=tpdfpar[iPID][1][imupar][ippar]*ppow;
        ppow*=momtmp;
      }
    }
    else {
      tcmupar[2*imupar]=gtpar[iPID][0][imupar]->Eval(momtmp);
      tcmupar[2*imupar+1]=gtpar[iPID][1][imupar]->Eval(momtmp);
    }
  }
}

void fiTQun_shared::DeletePtclProf(int iPID){
  
  if (!fRngTyp[iPID]) return;
  
  delete[] mombins[iPID];
  
  if (fUseFitCprof) {
    delete[] InParPtr[iPID];
    for (int i=0; i<3; i++) {
      delete[] InPar[iPID][i];
    }
  }
  else {
    for (int i=0; i<3; i++) {
      for (int j=0; j<nR0bin; j++) {
        for (int k=0; k<nth0bin; k++) {
          delete[] InTbl[iPID][i][j][k];
        }
        delete[] InTbl[iPID][i][j];
      }
      delete[] InTbl[iPID][i];
    }
  }
  
  delete gIniso[iPID][0];
  delete gIniso[iPID][1];
  
  delete gNphot[iPID];
  delete gSmax[iPID];
  delete gNphtmom[iPID];
  
  if (!fUseFittpdf) {
    for (int i=0; i<ntpdfmupar; i++) {
      delete gtpar[iPID][0][i];
      delete gtpar[iPID][1][i];
    }
  }
  
  fRngTyp[iPID]=0;
  
}

void fiTQun_shared::DeleteProfiles(){
  
  delete[] R0bins;
  delete[] th0bins;
  for (int iPID=0; iPID<nPID; iPID++) DeletePtclProf(iPID);
  
}

void fiTQun_shared::SaveMSFit(int idMRFit, //id of multi-ring fit seed
                              int nseg, //# of segments om fot
                              int* segIndx,  //index of each segmentin ParArr
                              int iFit,  //id of ms fit
                              int *TypeArr, //array of particle types
                              double ParArr[][fiTQun_shared::nSnglTrkParams], //array of track params
                              double totmu, //total charge
                              double nll, //negative log likelihood
                              int PCflg){
  

  fitqunms_.fqmsnfit=iFit+1;  //# of ms fit

  fitqunms_.fqmspcflg[iFit] = PCflg; 
  fitqunms_.fqmsnll[iFit] = nll;
  fitqunms_.fqmstotmu[iFit] = totmu;
  fitqunms_.fqmsnseg[iFit]  = nseg;
  fitqunms_.fqmsifit[iFit]  = idMRFit; //id of MR fit used to seed, 1 for stand-alone
  fitqunms_.fqmspid[iFit] = TypeArr[segIndx[0]];
  fitqunms_.fqmsimer[iFit] = segIndx[0];

  for (int i=0;i<nseg;i++){
    fitqunms_.fqmsmom[iFit][i] = ParArr[segIndx[i]][6];
    
    fitqunms_.fqmseloss[iFit][i] = ParArr[segIndx[i]][7];
    fitqunms_.fqmst0[iFit][i] = ParArr[segIndx[i]][3];
    
    //    fitqunms_.fqmspos[iFit][i][0] = ParArr[segIndx[i]][0];
    //    fitqunms_.fqmspos[iFit][i][1] = ParArr[segIndx[i]][1];
    //    fitqunms_.fqmspos[iFit][i][2] = ParArr[segIndx[i]][2];
    TVector3 tmppos = TfmDetToGlblCoord(TVector3(ParArr[segIndx[i]]),true,false);
    tmppos.GetXYZ(fitqunms_.fqmspos[iFit][i]);

    //    fitqunms_.fqmsdir[iFit][i][0] = sin(ParArr[segIndx[i]][4])*cos(ParArr[segIndx[i]][5]);
    //    fitqunms_.fqmsdir[iFit][i][1] = sin(ParArr[segIndx[i]][4])*sin(ParArr[segIndx[i]][5]);
    //    fitqunms_.fqmsdir[iFit][i][2] = cos(ParArr[segIndx[i]][4]);
    TVector3 orgdir;
    orgdir.SetMagThetaPhi(1.,ParArr[segIndx[i]][4],ParArr[segIndx[i]][5]);
    TVector3 tmpdir = TfmDetToGlblCoord(orgdir,false,false);
    tmpdir.GetXYZ(fitqunms_.fqmsdir[iFit][i]);
/**/
  }
 // std::cout<<"saved MS fit info in common block "<<std::endl;

}

void fiTQun_shared::SetTotqMaxAng(double ctmp) {
  
  if (ctmp<0.) {
    ctmp=-ctmp;
    fmskbkPMT=true;//not use PMTs in backward direction for 1R fit
    std::cout << "PMTs in track backward direction are not used for 1R fit." << std::endl;
  }
  
  if (ctmp>180.) {
    cosmintotq=-1.;
  }
  else {
    cosmintotq=cos(ctmp/180.*pi);
  }
  
}

void fiTQun_shared::SetActivePMTs(int iflgMode){
  double R=this->DetGeomR;
  double h=this->DetGeomL;
  double x,y,z;
  double xh,yh;
  TFile *fDrPMT;
  TString ofnmtmp;
  
  int PhotCovFrac=iflgMode%10;
  int TBflag=(iflgMode-PhotCovFrac)/10;
  
  if (PhotCovFrac!=1 && PhotCovFrac!=2 && PhotCovFrac!=4) return;
  if (!(TBflag>=0 && TBflag<4)) return;
  
  if (ofnmflg==1) {
    ofnmtmp=Form("%s_PMTinfo.root",ofnmstr);
  }
  else {
    ofnmtmp="PMTinfo.root";
  }
  
  chrgfact=(double)PhotCovFrac;
  std::cout << Form("Reduce photo coverage to 1/%d",PhotCovFrac) << std::endl;
  std::cout << Form("TBflag = %d",TBflag) << std::endl;
  
  //  TH2D hPMTs(Form("hPMTinfo_%d",ievent),"",150,-(pi*R)-35.,(pi*R)-35.,147,-(h+2.*R),(h+2.*R));
  TH2D hPMTs(Form("hPMTinfo"),"",150,-(pi*R)-35.,(pi*R)-35.,147,-(h+2.*R),(h+2.*R));
  for (int i=0; i<fMaxNPMT; i++) {
    x = this->PMTpos[i][0];
    y = this->PMTpos[i][1];
    z = this->PMTpos[i][2];
    
    if (z>=h) {//top
      if (TBflag==1 || TBflag==3) this->PMTtype[i]=-1;
      xh=-y-35.;
      yh=-x+h+R;
    }
    else if (z<=-h) {//bottom
      if (TBflag==2 || TBflag==3) this->PMTtype[i]=-1;
      xh=-y-35.;
      yh=x-h-R;
    }
    else {//side
      xh=R*acos(x/R);
      if (y>0) xh=-xh;
      yh=z;
    }
    
    int ixPMT=hPMTs.GetXaxis()->FindBin(xh);
    int iyPMT=hPMTs.GetYaxis()->FindBin(yh);
    
    if (PhotCovFrac==2) {// 1/2 PMT
      if ((ixPMT+iyPMT)%2) this->PMTtype[i]=-1;
    }
    else if (PhotCovFrac==4) {// 1/4 PMT
      if (ixPMT%2 || iyPMT%2) this->PMTtype[i]=-1;
    }
    
    hPMTs.Fill(xh,yh,this->PMTtype[i]);
  }
  
  fDrPMT = new TFile(ofnmtmp,"RECREATE");
  
  hPMTs.SetOption("COLZ");
  hPMTs.SetMinimum(0.);
  
  hPMTs.Write();
  
  delete fDrPMT;
  
}

void fiTQun_shared::Setofilenm(char *ofname){
  for (int i=0; i<256; i++) {
    ofnmstr[i]=ofname[i];
    fQEventManager::ofname[i]=ofname[i];
    if (ofname[i]=='\0') {
      break;
    }
  }
  ofnmflg=1;
}

void fiTQun_shared::SetPhi0(int iPhi0, double Phi0val){
  if (iPhi0 >= 0) {
    Phi0[iPhi0] = Phi0val;
    return;
  }
  for (int i=0; i<nPID; i++) {
    Phi0[i]=Phi0val;
  }
}

void fiTQun_shared::SetdEdx(int iPID, double val){
  dEdx[iPID]=val;
  std::cout << Form("%s setting dEdx=%f",fiTQun_shared::PIDnames[iPID].Data(),dEdx[iPID]) << std::endl;
}

int fiTQun_shared::ClrInfo(int iMod) {// clear all common block variables other than the subevent and pi0 cut variables
  
  int idum=0;
  if (iMod==0) {
    clrfqcmns_(&idum);
  }
  
  return 0;
}

void fiTQun_shared::SavePreFit(int itwnd, double* fitparams){
  
  //  fitquntwnd_.fqtwnd_prftpos[itwnd][0] = fitparams[0];
  //  fitquntwnd_.fqtwnd_prftpos[itwnd][1] = fitparams[1];
  //  fitquntwnd_.fqtwnd_prftpos[itwnd][2] = fitparams[2];
  TVector3 tmppos = TfmDetToGlblCoord(TVector3(fitparams),true,false);
  tmppos.GetXYZ(fitquntwnd_.fqtwnd_prftpos[itwnd]);
  
  fitquntwnd_.fqtwnd_prftt0 [itwnd]    = fitparams[3];
}