Index: /trunk/MagicSoft/Mars/Changelog
===================================================================
--- /trunk/MagicSoft/Mars/Changelog	(revision 1851)
+++ /trunk/MagicSoft/Mars/Changelog	(revision 1852)
@@ -1,3 +1,18 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/03/21: Thomas Bretz
+    * manalysis/MEnergyEstParam.[h,cc]
+      - Added StopMapping and Print functions.
+
+ 2003/03/21: Abelardo Moralejo
+
+    * mhist/MHMatrix.[h,cc]:
+      - Added third argument (a filter) to the second instantiation 
+	of the Fill procedure.
+
+    * macros/CT1EnergyEst.C:
+      - Example of the parameter calculation and use of the energy 
+	estimation method for CT1.
+
 
  2003/03/21: Thomas Bretz
Index: /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C
===================================================================
--- /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C	(revision 1852)
+++ /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C	(revision 1852)
@@ -0,0 +1,403 @@
+/* ======================================================================== *\
+!
+! *
+! * 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");
+}
Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc
===================================================================
--- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc	(revision 1851)
+++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc	(revision 1852)
@@ -177,5 +177,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the four coefficients for the estimation of the impact parameter.
+// Set the five coefficients for the estimation of the impact parameter.
 // Argument must ba a TArrayD of size 5.
 //
@@ -193,5 +193,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the four coefficients for the estimation of the energy.
+// Set the seven coefficients for the estimation of the energy.
 // Argument must ba a TArrayD of size 7.
 //
@@ -209,5 +209,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the four coefficients for the estimation of impact and energy.
+// Set the twelve coefficients for the estimation of impact and energy.
 // Argument must ba a TArrayD of size 12.
 //
@@ -261,4 +261,12 @@
 }
 
+void MEnergyEstParam::StopMapping()
+{
+    fMatrix = NULL; 
+    fPairs->Clear();
+    fHillasSrc->Clear();
+    fEnergy->Clear();
+}
+
 // --------------------------------------------------------------------------
 //
@@ -272,4 +280,12 @@
 
     AddToBranchList(hillas+".fDist");
+}
+
+void MEnergyEstParam::Print(Option_t *opt)
+{
+    for (int i=0; i<fA.GetSize(); i++)
+      *fLog << all << "fA[" << i << "]=" << fA[i] << endl;
+    for (int i=0; i<fB.GetSize(); i++)
+      *fLog << all << "fB[" << i << "]=" << fB[i] << endl;
 }
 
@@ -295,10 +311,10 @@
 
     /* MY PARAM */
-    const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*width;
-    const Double_t e2 = fB[2] + fB[5]*size*width;
+    //const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*width;
+    //const Double_t e2 = fB[2] + fB[5]*size*width;
 
     /* MARCOS */
-    //const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*length;
-    //const Double_t e2 = fB[2] + fB[5]*length;
+    const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*length;
+    const Double_t e2 = fB[2] + fB[5]*length;
 
     TIter NextH(fHillasSrc);
@@ -317,11 +333,11 @@
 
         /* MARCOS */
-        //const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm]
-        /*const*/// Double_t er = e0 * (e1 + e2*ir);      // [GeV]
+        const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm]
+        /*const*/ Double_t er = e0 * (e1 + e2*0/*ir*/);      // [GeV]
 
         /* MY PARAM */
         // if (width==0) return kCONTINUE;
-        const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm]
-        Double_t er = e0 * (e1 + e2*ir);      // [GeV]
+        //const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm]
+        //Double_t er = e0 * (e1 + e2*ir);      // [GeV]
 
         /* MKA */
@@ -339,5 +355,5 @@
         }
 
-        est->SetEnergy(er);
+	est->SetEnergy(er);
         est->SetImpact(ir);
         est->SetReadyToSave();
Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h
===================================================================
--- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h	(revision 1851)
+++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h	(revision 1852)
@@ -50,4 +50,5 @@
 
     void InitMapping(MHMatrix *mat);
+    void StopMapping();
 
     Int_t GetNumCoeff() const { return fA.GetSize()+fB.GetSize(); }
@@ -57,4 +58,6 @@
     void SetCoeffB(const TArrayD &arr);
 
+    void Print(Option_t *o=NULL);
+
     ClassDef(MEnergyEstParam, 0) // Task to estimate the energy
 };
Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.cc
===================================================================
--- /trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1851)
+++ /trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1852)
@@ -68,4 +68,6 @@
 #include "MData.h"
 #include "MDataArray.h"
+#include "MF.h"
+
 
 ClassImp(MHMatrix);
@@ -600,5 +602,5 @@
 // --------------------------------------------------------------------------
 //
-Bool_t MHMatrix::Fill(MParList *plist, MTask *read)
+Bool_t MHMatrix::Fill(MParList *plist, MTask *read, MF *filter)
 {
     //
@@ -614,5 +616,14 @@
 
     tlist.AddToList(read);
+
+
+    if (filter)
+      {
+	tlist.AddToList(filter);
+	fillh.SetFilter(filter);
+      }
+
     tlist.AddToList(&fillh);
+
 
     MEvtLoop evtloop;
Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.h
===================================================================
--- /trunk/MagicSoft/Mars/mhist/MHMatrix.h	(revision 1851)
+++ /trunk/MagicSoft/Mars/mhist/MHMatrix.h	(revision 1852)
@@ -22,4 +22,5 @@
 class MParList;
 class MDataArray;
+class MF;
 
 class MHMatrix : public MH
@@ -92,5 +93,5 @@
     Double_t operator[](Int_t col) { return fM(fRow, col); }
 
-    Bool_t Fill(MParList *plist, MTask *read);
+    Bool_t Fill(MParList *plist, MTask *read, MF *filter=0);
 
     TString GetDataMember() const;
Index: /trunk/MagicSoft/Mars/mhist/Makefile
===================================================================
--- /trunk/MagicSoft/Mars/mhist/Makefile	(revision 1851)
+++ /trunk/MagicSoft/Mars/mhist/Makefile	(revision 1852)
@@ -23,5 +23,5 @@
 #
 INCLUDES = -I. -I../mbase -I../mraw -I../manalysis -I../mmc \
-	   -I../mgui -I../mgeom -I../mdata
+	   -I../mgui -I../mgeom -I../mdata -I../mfilter
 
 #------------------------------------------------------------------------------