Changeset 7142

Timestamp:

06/10/05 13:10:09 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (2 diffs)
• datacenter/macros/plotdb.C (modified) (1 diff)
• mdata/MDataChain.cc (modified) (1 diff)
• mhbase/MBinning.cc (modified) (5 diffs)
• mhbase/MBinning.h (modified) (1 diff)
• mhflux/MAlphaFitter.cc (modified) (3 diffs)
• mhflux/MHAlpha.cc (modified) (2 diffs)
• mhflux/MHCollectionArea.cc (modified) (2 diffs)
• mhflux/MHDisp.cc (modified) (2 diffs)
• mhflux/MHEffectiveOnTime.cc (modified) (1 diff)
• mhflux/MHEnergyEst.cc (modified) (1 diff)
• mhflux/MMcSpectrumWeight.cc (modified) (5 diffs)
• mhflux/MMcSpectrumWeight.h (modified) (4 diffs)
• mimage/MHVsSize.cc (modified) (2 diffs)
• mjobs/MDataSet.h (modified) (1 diff)
• mjobs/MJCut.cc (modified) (2 diffs)
• mjobs/MJOptimize.h (modified) (1 diff)
• mpointing/MSrcPosCam.cc (modified) (4 diffs)
• mpointing/MSrcPosCam.h (modified) (1 diff)
• mpointing/Makefile (modified) (1 diff)
• mranforest/MRFEnergyEst.cc (modified) (9 diffs)
• mranforest/MRFEnergyEst.h (modified) (3 diffs)
• mranforest/MRanTree.cc (modified) (4 diffs)
• mranforest/MRanTree.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -21 +21 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2005/06/10 Daniela Dorner
+
+   * datacenter/macros/plotdb.C:
+     - set minimum and maximum also for ZA graph
+
+   * mdata/MDataChain.cc:
+     - fixed a possible crash to to NULL-acess in GetDataMember
+
+   * mhbase/MBinning.[h,cc]:
+     - added a new binning type "asin" which is used for the
+       ZA binning
+
+   * mhflux/MAlphaFitter.cc:
+     - allow polynom order 1
+     - improved result line in PaintResult
+
+   * mhflux/MHAlpha.cc:
+     - call PaintResult also in DrawAll
+
+   * mhflux/MHCollectionArea.cc:
+     - removede nonsense A_{eff} from plot
+
+   * mhflux/MHDisp.cc:
+     - fixed a possible crash if fSrcPos==NULL
+
+   * mhflux/MHAlpha.cc, mhflux/MHCollectionArea.cc, 
+     mhflux/MHEffectiveOnTime.cc, mhflux/MHEnergyEst.cc,
+     mjobs/MJCut.cc:
+     - replaced SetEdgesCos by new SetEdgesASin (set the correct binning)
+       the old binning was not well aligned with the MC binning
+
+   * mhflux/MMcSpectrumWeight.[h,cc]:
+     - added the possibility to set ZA-weights
+       (Could be improved calculating correst sine-weights)
+
+   * mimage/MHVsSize.cc:
+     - fixed. Conc1 was incorrectly scaled
+
+   * mjobs/MDataSet.h:
+     - added getter for TLists
+
+   * mjobs/MJCut.cc:
+     - in non-wobble mode the hcalc2 had been wrrornously added 
+       to the tasklist -> removed
+       
+   * mjobs/MJOptimize.h:
+     - added a comment for EnableTestTrain
+
+   * mpointing/MSrcPosCam.[h,cc]:
+     - implemented Add member function
+
+   * mranforest/MRFEnergyEst.[h,cc]:
+     - fixed a problem with the removal of the last columns
+     - implemented debug-mode
+     - interpolate energy in logarithmic grid
+     - removed Test function completely
+     - added ReadEnv
+
+   * mranforest/MRanTree.[h,cc]:
+     - replaced some TMatrixFRow by TMatrixFRow_const
+
+
+
  2005/06/08 Daniela Dorner
 

trunk/MagicSoft/Mars/NEWS

