#include <iostream>
#include "PDK_MuGamma.h"
#include "TMath.h"
#include "TFitter.h"
#include <TMinuit.h>
#include "fiTQun_shared.h"
#include "mcguts.h"
#include "apscndryC.h"
#include <fstream>

using namespace fiTQun_parameters;

PDK_MuGamma* PDK_MuGamma::static_PDK_MuGamma;

extern"C" {

  void mcguts(int *ntracks);
  void trginfo_(float *);
  void vcrdvccm_();

}

PDK_MuGamma::PDK_MuGamma() {
  std::cout << "Starting PDK_MuGamma Fitter." << std::endl;
  static_PDK_MuGamma = this;
  //GetTruePDK_MuGammaParams(tParams);
}

PDK_MuGamma::~PDK_MuGamma() {
  std::cout << "PDK routine ended." << std::endl;
}

//************************************************************************************************************************************
//void PDK_MuGamma::PDK_MuGammaWrapper(int& nDim, double* gout, double& result, double par[], int flg);
//************************************************************************************************************************************
void PDK_MuGamma::PDK_MuGammaWrapper(int& nDim, double* gout, double& result, double par[], int flg){

  double partmp[fiTQun_shared::nPDK_MuGammaParams];
  for (int i=0; i<fiTQun_shared::nPDK_MuGammaParams; i++)
    partmp[i]=par[i];  

  double lnLpnl=0.;

  if (par[6]<fiTQun_shared::Get()->pCherenkovThr(imu))
    lnLpnl+=pow((fiTQun_shared::Get()->pCherenkovThr(imu)-par[6]),2.)*10.;
  
  if (par[13]<fiTQun_shared::Get()->pCherenkovThr(ie))
    lnLpnl+=pow((fiTQun_shared::Get()->pCherenkovThr(ie)-par[13]),2.)*10.;

  int PCdumflg;
  
  //lnLpnl=0.;
  result = static_PDK_MuGamma->MuGammaNglogL(partmp,PCdumflg) + lnLpnl;
}
//************************************************************************************************************************************
//double PDK_MuGamma::RunMinimizer(double *X, int& PCflg, int *flgparfix, double *arStps){
//************************************************************************************************************************************
double PDK_MuGamma::RunMinimizer(double *X, int& PCflg, int *flgparfix, double *arStps){

  int nPar = fiTQun_shared::nPDK_MuGammaParams;
  //      0   1   2   3     4      5    6   7   8   9   10   11     12    13
  // X = {X1, Y1, Z1, T1, theta1, phi1, p1, X2, Y2, Z2, T2, theta2, phi2, p2}

  const double Vtxres=50.;
  const double t0res=10.;
  const double angres=TMath::Pi()/16.;
  const double momres1=X[6]*0.1+20.;
  const double momres2=X[13]*0.1+5.;

//  const double Vtxres=500.;
//  const double t0res=20.;
//  const double angres=TMath::Pi()/8.;
//  const double momres1=X[6]*0.2+20.;
//  const double momres2=X[13]*0.1+5.;

  const int nRes = 5;
  double stpsize[nRes]={Vtxres,t0res,angres,momres1,momres2};

  if (arStps!=NULL) {
    for (int i=0; i<nRes; i++) stpsize[i]=arStps[i];
  }

  TFitter* minimizer = new TFitter(nPar);
  if ( fiTQun_shared::Get()->IsQuiet() ) {
    double p1 = -1;
    minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);
  }

  minimizer->SetFCN(PDK_MuGammaWrapper);
  minimizer->SetParameter(0, "x1",     X[0],stpsize[0],0,0);
  minimizer->SetParameter(1, "y1",     X[1],stpsize[0],0,0);
  minimizer->SetParameter(2, "z1",     X[2],stpsize[0],0,0);
  minimizer->SetParameter(3, "t1",     X[3],stpsize[1],0,0);
  minimizer->SetParameter(4, "theta1", X[4],stpsize[2],0,0);
  minimizer->SetParameter(5, "phi1",   X[5],stpsize[2],0,0);
  minimizer->SetParameter(6, "p1",     X[6],stpsize[3],0,0);
  minimizer->SetParameter(7, "x2",     X[7],stpsize[0],0,0);
  minimizer->SetParameter(8, "y2",     X[8],stpsize[0],0,0);
  minimizer->SetParameter(9, "z2",     X[9],stpsize[0],0,0);
  minimizer->SetParameter(10,"t2",     X[10],stpsize[1],0,0);
  minimizer->SetParameter(11,"theta2", X[11],stpsize[2],0,0);
  minimizer->SetParameter(12,"phi2",   X[12],stpsize[2],0,0);
  minimizer->SetParameter(13,"p2",     X[13],stpsize[4],0,0);
  
  if (flgparfix!=NULL) {
    for (int i=0; i<nPar; i++) {
      if (flgparfix[i]==1) {
	std::cout << "Fixing parameter " << i << std::endl;
        minimizer->FixParameter(i);
      }
    }
  }

  if (fiTQun_shared::Get()->GetTimeflg() == false) minimizer->FixParameter(3);

  Double_t CmdOpts[2]={3000,0.1};
  int ierflg = minimizer->ExecuteCommand("SIMPLEX",CmdOpts,1);

  //      minimizer->ExecuteCommand("IMPROVE",0,0);

  for (int i=0; i < nPar; i++) { 
    X[i]=minimizer->GetParameter(i);
  }

  if(same3Vertex)
    for(int i = 0; i < 3; i++)
      X[i+7] = X[i];

  double pCthr1=fqshared->pCherenkovThr(imu);
  if (X[6]<pCthr1/2.) {
    X[6]=pCthr1/2.;//momentum lower bound (arbitrary choice)
  }
  double pCthr2=fqshared->pCherenkovThr(ie);
  if (X[13]<pCthr2/2.) {
    X[13]=pCthr2/2.;//momentum lower bound (arbitrary choice)
  }

  int nIter = minimizer->GetMinuit()->fNfcn;

  delete minimizer;

  double tmpnglnL = MuGammaNglogL(X,PCflg);

  if (ierflg!=0) {
    PCflg=-1-1000*nIter;//didn't converge!
    std::cout << "MINUIT did not converge!!! n_Iter=" << nIter << std::endl;
    //    return 1e100;
  }

  return tmpnglnL;

}

