/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz et al,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();

  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); //GetTaskList();

  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
  MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");

  eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));

  evtloop->Eventloop();

  f = eval->GetChisq();
}

//------------------------------------------------------------------------

void CT1EnergyEst(Float_t maxhadronness=1.)
{
  //  Bool_t evalenergy=kFALSE;
  //
  // Fill events into a MHMatrix
  //
  MParList parlist;
  MHMatrix matrix;

  Int_t col = matrix.AddColumn("MMcEvt.fEnergy");

  MEnergyEstParam eest("Hillas");
  eest.Add("HillasSrc");
  eest.InitMapping(&matrix);

  MReadTree read("Events");
  read.AddFile("MC_ON2.root", 200);
  read.DisableAutoScheme();

  TString hcut("MHadronness.fHadronness<");
  hcut += maxhadronness;

  MF filterhadrons(hcut);

  if (!matrix.Fill(&parlist, &read))
    return;

  //
  // Setup the tasklist used to evaluate the needed chisq
  //
  MTaskList tasklist;
  parlist.AddToList(&tasklist);

  MMatrixLoop loop(&matrix);


  MChisqEval eval;
  eval.SetY1(new MDataElement(&matrix, col));
  eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
  eval.SetOwner();

  tasklist.AddToList(&loop);
  tasklist.AddToList(&eest);
  tasklist.AddToList(&eval);

  MEvtLoop evtloop;
  evtloop.SetParList(&parlist);
 
  //
  // Be careful: This is not thread safe
  //
  TMinuit minuit(12);
  minuit.SetPrintLevel(-1);
  minuit.SetFCN(fcn);

  // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
  if (minuit.SetErrorDef(1))
    {
      cout << "SetErrorDef failed." << endl;
      return;
    }

  //
  // Set initial values
  //
  TArrayD fA(5);
  fA[0] =1;//4916.4;   */-2414.75;
  fA[1] =1;//149.549;  */   1134.28;
  fA[2] =1;//-558.209; */   132.932;
  fA[3] =1;//0.270725; */  0.292845;
  fA[4] =1;//107.001;  */   107.001;
  
  TArrayD fB(7);
  fB[0] = 1;//-8234.12; */ -4282.25;
  fB[1] = 1;// 23.2153; */    18.892;
  fB[2] = 1;// 0.416372;*/  0.193373;
  fB[3] = 1;// 332.42;  */  203.803;
  fB[4] = 1;// -0.701764;*/ -0.534876;
  fB[5] = 1;//-0.0131774;*/ -0.00789539;
  fB[6] = 1;//-0.162687;*/   0.111913;
  
  // Set starting values and step sizes for parameters
  for (Int_t i=0; i<fA.GetSize(); i++)
    {
      TString name = "fA[";
      name += i;
      name += "]";
      Double_t vinit = fA[i];
      Double_t step  = fabs(fA[i]/3);

      Double_t limlo = 0; // limlo=limup=0: no limits
      Double_t limup = 0; 

      Bool_t rc = minuit.DefineParameter(i, name, vinit, step, limlo, limup);
      if (!rc)
	continue;

      cout << "Error in defining parameter #" << i << endl;
      return;
    }

  for (Int_t i=0; i<fB.GetSize(); i++)
    {
      TString name = "fB[";
      name += i;
      name += "]";
      Double_t vinit = fB[i];
      Double_t step  = fabs(fB[i]/3);

      Double_t limlo = 0; // limlo=limup=0: no limits
      Double_t limup = 0;

      Bool_t rc = minuit.DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
      if (!rc)
	continue;

      cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
      return;
    }

  //
  // Setup globals used in FCN
  //
  minuit.SetObjectFit(&evtloop);

   
  TStopwatch clock;
  clock.Start();

  // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
  gLog.SetNullOutput(kTRUE);
  Bool_t rc = minuit.Migrad();
  gLog.SetNullOutput(kFALSE);

  if (rc)
    {
      cout << "Migrad failed." << endl;
      return;
    }

  cout << endl;
  clock.Stop();
  clock.Print();
  cout << endl;

  cout << "Resulting Chisq: " << minuit.fAmin << endl;

  for (Int_t i=0; i<fA.GetSize()+fB.GetSize(); i++)
    {
      Double_t val;
      Double_t er;

      if (!minuit.GetParameter(i, val, er))
        {
	  cout << "Error getting parameter #" << i << endl;
	  return;
        }

      cout << i << ":  " << val << "  +-  " << er << endl;
    }

  /*
    // Print results
    Double_t amin, edm, errdef;
    Int_t nvpar, nparx, icstat;
    minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
    minuit.mnprin(3, amin);
  */
  eest.Print();



  //---------------------------------------------------------------

  // Part 2: Now test how the energy estimation method works.
  //
  //
  // Create a empty Parameter List and an empty Task List
  // The tasklist is identified in the eventloop by its name
  //


  MTaskList tlist2;
  parlist.Replace(&tlist2);

  //
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //


  read->SetEventNum(0);

  //
  // Use this to change the binnign of the histograms to CT1-style
  //
  Bool_t usect1 = kTRUE;

  MH3 mh3e("MMcEvt.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
  MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
  MH3 mh3eo("MMcEvt.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
  MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");

  MH3 mh3e2("MEnergyEst.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
  MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
  MH3 mh3eo2("MEnergyEst.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
  MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");

  MH3 mhe("MMcEvt.fEnergy",     "MEnergyEst.fEnergy");
  MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");

  mh3e.SetName("HistEnergy");
  mh3i.SetName("HistImpact");
  mh3eo.SetName("HistEnergyOffset");
  mh3io.SetName("HistImpactOffset");

  mh3e2.SetName("HistEnergy");
  mh3i2.SetName("HistImpact");
  mh3eo2.SetName("HistEnergyOffset");
  mh3io2.SetName("HistImpactOffset");

  mhe.SetName("HistEE");
  mhi.SetName("HistII");

  MFillH hfille(&mh3e);
  MFillH hfilli(&mh3i);
  MFillH hfilleo(&mh3eo);
  MFillH hfillio(&mh3io);

  MFillH hfille2(&mh3e2);
  MFillH hfilli2(&mh3i2);
  MFillH hfilleo2(&mh3eo2);
  MFillH hfillio2(&mh3io2);

  MFillH hfillee(&mhe);
  MFillH hfillii(&mhi);

  MBinning binsex("BinningHistEnergyX");
  MBinning binsey("BinningHistEnergyY");
  MBinning binsix("BinningHistImpactX");
  MBinning binsiy("BinningHistImpactY");
  MBinning binseox("BinningHistEnergyOffsetX");
  MBinning binseoy("BinningHistEnergyOffsetY");
  MBinning binsiox("BinningHistImpactOffsetX");
  MBinning binsioy("BinningHistImpactOffsetY");
  MBinning binseex("BinningHistEEX");
  MBinning binsiix("BinningHistIIX");
  MBinning binseey("BinningHistEEY");
  MBinning binsiiy("BinningHistIIY");

  binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
  binsey.SetEdges(50, 0, usect1 ? 0.8 : 1.75);
  binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
  binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);

  binsix.SetEdges(50, 0, usect1 ? 275 : 300);
  binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75);
  binsiox.SetEdges(50, 0, usect1 ? 275 : 300);
  binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);

  binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
  binseey.SetEdgesLog(50, usect1 ? 300 : 1,  usect1 ? 50000 : 2e3);
  binsiix.SetEdges(50, 0, usect1 ? 275 : 300);
  binsiiy.SetEdges(50, 0, usect1 ? 275 : 150);

  parlist.AddToList(&binsex);
  parlist.AddToList(&binsey);
  parlist.AddToList(&binsix);
  parlist.AddToList(&binsiy);
  parlist.AddToList(&binseox);
  parlist.AddToList(&binseoy);
  parlist.AddToList(&binsiox);
  parlist.AddToList(&binsioy);
  parlist.AddToList(&binseex);
  parlist.AddToList(&binseey);
  parlist.AddToList(&binsiix);
  parlist.AddToList(&binsiiy);

  eest.StopMapping();
  eest.Add("HillasSrc");

  //
  //  Setup tasklists
  //
  tlist2.AddToList(&read);

  TString hcut2("MHadronness.fHadronness>");
  hcut2 += maxhadronness;
  MContinue cont(hcut2);

  tlist2.AddToList(&cont);

  tlist2.AddToList(&eest);
  // tlist2.AddToList(new MPrint("MMcEvt"));
  // tlist2.AddToList(new MPrint("MEnergyEst"));

  tlist2.AddToList(&hfille);
  tlist2.AddToList(&hfilli);
  tlist2.AddToList(&hfilleo);
  tlist2.AddToList(&hfillio);

  tlist2.AddToList(&hfille2);
  tlist2.AddToList(&hfilli2);
  tlist2.AddToList(&hfilleo2);
  tlist2.AddToList(&hfillio2);

  tlist2.AddToList(&hfillee);
  tlist2.AddToList(&hfillii);

  //
  // Create and setup the eventloop
  //
  MProgressBar bar;

  MEvtLoop evtloop2;
  evtloop2.SetProgressBar(&bar);
  evtloop2.SetParList(&parlist);

  //
  // Execute your analysis
  //
  if (!evtloop2.Eventloop())
    return;

  tlist2.PrintStatistics();

  const TString text = "\\sqrt{<y>}=%.0f%%";

  char txt[1000];

  TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
  c->Divide(2,2);
  c->cd(1);
  mh3i.DrawClone("PROFXnonew");
  sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
  TLatex *t = new TLatex(180, 0.15, txt);
  t->Draw();
  c->cd(2);
  mh3e.DrawClone("PROFXnonew");
  sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
  t = new TLatex(3.5, 0.6, txt);
  t->Draw();
  c->cd(3);
  mh3io.DrawClone("PROFXnonew");
  c->cd(4);
  mh3eo.DrawClone("PROFXnonew");

  c=new TCanvas("Est2", "Estimates vs. E_{est}");
  c->Divide(2,2);
  c->cd(1);
  mh3i2.DrawClone("PROFXnonew");
  c->cd(2);
  mh3e2.DrawClone("PROFXnonew");
  c->cd(3);
  mh3io2.DrawClone("PROFXnonew");
  c->cd(4);
  mh3eo2.DrawClone("PROFXnonew");

  mhe.DrawClone("PROFX");
  mhi.DrawClone("PROFX");
}
