Changeset 7130


Ignore:
Timestamp:
06/03/05 18:02:36 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
21 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r7129 r7130  
    2323
    2424 2005/06/03 Thomas Bretz
     25
     26   * mjobs/hilocalib_sp1.root, mjobs/hilocalib_sp1_mc.root:
     27     - updated
     28
     29   * callisto.rc:
     30     - MJPedestalY2.MaxEvents: 2000 replaced by 5000 as in
     31       callisto_Dec04Jan05.txt
     32
     33   * manalysis/MMultiDimDistCalc.h:
     34     - changes to layout
     35
     36   * mbadpixels/MBadPixelsCalc.cc:
     37     - improved sanity checks
     38
     39   * mbase/MEvtLoop.cc:
     40     - fixed a bug which could cause Eventloop to crahs if
     41       parlist was not initialized
     42
     43   * mdata/MDataArray.[h,cc]:
     44     - added copy constructor
     45
     46   * mhbase/MFillH.cc:
     47     - made sure that no constructor can crash due to NULL pointers
     48
     49   * mhbase/MHMatrix.[h,cc]:
     50     - first check in AddColumn if the column is available. Afterwards
     51       check whether it can be added
     52     - added new interfaced to single rows
     53
     54   * mhflux/MMcSpectrumWeight.cc:
     55     - slight change to screen output
     56
     57   * mjobs/MJPedestal.cc:
     58     - slight change to screen output
     59
     60   * mpedestal/MPedCalcPedRun.cc:
     61     - fixed a bug which caused MC files not to work treat them now
     62       as pedestal files (always)
     63
     64   * mranforest/MRFEnergyEst.[h,cc]:
     65     - improved the code and the interface
     66
     67   * mranforest/MRanForestGrow.[h,cc]:
     68     - derives now from MRead to be able to use the bar in the
     69       display
     70
     71   * mtools/MTFillMatrix.[h,cc]:
     72     - allow to fill a single matrix with all events
    2573
    2674   * datacenter/macros/plotdb.C:

  • trunk/MagicSoft/Mars/NEWS

    r7128 r7130  
    5959   - callisto: In the resource file callisto_Dec04Jan05.rc
    6060     MJPedestalY2.ExtractWinRight has been reduced from 4.0 to 2.0
     61
     62   - callisto: new Hi-/Lo-Gain intercalibration constants
     63     hilocalib_sp1.root and hilocalib_sp1_mc.root
     64
     65   - callisto: changed default for MJPedestalY2.MaxEvents
     66     from 2000 to 5000 like in callisto_Dec04Jan05.txt
    6167
    6268   - callisto: in MCalibrationChargeCalc the limit fgPheErrLowerLimit

  • trunk/MagicSoft/Mars/callisto.rc

    r7127 r7130  
    317317#MJPedestalY3.UseData: Yes
    318318MJPedestalY1.MaxEvents: 500
    319 MJPedestalY2.MaxEvents: 2000
     319MJPedestalY2.MaxEvents: 5000
    320320MJPedestalY3.MaxEvents: 500
    321321

  • trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h

    r6982 r7130  
    1919    TString fHadronnessName;  // Name of container storing hadronness
    2020
    21     MHMatrix    *fMGammas;     //! Gammas describing matrix
    22     MHMatrix    *fMHadrons;    //! Hadrons (non gammas) describing matrix
     21    MHMatrix    *fMGammas;    //! Gammas describing matrix
     22    MHMatrix    *fMHadrons;   //! Hadrons (non gammas) describing matrix
    2323
    2424    MParameterD *fHadronness; //! Output container for calculated hadroness
    2525
    26     MDataArray  *fData;        //! Used to store the MDataChains to get the event values
     26    MDataArray  *fData;       //! Used to store the MDataChains to get the event values
    2727
    28     void StreamPrimitive(ofstream &out) const;
     28    void  StreamPrimitive(ofstream &out) const;
    2929    Int_t PreProcess(MParList *plist);
    3030    Int_t Process();

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc

    r7126 r7130  
    139139Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const
    140140{
    141     if (!fGeomCam || !fPedPhotCam || !fBadPixels)
    142     {
    143         *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - One of the necessary container are not initialized..." << endl;
    144         return kFALSE;
    145     }
    146 
    147     const Bool_t checklo  = fPedestalLevelVarianceLo>0;
    148     const Bool_t checkhi  = fPedestalLevelVarianceHi>0;
     141    const Bool_t checklo = fPedestalLevelVarianceLo>0;
     142    const Bool_t checkhi = fPedestalLevelVarianceHi>0;
    149143
    150144    if (fPedestalLevel<=0 && !checklo && !checkhi)
    151145        return kTRUE;
    152146
     147    if (!fGeomCam || !fPedPhotCam || !fBadPixels)
     148    {
     149        *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - One of the necessary container are not initialized..." << endl;
     150        return kFALSE;
     151    }
     152
    153153    const Int_t entries = fPedPhotCam->GetSize();
    154154
    … …  
    175175
    176176    //if no pixel has a minimum signal, return
     177    Int_t counter=0;
    177178    for (int i=0; i<na; i++)
    178179    {
    179         if (npix[i]==0 || meanrms[i]==0)
    180         {
    181             //fErrors[1]++;          //no valid Pedestals Rms
    182             return kFALSE;
     180        if (npix[i]==0 || meanrms[i]==0) //no valid Pedestals Rms
     181        {
     182            counter++;
     183            continue;
    183184        }
    184185
    185186        meanrms[i] /= npix[i];
    186187        npix[i]=0;
     188    }
     189
     190    if (counter==na)
     191    {
     192        *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - No pixel seems to contain a valid pedestal RMS..." << endl;
     193        return kFALSE;
    187194    }
    188195
    … …  
    209216    MArrayD lolim1(na), lolim2(na); // Precalcualtion of limits
    210217    MArrayD uplim1(na), uplim2(na); // for speeed reasons
     218    counter = 0;
    211219    for (int i=0; i<na; i++)
    212220    {
    213221        if (npix[i]==0 || meanrms2[i]==0)
    214222        {
    215             //fErrors[1]++;          //no valid Pedestals Rms
    216             return kFALSE;
     223            counter++;
     224            continue;
    217225        }
    218226
    … …  
    233241            uplim2[i]   = meanrms2[i]+fPedestalLevelVarianceHi*varrms2[i];
    234242        }
     243    }
     244
     245    if (counter==na)
     246    {
     247        *fLog << err << "MBadPixelsCalc::CheckPedestalRms: ERROR - No pixel seems to contain a valid pedestal RMS anymore..." << endl;
     248        return kFALSE;
    235249    }
    236250

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

    r7091 r7130  
    594594    // If Process has ever been called print statistics
    595595    //
    596     if (fTaskList->GetNumExecutions()>0)
     596    if (fTaskList && fTaskList->GetNumExecutions()>0)
    597597        switch (printstat)
    598598        {

  • trunk/MagicSoft/Mars/mdata/MDataArray.cc

    r5692 r7130  
    5151// --------------------------------------------------------------------------
    5252//
    53 // Constructor
     53// Default Constructor
    5454//
    5555MDataArray::MDataArray(const char *name, const char *title)
    … …  
    6060    gROOT->GetListOfCleanups()->Add(&fList);
    6161    fList.SetBit(kMustCleanup);
     62}
     63
     64// --------------------------------------------------------------------------
     65//
     66// Copy Constructor
     67//
     68MDataArray::MDataArray(MDataArray &a, const char *name, const char *title)
     69{
     70    fName  = name  ? name  : gsDefName.Data();
     71    fTitle = title ? title : gsDefTitle.Data();
     72
     73    gROOT->GetListOfCleanups()->Add(&fList);
     74    fList.SetBit(kMustCleanup);
     75
     76    TIter Next(&a.fList);
     77    MData *data = NULL;
     78    while ((data=(MData*)Next()))
     79        AddEntry(data->GetRule());
    6280}
    6381

  • trunk/MagicSoft/Mars/mdata/MDataArray.h

    r5692 r7130  
    2828public:
    2929    MDataArray(const char *name=NULL, const char *title=NULL);
     30    MDataArray(MDataArray &a, const char *name=NULL, const char *title=NULL);
    3031
    3132    void AddEntry(const TString rule);

  • trunk/MagicSoft/Mars/mhbase/MFillH.cc

    r6978 r7130  
    203203    fHName = hist;
    204204    fParContainer = par;
    205     fParContainerName = par->GetName();
     205    if (par)
     206        fParContainerName = par->GetName();
    206207
    207208    AddToBranchList(Form("%s.*", (const char*)ExtractName(hist)));
    208     AddToBranchList(Form("%s.*", par->GetName()));
     209    if (par)
     210        AddToBranchList(Form("%s.*", par->GetName()));
    209211
    210212    if (!title)
    211         fTitle = Form("Fill %s from %s", fName.Data(), par->GetDescriptor().Data());
     213        fTitle = Form("Fill %s from %s", fName.Data(), par?par->GetDescriptor().Data():"NULL");
    212214}
    213215
    … …  
    226228
    227229    fH = hist;
    228     fHName = hist->GetName();
     230    if (hist)
     231        fHName = hist->GetName();
    229232    fParContainerName = par;
    230233
    231     AddToBranchList(fH->GetDataMember());
     234    if (fH)
     235        AddToBranchList(fH->GetDataMember());
    232236    if (par)
    233237        AddToBranchList(Form("%s.*", (const char*)ExtractName(par)));
    … …  
    236240        return;
    237241
    238     fTitle = Form("Fill %s", hist->GetDescriptor().Data());
     242    fTitle = Form("Fill %s", hist ? hist->GetDescriptor().Data() : "NULL");
    239243    if (!par)
    240244        return;
    … …  
    257261
    258262    fH = hist;
    259     fHName = hist->GetName();
     263    if (hist)
     264        fHName = hist->GetName();
    260265    fParContainer = par;
    261     fParContainerName = par->GetName();
    262 
    263     AddToBranchList(fH->GetDataMember());
    264     AddToBranchList(Form("%s.*", par->GetName()));
     266    if (par)
     267        fParContainerName = par->GetName();
     268
     269    if (fH)
     270        AddToBranchList(fH->GetDataMember());
     271    if (par)
     272        AddToBranchList(Form("%s.*", par->GetName()));
    265273
    266274    if (!title)
    267         fTitle = Form("Fill %s from %s", hist->GetDescriptor().Data(), par->GetDescriptor().Data());
     275        fTitle = Form("Fill %s from %s",
     276                      hist?hist->GetDescriptor().Data():"NULL",
     277                      par ? par->GetDescriptor().Data():"NULL");
    268278}
    269279

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

    r6890 r7130  
    143143Int_t MHMatrix::AddColumn(const char *rule)
    144144{
     145    const Int_t idx = fData ? fData->FindRule(rule) : -1;
     146    if (idx>=0)
     147        return idx;
     148
    145149    if (IsValid(fM))
    146150    {
    … …  
    160164        SetBit(kIsOwner);
    161165    }
    162 
    163     const Int_t idx = fData->FindRule(rule);
    164     if (idx>=0)
    165         return idx;
    166166
    167167    fData->AddEntry(rule);
    … …  
    12621262    return kTRUE;
    12631263}
     1264
     1265// --------------------------------------------------------------------------
     1266//
     1267// Returns the row pointed to by fNumRow into TVector v
     1268//
     1269void MHMatrix::GetRow(TVector &v) const
     1270{
     1271    Int_t ncols = fM.GetNcols();
     1272
     1273    v.ResizeTo(ncols);
     1274
     1275    while (ncols--)
     1276        v(ncols) = fM(fRow, ncols);
     1277}

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

    r6890 r7130  
    7373    void AddColumns(MDataArray *mat);
    7474
     75    //MDataArray *GetColumns() { return fData; }
     76    const MDataArray *GetColumns() const { return fData; }
    7577    MDataArray *GetColumns() { return fData; }
    7678
    … …  
    102104    Double_t operator[](Int_t col) { return fM(fRow, col); }
    103105
     106    void GetRow(TVector &v) const;
     107    void operator>>(TVector &v) const
     108    {
     109        GetRow(v);
     110    }
     111
    104112    Bool_t Fill(MParList *plist, MTask *read, MFilter *filter=NULL);
    105113

  • trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

    r7094 r7130  
    328328    *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
    329329    *fLog << " Simulated spectral slope: " << fOldSlope << endl;
    330     *fLog << " New slope:                " << fNewSlope << endl;
     330    *fLog << " New spectral slope:       " << fNewSlope << endl;
    331331    *fLog << " User normalization:       " << fNorm << endl;
    332332    *fLog << " Old Spectrum:     " << GetFormulaSpecOldX() << "   (I=" << GetSpecOldIntegral() << ")" << endl;

  • trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

    r7125 r7130  
    11221122    {
    11231123        fillpul.SetFilter(&fcalib);
    1124         tlist.AddToList(&decode);
     1124        tlist.AddToList(&decode);
    11251125        tlist.AddToList(&fcalib);
    11261126        tlist.AddToList(&fillpul);
    … …  
    12551255        if (!calc.CheckPedestalRms(fBadPixels, fPedestalCamOut))
    12561256        {
    1257             *fLog << err << "ERROR - Checking PedestalRMS via MBadPixelsCalc failed...." << endl;
     1257            *fLog << err << "ERROR - MBadPixelsCalc::CheckPedestalRms failed...." << endl;
    12581258            return kFALSE;
    12591259        }

  • trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

    r7126 r7130  
    242242  // is not the first anymore
    243243  //
    244   if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal)
    245     {
     244  switch (fRunHeader->GetRunType())
     245  {
     246  case MRawRunHeader::kRTPedestal:
     247  case MRawRunHeader::kRTMonteCarlo:
    246248      fIsFirstPedRun = kFALSE;
    247249      fIsNotPedRun   = kFALSE;
    248250      return kTRUE;
    249     }
    250  
    251   if (fRunHeader->GetRunType()==MRawRunHeader::kRTCalibration)
    252     {
    253       TString proj(fRunHeader->GetProjectName());
    254       proj.ToLower();
    255      
    256       // test if a continuous light run has been declared as calibration...
    257       if (proj.Contains("cl"))
    258         {
    259           fIsFirstPedRun = kFALSE;
    260           fIsNotPedRun   = kFALSE;
    261           return kTRUE;
    262         }
    263     }
    264  
    265 
     251 
     252  case MRawRunHeader::kRTCalibration:
     253      {
     254          TString proj(fRunHeader->GetProjectName());
     255          proj.ToLower();
     256
     257          // test if a continuous light run has been declared as calibration...
     258          if (proj.Contains("cl"))
     259          {
     260              fIsFirstPedRun = kFALSE;
     261              fIsNotPedRun   = kFALSE;
     262              return kTRUE;
     263          }
     264      }
     265  }
    266266
    267267  fIsNotPedRun = kTRUE;
     268
    268269  //
    269270  // If this is the first call to ReInit (before reading the first file)
    … …  
    464465}
    465466
    466 //-------------------------------------------------------------
     467//-----------------------------------------------------------------------
    467468//
    468469// Return if the pedestal bit was set from the calibration trigger box.
    … …  
    473474Bool_t MPedCalcPedRun::IsPedBitSet()
    474475{
     476    if (fRunHeader->GetRunNumber()<gkFirstRunWithFinalBits)
     477        return kFALSE;
     478
    475479    if (!fTrigPattern)
    476         return kFALSE;
    477  
    478     if (fRunHeader->GetRunNumber() < gkFirstRunWithFinalBits)
    479480        return kFALSE;
    480481

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

    r6932 r7130  
    5757#include "MParameters.h"
    5858
    59 
    6059ClassImp(MRFEnergyEst);
    6160
    … …  
    6766// --------------------------------------------------------------------------
    6867//
    69 //
    70 MRFEnergyEst::MRFEnergyEst(const char *name, const char *title):fNumTrees(-1)
    71 {
    72     //
    73     //   set the name and title of this object
    74     //
     68//  Default constructor. Set the name and title of this object
     69//
     70MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
     71    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fData(0), fEnergyEst(0),
     72    fTestMatrix(0)
     73{
    7574    fName  = name  ? name  : gsDefName.Data();
    7675    fTitle = title ? title : gsDefTitle.Data();
    … …  
    7978// --------------------------------------------------------------------------
    8079//
    81 // Delete the data chains
    82 //
    83 MRFEnergyEst::~MRFEnergyEst()
    84 {
    85     //    delete fData;
    86 }
    87 
    88 // --------------------------------------------------------------------------
    89 Int_t MRFEnergyEst::Train()
    90 {
    91     if(!fMatrixTrain)
    92     {
    93         *fLog << err << dbginf << "fMatrixTrain not found... aborting." << endl;
    94         return kFALSE;
    95     }
    96 
    97     if(!fMatrixTrain->GetColumns())
    98     {
    99         *fLog << err << dbginf << "fMatrixTrain does not contain rules... aborting." << endl;
    100         return kFALSE;
    101     }
    102 
    103     const Int_t ncols = (fMatrixTrain->GetM()).GetNcols();
    104     const Int_t nrows = (fMatrixTrain->GetM()).GetNrows();
    105 
    106     cout<<"ncols="<<ncols<<" nrows="<<nrows<<endl<<flush;
    107 
    108     if(ncols<=0 || nrows <=0)
    109     {
    110         *fLog << err << dbginf << "No. of columns or no. of rows of fMatrixTrain equal 0 ... aborting." << endl;
     80// Train the RF with the goven MHMatrix. The last column must contain the
     81// True energy.
     82//
     83Int_t MRFEnergyEst::Train(const MHMatrix &matrixtrain, const TArrayD &grid)
     84{
     85    gLog.Separator("MRFEnergyEst - Train");
     86
     87    if (!matrixtrain.GetColumns())
     88    {
     89        *fLog << err << "ERROR - MHMatrix does not contain rules... abort." << endl;
     90        return kFALSE;
     91    }
     92
     93    const Int_t ncols = matrixtrain.GetM().GetNcols();
     94    const Int_t nrows = matrixtrain.GetM().GetNrows();
     95    if (ncols<=0 || nrows <=0)
     96    {
     97        *fLog << err << "ERROR - No. of columns or no. of rows of matrixtrain equal 0 ... abort." << endl;
     98        return kFALSE;
     99    }
     100
     101    const Int_t nbins = grid.GetSize()-1;
     102    if (nbins<=0)
     103    {
     104        *fLog << err << "ERROR - Energy grid not vaild... abort." << endl;
    111105        return kFALSE;
    112106    }
    113107
    114108    // rules (= combination of image par) to be used for energy estimation
    115     MDataArray used_rules;
    116     TString energy_rule;
    117     for(Int_t i=0;i<ncols;i++)
    118     {
    119         MDataArray *rules=fMatrixTrain->GetColumns();
    120         MData &data=(*rules)[i];
    121 
    122         if(i<ncols-1)
    123             used_rules.AddEntry(data.GetRule());
    124         else
    125             energy_rule=data.GetRule();
    126     }
    127 
    128     if(!energy_rule.Contains("fEnergy"))
    129     {
    130         *fLog << err << dbginf << "Can not accept "<<energy_rule<<" as true energy ... aborting." << endl;
    131         return kFALSE;
    132     }
    133 
    134     TFile fileRF(fRFfileName,"recreate");
    135     if(!fileRF.IsOpen())
    136     {
    137         *fLog << err << dbginf << "File to store RFs could not be opened... aborting." << endl;
    138         return kFALSE;
    139     }
    140 
    141     TMatrix *mptr=(TMatrix*)&(fMatrixTrain->GetM());
    142     const Int_t nbins = fEnergyGrid.GetSize()-1;
    143     if(nbins<=0)
    144     {
    145         *fLog << err << dbginf << "Energy grid not vaild... aborting." << endl;
    146         return kFALSE;
    147     }
     109    TFile fileRF(fFileName, "recreate");
     110    if (!fileRF.IsOpen())
     111    {
     112        *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl;
     113        return kFALSE;
     114    }
     115
     116    MDataArray usedrules;
     117    for(Int_t i=0; i<ncols-3; i++) // -3 is important!!!
     118        usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());
     119
     120    const TMatrix &m = matrixtrain.GetM();
    148121
    149122    // training of RF for each energy bin
    150     Int_t numbins=0;
    151     for(Int_t ie=0;ie<nbins;ie++)
    152     {
    153         TMatrix mat1(nrows,ncols);
    154         TMatrix mat0(nrows,ncols);
     123    for (Int_t ie=0; ie<nbins; ie++)
     124    {
     125        TMatrix mat1(nrows, ncols);
     126        TMatrix mat0(nrows, ncols);
    155127
    156128        // prepare matrix for current energy bin
    … …  
    158130        Int_t irow0=0;
    159131
    160         for(Int_t j=0;j<nrows;j++)
     132        for (Int_t j=0; j<nrows; j++)
    161133        {
    162             Double_t energy=(*mptr)(j,ncols-1);
    163             if(energy>pow(10.,fEnergyGrid[ie]) && energy<=pow(10.,fEnergyGrid[ie+1]))
    164             {
    165                 TMatrixRow(mat1, irow1) = TMatrixRow(*mptr,j);
    166                 irow1++;
    167             }else{
    168                 TMatrixRow(mat0, irow0) = TMatrixRow(*mptr,j);
    169                 irow0++;
    170             }
     134            const Double_t energy = m(j,ncols-1);
     135
     136            if (energy>grid[ie] && energy<=grid[ie+1])
     137                TMatrixFRow(mat1, irow1++) = TMatrixFRow_const(m,j);
     138            else
     139                TMatrixFRow(mat0, irow0++) = TMatrixFRow_const(m,j);
    171140        }
    172         if(irow1*irow0==0)continue;
    173 
    174         *fLog << inf << dbginf << "Training RF for energy bin "<<ie<< endl;
     141
     142        const Bool_t valid = irow1*irow0>0;
     143
     144        if (!valid)
     145            *fLog << warn << "WARNING - Skipping";
     146        else
     147            *fLog << inf << "Training RF for";
     148
     149        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ")" << endl;
     150
     151        if (!valid)
     152            continue;
     153
     154        gLog.SetNullOutput(kTRUE);
    175155
    176156        // last column must be removed (true energy col.)
    177         mat1.ResizeTo(irow1,ncols-1);
    178         mat0.ResizeTo(irow0,ncols-1);
     157        mat1.ResizeTo(irow1, ncols-1);
     158        mat0.ResizeTo(irow0, ncols-1);
    179159
    180160        // create MHMatrix as input for RF
    181         MHMatrix matrix1(mat1,"MatrixHadrons");
    182         MHMatrix matrix0(mat0,"MatrixGammas");
    183 
    184         MDataArray *rules1=matrix1.GetColumns();
    185         MDataArray *rules0=matrix0.GetColumns();
    186         // rules of new matrices be re-set
    187         if(rules1)delete rules1; rules1=new MDataArray();
    188         if(rules0)delete rules0; rules0=new MDataArray();
    189 
    190         for(Int_t i=0;i<ncols-1;i++)
    191         {
    192             //MDataArray *rules=fMatrixTrain->GetColumns();
    193             //MData &data=(*rules)[i];
    194             MData &data=used_rules[i];
    195             rules1->AddEntry(data.GetRule());
    196             rules0->AddEntry(data.GetRule());
    197         }
     161        MHMatrix matrix1(mat1, "MatrixHadrons");
     162        MHMatrix matrix0(mat0, "MatrixGammas");
     163
     164        matrix1.AddColumns(&usedrules);
     165        matrix0.AddColumns(&usedrules);
    198166
    199167        // training of RF
     168        MTaskList tlist;
    200169        MParList plist;
    201         MTaskList tlist;
    202170        plist.AddToList(&tlist);
    203171        plist.AddToList(&matrix0);
    … …  
    212180   
    213181        MEvtLoop evtloop;
     182        evtloop.SetDisplay(fDisplay);
    214183        evtloop.SetParList(&plist);
    215184
    216         //gLog.SetNullOutput(kTRUE);
    217 
    218         if (!evtloop.Eventloop()) return kFALSE;
    219 
    220         //gLog.SetNullOutput(kFALSE);
     185        gLog.SetNullOutput(kFALSE);
     186
     187        if (!evtloop.Eventloop())
     188            return kFALSE;
    221189
    222190        // save whole forest
    223191        MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
    224         forest->SetTitle(Form("%f",0.5*(fEnergyGrid[ie]+fEnergyGrid[ie+1])));
    225         forest->Write(Form("%d",numbins));
    226         numbins++;
     192        const TString title = Form("%f", TMath::Sqrt(grid[ie]*grid[ie+1]));
     193        //const TString title = Form("%f", (grid[ie]+grid[ie+1])/2);
     194        forest->SetTitle(title);
     195        forest->Write(title);
    227196    }
    228197
    229198    // save rules
    230     used_rules.Write("rules");
     199    usedrules.Write("rules");
    231200
    232201    fileRF.Close();
    … …  
    234203    return kTRUE;
    235204}
    236 
    237 Int_t MRFEnergyEst::Test()
    238 {
    239     if(!fMatrixTest)
    240     {
    241         *fLog << err << dbginf << "fMatrixTest not found... aborting." << endl;
    242         return kFALSE;
    243     }
    244 
    245     const Int_t ncols = (fMatrixTest->GetM()).GetNcols();
    246     const Int_t nrows = (fMatrixTest->GetM()).GetNrows();
    247 
    248     if(ncols<=0 || nrows <=0)
    249     {
    250         *fLog << err << dbginf << "No. of columns or no. of rows of fMatrixTrain equal 0 ... aborting." << endl;
    251         return kFALSE;
    252     }
    253 
    254     TMatrix *mptr=(TMatrix*)&(fMatrixTest->GetM());
    255 
    256     if(!ReadForests())
     205/*
     206Int_t MRFEnergyEst::Test(const MHMatrix &matrixtest)
     207{
     208    gLog.Separator("MRFEnergyEst - Test");
     209
     210    const Int_t ncols = matrixtest.GetM().GetNcols();
     211    const Int_t nrows = matrixtest.GetM().GetNrows();
     212
     213    if (ncols<=0 || nrows <=0)
     214    {
     215        *fLog << err << dbginf << "No. of columns or no. of rows of matrixtrain equal 0 ... aborting." << endl;
     216        return kFALSE;
     217    }
     218
     219    if (!ReadForests())
    257220    {
    258221        *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
    … …  
    262225    const Int_t nbins=fEForests.GetSize();
    263226
    264     Double_t e_low = 100;
    265     Double_t e_up  = 0;
     227    Double_t elow =  FLT_MAX;
     228    Double_t eup  = -FLT_MAX;
    266229
    267230    for(Int_t j=0;j<nbins;j++)
    268231    {
    269         e_low = TMath::Min(atof((fEForests[j])->GetTitle()),e_low);
    270         e_up  = TMath::Max(atof((fEForests[j])->GetTitle()),e_up);
    271     }
    272 
    273     TH1F hres("hres","",1000,-10,10);
    274     TH2F henergy("henergy","",100,e_low,e_up,100,e_low,e_up);
    275 
     232        elow = TMath::Min(atof(fEForests[j]->GetTitle()), elow);
     233        eup  = TMath::Max(atof(fEForests[j]->GetTitle()), eup);
     234    }
     235
     236    TH1F hres("hres", "", 1000, -10, 10);
     237    TH2F henergy("henergy", "",100, elow, eup,100, elow, eup);
     238
     239    const TMatrix &m=matrixtest.GetM();
    276240    for(Int_t i=0;i<nrows;i++)
    277241    {
    278         Double_t e_true = (*mptr)(i,ncols-1);
    279         Double_t e_est = 0;
    280         //Double_t hmax  = 0;
     242        Double_t etrue = m(i,ncols-1);
     243        Double_t eest  = 0;
    281244        Double_t hsum  = 0;
    282245
    283246        for(Int_t j=0;j<nbins;j++)
    284247        {
    285             const TVector v=TMatrixRow(*mptr,i);
    286             Double_t h = ((MRanForest*) (fEForests[j]))->CalcHadroness(v);
    287             Double_t e = atof((fEForests[j])->GetTitle());
    288             /*if(h>=hmax)
    289             {
    290                 hmax=h;
    291                 e_est=pow(10.,e);
    292             }*/
    293             hsum+=h;
    294             e_est+=h*e;
     248            const TVector  v = TMatrixFRow_const(m,i);
     249
     250            const Double_t h = ((MRanForest*)fEForests[j])->CalcHadroness(v);
     251            const Double_t e = atof(fEForests[j]->GetTitle());
     252
     253            hsum += h;
     254            eest += h*e;
    295255        }
    296         e_est/=hsum;
    297         e_est=pow(10.,e_est);
    298 
    299         if(e_true>80.) hres.Fill((e_est-e_true)/e_true);
    300         henergy.Fill(log10(e_true),log10(e_est));
    301     }
    302 
    303     if(TestBit(kEnableGraphicalOutput))
    304     {
    305         if(gStyle) gStyle->SetOptFit(1);
    306         // show results
    307         TCanvas *c=MH::MakeDefCanvas("c","",300,900);
    308 
    309         c->Divide(1,3);
    310         c->cd(1);
    311         henergy.SetTitle("Estimated vs true energy");
    312         henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])");
    313         henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])");
    314         henergy.DrawCopy();
    315 
    316         c->cd(2);
    317 
    318         TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S");
    319         hptr->SetTitle("Estimated vs true energy - profile");
    320         hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])");
    321         hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])");
    322         hptr->DrawCopy();
    323 
    324         c->cd(3);
    325         hres.Fit("gaus");
    326         hres.SetTitle("Energy resolution for E_{true}>80Gev");
    327         hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})");
    328         hres.GetYaxis()->SetTitle("counts");
    329         hres.DrawCopy();
    330 
    331 
    332         c->GetPad(1)->SetGrid();
    333         c->GetPad(2)->SetGrid();
    334         c->GetPad(3)->SetGrid();
    335 
    336     }
    337 
    338     return kTRUE;
    339 }
    340 
     256        eest /= hsum;
     257        eest  = pow(10., eest);
     258
     259        //if (etrue>80.)
     260        //    hres.Fill((eest-etrue)/etrue);
     261
     262        henergy.Fill(log10(etrue),log10(eest));
     263    }
     264
     265    if(gStyle)
     266        gStyle->SetOptFit(1);
     267
     268    // show results
     269    TCanvas *c=MH::MakeDefCanvas("c","",300,900);
     270
     271    c->Divide(1,3);
     272    c->cd(1);
     273    henergy.SetTitle("Estimated vs true energy");
     274    henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])");
     275    henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])");
     276    henergy.DrawCopy();
     277
     278    c->cd(2);
     279    TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S");
     280    hptr->SetTitle("Estimated vs true energy - profile");
     281    hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])");
     282    hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])");
     283    hptr->DrawCopy();
     284
     285    c->cd(3);
     286    hres.Fit("gaus");
     287    hres.SetTitle("Energy resolution for E_{true}>80Gev");
     288    hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})");
     289    hres.GetYaxis()->SetTitle("counts");
     290    hres.DrawCopy();
     291
     292    c->GetPad(1)->SetGrid();
     293    c->GetPad(2)->SetGrid();
     294    c->GetPad(3)->SetGrid();
     295
     296    return kTRUE;
     297}
     298*/
    341299Int_t MRFEnergyEst::ReadForests(MParList *plist)
    342300{
    343     TFile fileRF(fRFfileName,"read");
    344     if(!fileRF.IsOpen())
     301    TFile fileRF(fFileName,"read");
     302    if (!fileRF.IsOpen())
    345303    {
    346304        *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
    … …  
    348306    }
    349307
    350     TList *list=(TList*)fileRF.GetListOfKeys();
    351     const Int_t n=list->GetSize()-1;// subtract 1 since 1 key belongs to MDataArray
    352 
    353     fEForests.Expand(n);
    354 
    355     MRanForest forest;
    356     for(Int_t i=0;i<n;i++)
    357     {
    358         forest.Read(Form("%d",i));
    359         MRanForest *curforest=(MRanForest*)forest.Clone();
    360         const char *energy=list->FindObject(Form("%d",i))->GetTitle();
    361         const char *bin   =list->FindObject(Form("%d",i))->GetName();
    362 
    363         curforest->SetTitle(energy);
    364         curforest->SetName(bin);
    365 
    366         fEForests[i]=curforest;
    367     }
    368 
    369     if(plist)
    370     {
    371         fData=(MDataArray*)plist->FindCreateObj("MDataArray");
     308    fEForests.Delete();
     309
     310    TIter Next(fileRF.GetListOfKeys());
     311    TObject *o=0;
     312    while ((o=Next()))
     313    {
     314        MRanForest *forest;
     315        fileRF.GetObject(o->GetName(), forest);
     316        if (!forest)
     317            continue;
     318
     319        forest->SetTitle(o->GetTitle());
     320        forest->SetBit(kCanDelete);
     321
     322        fEForests.Add(forest);
     323
     324    }
     325
     326    if (plist)
     327    {
     328        fData = (MDataArray*)plist->FindCreateObj("MDataArray");
    372329        fData->Read("rules");
    373330    }
     331
    374332    fileRF.Close();
    375333
    … …  
    381339    fEnergyEst = (MParameterD*)plist->FindCreateObj("MParameterD", "MEnergyEst");
    382340    if (!fEnergyEst)
    383     {
    384         *fLog << err << dbginf << "MEnergyEst [MParameterD] not found... aborting." << endl;
    385         return kFALSE;
    386     }
    387 
    388     if(!ReadForests(plist))
    389     {
    390         *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
    391         return kFALSE;
    392     }
     341        return kFALSE;
     342
     343    if (!ReadForests(plist))
     344    {
     345        *fLog << err << "Reading RFs failed... aborting." << endl;
     346        return kFALSE;
     347    }
     348
     349    if (fTestMatrix)
     350        return kTRUE;
    393351
    394352    if (!fData)
    395353    {
    396         *fLog << err << dbginf << "MDataArray not found... aborting." << endl;
     354        *fLog << err << "MDataArray not found... aborting." << endl;
    397355        return kFALSE;
    398356    }
    … …  
    400358    if (!fData->PreProcess(plist))
    401359    {
    402         *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
     360        *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl;
    403361        return kFALSE;
    404362    }
    … …  
    413371{
    414372    TVector event;
    415     *fData >> event;
     373    if (fTestMatrix)
     374        *fTestMatrix >> event;
     375    else
     376        *fData >> event;
    416377
    417378    Double_t eest = 0;
    418     //Double_t hmax  = 0;
    419     Double_t hsum  = 0;
    420        
    421     for(Int_t j=0;j<fEForests.GetSize();j++)
    422     {
    423         Double_t h = ((MRanForest*) (fEForests[j]))->CalcHadroness(event);
    424         Double_t e = atof((fEForests[j])->GetTitle());
    425         /*if(h>=hmax)
    426         {
    427             hmax=h;
    428             e_est=pow(10.,e);
    429         }*/
    430         hsum+=h;
    431         eest+=h*e;
    432     }
    433     eest/=hsum;
    434     eest=pow(10.,eest);
    435 
    436     fEnergyEst->SetVal(eest);
     379    Double_t hsum = 0;
     380
     381    TIter Next(&fEForests);
     382    MRanForest *rf = 0;
     383    while ((rf=(MRanForest*)Next()))
     384    {
     385        const Double_t h = rf->CalcHadroness(event);
     386        const Double_t e = atof(rf->GetTitle());
     387
     388        hsum += h;
     389        eest += h*e;
     390    }
     391
     392    fEnergyEst->SetVal(eest/hsum);
    437393    fEnergyEst->SetReadyToSave();
    438394

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

    r6932 r7130  
    66#endif
    77
    8 #ifndef ROOT_TArrayD
    9 #include <TArrayD.h>
    10 #endif
    11 
    128#ifndef ROOT_TObjArray
    139#include <TObjArray.h>
    1410#endif
     11
     12class TArrayD;
    1513
    1614class MHMatrix;
    … …  
    2119{
    2220private:
    23     Int_t fNumTrees;
    24     Int_t fNumTry;
    25     Int_t fNdSize;
     21    Int_t        fNumTrees;  // Training parameters
     22    Int_t        fNumTry;    // Training parameters
     23    Int_t        fNdSize;    // Training parameters
    2624
    27     TString fRFfileName;
    28     MHMatrix *fMatrixTrain;
    29     MHMatrix *fMatrixTest;
    30     TArrayD fEnergyGrid;
     25    TString      fFileName;
     26    TObjArray    fEForests;
    3127
    3228    MDataArray  *fData;       //! Used to store the MDataChains to get the event values
    33     TObjArray fEForests;
     29    MParameterD *fEnergyEst;  //! Used to storeestimated energy
    3430
    35     MParameterD *fEnergyEst;
     31    MHMatrix *fTestMatrix;
    3632
    3733    Int_t PreProcess(MParList *plist);
    … …  
    4238public:
    4339    MRFEnergyEst(const char *name=NULL, const char *title=NULL);
    44     ~MRFEnergyEst();
    4540
    46     void SetMatrixTrain(MHMatrix *mat) { fMatrixTrain = mat; }
    47     void SetMatrixTest( MHMatrix *mat) { fMatrixTest  = mat; }
    48     void SetFile(TString str) { fRFfileName = str; }
     41    void  SetFileName(TString str)     { fFileName = str; }
    4942
    50     void SetLogEnergyGrid(TArrayD &egrid) { fEnergyGrid = egrid ; }
     43    void  SetNumTrees(UShort_t n=-1)   { fNumTrees = n; }
     44    void  SetNdSize(UShort_t n=-1)     { fNdSize   = n; }
     45    void  SetNumTry(UShort_t n=-1)     { fNumTry   = n; }
    5146
    52     void SetNumTrees(UShort_t n=100) { fNumTrees = n; }
    53     void SetNdSize(UShort_t n=1)     { fNdSize   = n; }
    54     void SetNumTry(UShort_t n)       { fNumTry   = n; }
     47    void  SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; }
     48    void  InitMapping(MHMatrix *m=0)   { fTestMatrix=m; }
    5549
    56     Int_t Train();
    57     Int_t Test();
     50    Int_t Train(const MHMatrix &n, const TArrayD &grid);
    5851
    5952    ClassDef(MRFEnergyEst, 0) // Task

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

    r2207 r7130  
    4646using namespace std;
    4747
    48 static const TString gsDefName  = "MRanForestGrow";
     48static const TString gsDefName  = "MRead";
    4949static const TString gsDefTitle = "Tree Classification Loop 1/2";
    5050
    … …  
    5555//
    5656MRanForestGrow::MRanForestGrow(const char *name, const char *title)
    57     : fNumTrees(100),fNumTry(3),fNdSize(1)
    5857{
    5958    //
    … …  
    6261    fName  = name  ? name  : gsDefName.Data();
    6362    fTitle = title ? title : gsDefTitle.Data();
     63
     64    SetNumTrees();
     65    SetNumTry();
     66    SetNdSize();
    6467}
    6568
    … …  
    124127Int_t MRanForestGrow::Process()
    125128{
    126     Bool_t not_last=fRanForest->GrowForest();
     129    const Bool_t not_last=fRanForest->GrowForest();
     130
    127131    fRanTree->SetReadyToSave();
    128132

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

    r2207 r7130  
    22#define MARS_MRanForestGrow
    33
    4 /////////////////////////////////////////////////////////////////////////////
    5 //                                                                         //
    6 // MRanForestGrow                                                          //
    7 //                                                                         //
    8 // Task to grow (train) a random forest                                    //
    9 //                                                                         //
    10 /////////////////////////////////////////////////////////////////////////////
    11 
    12 #ifndef MARS_MTask
    13 #include "MTask.h"
     4#ifndef MARS_MRead
     5#include "MRead.h"
    146#endif
    157
    … …  
    1911class MRanTree;
    2012
    21 class MRanForestGrow : public MTask
     13class MRanForestGrow : public MRead
    2214{
    2315private:
    … …  
    3527    Int_t PostProcess();
    3628
     29    UInt_t  GetEntries()            { return fNumTrees; }
     30    TString GetFileName() const     { return "MRanForestGrow"; }
     31    TString GetFullFileName() const { return "MRanForestGrow"; }
     32
    3733public:
    3834    MRanForestGrow(const char *name=NULL, const char *title=NULL);
    3935
    40     void SetNumTrees(Int_t n){    fNumTrees=n;}
    41     void SetNumTry(Int_t n)  {    fNumTry=n; }
    42     void SetNdSize(Int_t n)  {    fNdSize=n; }
     36    void SetNumTrees(Int_t n=-1) { fNumTrees=n>0?n:100; }
     37    void SetNumTry(Int_t   n=-1) { fNumTry  =n>0?n:  3; }
     38    void SetNdSize(Int_t   n=-1) { fNdSize  =n>0?n:  1; }
    4339
    4440    ClassDef(MRanForestGrow, 0) // Task to grow a random forest

  • trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc

    r2744 r7130  
    109109Bool_t MTFillMatrix::CheckResult(MHMatrix *m, Int_t num) const
    110110{
    111     if (!m)
     111    if (!m || num==0)
    112112        return kTRUE;
    113113
    … …  
    190190    fLog->Separator(GetDescriptor());
    191191    *fLog << "Fill " << fDestMatrix1->GetDescriptor() << " with " << fNumDestEvents1 << endl;
    192     *fLog << "Fill " << fDestMatrix2->GetDescriptor() << " with " << fNumDestEvents2 << endl;
     192    if (fDestMatrix2)
     193        *fLog << "Fill " << fDestMatrix2->GetDescriptor() << " with " << fNumDestEvents2 << endl;
    193194    *fLog << "Distribution choosen ";
    194195    if (fReference && fReference->GetHist().GetEntries()>0)
    … …  
    210211    // FIXME: Merge MFEventSelector and MFEventSelector2
    211212    MFilter *selector=0;
    212     if (fReference)
    213     {
    214         // Case of a reference/nominal distribution
    215         // The events must be read before selection
    216         MFEventSelector2 *sel = new MFEventSelector2(*fReference);
    217         sel->SetNumMax(fNumDestEvents1+fNumDestEvents2);
    218         sel->SetInverted();
    219 
    220         selector = sel;
    221     }
    222     else
    223     {
    224         // Case of a random distribution
    225         // The events can be selected before reading
    226         MFEventSelector *sel = new MFEventSelector;
    227         sel->SetNumSelectEvts(fNumDestEvents1+fNumDestEvents2);
    228         fReader->SetSelector(sel);
    229 
    230         selector = sel;
     213    if (fNumDestEvents1>0 || fNumDestEvents2>0)
     214    {
     215        if (fReference)
     216        {
     217            // Case of a reference/nominal distribution
     218            // The events must be read before selection
     219            MFEventSelector2 *sel = new MFEventSelector2(*fReference);
     220            sel->SetNumMax(fNumDestEvents1+fNumDestEvents2);
     221            sel->SetInverted();
     222
     223            selector = sel;
     224        }
     225        else
     226        {
     227            // Case of a random distribution
     228            // The events can be selected before reading
     229            MFEventSelector *sel = new MFEventSelector;
     230            sel->SetNumSelectEvts(fNumDestEvents1+fNumDestEvents2);
     231            fReader->SetSelector(sel);
     232
     233            selector = sel;
     234        }
    231235    }
    232236
    … …  
    256260    MFillH fill1(fDestMatrix1);
    257261    MFillH fill2(fDestMatrix2);
    258     fill1.SetFilter(&split);
    259     fill2.SetFilter(&invsplit);
     262    if (selector)
     263    {
     264        fill1.SetFilter(&split);
     265        fill2.SetFilter(&invsplit);
     266    }
    260267
    261268    // entries in MTaskList
    262269    tlist.AddToList(fReader);        // Read events
    263     if (fReference)
     270    if (fReference && selector)
    264271        tlist.AddToList(&cont);      // select a sample of events
    265272    tlist.AddToList(&invsplit);      // process invsplit (which implicitly processes split)
    266     if (fDestMatrix1 && fNumDestEvents1>0)
     273    if (fDestMatrix1)
    267274        tlist.AddToList(&fill1);     // fill matrix 1
    268     if (fDestMatrix2 && fNumDestEvents2>0)
     275    if (fDestMatrix2)
    269276        tlist.AddToList(&fill2);     // fill matrix 2
    270277    if (fWriteFile1)
    … …  
    289296    const Bool_t rc = evtloop.Eventloop();
    290297
    291     // Print execution statistics of the tasklist
    292     if (rc)
    293         tlist.PrintStatistics();
    294 
    295     delete selector;
     298    if (selector)
     299        delete selector;
    296300
    297301    if (!rc)
Note: See TracChangeset for help on using the changeset viewer.