@@ -2 +2 @@
 
  *** Version <cvs>
+
+   - general: Fixed the ZA binning. It didn't correctyl fit the
+     MC binning
+
+   - ganymed: the Conc1 plot was incorrectly scaled in MHVsSize
 
    - mars: show muon parameters graphically
@@ -7 +12 @@
    - mars: now the file to open can be given as commandline
      argument
- 
+
+
 
  *** Version 0.9.3 (2005/06/03)

trunk/MagicSoft/Mars/datacenter/macros/plotdb.C

@@ -97 +97 @@
         g.SetNameTitle(name, Form("%s vs Time", name.Data()));
         g.SetMarkerStyle(kFullDotMedium);
+
+        TGraph g2;
+        g2.SetNameTitle(name, Form("%s vs <Zd>", name.Data()));
+        g2.SetMarkerStyle(kFullDotMedium);
+
         if (fmax>fmin)
         {
             g.SetMinimum(fmin);
             g.SetMaximum(fmax);
+            g2.SetMinimum(fmin);
+            g2.SetMaximum(fmax);
         }
-
-        TGraph g2;
-        g2.SetNameTitle(name, Form("%s vs <Zd>", name.Data()));
-        g2.SetMarkerStyle(kFullDotMedium);
 
         Int_t first = -1;

trunk/MagicSoft/Mars/mdata/MDataChain.cc

@@ -832 +832 @@
 TString MDataChain::GetDataMember() const
 {
-    return fMember->GetDataMember();
+    return fMember ? fMember->GetDataMember() : "";
 }
 
 void MDataChain::SetVariables(const TArrayD &arr)
 {
-    return fMember->SetVariables(arr);
+    if (fMember)
+        fMember->SetVariables(arr);
 }
 

trunk/MagicSoft/Mars/mhbase/MBinning.cc

@@ -261 +261 @@
 //  lo: lowest edge
 //  hi: highest edge
-//  type: "lin" <default>, "log", "cos" (without quotationmarks)
+//  type: "lin" <default>, "log", "cos", "asin" (without quotationmarks)
 //  title: Whatever the title might be
 //
@@ -296 +296 @@
         str.Remove(0, pos);
         str = str.Strip(TString::kBoth);
-        if (typ!=(TString)"lin" && typ!=(TString)"log" && typ!=(TString)"cos")
+        if (typ!=(TString)"lin" && typ!=(TString)"log" && typ!=(TString)"cos" && typ!=(TString)"asin")
         {
             *fLog << warn << GetDescriptor() << "::SetEdges: Type " << typ << " unknown... ignored." << endl;
@@ -384 +384 @@
         return;
     }
