Index: trunk/MagicSoft/Mars/macros/CT1EnergyEst.C
===================================================================
--- trunk/MagicSoft/Mars/macros/CT1EnergyEst.C	(revision 1852)
+++ trunk/MagicSoft/Mars/macros/CT1EnergyEst.C	(revision 1886)
@@ -16,5 +16,6 @@
 !
 !
-!   Author(s): Thomas Bretz et al,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   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
@@ -23,9 +24,12 @@
 \* ======================================================================== */
 
+//
+// 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();
+  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
 
   MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
@@ -41,7 +45,14 @@
 //------------------------------------------------------------------------
 
-void CT1EnergyEst(Float_t maxhadronness=1.)
+void CT1EnergyEst(Float_t maxhadronness=0.22, Float_t maxdist = 1.)
 {
-  //  Bool_t evalenergy=kFALSE;
+  //
+  // 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
@@ -50,5 +61,10 @@
   MHMatrix matrix;
 
-  Int_t col = matrix.AddColumn("MMcEvt.fEnergy");
+  //
+  // 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");
@@ -56,14 +72,27 @@
   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", 200);
+  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);
 
-  if (!matrix.Fill(&parlist, &read))
+  //
+  // Fill the matrix used later in the minimization loop by Minuit:
+  //
+  if (!matrix.Fill(&parlist, &read, &filterhadrons))
     return;
 
@@ -76,8 +105,9 @@
   MMatrixLoop loop(&matrix);
 
-
   MChisqEval eval;
-  eval.SetY1(new MDataElement(&matrix, col));
+
+  eval.SetY1(new MDataElement(&matrix, colenergy));
   eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
+
   eval.SetOwner();
 
@@ -97,4 +127,5 @@
 
   // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
+
   if (minuit.SetErrorDef(1))
     {
@@ -104,23 +135,27 @@
 
   //
-  // 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;
-  
+  // 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);
-  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;
-  
+
+  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++)
     {
@@ -166,13 +201,13 @@
   minuit.SetObjectFit(&evtloop);
 
-   
   TStopwatch clock;
   clock.Start();
 
-  // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
+  // Now ready for minimization step:
+
   gLog.SetNullOutput(kTRUE);
   Bool_t rc = minuit.Migrad();
   gLog.SetNullOutput(kFALSE);
-
+  
   if (rc)
     {
@@ -199,5 +234,6 @@
         }
 
-      cout << i << ":  " << val << "  +-  " << er << endl;
+      cout << i << ":  " << val << endl; 
+	//	"  +-  " << er << endl;
     }
 
@@ -208,11 +244,32 @@
     minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
     minuit.mnprin(3, amin);
+    eest.Print();
   */
-  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.
   //
@@ -222,7 +279,8 @@
   //
 
-
   MTaskList tlist2;
-  parlist.Replace(&tlist2);
+  MParList  parlist2;  
+
+  parlist2.AddToList(&tlist2);
 
   //
@@ -231,122 +289,86 @@
   //
 
-
-  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");
+  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");
-  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");
+  mh3ed.SetName("HistEnergyDelta");
+  mh3ed2.SetName("HistEnergyDelta");
+
+  MBinning binsedx("BinningHistEnergyDeltaX");
+  MBinning binsedy("BinningHistEnergyDeltaY");
   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");
+
+  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(&read);
-
+
+  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(&eest);
+  tlist2.AddToList(&eest2);
+
   // 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(&hfilled);
+  tlist2.AddToList(&hfilled2);
   tlist2.AddToList(&hfillee);
-  tlist2.AddToList(&hfillii);
 
   //
@@ -357,5 +379,5 @@
   MEvtLoop evtloop2;
   evtloop2.SetProgressBar(&bar);
-  evtloop2.SetParList(&parlist);
+  evtloop2.SetParList(&parlist2);
 
   //
@@ -365,39 +387,53 @@
     return;
 
-  tlist2.PrintStatistics();
-
-  const TString text = "\\sqrt{<y>}=%.0f%%";
-
-  char txt[1000];
-
-  TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
+  //
+  // 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);
-  mh3i.DrawClone("PROFXnonew");
-  sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
-  TLatex *t = new TLatex(180, 0.15, txt);
-  t->Draw();
+  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);
-  mh3e.DrawClone("PROFXnonew");
-  sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
-  t = new TLatex(3.5, 0.6, txt);
-  t->Draw();
+  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);
-  mh3io.DrawClone("PROFXnonew");
+  energy_3->SetLeftMargin(0.15);
+  energy_3->SetBottomMargin(0.12);
+  mh3ed2.DrawClone("PROFXnonew");
+
+  ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae");
+
   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");
+  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;
+
 }
