#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TString.h>
#include <Riostream.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TFile.h>
#include <TF1.h>
#include <TH1D.h>
#include <TH2.h>
#include <TStyle.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TStopwatch.h>
#include <TPaveStats.h>
#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>
#include <Math/QuantFuncMathMore.h>
#include <TRandom3.h>

#include <signal.h>

#include <gsl/gsl_multimin.h>

using namespace std;

// #############################################################################

Int_t giNEvents;
Int_t giGenCor;
TString gtMacroDir;
TString gtUnpackDir;
TString gtCluDir;
TString gtInputFilePath;
Int_t giNDigiFiles;
TString gtMCLinkFileName;
Int_t giNLinkFiles;
TString gtTofGeoVersion;
TString gtSetupName;
Bool_t gbIdealCalibration;
Bool_t gbHeavyIonCollisions;
Bool_t gbIdealClustering;
Bool_t gbSeparateHitFile;
Bool_t gbAnalysisMCQA;
Int_t giDut;
Int_t giMRef;
Int_t giBRef;
Int_t giSel2;
Double_t gdScaleFactor;
Bool_t gbSimData;
Int_t giCalibrationMode;
Bool_t gbFinalRun;
Double_t gdTimeInSpill;
Int_t giNTISBins;
Double_t gdExtChi2Lim;
Double_t gdExtChi2Lim2;
Int_t giExtSigMode;
Int_t giTISBin;
Double_t gdUpperFitLimit;

// #############################################################################

const Bool_t gbRunSimplex = kFALSE; // FIXME
const Bool_t gbUseStartingValues = kTRUE;
const Bool_t gbDebugMode = kFALSE;

Double_t gdStartingResidualXSigma = 1.0; // [cm]
Double_t gdStartingResidualYSigma = 1.0; // [cm]
Double_t gdStartingResidualTSigma = 0.1; // [ns]

const Double_t gdDebugResidualXSigma = 0.3;  // [cm]
const Double_t gdDebugResidualYSigma = 0.7;  // [cm]
const Double_t gdDebugResidualTSigma = 0.08; // [ns]
const Double_t gdDebugResidualSigmaScale = 2.;

// #############################################################################

TH1* gtOriginalXResidual(NULL);
TH1* gtOriginalYResidual(NULL);
TH1* gtOriginalTResidual(NULL);
TH1* gtOriginalChiSquare(NULL);
TF1* gtOriginalChiSquareFit(NULL);
TH1* gtFinalXResidual(NULL);
TH1* gtFinalYResidual(NULL);
TH1* gtFinalTResidual(NULL);
TH1* gtFinalChiSquare(NULL);
TF1* gtFinalChiSquareFit(NULL);

Bool_t gbFirstCall(kTRUE);
Bool_t gbLastCall(kFALSE);

Bool_t gIsInterrupted(kFALSE);

void handler_ctrlc(int)
{
  gIsInterrupted = kTRUE;
}


Double_t chisqpdf(Double_t* x, Double_t* p);
Double_t chisqcdf(Double_t* x, Double_t* p);

