Changeset 2296 for trunk


Ignore:
Timestamp:
07/29/03 13:18:57 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 added
28 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2295 r2296  
    111111      - add member function GetSignificance()
    112112
    113    * mhist/MHMatrix.cc
    114      - add MProgressBar in Fill()
     113    * mhist/MHMatrix.cc
     114      - add MProgressBar in Fill()
    115115
    116116    * mmontecarlo/MMcEnergyEst.h

  • trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc

    r2274 r2296  
    172172    for (UShort_t i=0; i<entries; i++)
    173173    {
     174        // FIXME: Loop over blinds instead of all pixels!!!
    174175        MCerPhotPix &pix = (*fEvt)[i];
    175176
    … …  
    203204                perr[i]  += evtpix->GetErrorPhot();
    204205                // FIXME: perr[i] ???
     206                num++;
    205207            }
    206             num++;
    207208        }
    208209

  • trunk/MagicSoft/Mars/mbase/MArgs.cc

    r2275 r2296  
    3838ClassImp(MArgsEntry);
    3939ClassImp(MArgs);
     40
     41using namespace std;
    4042
    4143void MArgsEntry::Print(const Option_t *o) const

  • trunk/MagicSoft/Mars/mbase/MDirIter.h

    r2196 r2296  
    3838            AddDirectory(o->GetName(), o->GetTitle());
    3939    }
    40     MDirIter(const char *dir, const char *filter="") : fNext(&fList), fDirPtr(NULL)
     40    MDirIter(const char *dir, const char *filter="", Int_t rec=0) : fNext(&fList), fDirPtr(NULL)
    4141    {
    4242        fList.SetOwner();
    43         AddDirectory(dir, filter);
     43        AddDirectory(dir, filter, rec);
    4444    }
    4545    ~MDirIter()

  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

    r2247 r2296  
    125125// Add this file as the last entry in the chain
    126126//
    127 void MCT1ReadPreProc::AddFile(const char *txt)
     127Int_t MCT1ReadPreProc::AddFile(const char *txt, int)
    128128{
    129129    const char *name = gSystem->ExpandPathName(txt);
    … …  
    135135    {
    136136        *fLog << warn << "WARNING - Problem reading header... ignored." << endl;
    137         return;
     137        return 0;
    138138    }
    139139
    … …  
    142142    {
    143143        *fLog << warn << "WARNING - File contains no data... ignored." << endl;
    144         return;
     144        return 0;
    145145    }
    146146
    … …  
    150150
    151151    fFileNames->AddLast(new TNamed(txt, ""));
     152    return 1;
    152153}
    153154

  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

    r2206 r2296  
    9292    ~MCT1ReadPreProc();
    9393
    94     void AddFile(const char *fname);
     94    Int_t AddFile(const char *fname, int i=0);
    9595
    9696    UInt_t GetEntries() { return fEntries; }

  • trunk/MagicSoft/Mars/mfileio/MRead.cc

    r2173 r2296  
    3939#include "MLog.h"
    4040#include "MLogManip.h"
     41
     42#include "MDirIter.h"
    4143
    4244ClassImp(MRead);
    … …  
    98100    return i!=0;
    99101}
     102
     103Int_t MRead::AddFiles(MDirIter &next)
     104{
     105    Int_t rc=0;
     106    TString n;
     107    while (!(n=next()).IsNull())
     108    {
     109        const Int_t num = AddFile(n);
     110        if (num>0)
     111            rc += num;
     112    }
     113    return rc;
     114}
     115

  • trunk/MagicSoft/Mars/mfileio/MRead.h

    r2117 r2296  
    77
    88class MFilter;
     9class MDirIter;
    910
    1011class MRead : public MTask
    … …  
    2223    MFilter *GetSelector() { return fSelector; }
    2324
    24     Int_t AddFile(const char *fname, Int_t entries=-1) { return -1; }
     25    virtual Int_t AddFile(const char *fname, Int_t entries=-1) { return -1; }
     26    Int_t AddFiles(MDirIter &next);
    2527
    2628    Bool_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);

  • trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc

    r2206 r2296  
    399399// Add this file as the last entry in the chain
    400400//
    401 void MReadRflFile::AddFile(const char *txt)
     401Int_t MReadRflFile::AddFile(const char *txt, int)
    402402{
    403403    const char *name = gSystem->ExpandPathName(txt);
    … …  
    424424*/
    425425    fFileNames->AddLast(new TNamed(txt, ""));
     426    return 1;
    426427}
    427428

  • trunk/MagicSoft/Mars/mfileio/MReadRflFile.h

    r2206 r2296  
    4848    ~MReadRflFile();
    4949
    50     void AddFile(const char *fname);
     50    Int_t AddFile(const char *fname, int i=0);
    5151
    5252    Bool_t Rewind() { fNumFile=0; return kTRUE; }

  • trunk/MagicSoft/Mars/mfileio/MReadTree.cc

    r2206 r2296  
    121121
    122122    if (fname)
    123         if (fChain->Add(fname)>0)
    124             SetBit(kChainWasChanged);
     123        AddFile(fname);
     124//        if (fChain->Add(fname)>0)
     125//            SetBit(kChainWasChanged);
    125126}
    126127
    … …  
    586587    if (!GetEntries())
    587588    {
    588         *fLog << warn << dbginf << "No entries found in file(s)" << endl;
     589        *fLog << warn << GetDescriptor() << ": No entries found in file(s)" << endl;
    589590        return kFALSE;
    590591    }

  • trunk/MagicSoft/Mars/mhist/MH.cc

    r2209 r2296  
    796796// Otherwise the present gPad is returned.
    797797//
    798 TVirtualPad *MH::GetNewPad(Option_t *opt)
    799 {
    800     TString str(opt);
    801 
    802     if (!str.Contains("nonew", TString::kIgnoreCase))
     798TVirtualPad *MH::GetNewPad(TString &opt)
     799{
     800    opt.ToLower();
     801
     802    if (!opt.Contains("nonew"))
    803803        return NULL;
     804
     805    opt.ReplaceAll("nonew", "");
    804806
    805807    return gPad;
    … …  
    813815TObject *MH::Clone(const char *name) const
    814816{
    815     Bool_t store = TH1::AddDirectoryStatus();
     817    const Bool_t store = TH1::AddDirectoryStatus();
     818
    816819    TH1::AddDirectory(kFALSE);
    817 
    818820    TObject *o = MParContainer::Clone(name);
    819 
    820821    TH1::AddDirectory(store);
    821822
    … …  
    831832TObject *MH::DrawClone(Option_t *opt, Int_t w, Int_t h) const
    832833{
    833     TVirtualPad *p = GetNewPad(opt);
     834    TString option(opt);
     835
     836    TVirtualPad *p = GetNewPad(option);
    834837    if (!p)
    835838        p = MakeDefCanvas(this, w, h);
    … …  
    839842    gROOT->SetSelectedPad(NULL);
    840843
    841     TObject *o = MParContainer::DrawClone(opt);
     844    TObject *o = MParContainer::DrawClone(option);
    842845    o->SetBit(kCanDelete);
    843846    return o;

  • trunk/MagicSoft/Mars/mhist/MH.h

    r2150 r2296  
    7575    }
    7676
    77     static TVirtualPad *GetNewPad(Option_t *opt);
     77    static TVirtualPad *GetNewPad(TString &opt);
    7878
    7979    static void FindGoodLimits(Int_t nbins, Int_t &newbins, Double_t &xmin, Double_t &xmax, Bool_t isInteger);

  • trunk/MagicSoft/Mars/mhist/MHHadronness.cc

    r2225 r2296  
    430430// --------------------------------------------------------------------------
    431431//
     432// Search the hadronness corresponding to a given hadron acceptance.
     433//
     434Double_t MHHadronness::GetHadronness(Double_t acchad) const
     435{
     436    for (int i=1; i<fIntPhness->GetNbinsX()+1; i++)
     437        if (fIntPhness->GetBinContent(i)>acchad)
     438            return fIntPhness->GetBinLowEdge(i);
     439
     440    return -1;
     441}
     442
     443// --------------------------------------------------------------------------
     444//
    432445// Print the corresponding Gammas Acceptance for a hadron acceptance of
    433446// 10%, 20%, 30%, 40% and 50%. Also the minimum distance to the optimum
    … …  
    448461
    449462    *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
    450     *fLog << "acc(hadron)  acc(gamma)  acc(g)/acc(h)" << endl <<endl;
    451 
    452     *fLog << "    0.005    " << Form("%6.3f", GetGammaAcceptance(0.005)) << "    " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
    453     *fLog << "    0.02     " << Form("%6.3f", GetGammaAcceptance(0.02))  << "    " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
    454     *fLog << "    0.05     " << Form("%6.3f", GetGammaAcceptance(0.05))  << "    " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
    455     *fLog << "    0.1      " << Form("%6.3f", GetGammaAcceptance(0.1 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
    456     *fLog << "    0.2      " << Form("%6.3f", GetGammaAcceptance(0.2 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
    457     *fLog << "    0.3      " << Form("%6.3f", GetGammaAcceptance(0.3 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
     463    *fLog << "acc(hadron) acc(gamma) acc(g)/acc(h)  h" << endl <<endl;
     464
     465    *fLog << "    0.005    " << Form("%6.3f", GetGammaAcceptance(0.005)) << "      " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << "      " << GetHadronness(0.005) << endl;
     466    *fLog << "    0.02     " << Form("%6.3f", GetGammaAcceptance(0.02))  << "      " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02)   << "      " << GetHadronness(0.02) << endl;
     467    *fLog << "    0.05     " << Form("%6.3f", GetGammaAcceptance(0.05))  << "      " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05)   << "      " << GetHadronness(0.05) << endl;
     468    *fLog << "    0.1      " << Form("%6.3f", GetGammaAcceptance(0.1 ))  << "      " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1)     << "      " << GetHadronness(0.1) << endl;
     469    *fLog << "    0.2      " << Form("%6.3f", GetGammaAcceptance(0.2 ))  << "      " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2)     << "      " << GetHadronness(0.2) << endl;
     470    *fLog << "    0.3      " << Form("%6.3f", GetGammaAcceptance(0.3 ))  << "      " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3)     << "      " << GetHadronness(0.3) << endl;
    458471    *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << "        0.1  " << endl;
    459472    *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << "        0.2  " << endl;

  • trunk/MagicSoft/Mars/mhist/MHHadronness.h

    r2225 r2296  
    3838    Double_t GetGammaAcceptance(Double_t acchad) const;
    3939    Double_t GetHadronAcceptance(Double_t accgam) const;
     40    Double_t GetHadronness(Double_t acchad) const;
    4041
    4142    TH1D *Getghness() const  { return fGhness; }

  • trunk/MagicSoft/Mars/mhist/MHMatrix.cc

    r2282 r2296  
    260260    return kTRUE;
    261261}
     262
    262263/*
    263264// --------------------------------------------------------------------------
    … …  
    503504    if (!fM2.IsValid())
    504505    {
     506        if (!fM.IsValid())
     507        {
     508            *fLog << err << "MHMatrix::CalcDist - ERROR: fM not valid." << endl;
     509            return -1;
     510        }
     511
    505512        const TMatrix *m = InvertPosDef();
    506513        if (!m)
    … …  
    625632    tlist.AddToList(&fillh);
    626633
    627     MProgressBar bar;
     634    //MProgressBar bar;
    628635    MEvtLoop evtloop;
    629636    evtloop.SetParList(plist);
    630     evtloop.SetProgressBar(&bar);
     637    //evtloop.SetProgressBar(&bar);
    631638
    632639    if (!evtloop.Eventloop())
    633640        return kFALSE;
    634641
    635     tlist.PrintStatistics(0, kTRUE);
     642    tlist.PrintStatistics();
    636643
    637644    plist->Remove(&tlist);
    … …  
    10971104void MHMatrix::ShuffleRows(UInt_t seed)
    10981105{
    1099   TRandom rnd(seed);
    1100 
    1101   for (Int_t irow = 0; irow < fNumRow; irow++)
    1102     {
    1103       Int_t jrow = rnd.Integer(fNumRow);
    1104 
    1105       for (Int_t icol = 0; icol < fM.GetNcols(); icol++)
    1106         {
    1107           Real_t tmp = fM(irow,icol);
    1108           fM(irow,icol) = fM(jrow,icol);
    1109           fM(jrow,icol) = tmp;
    1110         }
    1111     }
    1112 
    1113   *fLog << warn << this->GetName() << " : Attention! Matrix rows have been shuffled." << endl;
    1114 
    1115 }
     1106    TRandom rnd(seed);
     1107
     1108    TVector v(fM.GetNcols());
     1109    TVector tmp(fM.GetNcols());
     1110    for (Int_t irow = 0; irow<fNumRow; irow++)
     1111    {
     1112        const Int_t jrow = rnd.Integer(fNumRow);
     1113
     1114        v = TMatrixRow(fM, irow);
     1115        TMatrixRow(fM, irow) = tmp = TMatrixRow(fM, jrow);
     1116        TMatrixRow(fM, jrow) = v;
     1117    }
     1118
     1119    *fLog << warn << GetDescriptor() << ": Attention! Matrix rows have been shuffled." << endl;
     1120}
     1121
     1122void MHMatrix::ReduceRows(UInt_t num)
     1123{
     1124    if ((Int_t)num>=fM.GetNrows())
     1125    {
     1126        *fLog << warn << GetDescriptor() << ": Warning - " << num << " >= rows=" << fM.GetNrows() << endl;
     1127        return;
     1128    }
     1129
     1130    const TMatrix m(fM);
     1131    fM.ResizeTo(num, fM.GetNcols());
     1132
     1133    TVector tmp(fM.GetNcols());
     1134    for (UInt_t irow=0; irow<num; irow++)
     1135        TMatrixRow(fM, irow) = tmp = TMatrixRow(m, irow);
     1136}

  • trunk/MagicSoft/Mars/mhist/MHMatrix.h

    r2133 r2296  
    113113
    114114    void ShuffleRows(UInt_t seed);
     115    void ReduceRows(UInt_t num);
    115116
    116117    ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events

  • trunk/MagicSoft/Mars/mimage/ImageLinkDef.h

    r2026 r2296  
    66
    77#pragma link C++ class MImgCleanStd+;
     8#pragma link C++ class MImgCleanTGB+;
    89#pragma link C++ class MCameraSmooth+;
    910

  • trunk/MagicSoft/Mars/mimage/MCameraSmooth.h

    r2209 r2296  
    2929    MCameraSmooth(Byte_t cnt=1, const char *name=NULL, const char *title=NULL);
    3030
    31     void SetUseCetralPixel(Bool_t b=kTRUE) { fUseCentralPixel=kTRUE; }
     31    void SetUseCentralPixel(Bool_t b=kTRUE) { fUseCentralPixel=kTRUE; }
    3232
    3333    ClassDef(MCameraSmooth, 0) // task to smooth the camera contants

  • trunk/MagicSoft/Mars/mimage/Makefile

    r2179 r2296  
    2929
    3030SRCFILES = MImgCleanStd.cc \
     31           MImgCleanTGB.cc \
    3132           MCameraSmooth.cc \
    3233           MHillas.cc \

  • trunk/MagicSoft/Mars/mmain/MStatusDisplay.cc

    r2274 r2296  
    840840Bool_t MStatusDisplay::HasCanvas(const TCanvas *c) const
    841841{
     842    if (!c)
     843        return kFALSE;
     844
    842845    if (gROOT->IsBatch())
    843846        return (Bool_t)fBatch->FindObject(c);

  • trunk/MagicSoft/Mars/mranforest/MHRanForest.cc

    r2173 r2296  
    6868    fGraphSigma = new TGraph;
    6969    fGraphSigma->SetTitle("Evolution of Standard deviation of estimated hadronness in tree combination");
    70     fGraphSigma->SetMaximum(1);
    7170    fGraphSigma->SetMarkerStyle(kFullDotSmall);
    7271}
    … …  
    113112{
    114113    fNumEvent++;
     114
    115115    Double_t hest=0;
    116116    Double_t htrue=fMcEvt->GetPartId()==kGAMMA ? 0. : 1.;
    … …  
    124124
    125125        hest/=i+1;
    126         fSigma[i]+=(htrue-hest)*(htrue-hest);
     126
     127        const Double_t val = htrue-hest;
     128        fSigma[i] += val*val;
    127129    }
    128130
    … …  
    145147    for (Int_t i=0; i<n; i++)
    146148    {
    147         Stat_t ip = i+1.;
    148         fSigma[i] = TMath::Sqrt(fSigma[i]/Stat_t(fNumEvent));
    149         Stat_t ig = fSigma[i];
    150         max=TMath::Max(max,ig);
    151         min=TMath::Min(min,ig);
    152         fGraphSigma->SetPoint(i,ip,ig);
     149        fSigma[i] = TMath::Sqrt(fSigma[i]/fNumEvent);
     150
     151        const Stat_t ig = fSigma[i];
     152        if (ig>max) max = ig;
     153        if (ig<min) min = ig;
     154        fGraphSigma->SetPoint(i, i+1, ig);
    153155    }
     156
     157    // This is used in root>3.04/? so that SetMaximum/Minimum can succeed
     158    fGraphSigma->GetHistogram();
     159
    154160    fGraphSigma->SetMaximum(1.05*max);
    155161    fGraphSigma->SetMinimum(0.95*min);
    … …  
    180186        return;
    181187
    182     h->GetXaxis()->SetRangeUser(0, 1);
     188    //h->GetXaxis()->SetRangeUser(0, fNumEvent+1);
    183189    h->SetXTitle("No.of Trees");
    184190    h->SetYTitle("\\sigma of est.hadronness");

  • trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc

    r2173 r2296  
    6666    fGraphGini = new TGraph;
    6767    fGraphGini->SetTitle("Importance of RF-input parameters measured by mean Gini decrease");
    68     fGraphGini->SetMaximum(1);
    6968    fGraphGini->SetMarkerStyle(kFullDotSmall);
    7069}
    … …  
    116115Bool_t MHRanForestGini::Finalize()
    117116{
    118     Int_t n = fGini.GetSize();
     117    const Int_t n = fGini.GetSize();
    119118
    120119    fGraphGini->Set(n);
    121120
    122     Stat_t max=0.;
    123     Stat_t min=0.;
     121    Stat_t max=0;
     122    Stat_t min=0;
    124123    for (Int_t i=0; i<n; i++)
    125124    {
    126         fGini[i]/=fRanForest->GetNumTrees();
    127         fGini[i]/=fRanForest->GetNumData();
     125        fGini[i] /= fRanForest->GetNumTrees();
     126        fGini[i] /= fRanForest->GetNumData();
    128127
    129         Stat_t ip = i+1.;
    130         Stat_t ig = fGini[i];
     128        const Stat_t ip = i+1;
     129        const Stat_t ig = fGini[i];
    131130
    132         max=TMath::Max(max,ig);
    133         min=TMath::Min(min,ig);
     131        if (ig>max) max=ig;
     132        if (ig<min) min=ig;
    134133
    135134        fGraphGini->SetPoint(i,ip,ig);
    136135    }
     136
     137    // This is used in root>3.04/? so that SetMaximum/Minimum can succeed
     138    fGraphGini->GetHistogram();
     139
    137140    fGraphGini->SetMaximum(1.05*max);
    138141    fGraphGini->SetMinimum(0.95*min);
    … …  
    163166        return;
    164167
    165     h->GetXaxis()->SetRangeUser(0, 1);
     168    //h->GetXaxis()->SetRangeUser(0, fGini.GetSize()+1);
    166169    h->SetXTitle("No.of RF-input parameter");
    167170    h->SetYTitle("Mean decrease in Gini-index [a.u.]");

  • trunk/MagicSoft/Mars/mranforest/MRanForest.cc

    r2173 r2296  
    103103    Int_t ntree=0;
    104104
    105     TIter forest(fForest);
     105    TIter Next(fForest);
    106106
    107107    MRanTree *tree;
    108     while ((tree=(MRanTree*)forest.Next()))
     108    while ((tree=(MRanTree*)Next()))
    109109    {
    110110        fTreeHad[ntree]=tree->TreeHad(event);
    … …  
    112112        ntree++;
    113113    }
    114     return hadroness/Double_t(ntree);
     114    return hadroness/ntree;
    115115}
    116116
    … …  
    190190Bool_t MRanForest::GrowForest()
    191191{
    192     Int_t ninbag=0;
    193     TArrayI datsortinbag;
    194     TArrayF classpopw;
    195     TArrayI jinbag;
    196     TArrayF winbag;
    197 
    198     jinbag.Set(fNumData);
    199     winbag.Set(fNumData);
    200     classpopw.Set(2);
    201 
    202     TMatrix hadrons=fHadrons->GetM();
    203     TMatrix gammas=fGammas->GetM();
     192    if(!gRandom)
     193    {
     194        *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
     195        return kFALSE;
     196    }
    204197
    205198    fTreeNo++;
    206199
    207200    // initialize running output
    208     if(fTreeNo==1)
    209     {
    210         cout<<endl<<endl<<"1st col.: no. of tree"<<endl;
    211         cout<<"2nd col.: error in % (calulated using oob-data -> overestim. of error)"<<endl;
    212     }
    213 
    214     jinbag.Reset();
    215     classpopw.Reset();
    216     winbag.Reset();
     201    if (fTreeNo==1)
     202    {
     203        *fLog << inf << endl;
     204        *fLog << underline; // << "1st col        2nd col" << endl;
     205        *fLog << "no. of tree    error in % (calulated using oob-data -> overestim. of error)" << endl;
     206    }
     207
     208    TArrayF classpopw(2);
     209    TArrayI jinbag(fNumData); // Initialization includes filling with 0
     210    TArrayF winbag(fNumData); // Initialization includes filling with 0
    217211
    218212    // bootstrap aggregating (bagging) -> sampling with replacement:
    … …  
    221215    // {0,1,...,fNumData-1}, which is the set of the index numbers of
    222216    // all events in the training sample
    223 
    224     for (Int_t n=0;n<fNumData;n++)
    225     {
    226         if(!gRandom)
    227         {
    228             *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
    229             return kFALSE;
    230         }
    231 
    232         Int_t k=Int_t(fNumData*gRandom->Rndm());
     217    for (Int_t n=0; n<fNumData; n++)
     218    {
     219        const Int_t k = Int_t(gRandom->Rndm()*fNumData);
    233220
    234221        classpopw[fHadTrue[k]]+=fWeight[k];
    … …  
    241228    // In bagging procedure ca. 2/3 of all elements in the original
    242229    // training sample are used to build the in-bag data
    243     datsortinbag=fDataSort;
    244 
    245     ModifyDataSort(datsortinbag,ninbag,jinbag);
     230    TArrayI datsortinbag=fDataSort;
     231    Int_t ninbag=0;
     232
     233    ModifyDataSort(datsortinbag, ninbag, jinbag);
     234
     235    const TMatrix &hadrons=fHadrons->GetM();
     236    const TMatrix &gammas =fGammas->GetM();
    246237
    247238    // growing single tree
    … …  
    258249    // determined from oob-data is underestimated, but can still be taken as upper limit.
    259250
    260     TVector event(fNumDim);
    261     for(Int_t ievt=0;ievt<fNumHad;ievt++)
    262     {
    263         if(jinbag[ievt]>0)continue;
    264         event=TMatrixRow(hadrons,ievt);
    265         fHadEst[ievt]+=fRanTree->TreeHad(event);
     251    for (Int_t ievt=0;ievt<fNumHad;ievt++)
     252    {
     253        if (jinbag[ievt]>0)
     254            continue;
     255        fHadEst[ievt] += fRanTree->TreeHad(hadrons, ievt);
    266256        fNTimesOutBag[ievt]++;
    267257    }
    268     for(Int_t ievt=0;ievt<fNumGam;ievt++)
    269     {
    270         if(jinbag[fNumHad+ievt]>0)continue;
    271         event=TMatrixRow(gammas,ievt);
    272         fHadEst[fNumHad+ievt]+=fRanTree->TreeHad(event);
     258    for (Int_t ievt=0;ievt<fNumGam;ievt++)
     259    {
     260        if (jinbag[fNumHad+ievt]>0)
     261            continue;
     262        fHadEst[fNumHad+ievt] += fRanTree->TreeHad(gammas, ievt);
    273263        fNTimesOutBag[fNumHad+ievt]++;
    274264    }
    275265
    276266    Int_t n=0;
    277     fErr=0;
    278     for(Int_t ievt=0;ievt<fNumData;ievt++)
    279         if(fNTimesOutBag[ievt]!=0)
     267    Double_t ferr=0;
     268    for (Int_t ievt=0;ievt<fNumData;ievt++)
     269        if (fNTimesOutBag[ievt]!=0)
    280270        {
    281             fErr+=TMath::Power(fHadEst[ievt]/fNTimesOutBag[ievt]-fHadTrue[ievt],2.);
     271            const Double_t val = fHadEst[ievt]/fNTimesOutBag[ievt]-fHadTrue[ievt];
     272            ferr += val*val;
    282273            n++;
    283274        }
    284275
    285     fErr/=Float_t(n);
    286     fErr=TMath::Sqrt(fErr);
     276    ferr = TMath::Sqrt(ferr/n);
    287277
    288278    // give running output
    289     cout << setw(5) << fTreeNo << setw(15) << Form("%.2f",100.*fErr) << endl;
     279    *fLog << inf << setw(5) << fTreeNo << Form("%15.2f", ferr*100) << endl;
    290280
    291281    // adding tree to forest
    292282    AddTree();
    293283
    294     return(fTreeNo<fNumTrees);
     284    return fTreeNo<fNumTrees;
    295285}
    296286
    … …  
    315305    TArrayI isort(fNumData);
    316306
    317     TMatrix hadrons=fHadrons->GetM();
    318     TMatrix gammas=fGammas->GetM();
     307    const TMatrix &hadrons=fHadrons->GetM();
     308    const TMatrix &gammas=fGammas->GetM();
    319309
    320310    for (Int_t j=0;j<fNumHad;j++)
    … …  
    346336        for(Int_t n=0;n<fNumData-1;n++)
    347337        {
    348             Int_t n1=isort[n];
    349             Int_t n2=isort[n+1];
     338            const Int_t n1=isort[n];
     339            const Int_t n2=isort[n+1];
     340
    350341            fDataSort[mvar*fNumData+n]=n1;
    351342            if(n==0) fDataRang[mvar*fNumData+n1]=0;
    … …  
    359350        fDataSort[(mvar+1)*fNumData-1]=isort[fNumData-1];
    360351    }
    361 
    362     return;
    363 }
    364 
    365 void MRanForest::ModifyDataSort(TArrayI &datsortinbag,Int_t ninbag,TArrayI &jinbag)
     352}
     353
     354void MRanForest::ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag)
    366355{
    367356    ninbag=0;
    … …  
    394383        }
    395384    }
    396     return;
    397385}
    398386

  • trunk/MagicSoft/Mars/mranforest/MRanForest.h

    r2114 r2296  
    6565    // estimates for classification error of growing forest
    6666    TArrayD fTreeHad;
    67     Float_t fErr;
    6867
    6968protected:
    7069    // create and modify (->due to bagging) fDataSort
    7170    void CreateDataSort();
    72     void ModifyDataSort(TArrayI &datsortinbag,Int_t ninbag,TArrayI &jinbag);
     71    void ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag);
    7372
    7473public:

  • trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc

    r2207 r2296  
    8484// Needs:
    8585//  - MatrixGammas  [MHMatrix]
    86 //  - MatrixHadrons {MHMatrix]
     86//  - MatrixHadrons [MHMatrix]
    8787//  - MHadronness
    8888//  - all data containers used to build the matrixes

  • trunk/MagicSoft/Mars/mranforest/MRanTree.cc

    r2173 r2296  
    7676}
    7777
    78 void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,
     78void MRanTree::GrowTree(const TMatrix &mhad, const TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,
    7979                        TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag,
    8080                        TArrayF &winbag,TArrayF &weight)
    … …  
    8888    for (Int_t n=0;n<numdata;n++)
    8989        if(jinbag[n]==1) ninbag++;
    90 
    91     // weighted class populations after split
    92     TArrayF wl(2); // left node
    93     TArrayF wc(2);
    94     TArrayF wr(2); // right node
    95     TArrayI nc(2);
    9690
    9791    TArrayI bestsplit(nrnodes);
    … …  
    106100    TArrayI idmove(numdata);
    107101
    108     idmove.Reset();
    109 
    110102    fBestVar.Set(nrnodes);
    111103    fTreeMap1.Set(nrnodes);
    … …  
    123115    BuildTree(datasort,datarang,hadtrue,numdim,numdata,bestsplit,
    124116              bestsplitnext,nodepop,nodestart,tclasspop,nrnodes,
    125               idmove,ncase,parent,jinbag,iv,winbag,wr,wc,wl,ninbag);
     117              idmove,ncase,parent,jinbag,iv,winbag,ninbag);
    126118
    127119    // post processing, determine cut (or split) values fBestSplit
    128120    Int_t nhad=mhad.GetNrows();
    129121
    130     for(Int_t k=0;k<nrnodes;k++)
    131     {
    132         Int_t bsp=bestsplit[k];
    133         Int_t bspn=bestsplitnext[k];
    134         Int_t msp=fBestVar[k];
    135 
    136         if (GetNodeStatus(k)!=-1)
    137         {
    138             fBestSplit[k] = bsp<nhad ? mhad(bsp,msp):mgam(bsp-nhad,msp);
    139             fBestSplit[k] += bspn<nhad ? mhad(bspn,msp):mgam(bspn-nhad,msp);
    140             fBestSplit[k] /=2.;
    141         }
     122    for(Int_t k=0; k<nrnodes; k++)
     123    {
     124        if (GetNodeStatus(k)==-1)
     125            continue;
     126
     127        const Int_t &bsp =bestsplit[k];
     128        const Int_t &bspn=bestsplitnext[k];
     129        const Int_t &msp =fBestVar[k];
     130
     131        fBestSplit[k]  = bsp<nhad  ? mhad(bsp, msp):mgam(bsp-nhad, msp);
     132        fBestSplit[k] += bspn<nhad ? mhad(bspn,msp):mgam(bspn-nhad,msp);
     133        fBestSplit[k] /= 2;
    142134    }
    143135
    … …  
    152144                             Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop,
    153145                             Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase,
    154                              TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,
    155                              TArrayF &wc,TArrayF &wl,Int_t kbuild)
    156 {
     146                             TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t kbuild)
     147{
     148    if(!gRandom)
     149    {
     150        *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
     151        return kFALSE;
     152    }
     153
     154    // weighted class populations after split
     155    TArrayF wc(2);
     156    TArrayF wr(2); // right node
     157
    157158    // For the best split, msplit is the index of the variable (e.g Hillas par., zenith angle ,...)
    158159    // split on. decsplit is the decreae in impurity measured by Gini-index.
    … …  
    160161    // and nsplitnext is the case number of the next larger value of msplit.
    161162
    162     Int_t mvar,nc,nbestvar=0,jstat,k;
    163     Float_t crit,crit0,critmax,critvar=0;
     163    Int_t nc,nbestvar=0,k;
     164    Float_t crit;
    164165    Float_t rrn, rrd, rln, rld, u;
    165166
    … …  
    171172    for (Int_t j=0;j<2;j++)
    172173    {
    173           pno+=tclasspop[j]*tclasspop[j];
    174           pdo+=tclasspop[j];
    175     }
    176     crit0=pno/pdo;
    177     jstat=0;
     174        pno+=tclasspop[j]*tclasspop[j];
     175        pdo+=tclasspop[j];
     176    }
     177
     178    const Double_t crit0=pno/pdo;
     179    Int_t jstat=0;
    178180
    179181    // start main loop through variables to find best split,
    180182    // (Gini-index as criterium crit)
    181183
    182     critmax=-1.0e20;  // FIXME: Replace by a constant from limits.h
     184    Double_t critmax=-1.0e20;  // FIXME: Replace by a constant from limits.h
    183185
    184186    // random split selection, number of trials = fNumTry
    185     if(!gRandom)
    186     {
    187         *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
    188         return kFALSE;
    189     }
    190187    for(Int_t mt=0;mt<fNumTry;mt++)
    191188    {
    192         mvar=Int_t(mdim*gRandom->Rndm());
     189        Int_t mvar=Int_t(gRandom->Rndm()*mdim);
    193190
    194191        // Gini index = rrn/rrd+rln/rld
    … …  
    197194        rln=0;
    198195        rld=0;
    199         wl.Reset();
    200 
    201         for (Int_t j=0;j<2;j++)
    202         {
    203             wr[j]=tclasspop[j];
    204         }
    205 
    206         critvar=-1.0e20;
     196
     197        TArrayF wl(2); // left node
     198        wr = tclasspop;
     199
     200        Double_t critvar=-1.0e20;
    207201
    208202        for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++)
    … …  
    223217            if (datarang[mvar*numdata+nc]<datarang[mvar*numdata+datasort[mvar*numdata+nsp+1]])
    224218            {
    225                 if(TMath::Min(rrd,rld)>1.0e-5)
     219                if (TMath::Min(rrd,rld)>1.0e-5)
    226220                {
    227221                    crit=(rln/rld)+(rrn/rrd);
    … …  
    309303                         TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop,
    310304                         Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent,
    311                          TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc,
    312                          TArrayF &wl,Int_t ninbag)
     305                         TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag)
    313306{
    314307    // Buildtree consists of repeated calls to two void functions, FindBestSplit and MoveData.
    … …  
    326319    // returns.
    327320
    328     Int_t msplit,nbest,ndendl,nc,jstat,ndend,ndstart;
     321    Int_t msplit, nbest, ndendl;
    329322    Float_t decsplit=0;
    330     Float_t popt1,popt2,pp;
    331     TArrayF classpop;
    332     TArrayI nodestatus;
    333 
    334     nodestatus.Set(nrnodes);
    335     classpop.Set(2*nrnodes);
    336 
    337     nodestatus.Reset();
     323    TArrayF classpop(2*nrnodes);
     324    TArrayI nodestatus(nrnodes);
     325
    338326    nodestart.Reset();
    339327    nodepop.Reset();
    340     classpop.Reset();
    341 
    342328
    343329    for (Int_t j=0;j<2;j++)
    … …  
    357343          // initialize for next call to FindBestSplit
    358344
    359           ndstart=nodestart[kbuild];
    360           ndend=ndstart+nodepop[kbuild]-1;
     345          const Int_t ndstart=nodestart[kbuild];
     346          const Int_t ndend=ndstart+nodepop[kbuild]-1;
    361347          for (Int_t j=0;j<2;j++)
    362348            tclasspop[j]=classpop[j*nrnodes+kbuild];
    363349
    364           jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata,
    365                               ndstart,ndend,tclasspop,msplit,decsplit,
    366                               nbest,ncase,jinbag,iv,winbag,wr,wc,wl,
    367                               kbuild);
     350          const Int_t jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata,
     351                                          ndstart,ndend,tclasspop,msplit,decsplit,
     352                                          nbest,ncase,jinbag,iv,winbag,kbuild);
    368353
    369354          if(jstat==1) {
    … …  
    392377          for (Int_t n=ndstart;n<=ndendl;n++)
    393378          {
    394               nc=ncase[n];
    395               Int_t j=hadtrue[nc];
     379              const Int_t nc=ncase[n];
     380              const Int_t j=hadtrue[nc];
    396381              classpop[j*nrnodes+ncur+1]+=winbag[nc];
    397382          }
    … …  
    399384          for (Int_t n=ndendl+1;n<=ndend;n++)
    400385          {
    401               nc=ncase[n];
    402               Int_t j=hadtrue[nc];
     386              const Int_t nc=ncase[n];
     387              const Int_t j=hadtrue[nc];
    403388              classpop[j*nrnodes+ncur+2]+=winbag[nc];
    404389          }
    … …  
    410395          if (nodepop[ncur+1]<=fNdSize) nodestatus[ncur+1]=-1;
    411396          if (nodepop[ncur+2]<=fNdSize) nodestatus[ncur+2]=-1;
    412           popt1=0;
    413           popt2=0;
     397
     398          Double_t popt1=0;
     399          Double_t popt2=0;
    414400          for (Int_t j=0;j<2;j++)
    415401          {
    416             popt1+=classpop[j*nrnodes+ncur+1];
    417             popt2+=classpop[j*nrnodes+ncur+2];
     402              popt1+=classpop[j*nrnodes+ncur+1];
     403              popt2+=classpop[j*nrnodes+ncur+2];
    418404          }
    419405
    420406          for (Int_t j=0;j<2;j++)
    421407          {
    422             if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1;
    423             if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1;
     408              if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1;
     409              if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1;
    424410          }
    425411
    … …  
    446432        {
    447433            fNumEndNodes++;
    448             pp=0;
     434            Double_t pp=0;
    449435            for (Int_t j=0;j<2;j++)
    450436            {
    … …  
    482468
    483469        const Int_t m=fBestVar[kt];
    484 
    485470        kt = event(m)<=fBestSplit[kt] ? fTreeMap1[kt] : fTreeMap2[kt];
    486471    }
    487472
    488473    return fBestSplit[kt];
     474}
     475
     476Double_t MRanTree::TreeHad(const TMatrixRow &event)
     477{
     478    Int_t kt=0;
     479    // to optimize on storage space node status and node class
     480    // are coded into fBestVar:
     481    // status of node kt = TMath::Sign(1,fBestVar[kt])
     482    // class  of node kt = fBestVar[kt]+2 (class defined by larger
     483    //  node population, actually not used)
     484    // hadronness assigned to node kt = fBestSplit[kt]
     485
     486    for (Int_t k=0;k<fNumNodes;k++)
     487    {
     488        if (fBestVar[kt]<0)
     489            break;
     490
     491        const Int_t m=fBestVar[kt];
     492        kt = event(m)<=fBestSplit[kt] ? fTreeMap1[kt] : fTreeMap2[kt];
     493    }
     494
     495    return fBestSplit[kt];
     496}
     497
     498Double_t MRanTree::TreeHad(const TMatrix &m, Int_t ievt)
     499{
     500    return TreeHad(TMatrixRow(m, ievt));
    489501}
    490502

  • trunk/MagicSoft/Mars/mranforest/MRanTree.h

    r2114 r2296  
    1515
    1616class TMatrix;
     17class TMatrixRow;
    1718class TVector;
    1819class TRandom;
    … …  
    6061
    6162    // functions used in tree growing process
    62     void GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,
     63    void GrowTree(const TMatrix &mhad, const TMatrix &mgam,Int_t numdata,
     64                  Int_t numdim,TArrayI &hadtrue,
    6365                  TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag,
    6466                  TArrayF &winbag,TArrayF &weight);
    … …  
    6769                        Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop,
    6870                        Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase,
    69                         TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,
    70                         TArrayF &wc,TArrayF &wl,Int_t kbuild);
     71                        TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t kbuild);
    7172
    7273    void MoveData(TArrayI &datasort,Int_t mdim,Int_t numdata,Int_t ndstart,
    … …  
    7879                   TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop,
    7980                   Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent,
    80                    TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc,
    81                    TArrayF &wl,Int_t ninbag);
     81                   TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag);
    8282
    8383    Double_t TreeHad(const TVector &event);
     84    Double_t TreeHad(const TMatrixRow &event);
     85    Double_t TreeHad(const TMatrix &m, Int_t ievt);
    8486    Double_t TreeHad();
    8587
Note: See TracChangeset for help on using the changeset viewer.