Changeset 7413 for trunk/MagicSoft

Timestamp:

11/21/05 11:09:12 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft

Files:

• Cosy/tpoint/tpoint0509.txt (modified) (4 diffs)
• Cosy/tpoint/tpoint0510.txt (modified) (1 diff)
• Cosy/tpoint/tpoint0511.txt (added)
• Mars/Changelog (modified) (1 diff)
• Mars/Makefile (modified) (1 diff)
• Mars/NEWS (modified) (1 diff)
• Mars/ganymed.rc (modified) (5 diffs)
• Mars/macros/dohtml.C (modified) (1 diff)
• Mars/macros/rootlogon.C (modified) (1 diff)
• Mars/macros/train (added)
• Mars/macros/train/traindisp.C (added)
• Mars/macros/train/trainenergy.C (added)
• Mars/macros/train/trainseparation.C (added)
• Mars/mbase/MParList.cc (modified) (3 diffs)
• Mars/mbase/MParList.h (modified) (1 diff)
• Mars/mhbase/MBinning.cc (modified) (4 diffs)
• Mars/mhbase/MBinning.h (modified) (1 diff)
• Mars/mhbase/MH3.cc (modified) (1 diff)
• Mars/mhbase/MHMatrix.cc (modified) (1 diff)
• Mars/mhbase/MHMatrix.h (modified) (1 diff)
• Mars/mhist/MHHadronness.cc (modified) (1 diff)
• Mars/mjobs/MDataSet.cc (modified) (1 diff)
• Mars/mjobs/MDataSet.h (modified) (1 diff)
• Mars/mjobs/MJCut.cc (modified) (5 diffs)
• Mars/mjobs/MJCut.h (modified) (2 diffs)
• Mars/mjoptim/MJOptimizeDisp.cc (modified) (5 diffs)
• Mars/mranforest/MHRanForestGini.cc (modified) (7 diffs)
• Mars/mranforest/MHRanForestGini.h (modified) (2 diffs)
• Mars/mranforest/MRanForest.cc (modified) (14 diffs)
• Mars/mranforest/MRanForestCalc.cc (modified) (5 diffs)
• Mars/mranforest/MRanForestCalc.h (modified) (2 diffs)
• Mars/mranforest/MRanTree.h (modified) (3 diffs)
• Mars/mranforest/Makefile (modified) (1 diff)
• Mars/mranforest/RanForestLinkDef.h (modified) (1 diff)
• Mars/mtools/MTFillMatrix.cc (modified) (4 diffs)
• Mars/mtools/MTFillMatrix.h (modified) (3 diffs)

trunk/MagicSoft/Cosy/tpoint/tpoint0509.txt

@@ -1 +1 @@
 Magic Model  TPOINT data file
 : ALTAZ
-#
-# ***** After a refocussing of the mirror *****
-#
 
 ### 24.9.: tpoints always with the white paper screen
@@ -11 +8 @@
 1.742509 57.99808 144.2056 49.02492 0.945 60.71667 0.0005732658 0.007748751 53637.075763 216.1 2.654
 # Hamal,  mag=2.00, Zd=5
-139.6868 83.23968 281.6764 74.22004 2.119444 23.4625 -0.04332626 0.03920611 53637.11646 225.2 2.635
+##139.6868 83.23968 281.6764 74.22004 2.119444 23.4625 -0.04332626 0.03920611 53637.11646 225.2 2.635
 
 # --- 24.09.2005 --- 02:59:09.399
 # Hamal,  mag=2.00, Zd=5
-161.5844 84.46453 303.3778 75.43595 2.119444 23.4625 -0.009743542 0.1953913 53637.124414 225.2 2.665
+##161.5844 84.46453 303.3778 75.43595 2.119444 23.4625 -0.009743542 0.1953913 53637.124414 225.2 2.665
 # And-51, mag=3.57, Zd=21
 -21.2567 68.112 121.1535 59.26484 1.633333 48.62833 0.02363187 0.04614806 53637.142328 218.1 3.353
@@ -24 +21 @@
 # --- 25.09.2005 --- 03:41:32.305
 # Nothing in the runbook?
--35.99076 37.46571 106.0467 28.7287 22.48611 58.41528 0.1168588 0.4782401 53638.153846 0 6.58
+##-35.99076 37.46571 106.0467 28.7287 22.48611 58.41528 0.1168588 0.4782401 53638.153846 0 6.58
 -38.29146 35.14824 104.1974 26.30449 22.25056 57.04361 0.01046067 0.02796104 53638.15701 220.4 3.97
 ## seems to be wrongly identified
@@ -76 +73 @@
 112.9061  44.22773 255.4126 35.24166 7.655     5.225    -0.005348532  0.000207558 53639.239228 215.7 2.556
 
+####################################################################
+# Private communication with Adrian:
+# Als wir Ende September rausfanden, dass die Laser unzuverlaessig
+# sind, wurde da eine Quick-and-Dirty Kalibration gemacht (da kein
+# Experte auf La Palma war)
+####################################################################
+

trunk/MagicSoft/Cosy/tpoint/tpoint0510.txt

@@ -2 +2 @@
 : ALTAZ
 
+#
+# ***** New focussing after collabortaion meeting in Tenerife (19.10.?) *****
+#
+
 # --- 25.10.2005 --- 02:11:32.993
+# Diphda. Magn: 2.04 Zd:59.21 Az:222.64
 -137.4972 30.88874 5.015726 21.99765 0.7263889 -17.98667 -0.03588449 -0.01783158 53668.091354 237.2 2.802
+# Algenib. Magn:2.83 Zd:45.55 Az:263.21
 -96.64771 44.18443 45.84061 35.33164 0.2205556 15.18361 -0.0004117089 -0.002440412 53668.098033 235.3 3.151
+# Markab. Magn:2.49 Zd:65.51 Az:273.02
 -87.16777 27.87191 55.32808 19.00867 23.07944 15.20528 -0.008018477 0.009832391 53668.10233 237.3 3.096
+# Markab. Magn:2.49 Zd:63.7 Az:273.65
 -86.52195 26.6309 55.98938 17.74595 23.07944 15.20528 -0.02510043 -0.001211299 53668.106316 233.3 3.102