+    if (o.Contains("asin", TString::kIgnoreCase))
+    {
+        SetEdgesASin(nbins, lo, up);
+        return;
+    }
     if (o.Contains("cos", TString::kIgnoreCase))
     {
@@ -420 +425 @@
     const Axis_t ud = up/kRad2Deg;
 
-    const Double_t binsize = (cos(ld)-cos(ud))/nbins;
+    const Double_t cld = ld<0 ? cos(ld)-1 : 1-cos(ld);
+    const Double_t cud = ud<0 ? cos(ud)-1 : 1-cos(ud);
+
+    SetEdgesASin(nbins, ld, ud);
+    /*
+    const Double_t binsize = (cld-cud)/nbins;
     fEdges.Set(nbins+1);
     for (int i=0; i<=nbins; i++)
-        fEdges[i] = acos(cos(ld)-binsize*i)*kRad2Deg;
+    {
+        const Double_t a = cld-binsize*i;
+        fEdges[i] = a<0 ? -acos(1+a)*kRad2Deg : acos(1-a)*kRad2Deg;
+        cout << a << " " << fEdges[i] << endl;
+    }
+
+    fType = kIsCosinic;*/
+}
+
+// --------------------------------------------------------------------------
+//
+// Specify the number of bins <nbins> (not the number of edges), the
+// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
+//
+void MBinning::SetEdgesASin(const Int_t nbins, Axis_t lo, Axis_t up)
+{
+    const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
+    fEdges.Set(nbins+1);
+    for (int i=0; i<=nbins; i++)
+    {
+        const Double_t a = binsize*i + lo;
+        fEdges[i] = a<0 ? -acos(1+a)*kRad2Deg : acos(1-a)*kRad2Deg;
+    }
 
     fType = kIsCosinic;
@@ -543 +575 @@
     {
         type = GetEnvValue(env, prefix, "Type", "lin");
-        if (type!=(TString)"lin" && type!=(TString)"log" && type!=(TString)"cos")
+        if (type!=(TString)"lin" && type!=(TString)"log" && type!=(TString)"cos" && type!=(TString)"acos")
         {
             *fLog << warn << GetDescriptor() << "::ReadEnv - WARNING: Type is not lin, log nor cos... assuming lin." << endl;

trunk/MagicSoft/Mars/mhbase/MBinning.h

@@ -66 +66 @@
     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);
 
     Int_t FindLoEdge(Double_t val) const

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

@@ -117 +117 @@
 
     fFunc->FixParameter(1, 0);
-    fFunc->FixParameter(4, 0);
+    if (fPolynomOrder!=1)
+        fFunc->FixParameter(4, 0);
     fFunc->SetParLimits(2, 0, 90);
     fFunc->SetParLimits(3, -1, 1);
@@ -159 +160 @@
     fFunc->ReleaseParameter(2);
     fFunc->FixParameter(3, fFunc->GetParameter(3));
+    fFunc->FixParameter(4, fFunc->GetParameter(4));
     for (int i=5; i<fFunc->GetNpar(); i++)
         fFunc->FixParameter(i, fFunc->GetParameter(i));
@@ -297 +299 @@
     const Int_t    l1 = w<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(w));
     const Int_t    l2 = m<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(m));
-    const TString fmt = Form("\\sigma_{Li/Ma}=%%.1f  \\omega=%%.%df\\circ  E=%%d  B=%%d  (x<%%.%df)  (\\chi_{b}^{2}/ndf=%%.1f  \\chi_{s}^{2}/ndf=%%.1f  c_{0}=%%.1f)",
+    const TString fmt = Form("\\sigma_{L/M}=%%.1f  \\omega=%%.%df\\circ  E=%%d B=%%d  x<%%.%df  \\tilde\\chi_{b}=%%.1f  \\tilde\\chi_{s}=%%.1f  c=%%.1f  f=%%.2f",
                              l1<1?1:l1+1, l2<1?1:l2+1);
 
     TLatex text(x, y, Form(fmt.Data(), fSignificance, w, (int)fEventsExcess,
                            (int)fEventsBackground, m, fChiSqBg, fChiSqSignal,
-                           fCoefficients[3]));
+                           fCoefficients[3], fScaleFactor));
 
     text.SetBit(TLatex::kTextNDC);

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -139 +139 @@
     binsa.SetEdges(18, 0, 90);
     binse.SetEdgesLog(15, 10, 100000);
-    binst.SetEdgesCos(50, 0, 60);
+    binst.SetEdgesASin(51, -0.005, 0.505);
     binse.Apply(fHEnergy);
     binst.Apply(fHTheta);
@@ -818 +818 @@
 
         if (hof ? fit.Fit(*hon, *hof, alpha) : fit.Fit(*hon))