//************************************************************************************************************************************
//double PDK_MuGamma::SetPDK_MuGammaSeed(double* fitpars, int& PCflg, bool seedpi0mass, bool translatevertex);
//************************************************************************************************************************************
double PDK_MuGamma::SetPDK_MuGammaSeed(double* fitpars, int& PCflg){

  // first 7 values of fitpars should be filled with 1-ring mu-like fit result as input

  const int nPoints = 100;

  double theta[nPoints];
  double phi[nPoints];
  double magic = 3.6;

  theta[0]    = TMath::Pi();
  phi  [0]    = 0.;
  for (int ipt  = 1; ipt<nPoints-1; ipt++) {
    double h = -1.+2.*double(ipt)/double(nPoints-1);
    double f = phi[ipt-1]+magic/sqrt(double(nPoints))/sqrt(1.-h*h);
    theta [ipt] = acos(h);
    phi   [ipt] = fmod(f,2.*TMath::Pi());
  }
  theta[nPoints-1]  = 0.;
  phi  [nPoints-1]  = 0.;

  //      0   1   2   3     4      5    6   7   8   9   10   11     12    13
  // X = {X1, Y1, Z1, T1, theta1, phi1, p1, X2, Y2, Z2, T2, theta2, phi2, p2}
  // Now start setting/changing parameters

  fitpars[7] = fitpars[0]; 
  fitpars[8] = fitpars[1]; 
  fitpars[9] = fitpars[2]; 
  fitpars[10] = fitpars[3]-12;
  fitpars[13] = 6.;

  const int nLowest = 3;
  double lowestnglnL[nLowest+1];
  int iLowestnglnL[nLowest+1];
  int itLowestnglnL[nLowest+1];
  for (int il=nLowest-1; il>=0; il--) {
    lowestnglnL[il] = 1.e100;
    iLowestnglnL[il] = -1;
    itLowestnglnL[il] = -1;
  }

  const int ntPoints = 11;
  double toffsets[ntPoints] = {50., 40., 30., 25., 20., 15., 12.5, 10., 7.5, 5., 2.5};

  for (int itoff=0; itoff<ntPoints; itoff++) {

    fitpars[10] = fitpars[3] - toffsets[itoff];

    for (int ipt=0; ipt<nPoints; ipt++) {
      fitpars[11] = theta[ipt];
      fitpars[12] = phi[ipt];

      double nglnL = MuGammaNglogL(fitpars,PCflg);
      //      std::cout << "(t,theta,phi,nglnL) = (" << fitpars[10] << "," << fitpars[11] << "," << fitpars[12] << "," << nglnL << ")" << std::endl;
      for (int il=nLowest-1; il>=0; il--) {
	if (nglnL < lowestnglnL[il]) {
	  lowestnglnL[il+1] = lowestnglnL[il];
	  iLowestnglnL[il+1] = iLowestnglnL[il];
	  itLowestnglnL[il+1] = itLowestnglnL[il];
	  lowestnglnL[il] = nglnL;
	  iLowestnglnL[il] = ipt;
	  itLowestnglnL[il] = itoff;
	  //  std::cout << "!!!!!! moving (t,theta,phi,nglnL) = (" << fitpars[10] << "," << fitpars[11] << "," << fitpars[12] << "," << nglnL << ") to position " << il << std::endl;
	} 
	else 
	  break;

      }//for il
    }// for ipt

  } // loop over time - for (int itoff=0;

  for (int il=0; il<nLowest; il++) {
    std::cout << "!!!!!! Lowest lglnL[" << il << "] = " << lowestnglnL[il] << " (toffset=" << toffsets[itLowestnglnL[il]] << ",theta=" << theta[iLowestnglnL[il]] << ",phi=" << phi[iLowestnglnL[il]] << ")" << std::endl;
  }

  fitpars[10] = fitpars[3] - toffsets[itLowestnglnL[0]];
  fitpars[11] = theta[iLowestnglnL[0]]; 
  fitpars[12] = phi[iLowestnglnL[0]];

  return MuGammaNglogL(fitpars,PCflg);

}
////************************************************************************************************************************************
////void PDK_MuGamma::SaveTestPDK_MuGammaFit(double* fitparams, int istage, double totmu, double nll, int PCflg); 
////************************************************************************************************************************************
//void PDK_MuGamma::SaveTestPDK_MuGammaFit(double* fitparams, int istage, double totmu, double nll, int PCflg){ //I have to change a lot in fQROOTOut.cc and also fitqunout.h!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Maybe this fct shold go to fiTQun_shared........
//  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];
//  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]);
//  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]);
//  fitqunpi0test_.fqtestpi0mom2[fitqunpi0test_.fqtestnpi0] = fitparams[10];
//  fitqunpi0test_.fqtestpi0dconv2[fitqunpi0test_.fqtestnpi0] = fitparams[11];
//
//  TVector3 dir1vect(sin(fitparams[4])*cos(fitparams[5]),
//                    sin(fitparams[4])*sin(fitparams[5]),
//                    cos(fitparams[4]));
//  TVector3 dir2vect(sin(fitparams[8])*cos(fitparams[9]),
//                    sin(fitparams[8])*sin(fitparams[9]),
//                    cos(fitparams[8]));
//
//  fitqunpi0test_.fqtestpi0mass[fitqunpi0test_.fqtestnpi0] = sqrt( 2.*fitparams[6]*fitparams[10]
//								  *(1 - dir1vect.Dot(dir2vect)) );
//  fitqunpi0test_.fqtestpi0photangle[fitqunpi0test_.fqtestnpi0] = acos(dir1vect.Dot(dir2vect));
//  TVector3 momtot = fitparams[6]*dir1vect + fitparams[10]*dir2vect;
//  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++;
//
//}
//************************************************************************************************************************************
//double PDK_MuGamma::PDK_MuGammaGetTwoRngnglogL(double *X, int& PCflg);
//************************************************************************************************************************************
double PDK_MuGamma::MuGammaNglogL(double *X, int& PCflg){

  Resetmu();
  //      0   1   2   3     4      5    6   7   8   9   10   11     12    13  
  // X = {X1, Y1, Z1, T1, theta1, phi1, p1, X2, Y2, Z2, T2, theta2, phi2, p2}

  double X1st[8];
  for (int i=0; i<=6; i++) X1st[i]=X[i];
  X1st[7] = 0.;

  double X2nd[8];
  for (int i=0; i<=6; i++) X2nd[i]=X[i+7];
  X2nd[7] = 0.;

  if (same3Vertex) {
    for (int i=0; i<=2; i++) X2nd[i]=X[i];
  }

  int PCflg1=fiTQun::OneRing(imu,X1st,0);
  int PCflg2=fiTQun::OneRing(ie,X2nd,1);
  //int PCflg2 = 0;
  PCflg = abs(PCflg1)+2*abs(PCflg2);
  if ((PCflg1<0)||(PCflg2<0)) PCflg = -1;

  double nll = fiTQun::nglogL();

  return nll;

}
//************************************************************************************************************************************
//PDK_MuGamma::GetTruePDK_MuGammaParams(double* tParams); Get true information for seed = 1 in the DoFit function.
//************************************************************************************************************************************