Double_t chisq(const gsl_vector* v, void* params)
{
  cout<<"Calling chisq function..."<<endl;
  cout<<"         --------------------------------------------------------------------------------"<<endl;
  cout<<Form("         residual X sigma:  %10.3f cm", gsl_vector_get(v, 0))<<endl;
  cout<<Form("         residual Y sigma:  %10.3f cm", gsl_vector_get(v, 1))<<endl;
  cout<<Form("         residual T sigma:  %10.3f ns", gsl_vector_get(v, 2))<<endl;

  Double_t dResidualXSigma(gsl_vector_get(v, 0));
  Double_t dResidualYSigma(gsl_vector_get(v, 1));
  Double_t dResidualTSigma(gsl_vector_get(v, 2));

  // negative parameters make no sense and might lead to unforeseen consequences
  if( 0. >= dResidualXSigma || 0. >= dResidualYSigma || 0. >= dResidualTSigma )
  {
    return 100000000.;
  }

  if(!gbDebugMode)
  {

    gSystem->Exec(Form("root -b -q '%s/ana_hits.C(%d, %d, \"%s\", \"%s\", \"%s\", %d, \"%s\", %d, \"%s\", \"%s\", %d, %d, %d, %d, %d, "
                                                 "%d, %d, %d, %d, %f, %d, %d, %d, %f, %d, %f, %f, %d, %f, %f, %f, %d)'",
                       gtMacroDir.Data(), giNEvents, giGenCor, gtUnpackDir.Data(), gtCluDir.Data(), gtInputFilePath.Data(), giNDigiFiles,
                       gtMCLinkFileName.Data(), giNLinkFiles, gtTofGeoVersion.Data(), gtSetupName.Data(), gbIdealCalibration,
                       gbHeavyIonCollisions, gbIdealClustering, gbSeparateHitFile, gbAnalysisMCQA, giDut, giMRef, giBRef, giSel2,
                       gdScaleFactor, gbSimData, giCalibrationMode, gbFinalRun, gdTimeInSpill, giNTISBins, gdExtChi2Lim, gdExtChi2Lim2,
                       giExtSigMode, dResidualXSigma, dResidualYSigma, dResidualTSigma, giTISBin));

    TFile::Open("./ana_cluster.hst.root", "READ");
    if(!gFile)
    {
      cout<<"Could not open histogram file. Abort macro execution."<<endl;
      gIsInterrupted = kTRUE;
      return 100000000.;
    }
  }

  TH1* tChiSquare(NULL);
  TH1* tXResidual(NULL);
  TH1* tYResidual(NULL);
  TH1* tTResidual(NULL);

  if(!gbDebugMode)
  {
    if(2 == giExtSigMode)
    {
      tChiSquare = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllChi2Sel24"));
      tXResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDXSel24"));
      tYResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDYSel24"));
      tTResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDTSel24"));
    }
    else
    {
      if(giTISBin)
      {
        TH2* tChiSquareTIS = dynamic_cast<TH2*>(gFile->FindObjectAny("hAllSelHitTupleMatchChiSqTIS"));
        TH2* tXResidualTIS = dynamic_cast<TH2*>(gFile->FindObjectAny("hAllSelHitTupleMatchResidualXTIS"));
        TH2* tYResidualTIS = dynamic_cast<TH2*>(gFile->FindObjectAny("hAllSelHitTupleMatchResidualYTIS"));
        TH2* tTResidualTIS = dynamic_cast<TH2*>(gFile->FindObjectAny("hAllSelHitTupleMatchResidualTTIS"));

        tChiSquare = tChiSquareTIS->ProjectionY("_py", giTISBin, giTISBin);
        tXResidual = tXResidualTIS->ProjectionY("_py", giTISBin, giTISBin);
        tYResidual = tYResidualTIS->ProjectionY("_py", giTISBin, giTISBin);
        tTResidual = tTResidualTIS->ProjectionY("_py", giTISBin, giTISBin);
      }
      else
      {
        tChiSquare = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllChi204"));
        tXResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDX04"));
        tYResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDY04"));
        tTResidual = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllDT04"));
      }
    }
  }
  else
  {
    tXResidual = new TH1F("hResidualX", "", 401, -2., 2.);
    tYResidual = new TH1F("hResidualY", "", 601, -3., 3.);
    tTResidual = new TH1F("hResidualT", "", 101, -0.5, 0.5);
    tChiSquare = new TH1F("hChiSquare", "", 1000, 0., 100.);

    if(gRandom)
    {
      delete gRandom;
    }
    gRandom = new TRandom3(4357);

    for(Int_t iEvent = 0; iEvent < 1000000; iEvent++)
    {
      Double_t dResidualX = gRandom->Gaus(0., gdDebugResidualXSigma);
      Double_t dResidualY = gRandom->Gaus(0., gdDebugResidualYSigma);
      Double_t dResidualT = gRandom->Gaus(0., gdDebugResidualTSigma);
      Double_t dChiSquare = TMath::Power(dResidualX/dResidualXSigma, 2.) + TMath::Power(dResidualY/dResidualYSigma, 2.) + TMath::Power(dResidualT/dResidualTSigma, 2.);

      tXResidual->Fill(dResidualX);
      tYResidual->Fill(dResidualY);
      tTResidual->Fill(dResidualT);
      tChiSquare->Fill(dChiSquare);
    }
  }

  if(!gbDebugMode)
  {
    if(!tXResidual || !tYResidual || !tTResidual || !tChiSquare)
    {
      cout<<"Chi2 residual histograms not found in file. Abort macro execution."<<endl;
      gIsInterrupted = kTRUE;
      return 100000000.;
    }
  }


  TF1* tChiSquareFit = new TF1("tChiSquareFit", chisqpdf, 0., 15., 3);
  tChiSquareFit->FixParameter(0, 3.);
  tChiSquareFit->FixParameter(1, 1.);
  tChiSquareFit->SetParameter(2, tChiSquare->Integral("width"));

  tChiSquare->Smooth();
  TFitResultPtr tFitResult = tChiSquare->Fit(tChiSquareFit, "SQN", "", 0., gdUpperFitLimit);

  if(static_cast<Int_t>(tFitResult))
  {
    cout<<"Fit to chi2 histogram failed."<<endl;
    return 100000000.;
  }

  Double_t dFitChiSquare = tChiSquareFit->GetChisquare();
  Int_t iFitNDegreesOfFreedom = tChiSquareFit->GetNDF();
  Double_t dUpperFitLimitIntegralShare = tChiSquare->Integral(1, tChiSquare->FindBin(gdUpperFitLimit) - 1)/tChiSquare->GetEntries();

  if(gbFirstCall)
  {
    gtOriginalXResidual = dynamic_cast<TH1*>(tXResidual->Clone("tOriginalXResidual"));
    gtOriginalYResidual = dynamic_cast<TH1*>(tYResidual->Clone("tOriginalYResidual"));
    gtOriginalTResidual = dynamic_cast<TH1*>(tTResidual->Clone("tOriginalTResidual"));
    gtOriginalChiSquare = dynamic_cast<TH1*>(tChiSquare->Clone("gtOriginalChiSquare"));

    gtOriginalXResidual->SetDirectory(gROOT);
    gtOriginalYResidual->SetDirectory(gROOT);
    gtOriginalTResidual->SetDirectory(gROOT);
    gtOriginalChiSquare->SetDirectory(gROOT);

    gtOriginalChiSquareFit = dynamic_cast<TF1*>(tChiSquareFit->Clone("gtOriginalChiSquareFit"));

    gbFirstCall = kFALSE;
  }

  if(gbLastCall)
  {
    gtFinalXResidual = dynamic_cast<TH1*>(tXResidual->Clone("tFinalXResidual"));
    gtFinalYResidual = dynamic_cast<TH1*>(tYResidual->Clone("tFinalYResidual"));
    gtFinalTResidual = dynamic_cast<TH1*>(tTResidual->Clone("tFinalTResidual"));
    gtFinalChiSquare = dynamic_cast<TH1*>(tChiSquare->Clone("gtFinalChiSquare"));

    gtFinalXResidual->SetDirectory(gROOT);
    gtFinalYResidual->SetDirectory(gROOT);
    gtFinalTResidual->SetDirectory(gROOT);
    gtFinalChiSquare->SetDirectory(gROOT);

    gtFinalChiSquareFit = dynamic_cast<TF1*>(tChiSquareFit->Clone("gtFinalChiSquareFit"));

    gbLastCall = kFALSE;

    Double_t dOriginalFitChiSquare = gtOriginalChiSquareFit->GetChisquare();
    Double_t dOriginalUpperFitLimitIntegralShare = gtOriginalChiSquare->Integral(1, gtOriginalChiSquare->FindBin(gdUpperFitLimit) - 1)/gtOriginalChiSquare->GetEntries();

    cout<<endl;
    cout<<Form("final    residual X sigma:  %12.5f cm <- %12.5f cm", dResidualXSigma, gdStartingResidualXSigma)<<endl;
    cout<<Form("final    residual Y sigma:  %12.5f cm <- %12.5f cm", dResidualYSigma, gdStartingResidualYSigma)<<endl;
    cout<<Form("final    residual T sigma:  %12.5f ns <- %12.5f ns", dResidualTSigma, gdStartingResidualTSigma)<<endl;
    cout<<Form("final    fit chi square:    %12.5f    <- %12.5f", dFitChiSquare, dOriginalFitChiSquare)<<endl;
    cout<<Form("final    fit chi^2/ndf:     %12.5f    <- %12.5f", dFitChiSquare/iFitNDegreesOfFreedom, dOriginalFitChiSquare/iFitNDegreesOfFreedom)<<endl;
    cout<<Form("final    integral share:    %12.5f    <- %12.5f", dUpperFitLimitIntegralShare, dOriginalUpperFitLimitIntegralShare)<<endl;
    cout<<Form("final    relative share:    %12.5f    <- %12.5f", dUpperFitLimitIntegralShare/ROOT::Math::chisquared_cdf(gdUpperFitLimit, 3),
                                                                dOriginalUpperFitLimitIntegralShare/ROOT::Math::chisquared_cdf(gdUpperFitLimit, 3))<<endl;
    cout<<Form("         nominal share:     %12.5f", ROOT::Math::chisquared_cdf(gdUpperFitLimit, 3))<<endl;
    cout<<Form("         fit ndf:           %12d", iFitNDegreesOfFreedom)<<endl;
    cout<<endl;
  }


  cout<<Form("         fit chi square:    %12.5f", dFitChiSquare)<<endl;
  cout<<Form("         fit ndf:           %12d", iFitNDegreesOfFreedom)<<endl;
  cout<<Form("         fit chi^2/ndf:     %12.5f", dFitChiSquare/iFitNDegreesOfFreedom)<<endl;
  cout<<Form("         integral share:    %12.5f", dUpperFitLimitIntegralShare)<<endl;
  cout<<Form("         nominal share:     %12.5f", ROOT::Math::chisquared_cdf(gdUpperFitLimit, 3))<<endl;
  cout<<endl;

  if(!gbDebugMode)
  {
    gFile->Close();
    delete gFile;
  }
  else
  {
    delete tXResidual;
    delete tYResidual;
    delete tTResidual;
    delete tChiSquare;
  }

  delete tChiSquareFit;

  if(!gbDebugMode)
  {
    return TMath::Abs(1. - dFitChiSquare/iFitNDegreesOfFreedom);
  }
  else
  {
    return  TMath::Power(1.*(1. - dFitChiSquare/iFitNDegreesOfFreedom), 2.)
           +TMath::Power(100.*(ROOT::Math::chisquared_cdf(gdUpperFitLimit, 3) - dUpperFitLimitIntegralShare), 2.);
  }
}