-            *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
+        {
+            *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << " scale=" << fit.GetScaleFactor() << endl;
+            fit.PaintResult();
+        }
         /*
         if (fit.FitEnergy(fHist, fOffData, i, kTRUE))

trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

@@ -95 +95 @@
     MBinning binsa, binse, binst;
     binse.SetEdgesLog(15, 10, 1000000);
-    binst.SetEdgesCos(50, 0, 60);
+    binst.SetEdgesASin(51, -0.005, 0.505);
 
     binse.Apply(fHEnergy);
@@ -289 +289 @@
     if (TString(option)=="paint4")
     {
-        const TString txt = Form("A_{eff}=%.0fm^{2}  A_{abs}=%.0fm^{2}  r=%.0fm",
-                                 GetCollectionAreaEff(),
+        //const TString txt = Form("A_{eff}=%.0fm^{2}  A_{abs}=%.0fm^{2}  r=%.0fm",
+        //                         GetCollectionAreaEff(),
+        //                         GetCollectionAreaAbs(), fMcAreaRadius);
+        const TString txt = Form("A_{abs}=%.0fm^{2}  r=%.0fm",
                                  GetCollectionAreaAbs(), fMcAreaRadius);
 

trunk/MagicSoft/Mars/mhflux/MHDisp.cc

@@ -203 +203 @@
 
     // Now we can safly derotate both position...
-    TVector2 srcp(fSrcPos->GetXY());
+    TVector2 srcp;
+    if (fSrcPos)
+        srcp = fSrcPos->GetXY();
+
     if (rho!=0)
     {
@@ -211 +214 @@
     }
 
-    if (fSrcPos)
-    {
-        pos1 -= srcp*fMm2Deg;
-        pos2 -= srcp*fMm2Deg;
-    }
+    pos1 -= srcp*fMm2Deg;
+    pos2 -= srcp*fMm2Deg;
 
     fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);

trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc

@@ -410 +410 @@
     // setup binning
     MBinning btheta("BinningTheta");
-    btheta.SetEdgesCos(100, 0, 60);
+    btheta.SetEdgesASin(51, -0.005, 0.505);
 
     MBinning btime("BinningDeltaT");

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -96 +96 @@
     MBinning binsi, binse, binst, binsr;
     binse.SetEdgesLog(15, 10, 1000000);
-    binst.SetEdgesCos(50, 0, 60);
+    binst.SetEdgesASin(51, -0.005, 0.505);
     binsi.SetEdges(10, 0, 400);
     binsr.SetEdges(50, -1.3, 1.3);

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -72 +72 @@
 
 #include <TF1.h>
+#include <TH1.h>
 
 #include "MLog.h"
@@ -79 +80 @@
 #include "MParameters.h"
 
+#include "MPointingPos.h"
+
 #include "MMcEvt.hxx"
 #include "MMcCorsikaRunHeader.h"
@@ -106 +109 @@
     fAllowChange = kFALSE;
 
-    fFunc   = NULL;
-    fMcEvt  = NULL;
-    fWeight = NULL;
+    fFunc      = NULL;
+    fMcEvt     = NULL;
+    fWeight    = NULL;
+    fZdWeights = NULL;
+    fPointing  = NULL;
 }
 
@@ -150 +155 @@
     if (!fWeight)
         return kFALSE;
+
+    if (!fZdWeights)
+        return kTRUE;
+
+    fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
+    if (!fPointing)
+    {
+        *fLog << err << "MPointingPos not found... abort." << endl;
+        return kFALSE;
+    }
 
     return kTRUE;
@@ -362 +377 @@
     const Double_t e = fMcEvt->GetEnergy();
 
-    fWeight->SetVal(fFunc->Eval(e));
+    Double_t w = 1;
+
+    if (fZdWeights)
+    {
+        const Int_t i = fZdWeights->GetXaxis()->FindFixBin(fPointing->GetZd());
+        w = fZdWeights->GetBinContent(i);
+        cout << i << " " << w << " " << fPointing->GetZd() << endl;
+    }
+
+    fWeight->SetVal(fFunc->Eval(e)*w);
 
     return kTRUE;

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

@@ -7 +7 @@
 
 class TF1;
+class TH1;
 class MParList;
 class MMcEvt;
 class MParameterD;
+class MPointingPos;
 class MMcCorsikaRunHeader;
 
@@ -17 +19 @@
     const MMcEvt *fMcEvt;   // Pointer to the container with the MC energy
     MParameterD  *fWeight;  // Pointer to the output MWeight container
+    MPointingPos *fPointing;
 
     TString fNameWeight;    // Name of the MWeight container
@@ -22 +25 @@
 
     TF1 *fFunc;             // Function calculating the weights
+    TH1 *fZdWeights;        // Set additional ZA weights
 
     Double_t fOldSlope;     // Slope of energy spectrum generated with Corsika
@@ -59 +63 @@
     void SetEnergyRange(Double_t min=-2, Double_t max=-1) { fEnergyMin=min; fEnergyMax=max; }
     void SetOldSlope(Double_t s=-2.6) { fOldSlope=s; }
+    void SetZdWeights(TH1 *h=0) { fZdWeights = h; }
     Bool_t Set(const MMcCorsikaRunHeader &h);
 

trunk/MagicSoft/Mars/mimage/MHVsSize.cc

@@ -99 +99 @@
     binsl.SetEdges(   100,    0, 296.7/2);
     binsd.SetEdges(   100,    0, 600);
-    binsc.SetEdgesLog(100, 1e-5, 5e-3);
+    binsc.SetEdgesLog(100, 3e-3, 1);
     binsa.SetEdges(   100,    0, 445*45);
     binsm.SetEdges(   100, -445, 445);
@@ -269 +269 @@
     fWidth.Fill( fHillas->GetSize(), scale*fHillas->GetWidth(),      w);
     fDist.Fill(  fHillas->GetSize(), scale*src->GetDist(),           w);
-    fConc1.Fill( fHillas->GetSize(), scale*fNewImagePar->GetConc1(), w);
+    fConc1.Fill( fHillas->GetSize(), fNewImagePar->GetConc1(),       w);
     fArea.Fill(  fHillas->GetSize(), scale*scale*fHillas->GetArea(), w);
     fM3Long.Fill(fHillas->GetSize(), scale*fHillasExt->GetM3Long()*TMath::Sign(1.0f, src->GetCosDeltaAlpha()), w);

trunk/MagicSoft/Mars/mjobs/MDataSet.h

@@ -60 +60 @@
     static Int_t  AddFilesToChain(MDirIter &files, TChain &chain);
 
+    const TList &GetSequencesOn() const  { return fSequencesOn; }
+    const TList &GetSequencesOff() const { return fSequencesOff; }
+
     Bool_t AddFiles(MRead &read) const;
     Bool_t AddFilesOn(MRead &read) const;

trunk/MagicSoft/Mars/mjobs/MJCut.cc

@@ -448 +448 @@
 
     // Initialize default binnings
-    MBinning bins1(18, 0,  90,   "BinningAlpha",     "lin");
-    MBinning bins2(15, 10, 1e6 , "BinningEnergyEst", "log");
-    MBinning bins3(50, 0,  60,   "BinningTheta",     "cos");
+    MBinning bins1(18,  0,     90,    "BinningAlpha",     "lin");
+    MBinning bins2(15, 10,     1e6 ,  "BinningEnergyEst", "log");
+    MBinning bins3(51, -0.005, 0.505, "BinningTheta",     "asin");
     MBinning bins4("BinningFalseSource");
     MBinning bins5("BinningWidth");
@@ -559 +559 @@
     tlist2.AddToList(&scalc);
     tlist2.AddToList(&hcalc);
-    //if (fIsWobble)
+    if (fIsWobble)
         tlist2.AddToList(&hcalc2);
     //tlist2.AddToList(&taskenv1);

trunk/MagicSoft/Mars/mjobs/MJOptimize.h

@@ -111 +111 @@
     void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; }
     void SetTolerance(Float_t tol=0)  { fTolerance=tol; }
-    void EnableTestTrain(Int_t b=2)   { fTestTrain=b; }
+    void EnableTestTrain(Int_t b=1)   { fTestTrain=b; } // Use 1 and -1
 
     // Parameter access

trunk/MagicSoft/Mars/mpointing/MSrcPosCam.cc

@@ -62 +62 @@
 // Copy constructor, set default name and title
 //
-MSrcPosCam::MSrcPosCam(const MSrcPosCam &p) : fX(p.fX), fY(p,fY)
+MSrcPosCam::MSrcPosCam(const MSrcPosCam &p) : fX(p.fX), fY(p.fY)
 {
     fName  = gsDefName.Data();
@@ -69 +69 @@
 
 // -----------------------------------------------------------------------
+//
+// Print contents
 //
 void MSrcPosCam::Print(Option_t *) const
@@ -78 +80 @@
 }
 
+// -----------------------------------------------------------------------
+//
+// Set fX to v.X() and fY to v.Y()
+//
 void MSrcPosCam::SetXY(const TVector2 &v)
 {
@@ -84 +90 @@
 }
 
+// -----------------------------------------------------------------------
+//
+// Add v.X() to fX and v.Y() to fY
+//
+void MSrcPosCam::Add(const TVector2 &v)
+{
+    fX+=v.X();
+    fY+=v.Y();
+}
+
+// -----------------------------------------------------------------------
+//
+// Get TVector2(fX, fY)
+//
 TVector2 MSrcPosCam::GetXY() const
 {

trunk/MagicSoft/Mars/mpointing/MSrcPosCam.h

@@ -25 +25 @@
     void SetXY(const TVector2 &v);
 
+    void Add(const TVector2 &v);
+
     Float_t GetX() const             { return fX; }
     Float_t GetY() const             { return fY; }
     TVector2 GetXY() const;
+
 
     void Print(Option_t *opt=NULL) const;

trunk/MagicSoft/Mars/mpointing/Makefile

@@ -28 +28 @@
            MSrcPosCam.cc \
            MSrcPosCalc.cc \
+           MSrcPosCorrect.cc \
            MSrcPosFromModel.cc
 

trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc

@@ -69 +69 @@
 //
 MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
-    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fData(0), fEnergyEst(0),
-    fTestMatrix(0)
+    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fDebug(kFALSE),
+    fData(0), fEnergyEst(0), fTestMatrix(0)
 {
     fName  = name  ? name  : gsDefName.Data();
@@ -91 +91 @@
     }
 
-    const Int_t ncols = matrixtrain.GetM().GetNcols();
-    const Int_t nrows = matrixtrain.GetM().GetNrows();
+    const TMatrix &m = matrixtrain.GetM();
+
+    const Int_t ncols = m.GetNcols();
+    const Int_t nrows = m.GetNrows();
     if (ncols<=0 || nrows <=0)
     {
@@ -114 +116 @@
     }
 
+    const Int_t nobs = 3; // Number of obsolete columns
+
     MDataArray usedrules;
-    for(Int_t i=0; i<ncols-3; i++) // -3 is important!!!
+    for(Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
         usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());
-
-    const TMatrix &m = matrixtrain.GetM();
 
     // training of RF for each energy bin
@@ -140 +142 @@
         }
 
-        const Bool_t valid = irow1*irow0>0;
-
-        if (!valid)
+        const Bool_t invalid = irow1==0 || irow0==0;
+
+        if (invalid)
             *fLog << warn << "WARNING - Skipping";
         else
             *fLog << inf << "Training RF for";
 
-        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ")" << endl;
-
-        if (!valid)
+        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irow0 << " " << irow1 << endl;
+
+        if (invalid)
             continue;
 
-        gLog.SetNullOutput(kTRUE);
+        if (fDebug)
+            gLog.SetNullOutput(kTRUE);
 
         // last column must be removed (true energy col.)
-        mat1.ResizeTo(irow1, ncols-1);
-        mat0.ResizeTo(irow0, ncols-1);
+        mat1.ResizeTo(irow1, ncols-nobs);
+        mat0.ResizeTo(irow0, ncols-nobs);
 
         // create MHMatrix as input for RF
@@ -183 +186 @@
         evtloop.SetParList(&plist);
 
-        gLog.SetNullOutput(kFALSE);
+        if (fDebug)
+            gLog.SetNullOutput(kFALSE);
 
         if (!evtloop.Eventloop())
             return kFALSE;
 
+        const Double_t E = (log10(grid[ie])+log10(grid[ie+1]))/2;
+
         // save whole forest
         MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
-        const TString title = Form("%f", TMath::Sqrt(grid[ie]*grid[ie+1]));
-        //const TString title = Form("%f", (grid[ie]+grid[ie+1])/2);
+        const TString title = Form("%.10f", E);
         forest->SetTitle(title);
         forest->Write(title);
@@ -203 +208 @@
     return kTRUE;
 }
-/*
-Int_t MRFEnergyEst::Test(const MHMatrix &matrixtest)
-{
-    gLog.Separator("MRFEnergyEst - Test");
-
-    const Int_t ncols = matrixtest.GetM().GetNcols();
-    const Int_t nrows = matrixtest.GetM().GetNrows();
-
-    if (ncols<=0 || nrows <=0)
-    {
-        *fLog << err << dbginf << "No. of columns or no. of rows of matrixtrain equal 0 ... aborting." << endl;
-        return kFALSE;
-    }
-
-    if (!ReadForests())
-    {
-        *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
-        return kFALSE;
-    }
-
-    const Int_t nbins=fEForests.GetSize();
-
-    Double_t elow =  FLT_MAX;
-    Double_t eup  = -FLT_MAX;
-
-    for(Int_t j=0;j<nbins;j++)
-    {
-        elow = TMath::Min(atof(fEForests[j]->GetTitle()), elow);
-        eup  = TMath::Max(atof(fEForests[j]->GetTitle()), eup);
-    }
-
-    TH1F hres("hres", "", 1000, -10, 10);
-    TH2F henergy("henergy", "",100, elow, eup,100, elow, eup);
-
-    const TMatrix &m=matrixtest.GetM();
-    for(Int_t i=0;i<nrows;i++)
-    {
-        Double_t etrue = m(i,ncols-1);
-        Double_t eest  = 0;
-        Double_t hsum  = 0;
-
-        for(Int_t j=0;j<nbins;j++)
-        {
-            const TVector  v = TMatrixFRow_const(m,i);
-
-            const Double_t h = ((MRanForest*)fEForests[j])->CalcHadroness(v);
-            const Double_t e = atof(fEForests[j]->GetTitle());
-
-            hsum += h;
-            eest += h*e;
-        }
-        eest /= hsum;
-        eest  = pow(10., eest);
-
-        //if (etrue>80.)
-        //    hres.Fill((eest-etrue)/etrue);
-
-        henergy.Fill(log10(etrue),log10(eest));
-    }
-
-    if(gStyle)
-        gStyle->SetOptFit(1);
-
-    // show results
-    TCanvas *c=MH::MakeDefCanvas("c","",300,900);
-
-    c->Divide(1,3);
-    c->cd(1);
-    henergy.SetTitle("Estimated vs true energy");
-    henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])");
-    henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])");
-    henergy.DrawCopy();
-
-    c->cd(2);
-    TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S");
-    hptr->SetTitle("Estimated vs true energy - profile");
-    hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])");
-    hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])");
-    hptr->DrawCopy();
-
-    c->cd(3);
-    hres.Fit("gaus");
-    hres.SetTitle("Energy resolution for E_{true}>80Gev");
-    hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})");
-    hres.GetYaxis()->SetTitle("counts");
-    hres.DrawCopy();
-
-    c->GetPad(1)->SetGrid();
-    c->GetPad(2)->SetGrid();
-    c->GetPad(3)->SetGrid();
-
-    return kTRUE;
-}
-*/
-Int_t MRFEnergyEst::ReadForests(MParList *plist)
+
+Int_t MRFEnergyEst::ReadForests(MParList &plist)
 {
     TFile fileRF(fFileName,"read");
@@ -321 +233 @@
 
         fEForests.Add(forest);
-
-    }
-
-    if (plist)
-    {
-        fData = (MDataArray*)plist->FindCreateObj("MDataArray");
-        fData->Read("rules");
-    }
-
-    fileRF.Close();
+    }
+
+    if (fData->Read("rules")<=0)
+    {
+        *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
+        return kFALSE;
+    }
 
     return kTRUE;