void PDK_MuGamma::GetTruePDK_MuGammaParams(double* tParams){ //NEEDS A LOT OF CHANGES!!!!!!!!!! not sure how to do them.

  bool printStuff = true;
  
  int temp_track = -1;

  for (int ipar=0; ipar<fiTQun_shared::nPDK_MuGammaParams; ipar++) {
    tParams[ipar] = 0.;
  }

#ifndef NOSKLIBRARIES
  /* load mc_particle structure */
  int nmctrks;
  mcguts(&nmctrks);

  float trgofst;
  vcrdvccm_();
  trginfo_(&trgofst);

  if (printStuff) std::cout << "This event has " << nmctrks << " true tracks" << std::endl;
  for (int itrk=0; itrk<nmctrks; itrk++) {
    int parentid = mc_particle.parent[itrk];
    int inOutOrSecondary = mc_particle.origin[itrk];
    int pdgCode = mc_particle.type[itrk];
    float mass = mc_particle.mass[itrk];
    float momentum = mc_particle.momentum[itrk];
    float energy = mc_particle.energy[itrk];
    float beta = mc_particle.beta[itrk];
    float dir[3];
    float inivtx[3];
    float finvtx[3];
    for (int ix=0; ix<3; ix++) {
      dir[ix] = mc_particle.dir[itrk][ix];
      /* NUANCE outgoing particles get mislabelled as incoming, so if the
         initial_vertex is not set, use the the_vertex instead */
      inivtx[ix] = mc_particle.initial_vertex[itrk][ix]; // so this is broken for Nuance input files
      finvtx[ix] = mc_particle.final_vertex[itrk][ix];
    } float time = mc_particle.time[itrk];
    if (printStuff) {
      std::cout << "track number " << itrk << " has PDG code, momentum, parent = " << pdgCode << ", " << momentum << ", " << parentid << \
	std::endl;
      std::cout << "direction = (" << dir[0] << "," << dir[1] << "," << dir[2] << ")" << std::endl;
      std::cout << "inivtx = (" << inivtx[0] << "," << inivtx[1] << "," << inivtx[2] << ")" << std::endl;
      std::cout << "finvtx = (" << finvtx[0] << "," << finvtx[1] << "," << finvtx[2] << ")" << std::endl;
      std::cout << "time, trgofst, mass, beta = " << time << ", " << trgofst << ", " << mass << ", " << beta << std::endl;
    }

    if (pdgCode == 321 && parentid == 0)
      temp_track = itrk;

    if (pdgCode == -13 && parentid == temp_track) {
      tParams[0] = inivtx[0];
      tParams[1] = inivtx[1];
      tParams[2] = inivtx[2];
      tParams[3] = trgofst + time;
      tParams[4] = acos(dir[2]);
      tParams[5] = atan2(dir[1],dir[0]);
      tParams[6] = momentum;
    }

    if (pdgCode == 22 && parentid == 0) {
      tParams[7] = inivtx[0];
      tParams[8] = inivtx[1];
      tParams[9] = inivtx[2];
      if (printStuff) std::cout << "setting tParams[3] to " << trgofst + time << std::endl;
      tParams[10] = trgofst + time;
      tParams[11] = acos(dir[2]);
      tParams[12] = atan2(dir[1],dir[0]);
      tParams[13] = momentum;
    }

  }
#endif

}