Bool_t get_sigma_scale()
{
  cout<<"Calling get_sigma_scale function..."<<endl;
  cout<<"         --------------------------------------------------------------------------------"<<endl;

  gSystem->Exec(Form("root -b -q '%s/ana_hits.C(%d, %d, \"%s\", \"%s\", \"%s\", %d, \"%s\", %d, \"%s\", \"%s\", %d, %d, %d, %d, %d, "
                                               "%d, %d, %d, %d, %f, %d, %d, %d, %f, %d, %f, %f, %d, %f, %f, %f, %d)'",
                     gtMacroDir.Data(), giNEvents, giGenCor, gtUnpackDir.Data(), gtCluDir.Data(), gtInputFilePath.Data(), giNDigiFiles,
                     gtMCLinkFileName.Data(), giNLinkFiles, gtTofGeoVersion.Data(), gtSetupName.Data(), gbIdealCalibration,
                     gbHeavyIonCollisions, gbIdealClustering, gbSeparateHitFile, gbAnalysisMCQA, giDut, giMRef, giBRef, giSel2,
                     gdScaleFactor, gbSimData, giCalibrationMode, gbFinalRun, gdTimeInSpill, giNTISBins, gdExtChi2Lim, gdExtChi2Lim2,
                     0, 0., 0., 0., giTISBin));

  TFile::Open("./ana_cluster.hst.root", "READ");
  if(!gFile)
  {
    cout<<"Could not open histogram file. Abort macro execution."<<endl;
    return kFALSE;
  }


  TH1* tChiSquare(NULL);

  if(2 == giExtSigMode)
  {
    tChiSquare = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllChi2Sel24"));
  }
  else
  {
    if(giTISBin)
    {
      TH2* tChiSquareTIS = dynamic_cast<TH2F*>(gFile->FindObjectAny("hAllSelHitTupleMatchChiSqTIS"));

      if(tChiSquareTIS)
      {
        tChiSquare = tChiSquareTIS->ProjectionY("_py", giTISBin, giTISBin);
      }
    }
    else
    {
      tChiSquare = dynamic_cast<TH1*>(gFile->FindObjectAny("hAllChi204"));
    }
  }

  if(!tChiSquare)
  {
    cout<<"Chi2 histogram not found in input calibration file. Abort macro execution."<<endl;
    return kFALSE;
  }

  TF1* tChiSquareFit = new TF1("tChiSquareFit", chisqpdf, 0., 15., 3);
  tChiSquareFit->FixParameter(0, 3.);
  tChiSquareFit->SetParameter(1, 1.5);
  tChiSquareFit->SetParameter(2, tChiSquare->Integral("width"));

  tChiSquare->Smooth();
  TFitResultPtr tFitResult = tChiSquare->Fit(tChiSquareFit, "SQN", "", 0., gdUpperFitLimit);

  if(static_cast<Int_t>(tFitResult))
  {
    cout<<"Fit to chi2 histogram failed. Abort macro execution."<<endl;
    return kFALSE;
  }

//  Double_t dChi2SigmaScale = tChiSquareFit->GetParameter(1); // TODO: modified to make use also of negative (?) sigma scales
  Double_t dChi2SigmaScale = TMath::Abs(tChiSquareFit->GetParameter(1));
  Double_t dFitChiSquare = tChiSquareFit->GetChisquare();
  Int_t iFitNDegreesOfFreedom = tChiSquareFit->GetNDF();

  gFile->Close();
  delete gFile;

  delete tChiSquareFit;


  Double_t dResidualXSigma;
  Double_t dResidualYSigma;
  Double_t dResidualTSigma;

  TFile::Open("./calib_ana.cor_out.root", "READ");
  if(!gFile)
  {
    cout<<"Could not open output calibration file. Abort macro execution."<<endl;
    return kFALSE;
  }

  TH1* tHisto;

  if(2 == giExtSigMode)
  {
    tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelHitTupleResidualXYT_Width"));

    dResidualXSigma = tHisto->GetBinContent(1);
    dResidualYSigma = tHisto->GetBinContent(2);
    dResidualTSigma = tHisto->GetBinContent(3);
  }
  else
  {
    if(giTISBin)
    {
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualX_WidthTIS"));
      dResidualXSigma = tHisto->GetBinContent(1, giTISBin);
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualY_WidthTIS"));
      dResidualYSigma = tHisto->GetBinContent(1, giTISBin);
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualT_WidthTIS"));
      dResidualTSigma = tHisto->GetBinContent(1, giTISBin);
    }
    else
    {
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualX_Width"));
      dResidualXSigma = tHisto->GetBinContent(1);
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualY_Width"));
      dResidualYSigma = tHisto->GetBinContent(1);
      tHisto = dynamic_cast<TH1*>(gFile->FindObjectAny("hSelTypeNNResidualT_Width"));
      dResidualTSigma = tHisto->GetBinContent(1);
    }
  }

  gFile->Close();
  delete gFile;

  gdStartingResidualXSigma = dResidualXSigma/dChi2SigmaScale;
  gdStartingResidualYSigma = dResidualYSigma/dChi2SigmaScale;
  gdStartingResidualTSigma = dResidualTSigma/dChi2SigmaScale;

  cout<<endl;
  cout<<Form("initial  residual X sigma:  %12.5f cm", dResidualXSigma)<<endl;
  cout<<Form("initial  residual Y sigma:  %12.5f cm", dResidualYSigma)<<endl;
  cout<<Form("initial  residual T sigma:  %12.5f ns", dResidualTSigma)<<endl;
  cout<<Form("         chi^2 sigma scale: %12.5f", dChi2SigmaScale)<<endl;
  cout<<Form("         fit chi square:    %12.5f", dFitChiSquare)<<endl;
  cout<<Form("         fit ndf:           %12d", iFitNDegreesOfFreedom)<<endl;
  cout<<Form("         fit chi^2/ndf:     %12.5f", dFitChiSquare/iFitNDegreesOfFreedom)<<endl;
  cout<<Form("scaled   residual X sigma:  %12.5f cm", gdStartingResidualXSigma)<<endl;
  cout<<Form("scaled   residual Y sigma:  %12.5f cm", gdStartingResidualYSigma)<<endl;
  cout<<Form("scaled   residual T sigma:  %12.5f ns", gdStartingResidualTSigma)<<endl;
  cout<<endl;

  return kTRUE;
}
#endif