@@ -341 +250 @@
         return kFALSE;
 
-    if (!ReadForests(plist))
+    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;
 
-    if (!fData)
-    {
-        *fLog << err << "MDataArray not found... aborting." << endl;
-        return kFALSE;
-    }
+    fData->Print();
 
     if (!fData->PreProcess(plist))
@@ -387 +298 @@
 
         hsum += h;
-        eest += h*e;
-    }
-
-    fEnergyEst->SetVal(eest/hsum);
+        eest += e*h;
+    }
+
+    fEnergyEst->SetVal(pow(10, eest/hsum));
     fEnergyEst->SetReadyToSave();
 
     return kTRUE;
 }
+
+// --------------------------------------------------------------------------
+//
+//
+Int_t MRFEnergyEst::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    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));
+    }
+    return rc;
+}

trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.h

@@ -23 +23 @@
     Int_t        fNdSize;    // Training parameters
 
+    Bool_t       fDebug;     // Debugging of eventloop while training on/off
+
     TString      fFileName;
     TObjArray    fEForests;
@@ -34 +36 @@
     Int_t Process();
 
-    Int_t ReadForests(MParList *plist=NULL);
+    Int_t ReadForests(MParList &plist);
+
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
 public:
@@ -41 +45 @@
     void  SetFileName(TString str)     { fFileName = str; }
 
