source: trunk/MagicSoft/Mars/macros/CT1EnergyEst.C@ 1893

/* ======================================================================== *\
!
! *
! * 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,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Abelardo Moralejo, 03/2003 <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//
// fcn is the function to be minimized using TMinuit::Migrad
//
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=0.22, Float_t maxdist = 1.)
{
  //
  // Upper cut in Maxdist taken from D. Kranich's Ph D.
  // We have to convert maxdist from degrees to mm:
  //

  MGeomCamCT1 gct1;
  maxdist /= gct1.GetConvMm2Deg();

  //
  // Fill events into a MHMatrix
  //
  MParList parlist;
  MHMatrix matrix;

  //
  // colenergy and colimpact are the indexes of the matrix columns containing 
  // the MC energy and impact parameter.
  //
  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
  Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact");

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

  //
  // Here we read in the events for the training. Second argument of AddFile is
  // the number of events. With 2000 it is rather fast, but only the lowest zenith
  // angles will enter in the minimization.
  //
  MReadTree read("Events");
  read.AddFile("MC_ON2.root",2000);
  read.DisableAutoScheme();

  //
  // Filter to keep the most gamma-like events:
  //
  TString hcut("MHadronness.fHadronness<");
  hcut += maxhadronness;
  hcut += " && HillasSrc.fDist < ";
  hcut += maxdist;

  MF filterhadrons(hcut);

  //
  // Fill the matrix used later in the minimization loop by Minuit:
  //
  if (!matrix.Fill(&parlist, &read, &filterhadrons))
    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, colenergy));
  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 of the parameters (close enough to the final ones, taken
  // from previous runs of the procedure). Parameter fA[4] is not used in the default
  // energy estimation model (from D. Kranich).
  //
  TArrayD fA(5);  
  TArrayD fB(7);

  fA[0] =  21006.2;
  fA[1] = -43.2648;
  fA[2] = -690.403;
  fA[3] = -0.0428544;
  fA[4] =  1.;
  fB[0] = -3391.05;
  fB[1] =  136.58;
  fB[2] =  0.253807;
  fB[3] =  254.363;
  fB[4] =  61.0386;
  fB[5] = -0.0190984;
  fB[6] = -421695;

  //
  // 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:

  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 << endl; 
        //      "  +-  " << 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();
  */

  eest.StopMapping();


  //
  // Now write the parameters of the energy estimator to file:
  //
  TVector* EnergyParams = new TVector(12);
  for (Int_t i=0; i<eest.GetNumCoeff(); i++)
    EnergyParams(i) = eest.GetCoeff(i);
  TFile out("energyest.root", "recreate");
  EnergyParams->Write("EnergyParams");
  out.Close();

  //
  // End of Part 1 (estimation of the parameters)
  //


  ////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////


  //
  // 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;
  MParList  parlist2;  

  parlist2.AddToList(&tlist2);

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

  MReadTree read2("Events");
  read2.AddFile("MC_ON2.root");
  read2.DisableAutoScheme();

  //
  // Create some histograms to display the results:
  //

  MH3 mh3ed("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
  MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
  MH3 mhe("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)");

  MFillH hfilled(&mh3ed);
  MFillH hfilled2(&mh3ed2);
  MFillH hfillee(&mhe);

  mhe.SetName("HistEE");
  mh3ed.SetName("HistEnergyDelta");
  mh3ed2.SetName("HistEnergyDelta");

  MBinning binsedx("BinningHistEnergyDeltaX");
  MBinning binsedy("BinningHistEnergyDeltaY");
  MBinning binseex("BinningHistEEX");
  MBinning binseey("BinningHistEEY");

  binsedx.SetEdges(25, 2.5, 5.);
  binsedy.SetEdges(30, -1., 1.);
  binseex.SetEdges(25, 2.5, 5.);
  binseey.SetEdges(30, 2., 5.);

  parlist2.AddToList(&binsedx);
  parlist2.AddToList(&binsedy);
  parlist2.AddToList(&binseex);
  parlist2.AddToList(&binseey);

  //
  // Create energy estimator task, add necessary containers and 
  // initialize parameters read from file:
  //

  MEnergyEstParam eest2("Hillas");
  eest2.Add("HillasSrc");

  TFile enparam("energyest.root");
  TArrayD parA(5);
  TArrayD parB(7);

  for (Int_t i=0; i<parA.GetSize(); i++)
    parA[i] = EnergyParams(i);
  for (Int_t i=0; i<parB.GetSize(); i++)
    parB[i] = EnergyParams(i+parA.GetSize());

  eest2.SetCoeffA(parA);
  eest2.SetCoeffB(parB);

  enparam.Close();

  //
  //  Setup tasklists
  //

  tlist2.AddToList(&read2);

  //
  // Include filter on hadronness and maximum value of dist:
  //
  TString hcut2("MHadronness.fHadronness>");
  hcut2 += maxhadronness;
  hcut2 += "|| HillasSrc.fDist > ";
  hcut2 += maxdist;

  MContinue cont(hcut2);

  tlist2.AddToList(&cont);
  tlist2.AddToList(&eest2);

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

  tlist2.AddToList(&hfilled);
  tlist2.AddToList(&hfilled2);
  tlist2.AddToList(&hfillee);

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

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

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

  //
  // Plot results:
  //
  mhe.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
  mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)");
  mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2);

  TCanvas *c = new TCanvas("energy","CT1 Energy estimation");
  c->Divide(2,2);
  c->cd(1);
  energy_1->SetBottomMargin(0.12);
  mhe.DrawClone("PROFXnonew");

  mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
  mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
  mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2);
  mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5);
  c->cd(2);
  energy_2->SetLeftMargin(0.15);
  energy_2->SetBottomMargin(0.12);
  mh3ed.DrawClone("PROFXnonew");

  mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)");
  mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
  mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2);
  mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5);
  c->cd(3);
  energy_3->SetLeftMargin(0.15);
  energy_3->SetBottomMargin(0.12);
  mh3ed2.DrawClone("PROFXnonew");

  ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae");

  c->cd(4);
  energy_4->SetLeftMargin(0.1);
  energy_4->SetRightMargin(0.05);
  energy_4->SetBottomMargin(0.12);
  deltae.GetXaxis()->SetTitleOffset(1.2);
  deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}");
  deltae.Draw("e");
  
  TFile out2("energyest.root", "update");
  ((TH2*)mh3ed.GetHist())->Write("mh3ed");
  ((TH2*)mh3ed2.GetHist())->Write("mh3ed2");
  ((TH2*)mhe.GetHist())->Write("mhe");
  deltae.Write();
  out2.Close();

  return;

}