+

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2005/11/21 Thomas Bretz
+
+   * mjtrain/*, macros/train/*.C:
+     - added new class to train the random forest
+
+   * Makefile:
+     - added mjtrain
+
+   * macros/dohtml.C, macros/rootlogon.C:
+     - added mjoptim and mjtrain
+
+   * mbase/MParList.[h,cc]:
+     - fixed copy constructor (it was crashing due to the
+       non-initialized lists)
+
+   * mhbase/MBinning.[h,cc]:
+     - renamed SetEdges to SetEdgesLin (SetEdges now is just and
+       abbreviation)
+     - added a check for the number of bins to all SetEdges*
+
+   * mhbase/MH3.cc:
+     - enable grid
+
+   * mhbase/MHMatrix.[h,cc]:
+     - new member function Addcolumns taking a TCollection
+
+   * mhist/MHHadronness.cc:
+     - align hadronnes into [0,1]
+
+   * mjobs/MDataSet.[h,cc]:
+     - call SetOwner for fSequencesOn and fSequencesOff
+     - added Copy member function
+
+   * mjobs/MJCut.[h,cc]:
+     - added new TaskEnv "CalcDisp" for a disp estimator
+
+   * mjoptim/MJOptimizeDisp.cc:
+     - fixed typo in a title
+     - added a new histogram showing theta-sq versus size
+
+   * mranforest/MHRanForestGini.[h,cc]:
+     - revised to display more information
+     - exchanged pointers to TGraph by objects
+
+   * mranforest/MRanForest.[h,cc]:
+     - replaced some weired copy of the train matrix by a direct access
+     - revised output
+     - pipe error/resolution to tree
+
+   * mranforest/MRanForestCalc.[h,cc]:
+     - copied code from MRFEnergyEst
+     - allow to set a name for the output container
+     - set numtry to 0 (auto)
+     - set ndsize to 1 (there is no auto mode)
+     - introduced variable for number of obsolete parameters
+
+   * mranforest/MRanTree.h:
+     - new data member to store resolution/error
+
+   * mranforest/Makefile, mranforest/RanForestLinkDef.h:
+     - removed MRFEnergyEst 
+     - added MRanForestCalc
+
+   * mtools/MTFillMatrix.[h,cc]:
+     - added the possibility to use PreCuts befre filling the matrix
+     - added ReadEnv
+
+
+
  2005/11/18 Thomas Bretz
 

trunk/MagicSoft/Mars/Makefile

@@ -67 +67 @@
           mjobs \
           mjoptim \
+          mjtrain \
           mtools \
           mmuon

trunk/MagicSoft/Mars/NEWS

@@ -13 +13 @@
    - general: MTFillMatrix (the class to fill one or two MHMatrix from
      files) now allows adding a pre-cut like in the optimization. E.g. this
-     is useful to perform g/h-separation cuts before training the random forest.
-
-   - general: Updated the random forest classes to support also the
-     regression method implemented by Thomas H.
-
-   - general: added new tutorial macro how to train the random forest
-     for energy estimation (macros/optim/rfenergyest.C)
+     is useful to perform g/h-separation cuts before training the random
+     forest.
+
+   - RanForest:
+      + Updated the random forest classes to support also the
+        regression method implemented by Thomas H.
+      + added new tutorial macro how to train the random forest
+        for energy estimation (macros/optim/rfenergyest.C)
+      + new classes to train the random forest (still in development)
+        mjtrain/MJTrainEnergy, mjtrain/MJTrainDisp, mjtrain/MJTrainSeparation
+      + new tutorial macros for random forest training in macros/train
+        trainenergy.C, traindisp.C, trainseparation.C
+
+   - ganymed: In addition to the Hadronness calculator (CalcHadronness)
+     a new option was implemented to estimate Disp (CalcDisp)
 
    - ganymed: Implemented two new options which allow

trunk/MagicSoft/Mars/ganymed.rc

@@ -39 +39 @@
 BinningDist.Raw:         25    0    2
 BinningMaxDist.Raw:      25    0    2
+
+# -------------------------------------------------------------------------
+# Using the starguider for pointing correction. 
+# To switch it off use "MPointingDevCalc.MaxAbsDev: -1"
+# -------------------------------------------------------------------------
+#MPointingDevCalc.MaxAbsDev:   15
+#MPointingDevCalc.NumMinStars: 8
+#MPointingDevCalc.NsbLevel:    3.0
+#MPointingDevCalc.NsbMin:      30
+#MPointingDevCalc.NsbMax:      60
+
+# -------------------------------------------------------------------------
+# Setup the misfocussing correction
+# -------------------------------------------------------------------------
+#MSrcPosCorrect.Dx: 0
+#MSrcPosCorrect.Dy: 0
 
 # -------------------------------------------------------------------------
@@ -83 +99 @@
 Cut0.Condition: ({0} || {1}) && {2} && {3} && {4} && {5}
 Cut0.0: MImagePar.fNumSatPixelsHG < 1
-Cut0.1: MHillas.GetArea*(MGeomCam.fConvMm2Deg^2) > (0.002*MImagePar.fNumSatPixelsHG) + 0.025
+Cut0.1: MHillas.GetArea*(MGeomCam.fConvMm2Deg^2) > (0.003*MImagePar.fNumSatPixelsHG) + 0.0325
 Cut0.2: MNewImagePar.fNumUsedPixels>5
 Cut0.3: MNewImagePar.fLeakage1 < 0.3
@@ -90 +106 @@
 
 # ---------- SETUP FOR ONOFF-MODE -----------
+MAlphaFitter.ScaleMin: 0.137
+MAlphaFitter.ScaleMax: 0.640
 MAlphaFitter.BackgroundFitMin: 0.137
-MAlphaFitter.BackgroundFitMax: 0.64
+MAlphaFitter.BackgroundFitMax: 1.000
 MAlphaFitter.PolynomOrder: 1
-MAlphaFitter.ScaleMode: Background
+MAlphaFitter.ScaleMode: OffRegion
 MAlphaFitter.SignalFunction: ThetaSq
 
@@ -99 +117 @@
 Cut1.ThetaCut: None
 Cut1.Param0:  1.3245
-Cut1.Param1:  0.189
-Cut1.Param2:  0.230
-Cut1.Param3:  5.320
-Cut1.Param4:  0.100
-Cut1.Param5: -0.0636
+Cut1.Param1:  0.204
+Cut1.Param2:  0.228
+Cut1.Param3:  5.50
+Cut1.Param4:  0.088
+Cut1.Param5: -0.07
 Cut1.Param6:  8.2957
 Cut1.Param7:  0.8677
@@ -137 +155 @@
 
 # -------------------------------------------------------------------------
-# Use this to executa a task (eg to calc hadronness) after energy
-# before Cut1
+# Use this to executa a task (eg to calc hadronness) before Cut1
 # -------------------------------------------------------------------------
 # CalcHadronness: MRanForestCalc
-# CalcHadronness.File: /home/tbretz/Mars.cvs/RFspot3cmAll.root
+# CalcHadronness.File: rf-energy.root
+# CalcHadronness.NameOutput: Hadronness
 # CalcHadronness.Debug: No
+
+# -------------------------------------------------------------------------
+# Use this to executa a task (eg to calc disp) before Cut1
+# -------------------------------------------------------------------------
+# CalcDisp: MRanForestCalc
+# CalcDisp.File: rf-disp.root
+# CalcDisp.NameOutput: Disp
+# CalcDisp.Debug: No
+# Cut1.CalcDisp: No

trunk/MagicSoft/Mars/macros/dohtml.C

@@ -77 +77 @@
     sourcedir += "mimage:";
     sourcedir += "mjobs:";
+    sourcedir += "mjoptim:";
+    sourcedir += "mjtrain:";
     sourcedir += "mmain:";
     sourcedir += "mmc:";

trunk/MagicSoft/Mars/macros/rootlogon.C

@@ -144 +144 @@
     gInterpreter->AddIncludePath(dir+"mimage");
     gInterpreter->AddIncludePath(dir+"mjobs");
+    gInterpreter->AddIncludePath(dir+"mjoptim");
+    gInterpreter->AddIncludePath(dir+"mjtrain");
     gInterpreter->AddIncludePath(dir+"mmain");
     gInterpreter->AddIncludePath(dir+"mmc");

trunk/MagicSoft/Mars/mbase/MParList.cc

@@ -62 +62 @@
 // --------------------------------------------------------------------------
 //
-//  default constructor
 //  creates an empty list
 //
-MParList::MParList(const char *name, const char *title)
+void MParList::Init(const char *name, const char *title)
 {
     fName  = name  ? name  : gsDefName.Data();
@@ -83 +82 @@
 }
 
+
+// --------------------------------------------------------------------------
+//
+//  default constructor
+//  creates an empty list
+//
+MParList::MParList(const char *name, const char *title)
+{
+    Init(name, title);
+}
+
 // --------------------------------------------------------------------------
 //
@@ -90 +100 @@
 //  entries)
 //
-MParList::MParList(MParList &ts)
-{
+MParList::MParList(const MParList &ts, const char *name, const char *title)
+{
+    Init(name, title);
+
     fContainer->AddAll(ts.fContainer);
 }

trunk/MagicSoft/Mars/mbase/MParList.h

@@ -37 +37 @@
     void StreamPrimitive(ofstream &out) const;
 
+    void Init(const char *name, const char *title);
+
 public:
     enum { kDoNotReset = BIT(17), kIsProcessing = BIT(18) };
 
     MParList(const char *name=NULL, const char *title=NULL);
-    MParList(MParList &ts);
+    MParList(const MParList &ts, const char *name=NULL, const char *title=NULL);
 
     virtual ~MParList();

trunk/MagicSoft/Mars/mhbase/MBinning.cc

@@ -242 +242 @@
 // lowest <lo> and highest <up> Edge (of your histogram)
 //
-void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up)
-{
+void MBinning::SetEdgesLin(Int_t nbins, Axis_t lo, Axis_t up)
+{
+    if (nbins<=0)
+    {
+        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
+        nbins = 1;
+    }
+
+    if (up<lo)
+    {
+        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
+
+        const Axis_t dummy(lo);
+        lo = up;
+        up = dummy;
+    }
+
     const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
     fEdges.Set(nbins+1);
@@ -402 +417 @@
 // lowest <lo> and highest <up> Edge (of your histogram)
 //
-void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up)
+void MBinning::SetEdgesLog(Int_t nbins, Axis_t lo, Axis_t up)
 {
     // if (lo==0) ...
+    if (nbins<=0)
+    {
+        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
+        nbins = 1;
+    }
+
+    if (up<lo)
+    {
+        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
+
+        const Axis_t dummy(lo);
+        lo = up;
+        up = dummy;
+    }
 
     const Double_t binsize = log10(up/lo)/nbins;
@@ -419 +448 @@
 // lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
 //
-void MBinning::SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up)
-{
+void MBinning::SetEdgesCos(Int_t nbins, Axis_t lo, Axis_t up)
+{
+    if (nbins<=0)
+    {
+        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
+        nbins = 1;
+    }
+
+    if (up<lo)
+    {
+        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
+
+        const Axis_t dummy(lo);
+        lo = up;
+        up = dummy;
+    }
+
     // if (lo==0) ...
     const Axis_t ld = lo/kRad2Deg;
@@ -436 +480 @@
 // lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
 //
-void MBinning::SetEdgesASin(const Int_t nbins, Axis_t lo, Axis_t up)
-{
+void MBinning::SetEdgesASin(Int_t nbins, Axis_t lo, Axis_t up)
+{
+    if (nbins<=0)
+    {
+        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
+        nbins = 1;
+    }
+
+    if (up<lo)
+    {
+        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
+
+        const Axis_t dummy(lo);
+        lo = up;
+        up = dummy;
+    }
+
     const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
     fEdges.Set(nbins+1);

trunk/MagicSoft/Mars/mhbase/MBinning.h

@@ -62 +62 @@
     void SetEdges(const MBinning &bins) { SetEdges(bins.fEdges); fType = bins.fType; fTitle = bins.fTitle; }
     void SetEdges(const TH1 &h, const Char_t axis='x');
-    void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up);
+    void SetEdges(Int_t nbins, const Axis_t lo, const Axis_t up) { SetEdgesLin(nbins, lo, up); }
     void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up, const char *opt);
-    void SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up);
-    void SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up);
-    void SetEdgesASin(const Int_t nbins, Axis_t lo, Axis_t up);
+    void SetEdgesLin(Int_t nbins, Axis_t lo, Axis_t up);
+    void SetEdgesLog(Int_t nbins, Axis_t lo, Axis_t up);
+    void SetEdgesCos(Int_t nbins, Axis_t lo, Axis_t up);
+    void SetEdgesASin(Int_t nbins, Axis_t lo, Axis_t up);
 
     Int_t FindLoEdge(Double_t val) const

trunk/MagicSoft/Mars/mhbase/MH3.cc

@@ -561 +561 @@
     TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(fHist);
     pad->SetBorderMode(0);
+    pad->SetGridx();
+    pad->SetGridy();
 
     fHist->SetFillStyle(4000);

trunk/MagicSoft/Mars/mhbase/MHMatrix.cc

@@ -209 +209 @@
 }
 
+void MHMatrix::AddColumns(const TCollection &list)
+{
+    TIter Next(&list);
+    TObject *obj = 0;
+    while ((obj=Next()))
+        AddColumn(obj->GetName());
+}
+
+
 // --------------------------------------------------------------------------
 //

trunk/MagicSoft/Mars/mhbase/MHMatrix.h

@@ -73 +73 @@
     Int_t AddColumn(const char *name);
     void AddColumns(MDataArray *mat);
+    void AddColumns(const TCollection &list);
 
     //MDataArray *GetColumns() { return fData; }

trunk/MagicSoft/Mars/mhist/MHHadronness.cc

@@ -209 +209 @@
     const MParameterD &had = par ? *(MParameterD*)par : *fHadronness;
 
-    const Double_t h = had.GetVal();
+    const Double_t h = TMath::Min(TMath::Max(had.GetVal(), 0.), 1.);
 
     if (TMath::IsNaN(h))

trunk/MagicSoft/Mars/mjobs/MDataSet.cc

@@ -196 +196 @@
     fName  = fname;
 
+    fSequencesOn.SetOwner();
+    fSequencesOff.SetOwner();
+
     const char *expname = gSystem->ExpandPathName(fname);
 

trunk/MagicSoft/Mars/mjobs/MDataSet.h

@@ -54 +54 @@
 
 public:
-    MDataSet() : fNumAnalysis((UInt_t)-1) { }
+    MDataSet(Int_t num=(UInt_t)-1) : fNumAnalysis(num) { }
     MDataSet(const char *fname, TString sequences="", TString data="");
+
+    void Copy(TObject &obj) const
+    {
+        MDataSet &ds = (MDataSet&)obj;
+        ds.fNumAnalysis = fNumAnalysis;
+        ds.fNumSequencesOn = fNumSequencesOn;
+        ds.fNumSequencesOff = fNumSequencesOff;
+        ds.fNameSource = fNameSource;
+        ds.fCatalog = fCatalog;
+        ds.fComment = fComment;
+        ds.fIsWobbleMode = fIsWobbleMode;
+
+        TObject *o=0;
+
+        ds.fSequencesOn.Delete();
+        ds.fSequencesOff.Delete();
+
+        TIter NextOn(&fSequencesOn);
+        while ((o=NextOn()))
+            ds.fSequencesOn.Add(o->Clone());
+
+        TIter NextOff(&fSequencesOff);
+        while ((o=NextOff()))
+            ds.fSequencesOff.Add(o->Clone());
+    }
 
     // Getter

trunk/MagicSoft/Mars/mjobs/MJCut.cc

@@ -86 +86 @@
     : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
     fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kTRUE),
-    fNameHist("MHThetaSq"), fCalcHadronness(0)
+    fNameHist("MHThetaSq"), fCalcHadronness(0), fCalcDisp(0)
 {
     fName  = name  ? name  : "MJCut";
@@ -102 +102 @@
     if (fCalcHadronness)
         delete fCalcHadronness;
+    if (fCalcDisp)
+        delete fCalcDisp;
 }
 
@@ -160 +162 @@
         delete fCalcHadronness;
     fCalcHadronness = task ? (MTask*)task->Clone() : 0;
+}
+
+// --------------------------------------------------------------------------
+//
+// Setup a task calculating disp. The given task is cloned.
+//
+void MJCut::SetDispCalculator(const MTask *task)
+{
+    if (fCalcDisp)
+        delete fCalcDisp;
+    fCalcDisp = task ? (MTask*)task->Clone() : 0;
 }
 
@@ -540 +553 @@
     taskenv2.SetDefault(fCalcHadronness);
 
+    MTaskEnv taskenv3("CalcDisp");
+    taskenv3.SetDefault(fCalcDisp);
+
     MFillH fill1a("MHHillasOffPre  [MHHillas]",      "MHillas",      "FillHillasPre");
     MFillH fill2a("MHHillasOffPost [MHHillas]",      "MHillas",      "FillHillasPost");
@@ -579 +595 @@
     //tlist2.AddToList(&taskenv1);
     tlist2.AddToList(&taskenv2);
+    tlist2.AddToList(&taskenv3);
     tlist2.AddToList(&cont0);
     if (write0)

trunk/MagicSoft/Mars/mjobs/MJCut.h

@@ -33 +33 @@
     //MTask *fEstimateEnergy;
     MTask *fCalcHadronness;
+    MTask *fCalcDisp;
 
     TString  GetOutputFile(UInt_t num) const;
@@ -68 +69 @@
     //void SetEnergyEstimator(const MTask *task=0);
     void SetHadronnessCalculator(const MTask *task=0);
+    void SetDispCalculator(const MTask *task=0);
 
     ClassDef(MJCut, 0) // Standard program to perform g/h-separation cuts

trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc

@@ -71 +71 @@
 
 // histograms
+#include "MH3.h"
+#include "MBinning.h"
 #include "../mhflux/MAlphaFitter.h"
 #include "../mhflux/MHThetaSq.h"
@@ -91 +93 @@
 //
 // Read all events from file which do match rules and optimize
-// energy estimator.
+// disp estimator.
 //
 Bool_t MJOptimizeDisp::RunDisp(const char *fname, const char *rule, MTask *weights)
 {
-    fLog->Separator("Preparing Energy optimization");
+    fLog->Separator("Preparing Disp optimization");
 
     MParList parlist;
@@ -124 +126 @@
     const Int_t num1 = m.AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
     const Int_t num2 = m.AddColumn("MHillasSrc.fAlpha*kDegToRad");
+    const Int_t num3 = m.AddColumn("MHillas.fSize");
 
     MHThetaSq hist;
     hist.SkipHistTime();
     hist.SkipHistTheta();
-    hist.SkipHistEnergy();
-    hist.InitMapping(&m); 
+    //hist.SkipHistEnergy();
+    //hist.ForceUsingSize();
+    hist.InitMapping(&m, 1);
 
     MFDataMember filter("DataType.fVal", '>', 0.5);
@@ -171 +175 @@
     MMatrixLoop loop(&m);
 
+    const char *n3   = Form("M[%d]", num3);
+    MH3 hdisp(n3, "sqrt(ThetaSquared.fVal)");
+    hdisp.SetTitle("\\vartheta^{2} distribution vs. Size:Size [phe]:\\vartheta^{2} [\\circ^{2}]");
+
+    MBinning binsx(100, 10, 100000, "BinningMH3X", "log");
+    MBinning binsy(100, 0,  2,      "BinningMH3Y", "lin");
+
+    parlist.AddToList(&binsx);
+    parlist.AddToList(&binsy);
+
+    MFillH fillh2(&hdisp);
+    fillh2.SetDrawOption("blue profx");
+
     tasklist.AddToList(&loop);
     tasklist.AddToList(&calc1);
@@ -177 +194 @@
         tasklist.AddToList(weights);
     tasklist.AddToList(&fill);
+    tasklist.AddToList(&fillh2);
     tasklist.AddToList(&eval);
 

trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc

@@ -43 +43 @@
 #include "MRanTree.h"
 #include "MRanForest.h"
+#include "MDataArray.h"
 
 #include "MLog.h"
@@ -57 +58 @@
 //
 MHRanForestGini::MHRanForestGini(Int_t nbins, const char *name, const char *title)
+    : fRules(0.01, 0.01, 0.99, 0.99)
 {
     //
@@ -64 +66 @@
     fTitle = title ? title : "Measure of importance of Random Forest-input parameters";
 
-    fGraphGini = new TGraph;
-    fGraphGini->SetTitle("Importance of RF-input parameters measured by mean Gini decrease");
-    fGraphGini->SetMarkerStyle(kFullDotSmall);
-}
-
-// --------------------------------------------------------------------------
-//
-// Delete the histograms.
-//
-MHRanForestGini::~MHRanForestGini()
-{
-    delete fGraphGini;
+    fGraphGini.SetNameTitle("Gini", "Importance of RF-input parameters measured by mean Gini decrease");
+    fGraphGini.SetMarkerStyle(kFullDotMedium);
+
+    fGraphError.SetNameTitle("ResErr", "Resolution/Error versus train step");
+    fGraphError.SetMarkerStyle(kFullDotMedium);
+
+    fGraphNodes.SetNameTitle("Nodes", "Number of nodes versus train step");
+    fGraphNodes.SetMarkerStyle(kFullDotMedium);
+
+    fRules.SetTextAlign(13);
+    fRules.SetTextSize(0.05);
 }
 
@@ -104 +105 @@
 Bool_t MHRanForestGini::Fill(const MParContainer *par, const Stat_t w)
 {
+    MRanTree *t = fRanForest->GetCurTree();
+
     for (Int_t i=0;i<fRanForest->GetNumDim();i++)
-        fGini[i]+=fRanForest->GetCurTree()->GetGiniDec(i);
+        fGini[i] += t->GetGiniDec(i);
+
+    fGraphError.SetPoint(fGraphError.GetN(), GetNumExecutions(), t->GetError());
+    fGraphNodes.SetPoint(fGraphError.GetN(), GetNumExecutions(), t->GetNumEndNodes());
 
     return kTRUE;
@@ -117 +123 @@
     const Int_t n = fGini.GetSize();
 
-    fGraphGini->Set(n);
-
-    Stat_t max=0;
-    Stat_t min=0;
+    fGraphGini.Set(n);
+
     for (Int_t i=0; i<n; i++)
     {
-        fGini[i] /= fRanForest->GetNumTrees();
-        fGini[i] /= fRanForest->GetNumData();
-
-        const Stat_t ip = i+1;
-        const Stat_t ig = fGini[i];
-
-        if (ig>max) max=ig;
-        if (ig<min) min=ig;
-
-        fGraphGini->SetPoint(i,ip,ig);
-    }
-
-    // This is used in root>3.04/? so that SetMaximum/Minimum can succeed
-    fGraphGini->GetHistogram();
-
-    fGraphGini->SetMaximum(1.05*max);
-    fGraphGini->SetMinimum(0.95*min);
+        fGini[i] /= fRanForest->GetNumTrees()*fRanForest->GetNumData();
+        fGraphGini.SetPoint(i, i+1, fGini[i]);
+    }
+
+    fRules.AddText("");
+    const MDataArray &arr = *fRanForest->GetRules();
+    int i;
+    for (i=0; i<arr.GetNumEntries(); i++)
+    {
+        TString s;
+        s += i+1;
+        s += ") ";
+        s += arr.GetRule(i);
+        fRules.AddText(s);
+    }
+    for (; i<20; i++)
+        fRules.AddText("");
 
     return kTRUE;
@@ -150 +154 @@
 void MHRanForestGini::Draw(Option_t *)
 {
-    if (fGraphGini->GetN()==0)
+    if (fGraphGini.GetN()==0)
         return;
 
@@ -158 +162 @@
     AppendPad("");
 
-    fGraphGini->Draw("ALP");
-    pad->Modified();
-    pad->Update();
-
-    TH1 *h = fGraphGini->GetHistogram();
-    if (!h)
-        return;
-
-    //h->GetXaxis()->SetRangeUser(0, fGini.GetSize()+1);
-    h->SetXTitle("No.of RF-input parameter");
-    h->SetYTitle("Mean decrease in Gini-index [a.u.]");
-
-    pad->Modified();
-    pad->Update();
-}
+    pad->Divide(2,2);
+
+    pad->cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    fGraphGini.Draw("ALP");
+
+    TH1 *h = fGraphGini.GetHistogram();
+    if (h)
+    {
+        h->SetXTitle("No.of RF-input parameter");
+        h->SetYTitle("Mean decrease in Gini-index [au]");
+        h->GetXaxis()->SetNdivisions(10);
+    }
+
+    pad->cd(2);
+    gPad->SetBorderMode(0);
+    fGraphError.Draw("ALP");
+    h = fGraphError.GetHistogram();
+    if (h)
+    {
+        h->SetXTitle("Train step/Tree number");
+        h->SetYTitle("Error/Resolution");
+    }
+
+    pad->cd(3);
+    gPad->SetBorderMode(0);
+    fGraphNodes.Draw("ALP");
+    h = fGraphNodes.GetHistogram();
+    if (h)
+    {
+        h->SetXTitle("Train step/Tree number");
+        h->SetYTitle("Number of end nodes");
+    }
+
+    pad->cd(4);
+    gPad->SetBorderMode(0);
+    fRules.Draw();
+}

trunk/MagicSoft/Mars/mranforest/MHRanForestGini.h

@@ -9 +9 @@
 #include <TArrayF.h>
 #endif
+#ifndef ROOT_TGraph
+#include <TGraph.h>
+#endif
+#ifndef ROOT_TPaveText
+#include <TPaveText.h>
+#endif
 
-class TH1D;
-class TGraph;
 class MParList;
 class MRanForest;
@@ -22 +26 @@
 
     TArrayF fGini;           //!
-    TGraph *fGraphGini;      //->
+
+    TGraph  fGraphGini;
+    TGraph  fGraphError;
+    TGraph  fGraphNodes;
+
+    TPaveText fRules;
 
 public:
     MHRanForestGini(Int_t nbins=100, const char *name=NULL, const char *title=NULL);
-    ~MHRanForestGini();
-
-    TGraph *GetGraphGini() const  { return fGraphGini; }
 
     Bool_t SetupFill(const MParList *plist);

trunk/MagicSoft/Mars/mranforest/MRanForest.cc

@@ -166 +166 @@
         if(grid[i]>=grid[i+1])
         {
-            *fLog<<inf<<"Grid points must be in increasing order! Ignoring grid."<<endl;
+            *fLog<<warn<<"Grid points must be in increasing order! Ignoring grid."<<endl;
             return;
         }
@@ -225 +225 @@
             fClass[j] = int(fHadTrue[j]+0.5);
     }
-
-    return;
 }
 
@@ -272 +270 @@
 Bool_t MRanForest::AddTree(MRanTree *rantree=NULL)
 {
-    fRanTree = rantree ? rantree:fRanTree;
+    fRanTree = rantree ? rantree : fRanTree;
 
     if (!fRanTree) return kFALSE;
@@ -287 +285 @@
     // access matrix, copy last column (target) preliminarily
     // into fHadTrue
-    TMatrix mat_tmp = mat->GetM();
-    int dim         = mat_tmp.GetNcols();
-    int numdata     = mat_tmp.GetNrows();
-
-    fMatrix=new TMatrix(mat_tmp);
+    fMatrix = new TMatrix(mat->GetM());
+
+    int dim     = fMatrix->GetNcols()-1;
+    int numdata = fMatrix->GetNrows();
 
     fHadTrue.Set(numdata);
@@ -297 +294 @@
 
     for (Int_t j=0;j<numdata;j++)
-        fHadTrue[j] = (*fMatrix)(j,dim-1);
+        fHadTrue[j] = (*fMatrix)(j,dim);
 
     // remove last col
-    fMatrix->ResizeTo(numdata,dim-1);
-    dim=fMatrix->GetNcols();
+    fMatrix->ResizeTo(numdata,dim);
 
     //-------------------------------------------------------------------
@@ -308 +304 @@
     fClass.Reset(0);
 
-    if(fClassify) PrepareClasses();
+    if (fClassify)
+        PrepareClasses();
 
     //-------------------------------------------------------------------
     // allocating and initializing arrays
-    fHadEst.Set(numdata);       fHadEst.Reset(0);
-    fNTimesOutBag.Set(numdata); fNTimesOutBag.Reset(0);
-    fDataSort.Set(dim*numdata); fDataSort.Reset(0);
-    fDataRang.Set(dim*numdata); fDataRang.Reset(0);
+    fHadEst.Set(numdata);
+    fHadEst.Reset(0);
+
+    fNTimesOutBag.Set(numdata);
+    fNTimesOutBag.Reset(0);
+
+    fDataSort.Set(dim*numdata);
+    fDataSort.Reset(0);
+
+    fDataRang.Set(dim*numdata);
+    fDataRang.Reset(0);
 
     if(fWeight.GetSize()!=numdata)
@@ -326 +330 @@
     //-------------------------------------------------------------------
     // setup rules to be used for classification/regression
-    MDataArray *allrules=(MDataArray*)mat->GetColumns();
+    const MDataArray *allrules=(MDataArray*)mat->GetColumns();
     if(allrules==NULL)
     {
@@ -333 +337 @@
     }
 
-    fRules=new MDataArray(); fRules->Reset();
-    TString target_rule;
-
-    for(Int_t i=0;i<dim+1;i++)
-    {
-        MData &data=(*allrules)[i];
-        if(i<dim)
-            fRules->AddEntry(data.GetRule());
-        else
-            target_rule=data.GetRule();
-    }
-
-    *fLog << inf <<endl;
-    *fLog << inf <<"Setting up RF for training on target:"<<endl<<" "<<target_rule.Data()<<endl;
-    *fLog << inf <<"Following rules are used as input to RF:"<<endl;
-
-    for(Int_t i=0;i<dim;i++)
-    {
-        MData &data=(*fRules)[i];
-        *fLog<<inf<<" "<<i<<") "<<data.GetRule()<<endl<<flush;
-    }
-
-    *fLog << inf <<endl;
+    fRules = new MDataArray();
+    fRules->Reset();
+
+    const TString target_rule = (*allrules)[dim];
+    for (Int_t i=0;i<dim;i++)
+        fRules->AddEntry((*allrules)[i].GetRule());
+
+    *fLog << inf << endl;
+    *fLog << "Setting up RF for training on target:" << endl;
+    *fLog << " " << target_rule.Data() << endl;
+    *fLog << "Following rules are used as input to RF:" << endl;
+    for (Int_t i=0;i<dim;i++)
+        *fLog << " " << i << ") " << (*fRules)[i].GetRule() << endl;
+
+    *fLog << endl;
 
     //-------------------------------------------------------------------
     // prepare (sort) data for fast optimization algorithm
-    if(!CreateDataSort()) return kFALSE;
+    if (!CreateDataSort())
+        return kFALSE;
 
     //-------------------------------------------------------------------
     // access and init tree container
     fRanTree = (MRanTree*)plist->FindCreateObj("MRanTree");
-
     if(!fRanTree)
     {
@@ -371 +367 @@
     }
 
+    const Int_t tryest = TMath::Nint(TMath::Sqrt(dim)+0.5);
+
+    *fLog << inf << endl;
+    *fLog << "Following input for the tree growing are used:"<<endl;
+    *fLog << " Number of Trees : "<<fNumTrees<<endl;
+    *fLog << " Number of Trials: "<<(fNumTry==0?tryest:fNumTry)<<(fNumTry==0?" (auto)":"")<<endl;
+    *fLog << " Final Node size : "<<fNdSize<<endl;
+    *fLog << " Using Grid:       "<<(fGrid.GetSize()>0?"Yes":"No")<<endl;
+    *fLog << " Number of Events: "<<numdata<<endl;
+    *fLog << " Number of Params: "<<dim<<endl;
+
+    if(fNumTry==0)
+    {
+        fNumTry=tryest;
+        *fLog << inf << endl;
+        *fLog << "Set no. of trials to the recommended value of round("<< TMath::Sqrt(dim) <<") = ";
+        *fLog << fNumTry << endl;
+
+
+    }
+    fRanTree->SetNumTry(fNumTry);
     fRanTree->SetClassify(fClassify);
     fRanTree->SetNdSize(fNdSize);
 
-    if(fNumTry==0)
-    {
-        double ddim = double(dim);
-
-        fNumTry=int(sqrt(ddim)+0.5);
-        *fLog<<inf<<endl;
-        *fLog<<inf<<"Set no. of trials to the recommended value of round("<<sqrt(ddim)<<") = ";
-        *fLog<<inf<<fNumTry<<endl;
-
-    }
-    fRanTree->SetNumTry(fNumTry);
-
-    *fLog<<inf<<endl;
-    *fLog<<inf<<"Following settings for the tree growing are used:"<<endl;
-    *fLog<<inf<<" Number of Trees : "<<fNumTrees<<endl;
-    *fLog<<inf<<" Number of Trials: "<<fNumTry<<endl;
-    *fLog<<inf<<" Final Node size : "<<fNdSize<<endl;
-
     fTreeNo=0;
 
@@ -415 +414 @@
     if (fTreeNo==1)
     {
-        *fLog << inf << endl;
-        *fLog << underline;
+        *fLog << inf << endl << underline;
 
         if(calcResolution)
@@ -435 +433 @@
     TArrayF winbag(numdata); // Initialization includes filling with 0
 
-    float square=0; float mean=0;
+    float square=0;
+    float mean=0;
 
     for (Int_t n=0; n<numdata; n++)
@@ -483 +482 @@
     for (Int_t ievt=0;ievt<numdata;ievt++)
     {
-        if (jinbag[ievt]>0) continue;
+        if (jinbag[ievt]>0)
+            continue;
 
         fHadEst[ievt] +=fRanTree->TreeHad((*fMatrix), ievt);
@@ -491 +491 @@
 
     Int_t n=0;
-    double ferr=0;
+    Float_t ferr=0;
 
     for (Int_t ievt=0;ievt<numdata;ievt++)
@@ -508 +508 @@
     //-------------------------------------------------------------------
     // give running output
-    *fLog << inf << setw(5)  << fTreeNo;
-    *fLog << inf << setw(18) << fRanTree->GetNumEndNodes();
-    *fLog << inf << Form("%18.2f", ferr*100.);
-    *fLog << inf << endl;
+    *fLog << setw(5)  << fTreeNo;
+    *fLog << setw(18) << fRanTree->GetNumEndNodes();
+    *fLog << Form("%18.2f", ferr*100.);
+    *fLog << endl;
+
+    fRanTree->SetError(ferr);
 
     // adding tree to forest

trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Hengstebeck 3/2003 <mailto:hengsteb@alwa02.physik.uni-siegen.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2003
+!   Author(s): Thomas Hengstebeck 2/2005 <mailto:hengsteb@physik.hu-berlin.de>
+!   Author(s): Thomas Bretz 8/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -27 +28 @@
 //  MRanForestCalc
 //
-//  Calculates the hadroness of an event. It calculates a mean value of all
-//  classifications by the trees in a previously grown random forest.
-//
-//  To use only n trees for your calculation use:
-//  MRanForestCalc::SetUseNumTrees(n);
 //
 ////////////////////////////////////////////////////////////////////////////
@@ -38 +34 @@
 #include <TVector.h>
 
-#include "MHMatrix.h" // must be before MLogManip.h
-#include "MDataArray.h"
+#include "MHMatrix.h"
 
 #include "MLog.h"
 #include "MLogManip.h"
 
+#include "MData.h"
+#include "MDataArray.h"
+
+#include "MRanForest.h"
+#include "MParameters.h"
+
 #include "MParList.h"
-
-#include "MRanTree.h"
-#include "MRanForest.h"
-
-#include "MParameters.h"
-
+#include "MTaskList.h"
 #include "MEvtLoop.h"
-#include "MTaskList.h"
+#include "MRanForestGrow.h"
 #include "MFillH.h"
-#include "MStatusDisplay.h"
-#include "MRanForestGrow.h"
-#include "MRanForestFill.h"
-
-#include "MWriteRootFile.h"
-#include "MReadTree.h"
 
 ClassImp(MRanForestCalc);
@@ -65 +55 @@
 using namespace std;
 
-static const TString gsDefName  = "MRanForestCalc";
-static const TString gsDefTitle = "Tree Classification Loop 1/2";
-
-// --------------------------------------------------------------------------
-//
-// Setup histograms and the number of distances which are used for
-// avaraging in CalcDist
-//
+const TString MRanForestCalc::gsDefName    = "MRanForestCalc";
+const TString MRanForestCalc::gsDefTitle   = "RF for energy estimation";
+
+const TString MRanForestCalc::gsNameOutput = "RanForestOut";
+
 MRanForestCalc::MRanForestCalc(const char *name, const char *title)
-    : fNum(100), fHadronnessName("MHadronness"), fData(NULL), fRanForest(0), fRanTree(0)
-{
-    //
-    //   set the name and title of this object
-    //
+    : fDebug(kFALSE), fData(0), fRFOut(0),
+    fNumTrees(-1), fNumTry(-1), fNdSize(-1), fNumObsoleteVariables(1),
+    fTestMatrix(0), fEstimationMode(kMean)
+{
     fName  = name  ? name  : gsDefName.Data();
     fTitle = title ? title : gsDefTitle.Data();
-}
-
-// --------------------------------------------------------------------------
-//
-// Delete the data chains
-//
+
+    // FIXME:
+    fNumTrees = 100; //100
+    fNumTry   = 0;   //3   0 means: in MRanForest estimated best value will be calculated
+    fNdSize   = 1;   //1   
+}
+
 MRanForestCalc::~MRanForestCalc()
 {
-    if (!IsOwner())
-        return;
-
-    delete fRanForest;
-    delete fRanTree;
-}
-
-// --------------------------------------------------------------------------
-//
-// Needs:
-//  - MatrixGammas  [MHMatrix]
-//  - MatrixHadrons [MHMatrix]
-//  - MHadronness
-//  - all data containers used to build the matrixes
-//
-// The matrix object can be filles using MFillH. And must be of the same
-// number of columns (with the same meaning).
-//
+    fEForests.Delete();
+}
+
+Int_t MRanForestCalc::Train(const MHMatrix &matrixtrain, const TArrayD &grid, Int_t ver)
+{
+    gLog.Separator("MRanForestCalc - Train");
+
+    if (!matrixtrain.GetColumns())
+    {
+        *fLog << err << "ERROR - MHMatrix does not contain rules... abort." << endl;
+        return kFALSE;
+    }
+
+    const Int_t ncols = matrixtrain.GetM().GetNcols();
+    const Int_t nrows = matrixtrain.GetM().GetNrows();
+    if (ncols<=0 || nrows <=0)
+    {
+        *fLog << err << "ERROR - No. of columns or no. of rows of matrixtrain equal 0 ... abort." << endl;
+        return kFALSE;
+    }
+
+    // rules (= combination of image par) to be used for energy estimation
+    TFile fileRF(fFileName, "recreate");
+    if (!fileRF.IsOpen())
+    {
+        *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl;
+        return kFALSE;
+    }
+
+    const Int_t nobs = fNumObsoleteVariables; // Number of obsolete columns
+
+    const MDataArray &dcol = *matrixtrain.GetColumns();
+
+    MDataArray usedrules;
+    for (Int_t i=0; i<ncols; i++)
+        if (i<ncols-nobs)  // -3 is important!!!
+            usedrules.AddEntry(dcol[i].GetRule());
+        else
+            *fLog << inf << "Skipping " << dcol[i].GetRule() << " for training" << endl;
+
+    MDataArray rules(usedrules);
+    rules.AddEntry(ver<2?"Classification":dcol[ncols-1].GetRule());
+
+    // prepare matrix for current energy bin
+    TMatrix mat(matrixtrain.GetM());
+
+    // last column must be removed (true energy col.)
+    mat.ResizeTo(nrows, ncols-nobs+1);
+
+    if (fDebug)
+        gLog.SetNullOutput(kTRUE);
+
+    const Int_t nbins = ver>0 ? 1 : grid.GetSize()-1;
+    for (Int_t ie=0; ie<nbins; ie++)
+    {
+        switch (ver)
+        {
+        case 0: // Replace Energy Grid by classification
+            {
+                Int_t irows=0;
+                for (Int_t j=0; j<nrows; j++)
+                {
+                    const Double_t energy = matrixtrain.GetM()(j,ncols-1);
+                    const Bool_t   inside = energy>grid[ie] && energy<=grid[ie+1];
+
+                    mat(j, ncols-nobs) = inside ? 1 : 0;
+
+                    if (inside)
+                        irows++;
+                }
+                if (irows==0)
+                    *fLog << warn << "WARNING - Skipping";
+                else
+                    *fLog << inf << "Training RF for";
+
+                *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irows << "/" << nrows << endl;
+
+                if (irows==0)
+                    continue;
+            }
+            break;
+
+        case 1: // Use Energy as classifier
+        case 2:
+            for (Int_t j=0; j<nrows; j++)
+                mat(j, ncols-nobs) = matrixtrain.GetM()(j,ncols-1);
+            break;
+        }
+
+        MHMatrix matrix(mat, &rules, "MatrixTrain");
+
+        MParList plist;
+        MTaskList tlist;
+        plist.AddToList(&tlist);
+        plist.AddToList(&matrix);
+
+        MRanForest rf;
+        rf.SetNumTrees(fNumTrees);
+        rf.SetNumTry(fNumTry);
+        rf.SetNdSize(fNdSize);
+        rf.SetClassify(ver<2 ? 1 : 0);
+        if (ver==1)
+            rf.SetGrid(grid);
+
+        plist.AddToList(&rf);
+
+        MRanForestGrow rfgrow;
+        tlist.AddToList(&rfgrow);
+
+        MFillH fillh("MHRanForestGini");
+        tlist.AddToList(&fillh);
+
+        MEvtLoop evtloop;
+        evtloop.SetParList(&plist);
+        evtloop.SetDisplay(fDisplay);
+        evtloop.SetLogStream(fLog);
+
+        if (!evtloop.Eventloop())
+            return kFALSE;
+
+        if (fDebug)
+            gLog.SetNullOutput(kFALSE);
+
+        if (ver==0)
+        {
+            // Calculate bin center
+            const Double_t E = (TMath::Log10(grid[ie])+TMath::Log10(grid[ie+1]))/2;
+
+            // save whole forest
+            rf.SetUserVal(E);
+            rf.SetName(Form("%.10f", E));
+        }
+
+        rf.Write();
+    }
+
+    // save rules
+    usedrules.Write("rules");
+
+    return kTRUE;
+}
+
+Int_t MRanForestCalc::ReadForests(MParList &plist)
+{
+    TFile fileRF(fFileName, "read");
+    if (!fileRF.IsOpen())
+    {
+        *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
+        return kFALSE;
+    }
+
+    fEForests.Delete();
+
+    TIter Next(fileRF.GetListOfKeys());
+    TObject *o=0;
+    while ((o=Next()))
+    {
+        MRanForest *forest=0;
+        fileRF.GetObject(o->GetName(), forest);
+        if (!forest)
+            continue;
+
+        forest->SetUserVal(atof(o->GetName()));
+
+        fEForests.Add(forest);
+    }
+
+    // Maybe fEForests[0].fRules yould be used instead?
+
+    if (fData->Read("rules")<=0)
+    {
+        *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
+        return kFALSE;
+    }
+
+    return kTRUE;
+}
+
 Int_t MRanForestCalc::PreProcess(MParList *plist)
 {
-    if (!fRanForest)
-    {
-        fRanForest = (MRanForest*)plist->FindObject("MRanForest");
-        if (!fRanForest)
+    fRFOut = (MParameterD*)plist->FindCreateObj("MParameterD", fNameOutput);
+    if (!fRFOut)
+        return kFALSE;
+
+    fData = (MDataArray*)plist->FindCreateObj("MDataArray");
+    if (!fData)
+        return kFALSE;
+
+    if (!ReadForests(*plist))
+    {
+        *fLog << err << "Reading RFs failed... aborting." << endl;
+        return kFALSE;
+    }
+
+    *fLog << inf << "RF read from " << fFileName << endl;
+
+    if (fTestMatrix)
+        return kTRUE;
+
+    fData->Print();
+
+    if (!fData->PreProcess(plist))
+    {
+        *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl;
+        return kFALSE;
+    }
+
+    return kTRUE;
+}
+
+#include <TGraph.h>
+#include <TF1.h>
+Int_t MRanForestCalc::Process()
+{
+    TVector event;
+    if (fTestMatrix)
+        *fTestMatrix >> event;
+    else
+        *fData >> event;
+
+    // --------------- Single Tree RF -------------------
+    if (fEForests.GetEntries()==1)
+    {
+        MRanForest *rf = (MRanForest*)fEForests[0];
+        fRFOut->SetVal(rf->CalcHadroness(event));
+        fRFOut->SetReadyToSave();
+
+        return kTRUE;
+    }
+
+    // --------------- Multi Tree RF -------------------
+    static TF1 f1("f1", "gaus");
+
+    Double_t sume = 0;
+    Double_t sumh = 0;
+    Double_t maxh = 0;
+    Double_t maxe = 0;
+
+    Double_t max  = -1e10;
+    Double_t min  =  1e10;
+
+    TIter Next(&fEForests);
+    MRanForest *rf = 0;
+
+    TGraph g;
+    while ((rf=(MRanForest*)Next()))
+    {
+        const Double_t h = rf->CalcHadroness(event);
+        const Double_t e = rf->GetUserVal();
+
+        g.SetPoint(g.GetN(), e, h);
+
+        sume += e*h;
+        sumh += h;
+
+        if (h>maxh)
         {
-            *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
-            return kFALSE;
+            maxh = h;
+            maxe = e;
         }
-    }
-
-    if (!fRanTree)
-    {
-        fRanTree = (MRanTree*)plist->FindObject("MRanTree");
-        if (!fRanTree)
-        {
-            *fLog << err << dbginf << "MRanTree not found... aborting." << endl;
-            return kFALSE;
-        }
-    }
-
-    fData = fRanTree->GetRules();
-
-    if (!fData)
-    {
-        *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl;
-        return kFALSE;
-    }
-
-    if (!fData->PreProcess(plist))
-    {
-        *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
-        return kFALSE;
-    }
-
-    fHadroness = (MParameterD*)plist->FindCreateObj("MParameterD", fHadronnessName);
-    if (!fHadroness)
-        return kFALSE;
+        if (e>max)
+            max = e;
+        if (e<min)
+            min = e;
+    }
+
+    switch (fEstimationMode)
+    {
+    case kMean:
+        fRFOut->SetVal(pow(10, sume/sumh));
+        break;
+    case kMaximum:
+        fRFOut->SetVal(pow(10, maxe));
+        break;
+    case kFit:
+        f1.SetParameter(0, maxh);
+        f1.SetParameter(1, maxe);
+        f1.SetParameter(2, 0.125);
+        g.Fit(&f1, "Q0N");
+        fRFOut->SetVal(pow(10, f1.GetParameter(1)));
+        break;
+    }
+
+    fRFOut->SetReadyToSave();
 
     return kTRUE;
@@ -153 +367 @@
 //
 //
-Int_t MRanForestCalc::Process()
-{
-    // first copy the data from the data array to a vector event
-    TVector event;
-    *fData >> event;
-
-    Double_t hadroness=fRanForest->CalcHadroness(event);
-    fHadroness->SetVal(hadroness);
-
-    return kTRUE;
-}
-
-Bool_t MRanForestCalc::Grow(MHMatrix *matrixg,MHMatrix *matrixh,Int_t ntree,
-                            Int_t numtry,Int_t ndsize,const char* treefile,
-                            const char* treename,const char* contname,
-                            const char* hgininame)
-{
-
-    treename  = treename  ? treename  : "Tree";
-    contname  = contname  ? contname  : "MRanTree";
-    hgininame  = hgininame  ? hgininame  : "MHRanForestGini";
-
-    if (!matrixg->IsValid())
-    {
-        *fLog << err << dbginf << " MRanForestCalc::Grow - ERROR: matrixg not valid." << endl;
-        return kFALSE;
-    }
-    if(!matrixh->IsValid())
-    {
-        *fLog << err << dbginf << " MRanForestCalc::Grow - ERROR: matrixh not valid." << endl;
-        return kFALSE;
-    }
-
-    MEvtLoop run(GetName());
-    MTaskList tlist;
-    MParList plist;
-    plist.AddToList(&tlist);
-    plist.AddToList(matrixg);
-    plist.AddToList(matrixh);
-
-    // creating training task and setting parameters
-    MRanForestGrow rfgrow;
-    rfgrow.SetNumTrees(ntree); // number of trees
-    rfgrow.SetNumTry(numtry);  // number of trials in random split selection
-    rfgrow.SetNdSize(ndsize);  // limit for nodesize
-    tlist.AddToList(&rfgrow);
-
-    if(treefile){
-        MWriteRootFile rfwrite(treefile);
-        rfwrite.AddContainer(contname,treename);
-        tlist.AddToList(&rfwrite);
-    }
-
-    MFillH fillh(hgininame);
-    tlist.AddToList(&fillh);
-
-    run.SetParList(&plist);
-    
-    // Execute tree growing
-    if (!run.Eventloop())
-    {
-        *fLog << err << dbginf << "Evtloop in MRanForestCalc::Grow failed." << endl;
-        return kFALSE;
-    }
-    tlist.PrintStatistics(0, kTRUE);
-
-    if (TestBit(kEnableGraphicalOutput))
-        plist.FindObject(hgininame)->DrawClone();
-
-    return kTRUE;
-}
-
-Bool_t MRanForestCalc::Fill(Int_t ntree,const char* treefile,const char* treename)
-{
-    treefile  = treefile  ? treefile  : "RF.root";
-    treename  = treename  ? treename  : "Tree";
-
-    MParList  plist;
-
-    MTaskList tlist;
-    plist.AddToList(&tlist);
-
-    MReadTree read(treename,treefile);
-    read.DisableAutoScheme();
-
-    MRanForestFill rffill;
-    rffill.SetNumTrees(ntree);
-
-    tlist.AddToList(&read);
-    tlist.AddToList(&rffill);
-
-    MEvtLoop run(GetName());
-    run.SetParList(&plist);
-
-    //
-    // Execute tree reading
-    //
-    if (!run.Eventloop())
-    {
-        *fLog << err << dbginf << "Evtloop in MRanForestCalc::Fill failed." << endl;
-        return kFALSE;
-    }
-    tlist.PrintStatistics(0, kTRUE);
-
-    return kTRUE;
-}
-
 Int_t MRanForestCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
 {
-    if (!IsEnvDefined(env, prefix, "File",     print))
-        return kFALSE;
-
-    TString fname = GetEnvValue(env, prefix, "File", "RF.root");
-    TString tname = GetEnvValue(env, prefix, "Tree", "Tree");
-    const Int_t   num   = GetEnvValue(env, prefix, "NumTrees", 100);
-    const Bool_t  debug = GetEnvValue(env, prefix, "Debug", kFALSE);
-
-    fname.ReplaceAll("\015", "");
-    tname.ReplaceAll("\015", "");
-
-    *fLog << inf << dbginf << "Reading " << num << " trees " << tname << " from file " << fname << endl;
-
-    gLog.SetNullOutput(!debug);
-    MEvtLoop evtloop;
-    MParList  plist;
-    evtloop.SetParList(&plist);
-    MLog l;
-    l.SetNullOutput(!debug);
-    evtloop.SetLogStream(&l);
-    gLog.SetNullOutput(debug);
-
-    if (IsOwner())
-    {
-        delete fRanForest;
-        delete fRanTree;
-    }
-    fRanForest = new MRanForest;
-    fRanTree   = new MRanTree;
-    SetOwner();
-
-    plist.AddToList(fRanForest);
-    plist.AddToList(fRanTree);
-
-    MTaskList tlist;
-    plist.AddToList(&tlist);
-
-    MReadTree read(tname, fname);
-    read.DisableAutoScheme();
-
-    MRanForestFill rffill;
-    rffill.SetNumTrees(num);
-
-    tlist.AddToList(&read);
-    tlist.AddToList(&rffill);
-
-    if (!evtloop.Eventloop())
-    {
-        *fLog << err << "ERROR - Reading " << tname << " from file " << fname << endl;
-        return kERROR;
-    }
-
-    return kTRUE;
-}
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "FileName", print))
+    {
+        rc = kTRUE;
+        SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
+    }
+    if (IsEnvDefined(env, prefix, "Debug", print))
+    {
+        rc = kTRUE;
+        SetDebug(GetEnvValue(env, prefix, "Debug", fDebug));
+    }
+    if (IsEnvDefined(env, prefix, "NameOutput", print))
+    {
+        rc = kTRUE;
+        SetNameOutput(GetEnvValue(env, prefix, "NameOutput", fNameOutput));
+    }
+    if (IsEnvDefined(env, prefix, "EstimationMode", print))
+    {
+        TString txt = GetEnvValue(env, prefix, "EstimationMode", "");
+        txt = txt.Strip(TString::kBoth);
+        txt.ToLower();
+        if (txt==(TString)"mean")
+            fEstimationMode = kMean;
+        if (txt==(TString)"maximum")
+            fEstimationMode = kMaximum;
+        if (txt==(TString)"fit")
+            fEstimationMode = kFit;
+        rc = kTRUE;
+    }
+    return rc;
+}

trunk/MagicSoft/Mars/mranforest/MRanForestCalc.h

@@ -6 +6 @@
 #endif
 
-#ifndef MARS_MHMatrix
-#include "MHMatrix.h"
+#ifndef ROOT_TObjArray
+#include <TObjArray.h>
 #endif
 
-class MParList;
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+class MDataArray;
 class MParameterD;
-class MDataArray;
-class MRanTree;
-class MRanForest;
+class MHMatrix;
 
 class MRanForestCalc : public MTask
 {
+public:
+    enum EstimationMode_t
+    {
+        kMean,
+        kMaximum,
+        kFit
+    };
+
 private:
-    Int_t  fNum;              // number of trees used to compute hadronness
+    static const TString gsDefName;     //! Default Name
+    static const TString gsDefTitle;    //! Default Title
+    static const TString gsNameOutput;  //! Default Output name
 
-    TString fHadronnessName;  // Name of container storing hadronness
+    Bool_t       fDebug;      // Debugging of eventloop while training on/off
 
-    MParameterD *fHadroness;  //! Output container for calculated hadroness
+    TString      fFileName;   // File name to forest
+    TObjArray    fEForests;   // List of forests
+
+    TString      fNameOutput; // Name of output container
+
     MDataArray  *fData;       //! Used to store the MDataChains to get the event values
-    MRanForest  *fRanForest;
-    MRanTree    *fRanTree;
+    MParameterD *fRFOut;      //! Used to store result
 
+    Int_t        fNumTrees;   //! Training parameters
+    Int_t        fNumTry;     //! Training parameters
+    Int_t        fNdSize;     //! Training parameters
+
+    Int_t        fNumObsoleteVariables;
+
+    MHMatrix    *fTestMatrix; //! Test Matrix used in Process (together with MMatrixLoop)
+
+    EstimationMode_t fEstimationMode;
+
+private:
+    // MTask
     Int_t PreProcess(MParList *plist);
     Int_t Process();
 
-    enum { kIsOwner = BIT(14) };
+    // MRanForestCalc
+    Int_t ReadForests(MParList &plist);
 
-    Bool_t IsOwner() const { return TestBit(kIsOwner); }
-    void SetOwner() { SetBit(kIsOwner); }
+    // MParContainer
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+
+    // Train Interface
+    Int_t Train(const MHMatrix &n, const TArrayD &grid, Int_t ver=2);
 
 public:
@@ -40 +71 @@
     ~MRanForestCalc();
 
-    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+    // Setter for estimation
+    void SetFileName(TString filename)            { fFileName = filename; }
+    void SetEstimationMode(EstimationMode_t op)   { fEstimationMode = op; }
+    void SetNameOutput(TString name=gsNameOutput) { fNameOutput = name; }
 
-    void SetHadronnessName(const TString name) { fHadronnessName = name; }
-    TString GetHadronnessName() const { return fHadronnessName; }
+    // Setter for training
+    void SetNumTrees(UShort_t n=100) { fNumTrees = n; }
+    void SetNdSize(UShort_t n=5)     { fNdSize   = n; }
+    void SetNumTry(UShort_t n=0)     { fNumTry   = n; }
+    void SetDebug(Bool_t b=kTRUE)    { fDebug    = b; }
 
-    void SetUseNumTrees(UShort_t n=100) { fNum = n; }
+    void SetNumObsoleteVariables(Int_t n=1) { fNumObsoleteVariables = n; }
 
-    Bool_t Grow(MHMatrix *matrixg,MHMatrix *matrixh,Int_t ntree,
-                Int_t numtry,Int_t ndsize,const char* treefile=0,
-                const char* treename=0,const char* contname=0,
-                const char* hgininame=0);
+    // Train Interface
+    Int_t TrainMultiRF(const MHMatrix &n, const TArrayD &grid)
+    {
+        return Train(n, grid, 0);
+    }
+    Int_t TrainSingleRF(const MHMatrix &n, const TArrayD &grid=TArrayD())
+    {
+        return Train(n, grid, grid.GetSize()==0 ? 2 : 1);
+    }
 
-    Bool_t Fill(Int_t ntree,const char* treefile=0,const char* treename=0);
+    // Test Interface
+    void  SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; }
+    void  InitMapping(MHMatrix *m=0)   { fTestMatrix=m; }
 
-
-    ClassDef(MRanForestCalc, 0) // Task
+    ClassDef(MRanForestCalc, 0) // Task to calculate RF output and for RF training
 };
 

trunk/MagicSoft/Mars/mranforest/MRanTree.h

@@ -28 +28 @@
     Int_t fNumNodes;
     Int_t fNumEndNodes;
+    Float_t fError;
 
     TArrayI fBestVar;
@@ -70 +71 @@
     void SetNdSize(Int_t n);
     void SetNumTry(Int_t n);
+    void SetError(Float_t f) { fError = f; }
 
     Int_t GetNdSize() const { return fNdSize; }
@@ -75 +77 @@
     Int_t GetNumNodes()          const { return fNumNodes; }
     Int_t GetNumEndNodes()       const { return fNumEndNodes; }
+    Float_t GetError() const { return fError; }
 
     Int_t GetBestVar(Int_t i)    const { return fBestVar.At(i); }

trunk/MagicSoft/Mars/mranforest/Makefile

@@ -25 +25 @@
            MRanForest.cc \
            MRanForestGrow.cc \
+           MRanForestCalc.cc \
            MHRanForest.cc \
-           MHRanForestGini.cc \
-           MRFEnergyEst.cc
+           MHRanForestGini.cc
 
 ############################################################

trunk/MagicSoft/Mars/mranforest/RanForestLinkDef.h

@@ -8 +8 @@
 #pragma link C++ class MRanForest+;
 #pragma link C++ class MRanForestGrow+;
+#pragma link C++ class MRanForestCalc+;
 
 #pragma link C++ class MHRanForest+;
 #pragma link C++ class MHRanForestGini+;
 
-#pragma link C++ class MRFEnergyEst+;
-
 #endif

trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc

@@ -204 +204 @@
 }
 
+//------------------------------------------------------------------------
+//
+// Add all entries deriving from MFilter from list to PreCuts.
+// The ownership is not affected.
+//
+void MTFillMatrix::AddPreCuts(const TList &list)
+{
+    TIter Next(&list);
+    TObject *obj=0;
+    while ((obj=Next()))
+        if (obj->InheritsFrom(MFilter::Class()))
+            fPreCuts.Add(obj);
+}
 
 // --------------------------------------------------------------------------
@@ -209 +222 @@
 // Fill the matrix (FIXME: Flow diagram missing)
 //
-Bool_t MTFillMatrix::Process()
+Bool_t MTFillMatrix::Process(const MParList &parlist)
 {
     if (!fReader)
@@ -232 +245 @@
     // Create parameter list and task list, add tasklist to parlist
     //
-    MParList  plist;
+    parlist.Print();
+    MParList  plist(parlist);
     MTaskList tlist;
     plist.AddToList(&tlist);
@@ -372 +386 @@
     return WriteMatrix1(fname) && WriteMatrix2(fname);
 }
+
+Int_t MTFillMatrix::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "NumDestEvents1", print))
+    {
+        rc = kTRUE;
+        SetNumDestEvents1(GetEnvValue(env, prefix, "NumDestEvents1", fNumDestEvents1));
+    }
+    if (IsEnvDefined(env, prefix, "NumDestEvents2", print))
+    {
+        rc = kTRUE;
+        SetNumDestEvents2(GetEnvValue(env, prefix, "NumDestEvents2", fNumDestEvents2));
+    }
+    return rc;
+}

trunk/MagicSoft/Mars/mtools/MTFillMatrix.h

@@ -1 +1 @@
 #ifndef MARS_MTFillMatrix
 #define MARS_MTFillMatrix
+
+#ifndef MARS_MParList
+#include "MParList.h"
+#endif
 
 #ifndef MARS_MH3
@@ -27 +31 @@
 
     TList     fPreCuts;
+
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
     Bool_t CheckResult(MHMatrix *m, Int_t num) const;
@@ -72 +78 @@
     void AddPreCut(const char *rule);
     void AddPreCut(MFilter *f);
+    void AddPreCuts(const TList &list);
 
-    Bool_t Process();
+    Bool_t Process(const MParList &plist=MParList());
     Bool_t WriteMatrix1(const TString &fname) const { return WriteMatrix(fDestMatrix1, fname, 1); }
     Bool_t WriteMatrix2(const TString &fname) const { return WriteMatrix(fDestMatrix2, fname, 2); }