-    void  SetNumTrees(Int_t n=-1)   { fNumTrees = n; }
-    void  SetNdSize(Int_t n=-1)     { fNdSize   = n; }
-    void  SetNumTry(Int_t n=-1)     { fNumTry   = n; }
+    void  SetNumTrees(Int_t n=-1)      { fNumTrees = n; }
+    void  SetNdSize(Int_t n=-1)        { fNdSize   = n; }
+    void  SetNumTry(Int_t n=-1)        { fNumTry   = n; }
 
     void  SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; }
     void  InitMapping(MHMatrix *m=0)   { fTestMatrix=m; }
+
+    void SetDebug(Bool_t b=kTRUE)      { fDebug = b; }
 
     Int_t Train(const MHMatrix &n, const TArrayD &grid);

trunk/MagicSoft/Mars/mranforest/MRanTree.cc

@@ -447 +447 @@
     // are coded into fBestVar:
     // status of node kt = TMath::Sign(1,fBestVar[kt])
-    // class  of node kt = fBestVar[kt]+2 (class defined by larger
-    //  node population, actually not used)
+    // class  of node kt = fBestVar[kt]+2
+    //  (class defined by larger node population, actually not used)
     // hadronness assigned to node kt = fBestSplit[kt]
 
@@ -463 +463 @@
 }
 