//************************************************************************************************************************************
//double PDK_MuGamma::DoFit( int seedtype, double *X, int& PCflg, int defaultIndex, double *X0 );
//************************************************************************************************************************************
double PDK_MuGamma::DoFit( int ievent_, int itwnd_, int seedtype, int fittype ) {

  fiTQun_shared::nring = 2;
  
  ievent = ievent_;
  itwnd = itwnd_;

  int PCflg;

  double tmpctotqpnl=fiTQun_shared::Get()->GetTotqCoef(); // changed ctotqpnl to fiTQun_shared::
  fiTQun_shared::Get()->SetTotqCoef(0);//no total charge constraint                                    

  fiTQun_shared::Get()->SetFullScatflg(true); //Same as ctotqpnl bool var.

  //!!!!!!!!!!!!!!!!! seedtype values  CHANGE!!!!!!!!!!!!!!!!    these are from Pi0 DoFit      
  //  1 = seed with truth                                 
  //  2 = seed with fixed efit + scan with 50 MeV ring                                                              

  double tfitpars[fiTQun_shared::nPDK_MuGammaParams];
  double fitpars[fiTQun_shared::nPDK_MuGammaParams];
  int fixedpars[fiTQun_shared::nPDK_MuGammaParams];
  double nglnL;
  double totmu;

  GetTruePDK_MuGammaParams(tfitpars);
  if (tfitpars[3] != 0) {
    nglnL = MuGammaNglogL(tfitpars,PCflg);
    totmu=GetTotmu();
    std::cout << "(1) totmu = " << totmu << std::endl;
    //    fqshared->SaveTestPDK_MuGammaFit(tfitpars,0,totmu,nglnL,PCflg); //CHANGE!!!!!! probably from shared.cc
  } else {
    std::cout << "Could not get the true time of the muon in this event!" << std::endl;
  }

  //These are for many different seed types. For now I'm trying just one, so check after the comment BLOCK.
  //Only working seedtype currently.....seed=2 is for seed with fixed efit + scan with 50 MeV ring
  if(seedtype == 1){
    // seed with truth
    same3Vertex = true;
    for (int ipar=0; ipar<fiTQun_shared::nPDK_MuGammaParams; ipar++) {
      fitpars[ipar] = tfitpars[ipar];
    }
    
    /*    prints10 = 0;
    
    printstuff = true;
    
    /*
    fitpars[0] = -938.085;
    fitpars[1] = 858.936;
    fitpars[2] = 187.512;
    fitpars[3] = 950.;
    fitpars[4] = 0.731066;
    fitpars[5] = 0.60179;
    fitpars[6] = 235.532;    
    
    std::cout << "calling nglogL inside DoFit" << std::endl;
    */
    nglnL = MuGammaNglogL(fitpars,PCflg);
    //printstuff = false;
  }
  else if (seedtype == 2){
    // scan with 1-ring e-like fit fixed
    // start with 1-ring e-like fit result
    fqshared->GetDefaultSnglTrkFit(fitpars,itwnd,0,imu,totmu,nglnL,PCflg);
    nglnL = SetPDK_MuGammaSeed(fitpars,PCflg);
  }

  else if (seedtype == 3){
    double mufittime = 0;
    fqshared->GetDefaultSnglTrkFit(fitpars,itwnd,0,imu,totmu,nglnL,PCflg);
    mufittime = fitpars[3];
    same3Vertex = true;
    for (int ipar=0; ipar<fiTQun_shared::nPDK_MuGammaParams; ipar++) 
      fitpars[ipar] = tfitpars[ipar];
    fitpars[10] =  mufittime + fitpars[10] - fitpars[3];
    fitpars[3] = mufittime;
  }
  else
    std::cout << "ERROR: Unknown seed type:" << seedtype << " For now, this can only handle seed = 2. " << std::endl;

  totmu=GetTotmu();
  std::cout << "(2) totmu = " << totmu << std::endl;

  TString strBase;
  strBase=Form("%d_KNuMuG",ievent);

  //  if (fiTQun_shared::Get()->GetDrawflg() ) DrawPreddist(Form("%s_0",strBase.Data()));

  //  fqshared->SaveTestPDK_MuGammaFit(fitpars,1,totmu,nglnL,PCflg);

  std::cout << "###########Initial parameters: ";
  for (int i=0; i<fiTQun_shared::nPDK_MuGammaParams; i++) {
    std::cout << fitpars[i] << ", ";
  }

  std::cout << "###############" << std::endl;

  for (int ipar=0; ipar<fiTQun_shared::nPDK_MuGammaParams; ipar++) fixedpars[ipar] = 0;
  if (fittype == 0) { // Do full 14 parameter fit 
    same3Vertex = false;
  } else if (fittype == 1) { // Assume common 3-vertex (11 parameters)
    same3Vertex = true;
    for (int ipar=7; ipar<=9; ipar++) fixedpars[ipar] = 1;
  } else if (fittype == 2) { // Fit only the muon parameters
    same3Vertex = false;
    for (int ipar=7; ipar<14; ipar++) fixedpars[ipar] = 1;
    fitpars[13] = 0.01;
  }

  //  if (fiTQun_shared::Get()->GetDrawflg() ) DrawPreddist();

  // Do final fit
  nglnL = RunMinimizer(fitpars,PCflg,fixedpars); 
  totmu=GetTotmu();
  //  printstuff = false;

  //  if (fiTQun_shared::Get()->GetDrawflg() ) DrawPreddist(Form("%s_fit",strBase.Data()));

  std::cout << "###########        Fit result: ";
  for (int i=0; i<fiTQun_shared::nPDK_MuGammaParams; i++) {
    std::cout << fitpars[i] << ", ";
  }
  //pi0PhotAngle[1] = GetPhotAngle(pi0Params[1][4],pi0Params[1][5],pi0Params[1][8],pi0Params[1][9]);
  std::cout << "###############" << std::endl;
  std::cout << "###########True parameters: ";
  for (int i=0; i<fiTQun_shared::nPDK_MuGammaParams; i++) {
    std::cout << tfitpars[i] << ", ";
  }
  std::cout << "###############" << std::endl;
  std::cout << "nglnL = " << nglnL << std::endl;

  fqshared->SaveDefaultPDK_MuGammaFit(fitpars,0,totmu,nglnL,PCflg);

  fiTQun_shared::Get()->SetFullScatflg(false);
  fiTQun_shared::Get()->SetTotqCoef(tmpctotqpnl);

  return nglnL;

}