Changeset 7413 for trunk


Ignore:
Timestamp:
11/21/05 11:09:12 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft
Files:
5 added
31 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Cosy/tpoint/tpoint0509.txt

    r7391 r7413  
    11Magic Model  TPOINT data file
    22: ALTAZ
    3 #
    4 # ***** After a refocussing of the mirror *****
    5 #
    63
    74### 24.9.: tpoints always with the white paper screen
     
    1181.742509 57.99808 144.2056 49.02492 0.945 60.71667 0.0005732658 0.007748751 53637.075763 216.1 2.654
    129# Hamal,  mag=2.00, Zd=5
    13 139.6868 83.23968 281.6764 74.22004 2.119444 23.4625 -0.04332626 0.03920611 53637.11646 225.2 2.635
     10##139.6868 83.23968 281.6764 74.22004 2.119444 23.4625 -0.04332626 0.03920611 53637.11646 225.2 2.635
    1411
    1512# --- 24.09.2005 --- 02:59:09.399
    1613# Hamal,  mag=2.00, Zd=5
    17 161.5844 84.46453 303.3778 75.43595 2.119444 23.4625 -0.009743542 0.1953913 53637.124414 225.2 2.665
     14##161.5844 84.46453 303.3778 75.43595 2.119444 23.4625 -0.009743542 0.1953913 53637.124414 225.2 2.665
    1815# And-51, mag=3.57, Zd=21
    1916-21.2567 68.112 121.1535 59.26484 1.633333 48.62833 0.02363187 0.04614806 53637.142328 218.1 3.353
     
    2421# --- 25.09.2005 --- 03:41:32.305
    2522# Nothing in the runbook?
    26 -35.99076 37.46571 106.0467 28.7287 22.48611 58.41528 0.1168588 0.4782401 53638.153846 0 6.58
     23##-35.99076 37.46571 106.0467 28.7287 22.48611 58.41528 0.1168588 0.4782401 53638.153846 0 6.58
    2724-38.29146 35.14824 104.1974 26.30449 22.25056 57.04361 0.01046067 0.02796104 53638.15701 220.4 3.97
    2825## seems to be wrongly identified
     
    7673112.9061  44.22773 255.4126 35.24166 7.655     5.225    -0.005348532  0.000207558 53639.239228 215.7 2.556
    7774
     75####################################################################
     76# Private communication with Adrian:
     77# Als wir Ende September rausfanden, dass die Laser unzuverlaessig
     78# sind, wurde da eine Quick-and-Dirty Kalibration gemacht (da kein
     79# Experte auf La Palma war)
     80####################################################################
     81
  • trunk/MagicSoft/Cosy/tpoint/tpoint0510.txt

    r7391 r7413  
    22: ALTAZ
    33
     4#
     5# ***** New focussing after collabortaion meeting in Tenerife (19.10.?) *****
     6#
     7
    48# --- 25.10.2005 --- 02:11:32.993
     9# Diphda. Magn: 2.04 Zd:59.21 Az:222.64
    510-137.4972 30.88874 5.015726 21.99765 0.7263889 -17.98667 -0.03588449 -0.01783158 53668.091354 237.2 2.802
     11# Algenib. Magn:2.83 Zd:45.55 Az:263.21
    612-96.64771 44.18443 45.84061 35.33164 0.2205556 15.18361 -0.0004117089 -0.002440412 53668.098033 235.3 3.151
     13# Markab. Magn:2.49 Zd:65.51 Az:273.02
    714-87.16777 27.87191 55.32808 19.00867 23.07944 15.20528 -0.008018477 0.009832391 53668.10233 237.3 3.096
     15# Markab. Magn:2.49 Zd:63.7 Az:273.65
    816-86.52195 26.6309 55.98938 17.74595 23.07944 15.20528 -0.02510043 -0.001211299 53668.106316 233.3 3.102
     17
  • trunk/MagicSoft/Mars/Changelog

    r7409 r7413  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20 2005/11/21 Thomas Bretz
     21
     22   * mjtrain/*, macros/train/*.C:
     23     - added new class to train the random forest
     24
     25   * Makefile:
     26     - added mjtrain
     27
     28   * macros/dohtml.C, macros/rootlogon.C:
     29     - added mjoptim and mjtrain
     30
     31   * mbase/MParList.[h,cc]:
     32     - fixed copy constructor (it was crashing due to the
     33       non-initialized lists)
     34
     35   * mhbase/MBinning.[h,cc]:
     36     - renamed SetEdges to SetEdgesLin (SetEdges now is just and
     37       abbreviation)
     38     - added a check for the number of bins to all SetEdges*
     39
     40   * mhbase/MH3.cc:
     41     - enable grid
     42
     43   * mhbase/MHMatrix.[h,cc]:
     44     - new member function Addcolumns taking a TCollection
     45
     46   * mhist/MHHadronness.cc:
     47     - align hadronnes into [0,1]
     48
     49   * mjobs/MDataSet.[h,cc]:
     50     - call SetOwner for fSequencesOn and fSequencesOff
     51     - added Copy member function
     52
     53   * mjobs/MJCut.[h,cc]:
     54     - added new TaskEnv "CalcDisp" for a disp estimator
     55
     56   * mjoptim/MJOptimizeDisp.cc:
     57     - fixed typo in a title
     58     - added a new histogram showing theta-sq versus size
     59
     60   * mranforest/MHRanForestGini.[h,cc]:
     61     - revised to display more information
     62     - exchanged pointers to TGraph by objects
     63
     64   * mranforest/MRanForest.[h,cc]:
     65     - replaced some weired copy of the train matrix by a direct access
     66     - revised output
     67     - pipe error/resolution to tree
     68
     69   * mranforest/MRanForestCalc.[h,cc]:
     70     - copied code from MRFEnergyEst
     71     - allow to set a name for the output container
     72     - set numtry to 0 (auto)
     73     - set ndsize to 1 (there is no auto mode)
     74     - introduced variable for number of obsolete parameters
     75
     76   * mranforest/MRanTree.h:
     77     - new data member to store resolution/error
     78
     79   * mranforest/Makefile, mranforest/RanForestLinkDef.h:
     80     - removed MRFEnergyEst
     81     - added MRanForestCalc
     82
     83   * mtools/MTFillMatrix.[h,cc]:
     84     - added the possibility to use PreCuts befre filling the matrix
     85     - added ReadEnv
     86
     87
     88
    2089 2005/11/18 Thomas Bretz
    2190
  • trunk/MagicSoft/Mars/Makefile

    r7149 r7413  
    6767          mjobs \
    6868          mjoptim \
     69          mjtrain \
    6970          mtools \
    7071          mmuon
  • trunk/MagicSoft/Mars/NEWS

    r7409 r7413  
    1313   - general: MTFillMatrix (the class to fill one or two MHMatrix from
    1414     files) now allows adding a pre-cut like in the optimization. E.g. this
    15      is useful to perform g/h-separation cuts before training the random forest.
    16 
    17    - general: Updated the random forest classes to support also the
    18      regression method implemented by Thomas H.
    19 
    20    - general: added new tutorial macro how to train the random forest
    21      for energy estimation (macros/optim/rfenergyest.C)
     15     is useful to perform g/h-separation cuts before training the random
     16     forest.
     17
     18   - RanForest:
     19      + Updated the random forest classes to support also the
     20        regression method implemented by Thomas H.
     21      + added new tutorial macro how to train the random forest
     22        for energy estimation (macros/optim/rfenergyest.C)
     23      + new classes to train the random forest (still in development)
     24        mjtrain/MJTrainEnergy, mjtrain/MJTrainDisp, mjtrain/MJTrainSeparation
     25      + new tutorial macros for random forest training in macros/train
     26        trainenergy.C, traindisp.C, trainseparation.C
     27
     28   - ganymed: In addition to the Hadronness calculator (CalcHadronness)
     29     a new option was implemented to estimate Disp (CalcDisp)
    2230
    2331   - ganymed: Implemented two new options which allow
  • trunk/MagicSoft/Mars/ganymed.rc

    r7362 r7413  
    3939BinningDist.Raw:         25    0    2
    4040BinningMaxDist.Raw:      25    0    2
     41
     42# -------------------------------------------------------------------------
     43# Using the starguider for pointing correction.
     44# To switch it off use "MPointingDevCalc.MaxAbsDev: -1"
     45# -------------------------------------------------------------------------
     46#MPointingDevCalc.MaxAbsDev:   15
     47#MPointingDevCalc.NumMinStars: 8
     48#MPointingDevCalc.NsbLevel:    3.0
     49#MPointingDevCalc.NsbMin:      30
     50#MPointingDevCalc.NsbMax:      60
     51
     52# -------------------------------------------------------------------------
     53# Setup the misfocussing correction
     54# -------------------------------------------------------------------------
     55#MSrcPosCorrect.Dx: 0
     56#MSrcPosCorrect.Dy: 0
    4157
    4258# -------------------------------------------------------------------------
     
    8399Cut0.Condition: ({0} || {1}) && {2} && {3} && {4} && {5}
    84100Cut0.0: MImagePar.fNumSatPixelsHG < 1
    85 Cut0.1: MHillas.GetArea*(MGeomCam.fConvMm2Deg^2) > (0.002*MImagePar.fNumSatPixelsHG) + 0.025
     101Cut0.1: MHillas.GetArea*(MGeomCam.fConvMm2Deg^2) > (0.003*MImagePar.fNumSatPixelsHG) + 0.0325
    86102Cut0.2: MNewImagePar.fNumUsedPixels>5
    87103Cut0.3: MNewImagePar.fLeakage1 < 0.3
     
    90106
    91107# ---------- SETUP FOR ONOFF-MODE -----------
     108MAlphaFitter.ScaleMin: 0.137
     109MAlphaFitter.ScaleMax: 0.640
    92110MAlphaFitter.BackgroundFitMin: 0.137
    93 MAlphaFitter.BackgroundFitMax: 0.64
     111MAlphaFitter.BackgroundFitMax: 1.000
    94112MAlphaFitter.PolynomOrder: 1
    95 MAlphaFitter.ScaleMode: Background
     113MAlphaFitter.ScaleMode: OffRegion
    96114MAlphaFitter.SignalFunction: ThetaSq
    97115
     
    99117Cut1.ThetaCut: None
    100118Cut1.Param0:  1.3245
    101 Cut1.Param1:  0.189
    102 Cut1.Param2:  0.230
    103 Cut1.Param3:  5.320
    104 Cut1.Param4:  0.100
    105 Cut1.Param5: -0.0636
     119Cut1.Param1:  0.204
     120Cut1.Param2:  0.228
     121Cut1.Param3:  5.50
     122Cut1.Param4:  0.088
     123Cut1.Param5: -0.07
    106124Cut1.Param6:  8.2957
    107125Cut1.Param7:  0.8677
     
    137155
    138156# -------------------------------------------------------------------------
    139 # Use this to executa a task (eg to calc hadronness) after energy
    140 # before Cut1
     157# Use this to executa a task (eg to calc hadronness) before Cut1
    141158# -------------------------------------------------------------------------
    142159# CalcHadronness: MRanForestCalc
    143 # CalcHadronness.File: /home/tbretz/Mars.cvs/RFspot3cmAll.root
     160# CalcHadronness.File: rf-energy.root
     161# CalcHadronness.NameOutput: Hadronness
    144162# CalcHadronness.Debug: No
     163
     164# -------------------------------------------------------------------------
     165# Use this to executa a task (eg to calc disp) before Cut1
     166# -------------------------------------------------------------------------
     167# CalcDisp: MRanForestCalc
     168# CalcDisp.File: rf-disp.root
     169# CalcDisp.NameOutput: Disp
     170# CalcDisp.Debug: No
     171# Cut1.CalcDisp: No
  • trunk/MagicSoft/Mars/macros/dohtml.C

    r7301 r7413  
    7777    sourcedir += "mimage:";
    7878    sourcedir += "mjobs:";
     79    sourcedir += "mjoptim:";
     80    sourcedir += "mjtrain:";
    7981    sourcedir += "mmain:";
    8082    sourcedir += "mmc:";
  • trunk/MagicSoft/Mars/macros/rootlogon.C

    r5694 r7413  
    144144    gInterpreter->AddIncludePath(dir+"mimage");
    145145    gInterpreter->AddIncludePath(dir+"mjobs");
     146    gInterpreter->AddIncludePath(dir+"mjoptim");
     147    gInterpreter->AddIncludePath(dir+"mjtrain");
    146148    gInterpreter->AddIncludePath(dir+"mmain");
    147149    gInterpreter->AddIncludePath(dir+"mmc");
  • trunk/MagicSoft/Mars/mbase/MParList.cc

    r6915 r7413  
    6262// --------------------------------------------------------------------------
    6363//
    64 //  default constructor
    6564//  creates an empty list
    6665//
    67 MParList::MParList(const char *name, const char *title)
     66void MParList::Init(const char *name, const char *title)
    6867{
    6968    fName  = name  ? name  : gsDefName.Data();
     
    8382}
    8483
     84
     85// --------------------------------------------------------------------------
     86//
     87//  default constructor
     88//  creates an empty list
     89//
     90MParList::MParList(const char *name, const char *title)
     91{
     92    Init(name, title);
     93}
     94
    8595// --------------------------------------------------------------------------
    8696//
     
    90100//  entries)
    91101//
    92 MParList::MParList(MParList &ts)
    93 {
     102MParList::MParList(const MParList &ts, const char *name, const char *title)
     103{
     104    Init(name, title);
     105
    94106    fContainer->AddAll(ts.fContainer);
    95107}
  • trunk/MagicSoft/Mars/mbase/MParList.h

    r4828 r7413  
    3737    void StreamPrimitive(ofstream &out) const;
    3838
     39    void Init(const char *name, const char *title);
     40
    3941public:
    4042    enum { kDoNotReset = BIT(17), kIsProcessing = BIT(18) };
    4143
    4244    MParList(const char *name=NULL, const char *title=NULL);
    43     MParList(MParList &ts);
     45    MParList(const MParList &ts, const char *name=NULL, const char *title=NULL);
    4446
    4547    virtual ~MParList();
  • trunk/MagicSoft/Mars/mhbase/MBinning.cc

    r7151 r7413  
    242242// lowest <lo> and highest <up> Edge (of your histogram)
    243243//
    244 void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up)
    245 {
     244void MBinning::SetEdgesLin(Int_t nbins, Axis_t lo, Axis_t up)
     245{
     246    if (nbins<=0)
     247    {
     248        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
     249        nbins = 1;
     250    }
     251
     252    if (up<lo)
     253    {
     254        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
     255
     256        const Axis_t dummy(lo);
     257        lo = up;
     258        up = dummy;
     259    }
     260
    246261    const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
    247262    fEdges.Set(nbins+1);
     
    402417// lowest <lo> and highest <up> Edge (of your histogram)
    403418//
    404 void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up)
     419void MBinning::SetEdgesLog(Int_t nbins, Axis_t lo, Axis_t up)
    405420{
    406421    // if (lo==0) ...
     422    if (nbins<=0)
     423    {
     424        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
     425        nbins = 1;
     426    }
     427
     428    if (up<lo)
     429    {
     430        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
     431
     432        const Axis_t dummy(lo);
     433        lo = up;
     434        up = dummy;
     435    }
    407436
    408437    const Double_t binsize = log10(up/lo)/nbins;
     
    419448// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
    420449//
    421 void MBinning::SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up)
    422 {
     450void MBinning::SetEdgesCos(Int_t nbins, Axis_t lo, Axis_t up)
     451{
     452    if (nbins<=0)
     453    {
     454        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
     455        nbins = 1;
     456    }
     457
     458    if (up<lo)
     459    {
     460        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
     461
     462        const Axis_t dummy(lo);
     463        lo = up;
     464        up = dummy;
     465    }
     466
    423467    // if (lo==0) ...
    424468    const Axis_t ld = lo/kRad2Deg;
     
    436480// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
    437481//
    438 void MBinning::SetEdgesASin(const Int_t nbins, Axis_t lo, Axis_t up)
    439 {
     482void MBinning::SetEdgesASin(Int_t nbins, Axis_t lo, Axis_t up)
     483{
     484    if (nbins<=0)
     485    {
     486        *fLog << warn << "WARNING - Number of bins cannot be <= 0... reset to 1." << endl;
     487        nbins = 1;
     488    }
     489
     490    if (up<lo)
     491    {
     492        *fLog << warn << "WARNING - Upper edge must be greater than lower edge... exchanging." << endl;
     493
     494        const Axis_t dummy(lo);
     495        lo = up;
     496        up = dummy;
     497    }
     498
    440499    const Double_t binsize = nbins<=0 ? 0 : (up-lo)/nbins;
    441500    fEdges.Set(nbins+1);
  • trunk/MagicSoft/Mars/mhbase/MBinning.h

    r7142 r7413  
    6262    void SetEdges(const MBinning &bins) { SetEdges(bins.fEdges); fType = bins.fType; fTitle = bins.fTitle; }
    6363    void SetEdges(const TH1 &h, const Char_t axis='x');
    64     void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up);
     64    void SetEdges(Int_t nbins, const Axis_t lo, const Axis_t up) { SetEdgesLin(nbins, lo, up); }
    6565    void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up, const char *opt);
    66     void SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up);
    67     void SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up);
    68     void SetEdgesASin(const Int_t nbins, Axis_t lo, Axis_t up);
     66    void SetEdgesLin(Int_t nbins, Axis_t lo, Axis_t up);
     67    void SetEdgesLog(Int_t nbins, Axis_t lo, Axis_t up);
     68    void SetEdgesCos(Int_t nbins, Axis_t lo, Axis_t up);
     69    void SetEdgesASin(Int_t nbins, Axis_t lo, Axis_t up);
    6970
    7071    Int_t FindLoEdge(Double_t val) const
  • trunk/MagicSoft/Mars/mhbase/MH3.cc

    r6917 r7413  
    561561    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(fHist);
    562562    pad->SetBorderMode(0);
     563    pad->SetGridx();
     564    pad->SetGridy();
    563565
    564566    fHist->SetFillStyle(4000);
  • trunk/MagicSoft/Mars/mhbase/MHMatrix.cc

    r7399 r7413  
    209209}
    210210
     211void MHMatrix::AddColumns(const TCollection &list)
     212{
     213    TIter Next(&list);
     214    TObject *obj = 0;
     215    while ((obj=Next()))
     216        AddColumn(obj->GetName());
     217}
     218
     219
    211220// --------------------------------------------------------------------------
    212221//
  • trunk/MagicSoft/Mars/mhbase/MHMatrix.h

    r7399 r7413  
    7373    Int_t AddColumn(const char *name);
    7474    void AddColumns(MDataArray *mat);
     75    void AddColumns(const TCollection &list);
    7576
    7677    //MDataArray *GetColumns() { return fData; }
  • trunk/MagicSoft/Mars/mhist/MHHadronness.cc

    r6890 r7413  
    209209    const MParameterD &had = par ? *(MParameterD*)par : *fHadronness;
    210210
    211     const Double_t h = had.GetVal();
     211    const Double_t h = TMath::Min(TMath::Max(had.GetVal(), 0.), 1.);
    212212
    213213    if (TMath::IsNaN(h))
  • trunk/MagicSoft/Mars/mjobs/MDataSet.cc

    r7409 r7413  
    196196    fName  = fname;
    197197
     198    fSequencesOn.SetOwner();
     199    fSequencesOff.SetOwner();
     200
    198201    const char *expname = gSystem->ExpandPathName(fname);
    199202
  • trunk/MagicSoft/Mars/mjobs/MDataSet.h

    r7380 r7413  
    5454
    5555public:
    56     MDataSet() : fNumAnalysis((UInt_t)-1) { }
     56    MDataSet(Int_t num=(UInt_t)-1) : fNumAnalysis(num) { }
    5757    MDataSet(const char *fname, TString sequences="", TString data="");
     58
     59    void Copy(TObject &obj) const
     60    {
     61        MDataSet &ds = (MDataSet&)obj;
     62        ds.fNumAnalysis = fNumAnalysis;
     63        ds.fNumSequencesOn = fNumSequencesOn;
     64        ds.fNumSequencesOff = fNumSequencesOff;
     65        ds.fNameSource = fNameSource;
     66        ds.fCatalog = fCatalog;
     67        ds.fComment = fComment;
     68        ds.fIsWobbleMode = fIsWobbleMode;
     69
     70        TObject *o=0;
     71
     72        ds.fSequencesOn.Delete();
     73        ds.fSequencesOff.Delete();
     74
     75        TIter NextOn(&fSequencesOn);
     76        while ((o=NextOn()))
     77            ds.fSequencesOn.Add(o->Clone());
     78
     79        TIter NextOff(&fSequencesOff);
     80        while ((o=NextOff()))
     81            ds.fSequencesOff.Add(o->Clone());
     82    }
    5883
    5984    // Getter
  • trunk/MagicSoft/Mars/mjobs/MJCut.cc

    r7390 r7413  
    8686    : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
    8787    fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kTRUE),
    88     fNameHist("MHThetaSq"), fCalcHadronness(0)
     88    fNameHist("MHThetaSq"), fCalcHadronness(0), fCalcDisp(0)
    8989{
    9090    fName  = name  ? name  : "MJCut";
     
    102102    if (fCalcHadronness)
    103103        delete fCalcHadronness;
     104    if (fCalcDisp)
     105        delete fCalcDisp;
    104106}
    105107
     
    160162        delete fCalcHadronness;
    161163    fCalcHadronness = task ? (MTask*)task->Clone() : 0;
     164}
     165
     166// --------------------------------------------------------------------------
     167//
     168// Setup a task calculating disp. The given task is cloned.
     169//
     170void MJCut::SetDispCalculator(const MTask *task)
     171{
     172    if (fCalcDisp)
     173        delete fCalcDisp;
     174    fCalcDisp = task ? (MTask*)task->Clone() : 0;
    162175}
    163176
     
    540553    taskenv2.SetDefault(fCalcHadronness);
    541554
     555    MTaskEnv taskenv3("CalcDisp");
     556    taskenv3.SetDefault(fCalcDisp);
     557
    542558    MFillH fill1a("MHHillasOffPre  [MHHillas]",      "MHillas",      "FillHillasPre");
    543559    MFillH fill2a("MHHillasOffPost [MHHillas]",      "MHillas",      "FillHillasPost");
     
    579595    //tlist2.AddToList(&taskenv1);
    580596    tlist2.AddToList(&taskenv2);
     597    tlist2.AddToList(&taskenv3);
    581598    tlist2.AddToList(&cont0);
    582599    if (write0)
  • trunk/MagicSoft/Mars/mjobs/MJCut.h

    r7109 r7413  
    3333    //MTask *fEstimateEnergy;
    3434    MTask *fCalcHadronness;
     35    MTask *fCalcDisp;
    3536
    3637    TString  GetOutputFile(UInt_t num) const;
     
    6869    //void SetEnergyEstimator(const MTask *task=0);
    6970    void SetHadronnessCalculator(const MTask *task=0);
     71    void SetDispCalculator(const MTask *task=0);
    7072
    7173    ClassDef(MJCut, 0) // Standard program to perform g/h-separation cuts
  • trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc

    r7169 r7413  
    7171
    7272// histograms
     73#include "MH3.h"
     74#include "MBinning.h"
    7375#include "../mhflux/MAlphaFitter.h"
    7476#include "../mhflux/MHThetaSq.h"
     
    9193//
    9294// Read all events from file which do match rules and optimize
    93 // energy estimator.
     95// disp estimator.
    9496//
    9597Bool_t MJOptimizeDisp::RunDisp(const char *fname, const char *rule, MTask *weights)
    9698{
    97     fLog->Separator("Preparing Energy optimization");
     99    fLog->Separator("Preparing Disp optimization");
    98100
    99101    MParList parlist;
     
    124126    const Int_t num1 = m.AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
    125127    const Int_t num2 = m.AddColumn("MHillasSrc.fAlpha*kDegToRad");
     128    const Int_t num3 = m.AddColumn("MHillas.fSize");
    126129
    127130    MHThetaSq hist;
    128131    hist.SkipHistTime();
    129132    hist.SkipHistTheta();
    130     hist.SkipHistEnergy();
    131     hist.InitMapping(&m);
     133    //hist.SkipHistEnergy();
     134    //hist.ForceUsingSize();
     135    hist.InitMapping(&m, 1);
    132136
    133137    MFDataMember filter("DataType.fVal", '>', 0.5);
     
    171175    MMatrixLoop loop(&m);
    172176
     177    const char *n3   = Form("M[%d]", num3);
     178    MH3 hdisp(n3, "sqrt(ThetaSquared.fVal)");
     179    hdisp.SetTitle("\\vartheta^{2} distribution vs. Size:Size [phe]:\\vartheta^{2} [\\circ^{2}]");
     180
     181    MBinning binsx(100, 10, 100000, "BinningMH3X", "log");
     182    MBinning binsy(100, 0,  2,      "BinningMH3Y", "lin");
     183
     184    parlist.AddToList(&binsx);
     185    parlist.AddToList(&binsy);
     186
     187    MFillH fillh2(&hdisp);
     188    fillh2.SetDrawOption("blue profx");
     189
    173190    tasklist.AddToList(&loop);
    174191    tasklist.AddToList(&calc1);
     
    177194        tasklist.AddToList(weights);
    178195    tasklist.AddToList(&fill);
     196    tasklist.AddToList(&fillh2);
    179197    tasklist.AddToList(&eval);
    180198
  • trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc

    r2296 r7413  
    4343#include "MRanTree.h"
    4444#include "MRanForest.h"
     45#include "MDataArray.h"
    4546
    4647#include "MLog.h"
     
    5758//
    5859MHRanForestGini::MHRanForestGini(Int_t nbins, const char *name, const char *title)
     60    : fRules(0.01, 0.01, 0.99, 0.99)
    5961{
    6062    //
     
    6466    fTitle = title ? title : "Measure of importance of Random Forest-input parameters";
    6567
    66     fGraphGini = new TGraph;
    67     fGraphGini->SetTitle("Importance of RF-input parameters measured by mean Gini decrease");
    68     fGraphGini->SetMarkerStyle(kFullDotSmall);
    69 }
    70 
    71 // --------------------------------------------------------------------------
    72 //
    73 // Delete the histograms.
    74 //
    75 MHRanForestGini::~MHRanForestGini()
    76 {
    77     delete fGraphGini;
     68    fGraphGini.SetNameTitle("Gini", "Importance of RF-input parameters measured by mean Gini decrease");
     69    fGraphGini.SetMarkerStyle(kFullDotMedium);
     70
     71    fGraphError.SetNameTitle("ResErr", "Resolution/Error versus train step");
     72    fGraphError.SetMarkerStyle(kFullDotMedium);
     73
     74    fGraphNodes.SetNameTitle("Nodes", "Number of nodes versus train step");
     75    fGraphNodes.SetMarkerStyle(kFullDotMedium);
     76
     77    fRules.SetTextAlign(13);
     78    fRules.SetTextSize(0.05);
    7879}
    7980
     
    104105Bool_t MHRanForestGini::Fill(const MParContainer *par, const Stat_t w)
    105106{
     107    MRanTree *t = fRanForest->GetCurTree();
     108
    106109    for (Int_t i=0;i<fRanForest->GetNumDim();i++)
    107         fGini[i]+=fRanForest->GetCurTree()->GetGiniDec(i);
     110        fGini[i] += t->GetGiniDec(i);
     111
     112    fGraphError.SetPoint(fGraphError.GetN(), GetNumExecutions(), t->GetError());
     113    fGraphNodes.SetPoint(fGraphError.GetN(), GetNumExecutions(), t->GetNumEndNodes());
    108114
    109115    return kTRUE;
     
    117123    const Int_t n = fGini.GetSize();
    118124
    119     fGraphGini->Set(n);
    120 
    121     Stat_t max=0;
    122     Stat_t min=0;
     125    fGraphGini.Set(n);
     126
    123127    for (Int_t i=0; i<n; i++)
    124128    {
    125         fGini[i] /= fRanForest->GetNumTrees();
    126         fGini[i] /= fRanForest->GetNumData();
    127 
    128         const Stat_t ip = i+1;
    129         const Stat_t ig = fGini[i];
    130 
    131         if (ig>max) max=ig;
    132         if (ig<min) min=ig;
    133 
    134         fGraphGini->SetPoint(i,ip,ig);
    135     }
    136 
    137     // This is used in root>3.04/? so that SetMaximum/Minimum can succeed
    138     fGraphGini->GetHistogram();
    139 
    140     fGraphGini->SetMaximum(1.05*max);
    141     fGraphGini->SetMinimum(0.95*min);
     129        fGini[i] /= fRanForest->GetNumTrees()*fRanForest->GetNumData();
     130        fGraphGini.SetPoint(i, i+1, fGini[i]);
     131    }
     132
     133    fRules.AddText("");
     134    const MDataArray &arr = *fRanForest->GetRules();
     135    int i;
     136    for (i=0; i<arr.GetNumEntries(); i++)
     137    {
     138        TString s;
     139        s += i+1;
     140        s += ") ";
     141        s += arr.GetRule(i);
     142        fRules.AddText(s);
     143    }
     144    for (; i<20; i++)
     145        fRules.AddText("");
    142146
    143147    return kTRUE;
     
    150154void MHRanForestGini::Draw(Option_t *)
    151155{
    152     if (fGraphGini->GetN()==0)
     156    if (fGraphGini.GetN()==0)
    153157        return;
    154158
     
    158162    AppendPad("");
    159163
    160     fGraphGini->Draw("ALP");
    161     pad->Modified();
    162     pad->Update();
    163 
    164     TH1 *h = fGraphGini->GetHistogram();
    165     if (!h)
    166         return;
    167 
    168     //h->GetXaxis()->SetRangeUser(0, fGini.GetSize()+1);
    169     h->SetXTitle("No.of RF-input parameter");
    170     h->SetYTitle("Mean decrease in Gini-index [a.u.]");
    171 
    172     pad->Modified();
    173     pad->Update();
    174 }
     164    pad->Divide(2,2);
     165
     166    pad->cd(1);
     167    gPad->SetBorderMode(0);
     168    gPad->SetGridx();
     169    gPad->SetGridy();
     170    fGraphGini.Draw("ALP");
     171
     172    TH1 *h = fGraphGini.GetHistogram();
     173    if (h)
     174    {
     175        h->SetXTitle("No.of RF-input parameter");
     176        h->SetYTitle("Mean decrease in Gini-index [au]");
     177        h->GetXaxis()->SetNdivisions(10);
     178    }
     179
     180    pad->cd(2);
     181    gPad->SetBorderMode(0);
     182    fGraphError.Draw("ALP");
     183    h = fGraphError.GetHistogram();
     184    if (h)
     185    {
     186        h->SetXTitle("Train step/Tree number");
     187        h->SetYTitle("Error/Resolution");
     188    }
     189
     190    pad->cd(3);
     191    gPad->SetBorderMode(0);
     192    fGraphNodes.Draw("ALP");
     193    h = fGraphNodes.GetHistogram();
     194    if (h)
     195    {
     196        h->SetXTitle("Train step/Tree number");
     197        h->SetYTitle("Number of end nodes");
     198    }
     199
     200    pad->cd(4);
     201    gPad->SetBorderMode(0);
     202    fRules.Draw();
     203}
  • trunk/MagicSoft/Mars/mranforest/MHRanForestGini.h

    r2173 r7413  
    99#include <TArrayF.h>
    1010#endif
     11#ifndef ROOT_TGraph
     12#include <TGraph.h>
     13#endif
     14#ifndef ROOT_TPaveText
     15#include <TPaveText.h>
     16#endif
    1117
    12 class TH1D;
    13 class TGraph;
    1418class MParList;
    1519class MRanForest;
     
    2226
    2327    TArrayF fGini;           //!
    24     TGraph *fGraphGini;      //->
     28
     29    TGraph  fGraphGini;
     30    TGraph  fGraphError;
     31    TGraph  fGraphNodes;
     32
     33    TPaveText fRules;
    2534
    2635public:
    2736    MHRanForestGini(Int_t nbins=100, const char *name=NULL, const char *title=NULL);
    28     ~MHRanForestGini();
    29 
    30     TGraph *GetGraphGini() const  { return fGraphGini; }
    3137
    3238    Bool_t SetupFill(const MParList *plist);
  • trunk/MagicSoft/Mars/mranforest/MRanForest.cc

    r7409 r7413  
    166166        if(grid[i]>=grid[i+1])
    167167        {
    168             *fLog<<inf<<"Grid points must be in increasing order! Ignoring grid."<<endl;
     168            *fLog<<warn<<"Grid points must be in increasing order! Ignoring grid."<<endl;
    169169            return;
    170170        }
     
    225225            fClass[j] = int(fHadTrue[j]+0.5);
    226226    }
    227 
    228     return;
    229227}
    230228
     
    272270Bool_t MRanForest::AddTree(MRanTree *rantree=NULL)
    273271{
    274     fRanTree = rantree ? rantree:fRanTree;
     272    fRanTree = rantree ? rantree : fRanTree;
    275273
    276274    if (!fRanTree) return kFALSE;
     
    287285    // access matrix, copy last column (target) preliminarily
    288286    // into fHadTrue
    289     TMatrix mat_tmp = mat->GetM();
    290     int dim         = mat_tmp.GetNcols();
    291     int numdata     = mat_tmp.GetNrows();
    292 
    293     fMatrix=new TMatrix(mat_tmp);
     287    fMatrix = new TMatrix(mat->GetM());
     288
     289    int dim     = fMatrix->GetNcols()-1;
     290    int numdata = fMatrix->GetNrows();
    294291
    295292    fHadTrue.Set(numdata);
     
    297294
    298295    for (Int_t j=0;j<numdata;j++)
    299         fHadTrue[j] = (*fMatrix)(j,dim-1);
     296        fHadTrue[j] = (*fMatrix)(j,dim);
    300297
    301298    // remove last col
    302     fMatrix->ResizeTo(numdata,dim-1);
    303     dim=fMatrix->GetNcols();
     299    fMatrix->ResizeTo(numdata,dim);
    304300
    305301    //-------------------------------------------------------------------
     
    308304    fClass.Reset(0);
    309305
    310     if(fClassify) PrepareClasses();
     306    if (fClassify)
     307        PrepareClasses();
    311308
    312309    //-------------------------------------------------------------------
    313310    // allocating and initializing arrays
    314     fHadEst.Set(numdata);       fHadEst.Reset(0);
    315     fNTimesOutBag.Set(numdata); fNTimesOutBag.Reset(0);
    316     fDataSort.Set(dim*numdata); fDataSort.Reset(0);
    317     fDataRang.Set(dim*numdata); fDataRang.Reset(0);
     311    fHadEst.Set(numdata);
     312    fHadEst.Reset(0);
     313
     314    fNTimesOutBag.Set(numdata);
     315    fNTimesOutBag.Reset(0);
     316
     317    fDataSort.Set(dim*numdata);
     318    fDataSort.Reset(0);
     319
     320    fDataRang.Set(dim*numdata);
     321    fDataRang.Reset(0);
    318322
    319323    if(fWeight.GetSize()!=numdata)
     
    326330    //-------------------------------------------------------------------
    327331    // setup rules to be used for classification/regression
    328     MDataArray *allrules=(MDataArray*)mat->GetColumns();
     332    const MDataArray *allrules=(MDataArray*)mat->GetColumns();
    329333    if(allrules==NULL)
    330334    {
     
    333337    }
    334338
    335     fRules=new MDataArray(); fRules->Reset();
    336     TString target_rule;
    337 
    338     for(Int_t i=0;i<dim+1;i++)
    339     {
    340         MData &data=(*allrules)[i];
    341         if(i<dim)
    342             fRules->AddEntry(data.GetRule());
    343         else
    344             target_rule=data.GetRule();
    345     }
    346 
    347     *fLog << inf <<endl;
    348     *fLog << inf <<"Setting up RF for training on target:"<<endl<<" "<<target_rule.Data()<<endl;
    349     *fLog << inf <<"Following rules are used as input to RF:"<<endl;
    350 
    351     for(Int_t i=0;i<dim;i++)
    352     {
    353         MData &data=(*fRules)[i];
    354         *fLog<<inf<<" "<<i<<") "<<data.GetRule()<<endl<<flush;
    355     }
    356 
    357     *fLog << inf <<endl;
     339    fRules = new MDataArray();
     340    fRules->Reset();
     341
     342    const TString target_rule = (*allrules)[dim];
     343    for (Int_t i=0;i<dim;i++)
     344        fRules->AddEntry((*allrules)[i].GetRule());
     345
     346    *fLog << inf << endl;
     347    *fLog << "Setting up RF for training on target:" << endl;
     348    *fLog << " " << target_rule.Data() << endl;
     349    *fLog << "Following rules are used as input to RF:" << endl;
     350    for (Int_t i=0;i<dim;i++)
     351        *fLog << " " << i << ") " << (*fRules)[i].GetRule() << endl;
     352
     353    *fLog << endl;
    358354
    359355    //-------------------------------------------------------------------
    360356    // prepare (sort) data for fast optimization algorithm
    361     if(!CreateDataSort()) return kFALSE;
     357    if (!CreateDataSort())
     358        return kFALSE;
    362359
    363360    //-------------------------------------------------------------------
    364361    // access and init tree container
    365362    fRanTree = (MRanTree*)plist->FindCreateObj("MRanTree");
    366 
    367363    if(!fRanTree)
    368364    {
     
    371367    }
    372368
     369    const Int_t tryest = TMath::Nint(TMath::Sqrt(dim)+0.5);
     370
     371    *fLog << inf << endl;
     372    *fLog << "Following input for the tree growing are used:"<<endl;
     373    *fLog << " Number of Trees : "<<fNumTrees<<endl;
     374    *fLog << " Number of Trials: "<<(fNumTry==0?tryest:fNumTry)<<(fNumTry==0?" (auto)":"")<<endl;
     375    *fLog << " Final Node size : "<<fNdSize<<endl;
     376    *fLog << " Using Grid:       "<<(fGrid.GetSize()>0?"Yes":"No")<<endl;
     377    *fLog << " Number of Events: "<<numdata<<endl;
     378    *fLog << " Number of Params: "<<dim<<endl;
     379
     380    if(fNumTry==0)
     381    {
     382        fNumTry=tryest;
     383        *fLog << inf << endl;
     384        *fLog << "Set no. of trials to the recommended value of round("<< TMath::Sqrt(dim) <<") = ";
     385        *fLog << fNumTry << endl;
     386
     387
     388    }
     389    fRanTree->SetNumTry(fNumTry);
    373390    fRanTree->SetClassify(fClassify);
    374391    fRanTree->SetNdSize(fNdSize);
    375392
    376     if(fNumTry==0)
    377     {
    378         double ddim = double(dim);
    379 
    380         fNumTry=int(sqrt(ddim)+0.5);
    381         *fLog<<inf<<endl;
    382         *fLog<<inf<<"Set no. of trials to the recommended value of round("<<sqrt(ddim)<<") = ";
    383         *fLog<<inf<<fNumTry<<endl;
    384 
    385     }
    386     fRanTree->SetNumTry(fNumTry);
    387 
    388     *fLog<<inf<<endl;
    389     *fLog<<inf<<"Following settings for the tree growing are used:"<<endl;
    390     *fLog<<inf<<" Number of Trees : "<<fNumTrees<<endl;
    391     *fLog<<inf<<" Number of Trials: "<<fNumTry<<endl;
    392     *fLog<<inf<<" Final Node size : "<<fNdSize<<endl;
    393 
    394393    fTreeNo=0;
    395394
     
    415414    if (fTreeNo==1)
    416415    {
    417         *fLog << inf << endl;
    418         *fLog << underline;
     416        *fLog << inf << endl << underline;
    419417
    420418        if(calcResolution)
     
    435433    TArrayF winbag(numdata); // Initialization includes filling with 0
    436434
    437     float square=0; float mean=0;
     435    float square=0;
     436    float mean=0;
    438437
    439438    for (Int_t n=0; n<numdata; n++)
     
    483482    for (Int_t ievt=0;ievt<numdata;ievt++)
    484483    {
    485         if (jinbag[ievt]>0) continue;
     484        if (jinbag[ievt]>0)
     485            continue;
    486486
    487487        fHadEst[ievt] +=fRanTree->TreeHad((*fMatrix), ievt);
     
    491491
    492492    Int_t n=0;
    493     double ferr=0;
     493    Float_t ferr=0;
    494494
    495495    for (Int_t ievt=0;ievt<numdata;ievt++)
     
    508508    //-------------------------------------------------------------------
    509509    // give running output
    510     *fLog << inf << setw(5)  << fTreeNo;
    511     *fLog << inf << setw(18) << fRanTree->GetNumEndNodes();
    512     *fLog << inf << Form("%18.2f", ferr*100.);
    513     *fLog << inf << endl;
     510    *fLog << setw(5)  << fTreeNo;
     511    *fLog << setw(18) << fRanTree->GetNumEndNodes();
     512    *fLog << Form("%18.2f", ferr*100.);
     513    *fLog << endl;
     514
     515    fRanTree->SetError(ferr);
    514516
    515517    // adding tree to forest
  • trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc

    r6893 r7413  
    1616!
    1717!
    18 !   Author(s): Thomas Hengstebeck 3/2003 <mailto:hengsteb@alwa02.physik.uni-siegen.de>
    19 !
    20 !   Copyright: MAGIC Software Development, 2000-2003
     18!   Author(s): Thomas Hengstebeck 2/2005 <mailto:hengsteb@physik.hu-berlin.de>
     19!   Author(s): Thomas Bretz 8/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
     20!
     21!   Copyright: MAGIC Software Development, 2000-2005
    2122!
    2223!
     
    2728//  MRanForestCalc
    2829//
    29 //  Calculates the hadroness of an event. It calculates a mean value of all
    30 //  classifications by the trees in a previously grown random forest.
    31 //
    32 //  To use only n trees for your calculation use:
    33 //  MRanForestCalc::SetUseNumTrees(n);
    3430//
    3531////////////////////////////////////////////////////////////////////////////
     
    3834#include <TVector.h>
    3935
    40 #include "MHMatrix.h" // must be before MLogManip.h
    41 #include "MDataArray.h"
     36#include "MHMatrix.h"
    4237
    4338#include "MLog.h"
    4439#include "MLogManip.h"
    4540
     41#include "MData.h"
     42#include "MDataArray.h"
     43
     44#include "MRanForest.h"
     45#include "MParameters.h"
     46
    4647#include "MParList.h"
    47 
    48 #include "MRanTree.h"
    49 #include "MRanForest.h"
    50 
    51 #include "MParameters.h"
    52 
     48#include "MTaskList.h"
    5349#include "MEvtLoop.h"
    54 #include "MTaskList.h"
     50#include "MRanForestGrow.h"
    5551#include "MFillH.h"
    56 #include "MStatusDisplay.h"
    57 #include "MRanForestGrow.h"
    58 #include "MRanForestFill.h"
    59 
    60 #include "MWriteRootFile.h"
    61 #include "MReadTree.h"
    6252
    6353ClassImp(MRanForestCalc);
     
    6555using namespace std;
    6656
    67 static const TString gsDefName  = "MRanForestCalc";
    68 static const TString gsDefTitle = "Tree Classification Loop 1/2";
    69 
    70 // --------------------------------------------------------------------------
    71 //
    72 // Setup histograms and the number of distances which are used for
    73 // avaraging in CalcDist
    74 //
     57const TString MRanForestCalc::gsDefName    = "MRanForestCalc";
     58const TString MRanForestCalc::gsDefTitle   = "RF for energy estimation";
     59
     60const TString MRanForestCalc::gsNameOutput = "RanForestOut";
     61
    7562MRanForestCalc::MRanForestCalc(const char *name, const char *title)
    76     : fNum(100), fHadronnessName("MHadronness"), fData(NULL), fRanForest(0), fRanTree(0)
    77 {
    78     //
    79     //   set the name and title of this object
    80     //
     63    : fDebug(kFALSE), fData(0), fRFOut(0),
     64    fNumTrees(-1), fNumTry(-1), fNdSize(-1), fNumObsoleteVariables(1),
     65    fTestMatrix(0), fEstimationMode(kMean)
     66{
    8167    fName  = name  ? name  : gsDefName.Data();
    8268    fTitle = title ? title : gsDefTitle.Data();
    83 }
    84 
    85 // --------------------------------------------------------------------------
    86 //
    87 // Delete the data chains
    88 //
     69
     70    // FIXME:
     71    fNumTrees = 100; //100
     72    fNumTry   = 0;   //3   0 means: in MRanForest estimated best value will be calculated
     73    fNdSize   = 1;   //1   
     74}
     75
    8976MRanForestCalc::~MRanForestCalc()
    9077{
    91     if (!IsOwner())
    92         return;
    93 
    94     delete fRanForest;
    95     delete fRanTree;
    96 }
    97 
    98 // --------------------------------------------------------------------------
    99 //
    100 // Needs:
    101 //  - MatrixGammas  [MHMatrix]
    102 //  - MatrixHadrons [MHMatrix]
    103 //  - MHadronness
    104 //  - all data containers used to build the matrixes
    105 //
    106 // The matrix object can be filles using MFillH. And must be of the same
    107 // number of columns (with the same meaning).
    108 //
     78    fEForests.Delete();
     79}
     80
     81Int_t MRanForestCalc::Train(const MHMatrix &matrixtrain, const TArrayD &grid, Int_t ver)
     82{
     83    gLog.Separator("MRanForestCalc - Train");
     84
     85    if (!matrixtrain.GetColumns())
     86    {
     87        *fLog << err << "ERROR - MHMatrix does not contain rules... abort." << endl;
     88        return kFALSE;
     89    }
     90
     91    const Int_t ncols = matrixtrain.GetM().GetNcols();
     92    const Int_t nrows = matrixtrain.GetM().GetNrows();
     93    if (ncols<=0 || nrows <=0)
     94    {
     95        *fLog << err << "ERROR - No. of columns or no. of rows of matrixtrain equal 0 ... abort." << endl;
     96        return kFALSE;
     97    }
     98
     99    // rules (= combination of image par) to be used for energy estimation
     100    TFile fileRF(fFileName, "recreate");
     101    if (!fileRF.IsOpen())
     102    {
     103        *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl;
     104        return kFALSE;
     105    }
     106
     107    const Int_t nobs = fNumObsoleteVariables; // Number of obsolete columns
     108
     109    const MDataArray &dcol = *matrixtrain.GetColumns();
     110
     111    MDataArray usedrules;
     112    for (Int_t i=0; i<ncols; i++)
     113        if (i<ncols-nobs)  // -3 is important!!!
     114            usedrules.AddEntry(dcol[i].GetRule());
     115        else
     116            *fLog << inf << "Skipping " << dcol[i].GetRule() << " for training" << endl;
     117
     118    MDataArray rules(usedrules);
     119    rules.AddEntry(ver<2?"Classification":dcol[ncols-1].GetRule());
     120
     121    // prepare matrix for current energy bin
     122    TMatrix mat(matrixtrain.GetM());
     123
     124    // last column must be removed (true energy col.)
     125    mat.ResizeTo(nrows, ncols-nobs+1);
     126
     127    if (fDebug)
     128        gLog.SetNullOutput(kTRUE);
     129
     130    const Int_t nbins = ver>0 ? 1 : grid.GetSize()-1;
     131    for (Int_t ie=0; ie<nbins; ie++)
     132    {
     133        switch (ver)
     134        {
     135        case 0: // Replace Energy Grid by classification
     136            {
     137                Int_t irows=0;
     138                for (Int_t j=0; j<nrows; j++)
     139                {
     140                    const Double_t energy = matrixtrain.GetM()(j,ncols-1);
     141                    const Bool_t   inside = energy>grid[ie] && energy<=grid[ie+1];
     142
     143                    mat(j, ncols-nobs) = inside ? 1 : 0;
     144
     145                    if (inside)
     146                        irows++;
     147                }
     148                if (irows==0)
     149                    *fLog << warn << "WARNING - Skipping";
     150                else
     151                    *fLog << inf << "Training RF for";
     152
     153                *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irows << "/" << nrows << endl;
     154
     155                if (irows==0)
     156                    continue;
     157            }
     158            break;
     159
     160        case 1: // Use Energy as classifier
     161        case 2:
     162            for (Int_t j=0; j<nrows; j++)
     163                mat(j, ncols-nobs) = matrixtrain.GetM()(j,ncols-1);
     164            break;
     165        }
     166
     167        MHMatrix matrix(mat, &rules, "MatrixTrain");
     168
     169        MParList plist;
     170        MTaskList tlist;
     171        plist.AddToList(&tlist);
     172        plist.AddToList(&matrix);
     173
     174        MRanForest rf;
     175        rf.SetNumTrees(fNumTrees);
     176        rf.SetNumTry(fNumTry);
     177        rf.SetNdSize(fNdSize);
     178        rf.SetClassify(ver<2 ? 1 : 0);
     179        if (ver==1)
     180            rf.SetGrid(grid);
     181
     182        plist.AddToList(&rf);
     183
     184        MRanForestGrow rfgrow;
     185        tlist.AddToList(&rfgrow);
     186
     187        MFillH fillh("MHRanForestGini");
     188        tlist.AddToList(&fillh);
     189
     190        MEvtLoop evtloop;
     191        evtloop.SetParList(&plist);
     192        evtloop.SetDisplay(fDisplay);
     193        evtloop.SetLogStream(fLog);
     194
     195        if (!evtloop.Eventloop())
     196            return kFALSE;
     197
     198        if (fDebug)
     199            gLog.SetNullOutput(kFALSE);
     200
     201        if (ver==0)
     202        {
     203            // Calculate bin center
     204            const Double_t E = (TMath::Log10(grid[ie])+TMath::Log10(grid[ie+1]))/2;
     205
     206            // save whole forest
     207            rf.SetUserVal(E);
     208            rf.SetName(Form("%.10f", E));
     209        }
     210
     211        rf.Write();
     212    }
     213
     214    // save rules
     215    usedrules.Write("rules");
     216
     217    return kTRUE;
     218}
     219
     220Int_t MRanForestCalc::ReadForests(MParList &plist)
     221{
     222    TFile fileRF(fFileName, "read");
     223    if (!fileRF.IsOpen())
     224    {
     225        *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
     226        return kFALSE;
     227    }
     228
     229    fEForests.Delete();
     230
     231    TIter Next(fileRF.GetListOfKeys());
     232    TObject *o=0;
     233    while ((o=Next()))
     234    {
     235        MRanForest *forest=0;
     236        fileRF.GetObject(o->GetName(), forest);
     237        if (!forest)
     238            continue;
     239
     240        forest->SetUserVal(atof(o->GetName()));
     241
     242        fEForests.Add(forest);
     243    }
     244
     245    // Maybe fEForests[0].fRules yould be used instead?
     246
     247    if (fData->Read("rules")<=0)
     248    {
     249        *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
     250        return kFALSE;
     251    }
     252
     253    return kTRUE;
     254}
     255
    109256Int_t MRanForestCalc::PreProcess(MParList *plist)
    110257{
    111     if (!fRanForest)
    112     {
    113         fRanForest = (MRanForest*)plist->FindObject("MRanForest");
    114         if (!fRanForest)
     258    fRFOut = (MParameterD*)plist->FindCreateObj("MParameterD", fNameOutput);
     259    if (!fRFOut)
     260        return kFALSE;
     261
     262    fData = (MDataArray*)plist->FindCreateObj("MDataArray");
     263    if (!fData)
     264        return kFALSE;
     265
     266    if (!ReadForests(*plist))
     267    {
     268        *fLog << err << "Reading RFs failed... aborting." << endl;
     269        return kFALSE;
     270    }
     271
     272    *fLog << inf << "RF read from " << fFileName << endl;
     273
     274    if (fTestMatrix)
     275        return kTRUE;
     276
     277    fData->Print();
     278
     279    if (!fData->PreProcess(plist))
     280    {
     281        *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl;
     282        return kFALSE;
     283    }
     284
     285    return kTRUE;
     286}
     287
     288#include <TGraph.h>
     289#include <TF1.h>
     290Int_t MRanForestCalc::Process()
     291{
     292    TVector event;
     293    if (fTestMatrix)
     294        *fTestMatrix >> event;
     295    else
     296        *fData >> event;
     297
     298    // --------------- Single Tree RF -------------------
     299    if (fEForests.GetEntries()==1)
     300    {
     301        MRanForest *rf = (MRanForest*)fEForests[0];
     302        fRFOut->SetVal(rf->CalcHadroness(event));
     303        fRFOut->SetReadyToSave();
     304
     305        return kTRUE;
     306    }
     307
     308    // --------------- Multi Tree RF -------------------
     309    static TF1 f1("f1", "gaus");
     310
     311    Double_t sume = 0;
     312    Double_t sumh = 0;
     313    Double_t maxh = 0;
     314    Double_t maxe = 0;
     315
     316    Double_t max  = -1e10;
     317    Double_t min  =  1e10;
     318
     319    TIter Next(&fEForests);
     320    MRanForest *rf = 0;
     321
     322    TGraph g;
     323    while ((rf=(MRanForest*)Next()))
     324    {
     325        const Double_t h = rf->CalcHadroness(event);
     326        const Double_t e = rf->GetUserVal();
     327
     328        g.SetPoint(g.GetN(), e, h);
     329
     330        sume += e*h;
     331        sumh += h;
     332
     333        if (h>maxh)
    115334        {
    116             *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
    117             return kFALSE;
     335            maxh = h;
     336            maxe = e;
    118337        }
    119     }
    120 
    121     if (!fRanTree)
    122     {
    123         fRanTree = (MRanTree*)plist->FindObject("MRanTree");
    124         if (!fRanTree)
    125         {
    126             *fLog << err << dbginf << "MRanTree not found... aborting." << endl;
    127             return kFALSE;
    128         }
    129     }
    130 
    131     fData = fRanTree->GetRules();
    132 
    133     if (!fData)
    134     {
    135         *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl;
    136         return kFALSE;
    137     }
    138 
    139     if (!fData->PreProcess(plist))
    140     {
    141         *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
    142         return kFALSE;
    143     }
    144 
    145     fHadroness = (MParameterD*)plist->FindCreateObj("MParameterD", fHadronnessName);
    146     if (!fHadroness)
    147         return kFALSE;
     338        if (e>max)
     339            max = e;
     340        if (e<min)
     341            min = e;
     342    }
     343
     344    switch (fEstimationMode)
     345    {
     346    case kMean:
     347        fRFOut->SetVal(pow(10, sume/sumh));
     348        break;
     349    case kMaximum:
     350        fRFOut->SetVal(pow(10, maxe));
     351        break;
     352    case kFit:
     353        f1.SetParameter(0, maxh);
     354        f1.SetParameter(1, maxe);
     355        f1.SetParameter(2, 0.125);
     356        g.Fit(&f1, "Q0N");
     357        fRFOut->SetVal(pow(10, f1.GetParameter(1)));
     358        break;
     359    }
     360
     361    fRFOut->SetReadyToSave();
    148362
    149363    return kTRUE;
     
    153367//
    154368//
    155 Int_t MRanForestCalc::Process()
    156 {
    157     // first copy the data from the data array to a vector event
    158     TVector event;
    159     *fData >> event;
    160 
    161     Double_t hadroness=fRanForest->CalcHadroness(event);
    162     fHadroness->SetVal(hadroness);
    163 
    164     return kTRUE;
    165 }
    166 
    167 Bool_t MRanForestCalc::Grow(MHMatrix *matrixg,MHMatrix *matrixh,Int_t ntree,
    168                             Int_t numtry,Int_t ndsize,const char* treefile,
    169                             const char* treename,const char* contname,
    170                             const char* hgininame)
    171 {
    172 
    173     treename  = treename  ? treename  : "Tree";
    174     contname  = contname  ? contname  : "MRanTree";
    175     hgininame  = hgininame  ? hgininame  : "MHRanForestGini";
    176 
    177     if (!matrixg->IsValid())
    178     {
    179         *fLog << err << dbginf << " MRanForestCalc::Grow - ERROR: matrixg not valid." << endl;
    180         return kFALSE;
    181     }
    182     if(!matrixh->IsValid())
    183     {
    184         *fLog << err << dbginf << " MRanForestCalc::Grow - ERROR: matrixh not valid." << endl;
    185         return kFALSE;
    186     }
    187 
    188     MEvtLoop run(GetName());
    189     MTaskList tlist;
    190     MParList plist;
    191     plist.AddToList(&tlist);
    192     plist.AddToList(matrixg);
    193     plist.AddToList(matrixh);
    194 
    195     // creating training task and setting parameters
    196     MRanForestGrow rfgrow;
    197     rfgrow.SetNumTrees(ntree); // number of trees
    198     rfgrow.SetNumTry(numtry);  // number of trials in random split selection
    199     rfgrow.SetNdSize(ndsize);  // limit for nodesize
    200     tlist.AddToList(&rfgrow);
    201 
    202     if(treefile){
    203         MWriteRootFile rfwrite(treefile);
    204         rfwrite.AddContainer(contname,treename);
    205         tlist.AddToList(&rfwrite);
    206     }
    207 
    208     MFillH fillh(hgininame);
    209     tlist.AddToList(&fillh);
    210 
    211     run.SetParList(&plist);
    212    
    213     // Execute tree growing
    214     if (!run.Eventloop())
    215     {
    216         *fLog << err << dbginf << "Evtloop in MRanForestCalc::Grow failed." << endl;
    217         return kFALSE;
    218     }
    219     tlist.PrintStatistics(0, kTRUE);
    220 
    221     if (TestBit(kEnableGraphicalOutput))
    222         plist.FindObject(hgininame)->DrawClone();
    223 
    224     return kTRUE;
    225 }
    226 
    227 Bool_t MRanForestCalc::Fill(Int_t ntree,const char* treefile,const char* treename)
    228 {
    229     treefile  = treefile  ? treefile  : "RF.root";
    230     treename  = treename  ? treename  : "Tree";
    231 
    232     MParList  plist;
    233 
    234     MTaskList tlist;
    235     plist.AddToList(&tlist);
    236 
    237     MReadTree read(treename,treefile);
    238     read.DisableAutoScheme();
    239 
    240     MRanForestFill rffill;
    241     rffill.SetNumTrees(ntree);
    242 
    243     tlist.AddToList(&read);
    244     tlist.AddToList(&rffill);
    245 
    246     MEvtLoop run(GetName());
    247     run.SetParList(&plist);
    248 
    249     //
    250     // Execute tree reading
    251     //
    252     if (!run.Eventloop())
    253     {
    254         *fLog << err << dbginf << "Evtloop in MRanForestCalc::Fill failed." << endl;
    255         return kFALSE;
    256     }
    257     tlist.PrintStatistics(0, kTRUE);
    258 
    259     return kTRUE;
    260 }
    261 
    262369Int_t MRanForestCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    263370{
    264     if (!IsEnvDefined(env, prefix, "File",     print))
    265         return kFALSE;
    266 
    267     TString fname = GetEnvValue(env, prefix, "File", "RF.root");
    268     TString tname = GetEnvValue(env, prefix, "Tree", "Tree");
    269     const Int_t   num   = GetEnvValue(env, prefix, "NumTrees", 100);
    270     const Bool_t  debug = GetEnvValue(env, prefix, "Debug", kFALSE);
    271 
    272     fname.ReplaceAll("\015", "");
    273     tname.ReplaceAll("\015", "");
    274 
    275     *fLog << inf << dbginf << "Reading " << num << " trees " << tname << " from file " << fname << endl;
    276 
    277     gLog.SetNullOutput(!debug);
    278     MEvtLoop evtloop;
    279     MParList  plist;
    280     evtloop.SetParList(&plist);
    281     MLog l;
    282     l.SetNullOutput(!debug);
    283     evtloop.SetLogStream(&l);
    284     gLog.SetNullOutput(debug);
    285 
    286     if (IsOwner())
    287     {
    288         delete fRanForest;
    289         delete fRanTree;
    290     }
    291     fRanForest = new MRanForest;
    292     fRanTree   = new MRanTree;
    293     SetOwner();
    294 
    295     plist.AddToList(fRanForest);
    296     plist.AddToList(fRanTree);
    297 
    298     MTaskList tlist;
    299     plist.AddToList(&tlist);
    300 
    301     MReadTree read(tname, fname);
    302     read.DisableAutoScheme();
    303 
    304     MRanForestFill rffill;
    305     rffill.SetNumTrees(num);
    306 
    307     tlist.AddToList(&read);
    308     tlist.AddToList(&rffill);
    309 
    310     if (!evtloop.Eventloop())
    311     {
    312         *fLog << err << "ERROR - Reading " << tname << " from file " << fname << endl;
    313         return kERROR;
    314     }
    315 
    316     return kTRUE;
    317 }
     371    Bool_t rc = kFALSE;
     372    if (IsEnvDefined(env, prefix, "FileName", print))
     373    {
     374        rc = kTRUE;
     375        SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
     376    }
     377    if (IsEnvDefined(env, prefix, "Debug", print))
     378    {
     379        rc = kTRUE;
     380        SetDebug(GetEnvValue(env, prefix, "Debug", fDebug));
     381    }
     382    if (IsEnvDefined(env, prefix, "NameOutput", print))
     383    {
     384        rc = kTRUE;
     385        SetNameOutput(GetEnvValue(env, prefix, "NameOutput", fNameOutput));
     386    }
     387    if (IsEnvDefined(env, prefix, "EstimationMode", print))
     388    {
     389        TString txt = GetEnvValue(env, prefix, "EstimationMode", "");
     390        txt = txt.Strip(TString::kBoth);
     391        txt.ToLower();
     392        if (txt==(TString)"mean")
     393            fEstimationMode = kMean;
     394        if (txt==(TString)"maximum")
     395            fEstimationMode = kMaximum;
     396        if (txt==(TString)"fit")
     397            fEstimationMode = kFit;
     398        rc = kTRUE;
     399    }
     400    return rc;
     401}
  • trunk/MagicSoft/Mars/mranforest/MRanForestCalc.h

    r6893 r7413  
    66#endif
    77
    8 #ifndef MARS_MHMatrix
    9 #include "MHMatrix.h"
     8#ifndef ROOT_TObjArray
     9#include <TObjArray.h>
    1010#endif
    1111
    12 class MParList;
     12#ifndef ROOT_TArrayD
     13#include <TArrayD.h>
     14#endif
     15
     16class MDataArray;
    1317class MParameterD;
    14 class MDataArray;
    15 class MRanTree;
    16 class MRanForest;
     18class MHMatrix;
    1719
    1820class MRanForestCalc : public MTask
    1921{
     22public:
     23    enum EstimationMode_t
     24    {
     25        kMean,
     26        kMaximum,
     27        kFit
     28    };
     29
    2030private:
    21     Int_t  fNum;              // number of trees used to compute hadronness
     31    static const TString gsDefName;     //! Default Name
     32    static const TString gsDefTitle;    //! Default Title
     33    static const TString gsNameOutput;  //! Default Output name
    2234
    23     TString fHadronnessName;  // Name of container storing hadronness
     35    Bool_t       fDebug;      // Debugging of eventloop while training on/off
    2436
    25     MParameterD *fHadroness;  //! Output container for calculated hadroness
     37    TString      fFileName;   // File name to forest
     38    TObjArray    fEForests;   // List of forests
     39
     40    TString      fNameOutput; // Name of output container
     41
    2642    MDataArray  *fData;       //! Used to store the MDataChains to get the event values
    27     MRanForest  *fRanForest;
    28     MRanTree    *fRanTree;
     43    MParameterD *fRFOut;      //! Used to store result
    2944
     45    Int_t        fNumTrees;   //! Training parameters
     46    Int_t        fNumTry;     //! Training parameters
     47    Int_t        fNdSize;     //! Training parameters
     48
     49    Int_t        fNumObsoleteVariables;
     50
     51    MHMatrix    *fTestMatrix; //! Test Matrix used in Process (together with MMatrixLoop)
     52
     53    EstimationMode_t fEstimationMode;
     54
     55private:
     56    // MTask
    3057    Int_t PreProcess(MParList *plist);
    3158    Int_t Process();
    3259
    33     enum { kIsOwner = BIT(14) };
     60    // MRanForestCalc
     61    Int_t ReadForests(MParList &plist);
    3462
    35     Bool_t IsOwner() const { return TestBit(kIsOwner); }
    36     void SetOwner() { SetBit(kIsOwner); }
     63    // MParContainer
     64    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
     65
     66    // Train Interface
     67    Int_t Train(const MHMatrix &n, const TArrayD &grid, Int_t ver=2);
    3768
    3869public:
     
    4071    ~MRanForestCalc();
    4172
    42     Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
     73    // Setter for estimation
     74    void SetFileName(TString filename)            { fFileName = filename; }
     75    void SetEstimationMode(EstimationMode_t op)   { fEstimationMode = op; }
     76    void SetNameOutput(TString name=gsNameOutput) { fNameOutput = name; }
    4377
    44     void SetHadronnessName(const TString name) { fHadronnessName = name; }
    45     TString GetHadronnessName() const { return fHadronnessName; }
     78    // Setter for training
     79    void SetNumTrees(UShort_t n=100) { fNumTrees = n; }
     80    void SetNdSize(UShort_t n=5)     { fNdSize   = n; }
     81    void SetNumTry(UShort_t n=0)     { fNumTry   = n; }
     82    void SetDebug(Bool_t b=kTRUE)    { fDebug    = b; }
    4683
    47     void SetUseNumTrees(UShort_t n=100) { fNum = n; }
     84    void SetNumObsoleteVariables(Int_t n=1) { fNumObsoleteVariables = n; }
    4885
    49     Bool_t Grow(MHMatrix *matrixg,MHMatrix *matrixh,Int_t ntree,
    50                 Int_t numtry,Int_t ndsize,const char* treefile=0,
    51                 const char* treename=0,const char* contname=0,
    52                 const char* hgininame=0);
     86    // Train Interface
     87    Int_t TrainMultiRF(const MHMatrix &n, const TArrayD &grid)
     88    {
     89        return Train(n, grid, 0);
     90    }
     91    Int_t TrainSingleRF(const MHMatrix &n, const TArrayD &grid=TArrayD())
     92    {
     93        return Train(n, grid, grid.GetSize()==0 ? 2 : 1);
     94    }
    5395
    54     Bool_t Fill(Int_t ntree,const char* treefile=0,const char* treename=0);
     96    // Test Interface
     97    void  SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; }
     98    void  InitMapping(MHMatrix *m=0)   { fTestMatrix=m; }
    5599
    56 
    57     ClassDef(MRanForestCalc, 0) // Task
     100    ClassDef(MRanForestCalc, 0) // Task to calculate RF output and for RF training
    58101};
    59102
  • trunk/MagicSoft/Mars/mranforest/MRanTree.h

    r7396 r7413  
    2828    Int_t fNumNodes;
    2929    Int_t fNumEndNodes;
     30    Float_t fError;
    3031
    3132    TArrayI fBestVar;
     
    7071    void SetNdSize(Int_t n);
    7172    void SetNumTry(Int_t n);
     73    void SetError(Float_t f) { fError = f; }
    7274
    7375    Int_t GetNdSize() const { return fNdSize; }
     
    7577    Int_t GetNumNodes()          const { return fNumNodes; }
    7678    Int_t GetNumEndNodes()       const { return fNumEndNodes; }
     79    Float_t GetError() const { return fError; }
    7780
    7881    Int_t GetBestVar(Int_t i)    const { return fBestVar.At(i); }
  • trunk/MagicSoft/Mars/mranforest/Makefile

    r7396 r7413  
    2525           MRanForest.cc \
    2626           MRanForestGrow.cc \
     27           MRanForestCalc.cc \
    2728           MHRanForest.cc \
    28            MHRanForestGini.cc \
    29            MRFEnergyEst.cc
     29           MHRanForestGini.cc
    3030
    3131############################################################
  • trunk/MagicSoft/Mars/mranforest/RanForestLinkDef.h

    r7396 r7413  
    88#pragma link C++ class MRanForest+;
    99#pragma link C++ class MRanForestGrow+;
     10#pragma link C++ class MRanForestCalc+;
    1011
    1112#pragma link C++ class MHRanForest+;
    1213#pragma link C++ class MHRanForestGini+;
    1314
    14 #pragma link C++ class MRFEnergyEst+;
    15 
    1615#endif
  • trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc

    r7401 r7413  
    204204}
    205205
     206//------------------------------------------------------------------------
     207//
     208// Add all entries deriving from MFilter from list to PreCuts.
     209// The ownership is not affected.
     210//
     211void MTFillMatrix::AddPreCuts(const TList &list)
     212{
     213    TIter Next(&list);
     214    TObject *obj=0;
     215    while ((obj=Next()))
     216        if (obj->InheritsFrom(MFilter::Class()))
     217            fPreCuts.Add(obj);
     218}
    206219
    207220// --------------------------------------------------------------------------
     
    209222// Fill the matrix (FIXME: Flow diagram missing)
    210223//
    211 Bool_t MTFillMatrix::Process()
     224Bool_t MTFillMatrix::Process(const MParList &parlist)
    212225{
    213226    if (!fReader)
     
    232245    // Create parameter list and task list, add tasklist to parlist
    233246    //
    234     MParList  plist;
     247    parlist.Print();
     248    MParList  plist(parlist);
    235249    MTaskList tlist;
    236250    plist.AddToList(&tlist);
     
    372386    return WriteMatrix1(fname) && WriteMatrix2(fname);
    373387}
     388
     389Int_t MTFillMatrix::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
     390{
     391    Bool_t rc = kFALSE;
     392    if (IsEnvDefined(env, prefix, "NumDestEvents1", print))
     393    {
     394        rc = kTRUE;
     395        SetNumDestEvents1(GetEnvValue(env, prefix, "NumDestEvents1", fNumDestEvents1));
     396    }
     397    if (IsEnvDefined(env, prefix, "NumDestEvents2", print))
     398    {
     399        rc = kTRUE;
     400        SetNumDestEvents2(GetEnvValue(env, prefix, "NumDestEvents2", fNumDestEvents2));
     401    }
     402    return rc;
     403}
  • trunk/MagicSoft/Mars/mtools/MTFillMatrix.h

    r7402 r7413  
    11#ifndef MARS_MTFillMatrix
    22#define MARS_MTFillMatrix
     3
     4#ifndef MARS_MParList
     5#include "MParList.h"
     6#endif
    37
    48#ifndef MARS_MH3
     
    2731
    2832    TList     fPreCuts;
     33
     34    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    2935
    3036    Bool_t CheckResult(MHMatrix *m, Int_t num) const;
     
    7278    void AddPreCut(const char *rule);
    7379    void AddPreCut(MFilter *f);
     80    void AddPreCuts(const TList &list);
    7481
    75     Bool_t Process();
     82    Bool_t Process(const MParList &plist=MParList());
    7683    Bool_t WriteMatrix1(const TString &fname) const { return WriteMatrix(fDestMatrix1, fname, 1); }
    7784    Bool_t WriteMatrix2(const TString &fname) const { return WriteMatrix(fDestMatrix2, fname, 2); }
Note: See TracChangeset for help on using the changeset viewer.