-Double_t MRanTree::TreeHad(const TMatrixRow &event)
+Double_t MRanTree::TreeHad(const TMatrixFRow_const &event)
 {
     Int_t kt=0;
@@ -504 +504 @@
 Bool_t MRanTree::AsciiWrite(ostream &out) const
 {
-    TString str;
+    out.width(5);
+    out << fNumNodes << endl;
+
     Int_t k;
-
-    out.width(5);out<<fNumNodes<<endl;
-
-    for (k=0;k<fNumNodes;k++)
-    {
-        str=Form("%f",GetBestSplit(k));
+    for (k=0; k<fNumNodes; k++)
+    {
+        TString str=Form("%f", GetBestSplit(k));
 
         out.width(5);  out << k;
@@ -521 +520 @@
         out.width(5);  out << GetNodeClass(k);
     }
-    out<<endl;
-
-    return k==fNumNodes;
-}
+    out << endl;
+
+    return kTRUE;
+}

trunk/MagicSoft/Mars/mranforest/MRanTree.h

@@ -16 +16 @@
 class TMatrix;
 class TMatrixRow;
+class TMatrixFRow_const;
 class TVector;
 class TRandom;
@@ -84 +85 @@
 
     Double_t TreeHad(const TVector &event);
-    Double_t TreeHad(const TMatrixRow &event);
+    Double_t TreeHad(const TMatrixFRow_const &event);
     Double_t TreeHad(const TMatrix &m, Int_t ievt);
     Double_t TreeHad();