source: trunk/MagicSoft/Mars/macros/CT1EEst.C@ 1978

/* ======================================================================== *\
!
! *
! * 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>
!              Wolfgang Wittek,   04/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */
//
//------------------------------------------------------------------------
//
// fcn calculates 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 CT1EgyEst()
{
  TString inPath        = "~magican/ct1test/wittek/marsoutput/optionA/";
  TString fileNameIn    = "MC2.root";

  TString outPath       = "~magican/ct1test/wittek/marsoutput/optionA/";
  TString energyParName = "energyest.root";
 
  TString hilName    = "MHillas";
  TString hilSrcName = "MHillasSrc";

  TString hadronnessName = "MHadronness";

  Int_t howMany = 2200;

  Float_t maxhadronness = 0.4;
  Float_t maxalpha      = 20.0;

  CT1EEst(inPath,  fileNameIn, 
            outPath, energyParName, 
            hilName, hilSrcName,    hadronnessName, 
            howMany, maxhadronness, maxalpha);
}



//------------------------------------------------------------------------
//
// Arguments of CT1EEst :
// ------------------------
//
// inPath,  fileNameIn       path and name of input file (MC gamma events)
// outPath, energyParName    path and name of file containing the parameters
//                           of the energy estimator
// hilName, hilSrcName       names of "MHillas" and "MHillasSrc" containers
// hadronnessName            name of container holding the hadronness
// howMany                   no.of events to be read from input file
// maxhadronness             maximum value of hadronness to be accepted
// maxalpha                  maximum value of |ALPHA| to be accepted
//
void CT1EEst(TString inPath,  TString fileNameIn, 
               TString outPath, TString energyParName,
               TString hilName, TString hilSrcName, TString hadronnessName,
               Int_t howMany, 
               Float_t maxhadronness, Float_t maxalpha)
{
  cout << "================================================================"
       << endl;
  cout << "Start of energy estimation part" << endl;
  cout << "" << endl;
  cout << "CT1EEst input values : " << endl;
  cout << "---------------------- " << endl;
  cout << "inPath, fileNameIn = " 
       <<  inPath << ",  " << fileNameIn << endl;
  cout << "outPath, energyParName = " 
       <<  outPath << ",  " << energyParName << endl;
  cout << "hilName, hilSrcName, hadronnessName = " << hilName << ",  "
       << hilSrcName << ",  " << hadronnessName << endl;
  cout << "howMany, maxhadronness, maxalpha = " << howMany << ",  "
       << maxhadronness << ",  " << maxalpha << endl;
  cout << "" << endl;


  // convert from [deg] to [mm]
  //
  MGeomCamCT1 gct1;
  //maxdist /= gct1.GetConvMm2Deg();


  //==========================================================================
  // Start of Part 1 (determination of the parameters of the energy estimator)
  //

  //-----------------------------------------------------------------------
  // Fill events into a MHMatrix, 
  // to be used later in the minimization by MINUIT
  //

  MParList parlist;

  TString inputfile(outPath);
  inputfile += fileNameIn;

  MFEventSelector selector;
  selector.SetNumSelectEvts(howMany);

  cout << "Read events from file '" << inputfile << "'" << endl;    
  MReadTree read("Events");
  read.AddFile(inputfile);
  read.DisableAutoScheme();
  read.SetSelector(&selector);

  MFCT1SelFinal filterhadrons;
  filterhadrons.SetCuts(maxhadronness, maxalpha);
  filterhadrons.SetHadronnessName(hadronnessName);
  filterhadrons.SetInverted();

  cout << "Define columns of matrix" << endl;
  MHMatrix matrix;
  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
  Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact");

  if (colenergy < 0  ||  colimpact < 0)
  {
    cout << "colenergy, colimpact = " << colenergy << ",  " 
         << colimpact << endl;
    return;
  }

  MEnergyEstParam eest(hilName);
  eest.Add(hilSrcName);
  eest.InitMapping(&matrix);

  cout << "--------------------------------------" << endl;
  cout << "Fill events into the matrix" << endl;
  if ( !matrix.Fill(&parlist, &read, &filterhadrons) )
    return;
  cout << "Matrix was filled with " << matrix.GetNumRows() 
       << " events" << endl;  

  //-----------------------------------------------------------------------
  //
  // Setup the eventloop which will be executed in the fcn of MINUIT 
  //
  cout << "--------------------------------------" << endl;
  cout << "Setup event loop to be used in 'fcn'" << endl;

  MParList  parlist;
  MTaskList tasklist;

  MMatrixLoop loop(&matrix);


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

  //********************************
  // Entries in MParList

  parlist.AddToList(&tasklist);

  //********************************
  // Entries in MTaskList

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

  //********************************

  cout << "Event loop was setup" << endl; 
  MEvtLoop evtloop;
  evtloop.SetParList(&parlist);



  //..........................................................................

  //----------   Start of minimization part   -------------------- 
  //
  // Do the minimization with MINUIT
  //
  // Be careful: This is not thread safe
  //
  cout << "--------------------------------------" << endl;
  cout << "Start minimization" << endl;

  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();
  cout << "Minimization finished" << endl;

  //----------   End of minimization part   -------------------- 

  //..........................................................................


  //
  // Now write the parameters of the energy estimator to a file
  //
  TVector* EnergyParams = new TVector(12);

  for (Int_t i=0; i<eest.GetNumCoeff(); i++)
    EnergyParams(i) = eest.GetCoeff(i);

  TString paramout(outPath);
  paramout += energyParName;

  TFile outfile(paramout, "RECREATE");
  EnergyParams->Write("EnergyParams");
  outfile.Close();
 
  cout << "--------------------------------------" << endl;
  cout << "Write parameters of energy estimator onto file '" 
       << paramout << endl;
  cout << "--------------------------------------" << endl;
  //
  // End of Part 1 (determination of the parameters of the energy estimator)
  //=========================================================================


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


  //==========================================================================
  // Start of Part 2 (test quality of energy estimation)
  //
  //

  //--------------------------------------------
  // Read the parameters of the energy estimator 
  //

  TString paramin(outPath);
  paramin += energyParName;

  cout << "--------------------------------------" << endl;
  cout << "Read parameters of energy estimator from file '" 
       << paramin << endl;
  TFile enparam(paramin);

  //=======================================================================
  // Setup the event loop
  //
  cout << "--------------------------------------" << endl;
  cout << "Setup event loop for checking quality of energy estimation" << endl;


  MTaskList tlist2;
  MParList  parlist2;  

  //-----------------------------------------------
  // Read events
  //

  TString inputfile(inPath);
  inputfile += fileNameIn;

  cout << "Read events from file '" << inputfile << "'" << endl;    
  MReadTree read2("Events");
  read2.AddFile(inputfile);
  read2.DisableAutoScheme();


  //-----------------------------------------------
  // Select events
  //
  cout << "Select events with hadronness < " << maxhadronness 
       << " and |alpha| < " << maxalpha << endl; 
  MFCT1SelFinal hcut2;
  hcut2.SetCuts(maxhadronness, maxalpha);
  hcut2.SetHadronnessName(hadronnessName);

  MContinue cont(&hcut2);


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

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

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

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

  binsedx.SetEdges(35, 2.0, 5.5);
  binsedy.SetEdges(40,-1.5, 2.5);

  binseex.SetEdges(35, 2.0, 5.5);
  binseey.SetEdges(35, 2.0, 5.5);

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


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

  MEnergyEstParam eest2(hilName);
  eest2.Add(hilSrcName);

  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();



  //********************************
  // Entries in MParList

  parlist2.AddToList(&tlist2);

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


  //********************************
  // Entries in MTaskList

  tlist2.AddToList(&read2);

  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);

  //********************************

  cout << "Event loop was setup" << endl; 
  MProgressBar bar;

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

  if (!evtloop2.Eventloop())
    return;

  TH1& histed = mh3ed.GetHist();
  Float_t  meany = histed.GetMean(2);
  Float_t  RMSy  = histed.GetRMS(2);
  Float_t resolution = sqrt( meany*meany + RMSy*RMSy );
  cout << "" << endl;
  cout << "With  y = (E_EST -E_MC)/E_MC one obtains :" << endl;
  cout << "      Average y = "        << histed.GetMean(2) 
       << ",  RMS = "                 << histed.GetRMS(2)
       << ",  Average resolution = "  << resolution << endl;
  cout << "" << endl;

  //=======================================================================

  //
  // Plot results:
  //
  cout << "--------------------------------------" << endl;
  cout << "Plot the results" << endl;

  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");
  
  //-------------------------------------------
  //
  // include the histograms 
  // into the root file containing the parameters of the energy estimator
  //

  
  TString paramout(outPath);
  paramout += energyParName;

  TFile out2(paramout, "UPDATE");
  ((TH2*)mh3ed.GetHist())->Write("mh3ed");
  ((TH2*)mh3ed2.GetHist())->Write("mh3ed2");
  ((TH2*)mhe.GetHist())->Write("mhe");
  deltae.Write();
  out2.Close();

  cout << "Quality histograms were added onto the file '" << paramout << endl;

  cout << "" << endl;
  cout << "End of energy estimation part" << endl;
  cout << "================================================================"
       << endl;
  //
  // End of Part 2 (test quality of energy estimation)
  //==========================================================================


  return;

}
//============================================================================