void find_chi2_weights(Int_t iNEvents, Int_t iGenCor, const TString& tUnpackDir, const TString& tCluDir,
                       const TString& tInputFilePath, Int_t iNDigiFiles, const TString& tMCLinkFileName, Int_t iNLinkFiles,
                       const TString& tTofGeoVersion, const TString& tSetupName, Bool_t bIdealCalibration,
                       Bool_t bHeavyIonCollisions, Bool_t bIdealClustering, Bool_t bSeparateHitFile, Bool_t bAnalysisMCQA,
                       Int_t iDut, Int_t iMRef, Int_t iBRef, Int_t iSel2, Double_t dScaleFactor, Bool_t bSimData,
                       Int_t iCalibrationMode, Bool_t bFinalRun, Double_t dTimeInSpill, Int_t iNTISBins,
                       Double_t dExtChi2Lim, Double_t dExtChi2Lim2, Int_t iExtSigMode, const TString& tMacroDir,
                       Int_t iTISBin = 0)
{
#if defined(__CINT__) && !defined(__MAKECINT__)
  cout<<endl;
  cout<<"Macro needs to be compiled using ACLiC!"<<endl;
#else
  cout<<endl;

  gSystem->IgnoreInterrupt(kTRUE);
  gIsInterrupted = kFALSE;

  gErrorIgnoreLevel = kWarning;
  gStyle->SetOptStat("emr");


  // run parameter initialization
  giNEvents = iNEvents;
  giGenCor = iGenCor;
  gtMacroDir = tMacroDir;
  gtUnpackDir = tUnpackDir;
  gtCluDir = tCluDir;
  gtInputFilePath = tInputFilePath;
  giNDigiFiles = iNDigiFiles;
  gtMCLinkFileName = tMCLinkFileName;
  giNLinkFiles = iNLinkFiles;
  gtTofGeoVersion = tTofGeoVersion;
  gtSetupName = tSetupName;
  gbIdealCalibration = bIdealCalibration;
  gbHeavyIonCollisions = bHeavyIonCollisions;
  gbIdealClustering = bIdealClustering;
  gbSeparateHitFile = bSeparateHitFile;
  gbAnalysisMCQA = bAnalysisMCQA;
  giDut = iDut;
  giMRef = iMRef;
  giBRef = iBRef;
  giSel2 = iSel2;
  gdScaleFactor = dScaleFactor;
  gbSimData = bSimData;
  giCalibrationMode = iCalibrationMode;
  gbFinalRun = bFinalRun;
  gdTimeInSpill = dTimeInSpill;
  giNTISBins = iNTISBins;
  gdExtChi2Lim = dExtChi2Lim;
  gdExtChi2Lim2 = dExtChi2Lim2;
  giExtSigMode = iExtSigMode;
  giTISBin = iTISBin;

  if(2 == giExtSigMode)
  {
    gdUpperFitLimit = gdExtChi2Lim2;
    gdExtChi2Lim2 = 1000000.;
  }
  else if(1 == giExtSigMode)
  {
    gdUpperFitLimit = gdExtChi2Lim;
    gdExtChi2Lim = 1000000.;

  }
  else
  {
    cout<<Form("Invalid ExtSigMode %d requested. Abort macro execution.", giExtSigMode)<<endl;
    return;
  }

  if(15. < gdUpperFitLimit)
  {
    cout<<"Invalid chi2 limit > 15 requested. Abort macro execution."<<endl;
    return;
  }


  // fit parameter initialization
  // 0: residual X sigma
  // 1: residual Y sigma
  // 2: residual T sigma
  Int_t iNbVariables(3);

  gsl_vector* tVariableValues = gsl_vector_alloc(iNbVariables);
  gsl_vector* tVariableSteps  = gsl_vector_alloc(iNbVariables);

  // residual X sigma
  gsl_vector_set(tVariableValues, 0, 1.0);
  gsl_vector_set(tVariableSteps, 0, -0.05);

  // residual Y sigma
  gsl_vector_set(tVariableValues, 1, 1.0);
  gsl_vector_set(tVariableSteps, 1, -0.05);

  // residual T sigma
  gsl_vector_set(tVariableValues, 2, 0.1);
  gsl_vector_set(tVariableSteps, 2, -0.005);

  // final values from previous macro execution?
  if(gbUseStartingValues)
  {
    if(gbDebugMode)
    {
      gdStartingResidualXSigma = gdDebugResidualSigmaScale*gdDebugResidualXSigma;
      gdStartingResidualYSigma = gdDebugResidualSigmaScale*gdDebugResidualYSigma;
      gdStartingResidualTSigma = gdDebugResidualSigmaScale*gdDebugResidualTSigma;
    }
    else
    {
      if(!get_sigma_scale())
      {
        cout<<"Failed to determine the chi2 sigma scale. Abort macro execution."<<endl;
        return;
      }

      gSystem->Exec("mkdir -p -v ./plots/InitialSigmaPlots");
      gSystem->Exec("mv -v ./*.pdf ./plots/InitialSigmaPlots/");
    }

    gsl_vector_set(tVariableValues, 0, gdStartingResidualXSigma);
    gsl_vector_set(tVariableValues, 1, gdStartingResidualYSigma);
    gsl_vector_set(tVariableValues, 2, gdStartingResidualTSigma);
  }


  Int_t iIteration = 0;
  Int_t iMinimizationStatus;
  Double_t iSimplexSize;

  Double_t dDummyParameters[1] = {1.0};

  gsl_multimin_function chisq_func;
  chisq_func.n = iNbVariables;
  chisq_func.f = chisq;
  chisq_func.params = dDummyParameters;

  TStopwatch timer;
  timer.Start();

  const gsl_multimin_fminimizer_type* tMinimizerType = gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer* tMinimizer = gsl_multimin_fminimizer_alloc(tMinimizerType, iNbVariables);


  if(gbRunSimplex)
  {
    gsl_multimin_fminimizer_set(tMinimizer, &chisq_func, tVariableValues, tVariableSteps);

    do // beginning of minimization iteration loop
    {
      signal(SIGINT, handler_ctrlc);

      if(gIsInterrupted)
      {
        cout<<endl;
        cout<<"CTRL+C received. Skipping the remaining minimization iterations."<<endl;
        cout<<endl;
        break;
      }

      iIteration++;

      iMinimizationStatus = gsl_multimin_fminimizer_iterate(tMinimizer);

      if(iMinimizationStatus)
      {
        break;
      }

      iSimplexSize = gsl_multimin_fminimizer_size(tMinimizer);
      iMinimizationStatus = gsl_multimin_test_size(iSimplexSize, 1e-4);//1e-3); //1e-2 before // TODO

      if(GSL_SUCCESS == iMinimizationStatus)
      {
        cout<<"converged to minimum at"<<endl;
      }

      cout<<Form("%5d  chisq = %7.3f, simplex size = %7.3f", iIteration, tMinimizer->fval, iSimplexSize)<<endl;
      cout<<"       --------------------------------------------------------------------------------"<<endl;

    } // end of minimization iteration loop
    while(GSL_CONTINUE == iMinimizationStatus && 100 > iIteration);

    cout<<"       --------------------------------------------------------------------------------"<<endl;
    gbLastCall = kTRUE;
    chisq(tMinimizer->x, dDummyParameters);
  }
  else
  {
    cout<<"       --------------------------------------------------------------------------------"<<endl;
    gbLastCall = kTRUE;
    chisq(tVariableValues, dDummyParameters);
  }

  cout<<"       --------------------------------------------------------------------------------"<<endl;
  cout<<" final parameter values                                                                "<<endl;
  cout<<"       --------------------------------------------------------------------------------"<<endl;

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  cout<<endl;
  cout<<" real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout<<endl;


  TCanvas* tCanvas = new TCanvas("tCanvas", "", 0., 0., 1400., 700.);
  tCanvas->Divide(4, 2);

  TH1* tHisto(NULL);
  TF1* tFunc(NULL);

  tCanvas->cd(1);
  tHisto = gtOriginalXResidual;
  if(tHisto)
  {
    tHisto->SetTitle("original;#Deltax [cm];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-5., 5.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(2);
  tHisto = gtOriginalYResidual;
  if(tHisto)
  {
    tHisto->SetTitle("original;#Deltay [cm];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-5., 5.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(3);
  tHisto = gtOriginalTResidual;
  if(tHisto)
  {
    tHisto->SetTitle("original;#Deltat [ns];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-1., 1.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(4);
  tHisto = gtOriginalChiSquare;
  if(tHisto)
  {
    tHisto->SetTitle("original;#chi^{2} [];;");
    tHisto->Draw();
//    gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(0., 15.);
    gPad->Modified();
    gPad->Update();

    tFunc = gtOriginalChiSquareFit;
    if(tFunc)
    {
      tFunc->SetNpx(1000000);
      tFunc->SetLineWidth(2);
      tFunc->SetRange(0., 15.);
      tFunc->Draw("SAME");
    }
  }

  tCanvas->cd(5);
  tHisto = gtFinalXResidual;
  if(tHisto)
  {
    tHisto->SetTitle("final;#Deltax [cm];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-5., 5.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(6);
  tHisto = gtFinalYResidual;
  if(tHisto)
  {
    tHisto->SetTitle("final;#Deltay [cm];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-5., 5.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(7);
  tHisto = gtFinalTResidual;
  if(tHisto)
  {
    tHisto->SetTitle("final;#Deltat [ns];;");
    tHisto->Draw();
//  gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(-1., 1.);
    gPad->Modified();
    gPad->Update();
  }

  tCanvas->cd(8);
  tHisto = gtFinalChiSquare;
  if(tHisto)
  {
    tHisto->SetTitle("final;#chi^{2} [];;");
    tHisto->Draw();
//    gPad->SetLogy();
    tHisto->GetXaxis()->SetRangeUser(0., 15.);
    gPad->Modified();
    gPad->Update();

    tFunc = gtFinalChiSquareFit;
    if(tFunc)
    {
      tFunc->SetNpx(1000000);
      tFunc->SetLineWidth(2);
      tFunc->SetRange(0., 15.);
      tFunc->Draw("SAME");
    }
  }

  tCanvas->SaveAs("optimal_chi2_sigmas.pdf");


  gSystem->Exec("mkdir -p -v ./plots/FinalSigmaPlots");
  gSystem->Exec("mv -v ./*.pdf ./plots/FinalSigmaPlots/");


  gsl_vector_free(tVariableValues);
  gsl_vector_free(tVariableSteps);
  gsl_multimin_fminimizer_free(tMinimizer);

  gSystem->IgnoreInterrupt(kFALSE);
#endif
}


Double_t chisqpdf(Double_t* x, Double_t* p)
{
  // p[0] : ndof
  // p[1] : sigma scale
  // p[2] : normalization

  return p[2]*TMath::Power(p[1], 2.)*ROOT::Math::chisquared_pdf(TMath::Power(p[1], 2.)*x[0], p[0]);
}

Double_t chisqcdf(Double_t* x, Double_t* p)
{
  // p[0] : ndof

  return ROOT::Math::chisquared_cdf(x[0], p[0]);
}