Changeset 1489 for trunk


Ignore:
Timestamp:
08/08/02 11:58:59 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r1487 r1489  
    11                                                                  -*-*- END -*-*-
     2
     3 2002/08/08: Thomas Bretz
     4
     5   * manalysis/MHillasSrc.cc:
     6     - use double dist instead of single fDist for calculation
     7
     8   * manalysis/MMultiDimDistCalc.[h,cc]:
     9     - added support for the kernel method
     10     - added stream primitive
     11     - changed version number to 1
     12     - adapted to new MHMatrix (using MDataArray)
     13
     14   * mdata/MDataArray.[h,cc]:
     15     - added
     16
     17   * mdata/DataLinkDef.h, madata/Makefile:
     18     - added MDataArray
     19
     20   * mfileio/MWriteRootFile.cc:
     21     - fixed some bugs in StreamPrimitive
     22     - StreamPrimtive doesn't write the default name/title anymore
     23
     24   * mhist/MHMatrix.[h,cc]:
     25     - replaced the Arrays for the rules by a MDataArray
     26     - implemented StreamPrimitive
     27     - implement the use of the kernel function for num<0
     28     - multiply fM2 by nevts-1
     29     - added sanity check in case of dists[i]<0
     30
     31   * mhist/MHHillas.[h,cc]:
     32     - added fUsedPix, fCorePix
     33     - added fUsedPix, fCorePix to plots
     34     - changed layout of plots
     35     - changed name and title of MakeDefCanvas
     36
     37   * mhist/MHHillasSrc.[h,cc]:
     38     - changed plot of Alpha from fabs(fAlpha) to fAlpha
     39     - changed name and title of MakeDefCanvas
     40
     41   * mhist/MHillasExt.[h,cc]:
     42     - changed layout of plots
     43     - changed name and title of MakeDefCanvas
     44
     45   * mbase/MTask.cc:
     46     - fixed wrong streaming of filter name
     47
     48   * macros/starplot.C:
     49     - added
     50
     51   * macros/dohtml.C:
     52     - added starplot.C
     53
     54
    255
    356 2002/08/07: Thomas Bretz

  • trunk/MagicSoft/Mars/macros/dohtml.C

    r1487 r1489  
    5252    html.Convert("trigrate.C",     "MARS - Calculate Trigger Rate from a MC root file");
    5353    html.Convert("star.C",         "MARS - (St)andard (A)nalysis and (R)econstruction");
     54    html.Convert("starplot.C",     "MARS - Plot parameters from file created with star.C");
    5455    html.Convert("comprob.C",      "MARS - Calculation of composite probabilities for G/H-Seperation");
    5556    html.Convert("multidimdist.C", "MARS - Calculation of multidimensional distances for G/H-Seperation");

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

    r1488 r1489  
    6969{
    7070    fName  = name  ? name  : "MHillasSrc";
    71     fTitle = title ? title : "parameters depending in source position";
     71    fTitle = title ? title : "Parameters depending in source position";
    7272}
    7373

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

    r1488 r1489  
    5252#include "MParList.h"
    5353#include "MDataChain.h"
     54#include "MDataArray.h"
    5455
    5556#include "MHadroness.h"
    … …  
    6566//
    6667MMultiDimDistCalc::MMultiDimDistCalc(const char *name, const char *title)
    67     : fNum(0), fUseKernel(kFALSE)
     68    : fNum(0), fUseKernel(kFALSE), fData(NULL)
    6869{
    6970    //
    … …  
    7374    fTitle = title ? title : gsDefTitle.Data();
    7475
    75     fData = new TList;
    76     fData->SetOwner();
     76    /*
     77     fData = new TList;
     78     fData->SetOwner();
     79     */
    7780}
    7881
    … …  
    8386MMultiDimDistCalc::~MMultiDimDistCalc()
    8487{
    85     delete fData;
     88    //    delete fData;
    8689}
    8790
    … …  
    119122    }
    120123
     124    /*
    121125    TIter Next(fMGammas->GetRules());
    122126    TObject *data=NULL;
    … …  
    130134        }
    131135        fData->Add(chain);
     136    }
     137    */
     138    fData = fMGammas->GetColumns();
     139    if (!fData)
     140    {
     141        *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl;
     142        return kFALSE;
     143    }
     144
     145    if (!fData->PreProcess(plist))
     146    {
     147        *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
     148        return kFALSE;
    132149    }
    133150
    … …  
    156173    TVector event(ncols);
    157174
     175    for (int i=0; i<fData->GetNumEntries(); i++)
     176        event(i) = (*fData)(i);
     177
     178    /*
    158179    Int_t n=0;
    159180    TIter Next(fData);
    … …  
    161182    while ((data=(MData*)Next()))
    162183        event(n++) = data->GetValue();
     184        */
    163185
    164186    Double_t numg = fNum;

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

    r1488 r1489  
    99class MParList;
    1010class MHadroness;
     11class MDataArray;
    1112
    1213class MMultiDimDistCalc : public MTask
    … …  
    2122    MHadroness *fHadroness; //! Output container for calculated hadroness
    2223
    23     TList *fData;           //! Used to store the MDataChains to get the event values
     24    MDataArray *fData;      //! Used to store the MDataChains to get the event values
    2425
    2526    void StreamPrimitive(ofstream &out) const;

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

    r1479 r1489  
    281281
    282282     fFilter->SavePrimitive(out);
    283      out << "   " << ToLower(fName) << ".SetFilter(&" << ToLower(fFilter->GetName()) <<");" << endl;
    284283     */
    285 }
     284     out << "   " << GetUniqueName() << ".SetFilter(&" << fFilter->GetUniqueName() <<");" << endl;
     285}

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

    r1487 r1489  
    6767    fTitle = title ? title : "Source independent image parameters";
    6868
    69     fLength = new TH1F("Length", "Length of Ellipse",               100,   0, 296.7);
    70     fWidth  = new TH1F("Width",  "Width of Ellipse",                100,   0, 296.7);
    71     fDistC  = new TH1F("DistC",  "Distance from center of camera",  100,   0, 445);
    72     fDelta  = new TH1F("Delta",  "Angle (Main axis - x-axis)",      101, -90,  90);
     69    fLength  = new TH1F("Length",  "Length of Ellipse",               100,   0, 296.7);
     70    fWidth   = new TH1F("Width",   "Width of Ellipse",                100,   0, 296.7);
     71    fDistC   = new TH1F("DistC",   "Distance from center of camera",  100,   0, 445);
     72    fDelta   = new TH1F("Delta",   "Angle (Main axis - x-axis)",      101, -90,  90);
     73    fUsedPix = new TH1F("UsedPix", "Number of used pixels",           150,   0, 150);
     74    fCorePix = new TH1F("CorePix", "Number of core pixels",           150,   0, 150);
    7375
    7476    fLength->SetDirectory(NULL);
    … …  
    7678    fDistC->SetDirectory(NULL);
    7779    fDelta->SetDirectory(NULL);
     80    fUsedPix->SetDirectory(NULL);
     81    fCorePix->SetDirectory(NULL);
    7882
    7983    fLength->SetXTitle("Length [mm]");
    … …  
    8185    fDistC->SetXTitle("Distance [mm]");
    8286    fDelta->SetXTitle("Delta [\\circ]");
     87    fUsedPix->SetXTitle("Number of Pixels");
     88    fCorePix->SetXTitle("Number of Pixels");
    8389
    8490    fLength->SetYTitle("Counts");
    … …  
    8692    fDistC->SetYTitle("Counts");
    8793    fDelta->SetYTitle("Counts");
     94    fUsedPix->SetYTitle("Counts");
     95    fCorePix->SetYTitle("Counts");
    8896
    8997    MBinning bins;
    … …  
    122130    delete fSize;
    123131    delete fCenter;
     132
     133    delete fUsedPix;
     134    delete fCorePix;
    124135}
    125136
    … …  
    150161    ApplyBinning(*plist, "Delta",  fDelta);
    151162    ApplyBinning(*plist, "Size",   fSize);
     163    ApplyBinning(*plist, "Pixels", fUsedPix);
     164    ApplyBinning(*plist, "Pixels", fCorePix);
    152165
    153166    const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
    … …  
    235248    const Double_t scale = fUseMmScale ? 1 : fMm2Deg;
    236249
    237     fLength->Fill(scale*h.GetLength());
    238     fWidth ->Fill(scale*h.GetWidth());
    239     fDistC ->Fill(scale*d);
    240     fCenter->Fill(scale*h.GetMeanX(), scale*h.GetMeanY());
    241     fDelta ->Fill(kRad2Deg*h.GetDelta());
    242     fSize  ->Fill(h.GetSize());
     250    fLength ->Fill(scale*h.GetLength());
     251    fWidth  ->Fill(scale*h.GetWidth());
     252    fDistC  ->Fill(scale*d);
     253    fCenter ->Fill(scale*h.GetMeanX(), scale*h.GetMeanY());
     254    fDelta  ->Fill(kRad2Deg*h.GetDelta());
     255    fSize   ->Fill(h.GetSize());
     256    fUsedPix->Fill(h.GetNumUsedPixels());
     257    fCorePix->Fill(h.GetNumCorePixels());
    243258
    244259    return kTRUE;
    … …  
    269284TObject *MHHillas::DrawClone(Option_t *opt) const
    270285{
    271     TCanvas *c = MakeDefCanvas("Hillas", fTitle, 720, 810);
     286    TCanvas *c = MakeDefCanvas(this, 720, 810);
    272287    c->Divide(2,3);
    273288
    … …  
    275290
    276291    c->cd(1);
    277     fLength->DrawCopy();
     292    fWidth->DrawCopy();
     293    fLength->SetLineColor(kBlue);
     294    fLength->DrawCopy("same");
    278295
    279296    c->cd(2);
    280     fWidth->DrawCopy();
    281 
    282     c->cd(3);
    283297    gPad->SetLogx();
    284298    fSize->DrawCopy();
     299
     300    c->cd(3);
     301    fCorePix->SetLineColor(kRed);
     302    fCorePix->DrawCopy();
     303    fUsedPix->SetLineColor(kGreen);
     304    fUsedPix->DrawCopy("same");
    285305
    286306    c->cd(4);
    … …  
    309329{
    310330    if (!gPad)
    311         MakeDefCanvas("Hillas", fTitle, 720, 810);
     331        MakeDefCanvas(this, 720, 810);
    312332
    313333    gPad->Divide(2,3);
    314334
    315335    gPad->cd(1);
    316     fLength->Draw();
     336    fWidth->Draw();
     337    fLength->SetLineColor(kBlue);
     338    fLength->Draw("same");
    317339
    318340    gPad->cd(2);
    319     fWidth->Draw();
    320 
    321     gPad->cd(3);
    322341    gPad->SetLogx();
    323342    fSize->Draw();
     343
     344    gPad->cd(3);
     345    fCorePix->SetLineColor(kRed);
     346    fCorePix->Draw();
     347    fUsedPix->SetLineColor(kGreen);
     348    fUsedPix->Draw("same");
    324349
    325350    gPad->cd(4);

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

    r1463 r1489  
    2222    TH1F *fSize;   //->
    2323    TH2F *fCenter; //->
     24
     25    TH1F *fUsedPix; //->
     26    TH1F *fCorePix; //->
    2427
    2528    void SetColors() const;

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

    r1465 r1489  
    6262    //
    6363    fName  = name  ? name  : "MHHillasExt";
    64     fTitle = title ? title : "Container for Hillas (ext) histograms";
     64    fTitle = title ? title : "Container for extended Hillas histograms";
    6565
    6666    //
    … …  
    6969    // connect all the histogram with the container fHist
    7070    //
     71    fHConc1.SetLineColor(kBlue);
     72    fHConc.SetFillStyle(0);
     73
    7174    fHConc.SetDirectory(NULL);
    7275    fHConc1.SetDirectory(NULL);
    … …  
    220223TObject *MHHillasExt::DrawClone(Option_t *opt) const
    221224{
    222     TCanvas &c = *MakeDefCanvas("Hillas", "Histograms of Hillas Parameters",
    223                                 720, 810);
    224     c.Divide(2, 3);
     225    TCanvas &c = *MakeDefCanvas(this, 720, 540);
     226    c.Divide(2, 2);
    225227
    226228    gROOT->SetSelectedPad(NULL);
    227229
    228     //
    229     // This is necessary to get the expected bahviour of DrawClone
    230     //
     230
    231231    c.cd(1);
    232     ((TH1F&)fHConc).DrawCopy();
     232    ((TH1F)fHConc1).DrawCopy();
     233    ((TH1F)fHConc).DrawCopy("same");
    233234
    234235    c.cd(2);
    235     ((TH1F&)fHConc1).DrawCopy();
    236 
    237     c.cd(5);
    238     ((TH1F&)fHAsym).DrawCopy();
     236    ((TH1F)fHAsym).DrawCopy();
    239237
    240238    c.cd(3);
    241     ((TH1F&)fHM3Long).DrawCopy();
     239    ((TH1F)fHM3Long).DrawCopy();
    242240
    243241    c.cd(4);
    244     ((TH1F&)fHM3Trans).DrawCopy();
     242    ((TH1F)fHM3Trans).DrawCopy();
    245243
    246244    c.Modified();
    … …  
    261259{
    262260    if (!gPad)
    263         MakeDefCanvas("Hillas", "Histograms of Hillas Parameters",
    264                       720, 810);
    265 
    266     gPad->Divide(2, 3);
     261        MakeDefCanvas(this, 720, 540);
     262
     263    gPad->Divide(2, 2);
    267264
    268265    gPad->cd(1);
    269     fHConc.DrawCopy();
     266    fHConc1.SetLineColor(kBlue);
     267    fHConc1.Draw();
     268    fHConc.Draw("same");
    270269
    271270    gPad->cd(2);
    272     fHConc1.DrawCopy();
    273 
    274     gPad->cd(5);
    275     fHAsym.DrawCopy();
     271    fHAsym.Draw();
    276272
    277273    gPad->cd(3);
    278     fHM3Long.DrawCopy();
     274    fHM3Long.Draw();
    279275
    280276    gPad->cd(4);
    281     fHM3Trans.DrawCopy();
     277    fHM3Trans.Draw();
    282278
    283279    gPad->Modified();

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

    r1465 r1489  
    6868    // connect all the histogram with the container fHist
    6969    //
    70     fAlpha    = new TH1F("Alpha",    "Alpha of Ellipse",             90,    0,  90);
     70    fAlpha    = new TH1F("Alpha",    "Alpha of Ellipse",            181,  -90,  90);
    7171    fDist     = new TH1F("Dist",     "Dist of Ellipse",             100,    0, 445);
    7272    fHeadTail = new TH1F("HeadTail", "HeadTail of Ellipse",         101, -445, 445);
    … …  
    8181    fDist->SetXTitle("Dist [mm]");
    8282    fHeadTail->SetXTitle("Head-Tail [mm]");
    83     fCosDA->SetXTitle("cos(\\delta,\\alpha) [mm]");
     83    fCosDA->SetXTitle("cos(\\delta,\\alpha)");
    8484
    8585    fAlpha->SetYTitle("Counts");
    … …  
    138138    const MHillasSrc &h = *(MHillasSrc*)par;
    139139
    140     fAlpha   ->Fill(fabs(h.GetAlpha()));
     140    fAlpha   ->Fill(h.GetAlpha());
    141141    fDist    ->Fill(fUseMmScale ? h.GetDist() : fMm2Deg*h.GetDist());
    142142    fHeadTail->Fill(fUseMmScale ? h.GetHeadTail() : fMm2Deg*h.GetHeadTail());
    … …  
    211211TObject *MHHillasSrc::DrawClone(Option_t *opt) const
    212212{
    213     TCanvas *c = MakeDefCanvas("HillasSrc", "Histograms of Source dependant Parameters",
    214                                700, 500);
     213    TCanvas *c = MakeDefCanvas(this, 700, 500);
    215214    c->Divide(2, 2);
    216215
    … …  
    250249{
    251250    if (!gPad)
    252         MakeDefCanvas("HillasSrc", "Histograms of Src dependant Parameters",
    253                       700, 500);
     251        MakeDefCanvas(this, 700, 500);
    254252
    255253    // FIXME: Display Source position

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

    r1355 r1489  
    4545#include "MHMatrix.h"
    4646
     47#include <fstream.h>
     48
    4749#include <TList.h>
    4850#include <TArrayD.h>
    … …  
    5456#include "MParList.h"
    5557
    56 #include "MDataChain.h"
     58#include "MData.h"
     59#include "MDataArray.h"
    5760
    5861ClassImp(MHMatrix);
    5962
    60 // --------------------------------------------------------------------------
    61 //
    62 // Default Constructor
     63static const TString gsDefName  = "MHMatrix";
     64static const TString gsDefTitle = "Multidimensional Matrix";
     65
     66// --------------------------------------------------------------------------
     67//
     68//  Default Constructor
    6369//
    6470MHMatrix::MHMatrix(const char *name, const char *title)
    65     : fNumRow(0)
    66 {
    67     fName  = name  ? name  : "MHMatrix";
    68     fTitle = title ? title : "Multidimensional Matrix";
    69 
    70     fData = new TList;
    71     fData->SetOwner();
    72 
    73     fRules = new TList;
    74     fRules->SetOwner();
    75 }
    76 
    77 // --------------------------------------------------------------------------
    78 //
    79 // Destructor
     71    : fNumRow(0), fData(NULL)
     72{
     73    fName  = name  ? name  : gsDefName.Data();
     74    fTitle = title ? title : gsDefTitle.Data();
     75}
     76
     77// --------------------------------------------------------------------------
     78//
     79//  Constructor. Initializes the columns of the matrix with the entries
     80//  from a MDataArray
     81//
     82MHMatrix::MHMatrix(MDataArray *mat, const char *name, const char *title)
     83    : fNumRow(0), fData(mat)
     84{
     85    fName  = name  ? name  : gsDefName.Data();
     86    fTitle = title ? title : gsDefTitle.Data();
     87}
     88
     89// --------------------------------------------------------------------------
     90//
     91//  Destructor. Does not deleted a user given MDataArray, except IsOwner
     92//  was called.
    8093//
    8194MHMatrix::~MHMatrix()
    8295{
    83     delete fData;
    84     delete fRules;
     96    if (TestBit(kIsOwner) && fData)
     97        delete fData;
    8598}
    8699
    … …  
    99112    }
    100113
    101     MDataChain &chain = *new MDataChain(rule);
    102     if (!chain.IsValid())
    103     {
    104         *fLog << err << "Error - Rule cannot be translated... ignored." << endl;
    105         delete &chain;
     114    if (!fData)
     115    {
     116        fData = new MDataArray;
     117        SetBit(kIsOwner);
     118    }
     119
     120    fData->AddEntry(rule);
     121}
     122
     123void MHMatrix::AddColumns(MDataArray *matrix)
     124{
     125    if (fM.IsValid())
     126    {
     127        *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl;
    106128        return;
    107129    }
    108     fData->Add(&chain);
    109 
    110     TNamed *name = new TNamed(rule, rule); // Fimxe, in 3.02/07 the title can't be "", why?
    111     fRules->Add(name);
    112 }
    113 
    114 void MHMatrix::AddColumns(const MHMatrix *matrix)
    115 {
    116     TIter Next(matrix->fRules);
    117     TObject *rule=NULL;
    118     while ((rule=Next()))
    119         AddColumn(rule->GetName());
     130
     131    if (fData)
     132        *fLog << warn << "Warning - columns already added... replacing." << endl;
     133
     134    if (fData && TestBit(kIsOwner))
     135    {
     136        delete fData;
     137        ResetBit(kIsOwner);
     138    }
     139
     140    fData = matrix;
    120141}
    121142
    … …  
    127148Bool_t MHMatrix::SetupFill(const MParList *plist)
    128149{
    129     if (fData->GetSize()==0)
    130     {
    131         *fLog << err << "Error - No Column specified... aborting." << endl;
     150    if (!fData)
     151    {
     152        *fLog << err << "Error - No Columns initialized... aborting." << endl;
    132153        return kFALSE;
    133154    }
    134155
    135     TIter Next(fData);
    136     MData *data = NULL;
    137     while ((data=(MData*)Next()))
    138         if (!data->PreProcess(plist))
    139             return kFALSE;
    140 
    141     return kTRUE;
     156    return fData->PreProcess(plist);
    142157}
    143158
    … …  
    155170    if (!fM.IsValid())
    156171    {
    157         fM.ResizeTo(1, fData->GetSize());
     172        fM.ResizeTo(1, fData->GetNumEntries());
    158173        return;
    159174    }
    … …  
    161176    TMatrix m(fM);
    162177
    163     fM.ResizeTo(fM.GetNrows()*2, fData->GetSize());
     178    fM.ResizeTo(fM.GetNrows()*2, fData->GetNumEntries());
    164179
    165180    for (int x=0; x<m.GetNcols(); x++)
    … …  
    176191    AddRow();
    177192
    178     int col=0;
    179     TIter Next(fData);
    180     MData *data = NULL;
    181     while ((data=(MData*)Next()))
    182         fM(fNumRow-1, col++) = data->GetValue();
     193    for (int col=0; col<fData->GetNumEntries(); col++)
     194        fM(fNumRow-1, col) = (*fData)(col);
    183195
    184196    return kTRUE;
    … …  
    195207    // It's not a fatal error so we don't need to stop PostProcessing...
    196208    //
    197     if (fData->GetSize()<1 || fNumRow<1)
     209    if (fData->GetNumEntries()==0 || fNumRow<1)
    198210        return kTRUE;
    199211
    200212    TMatrix m(fM);
    201213
    202     fM.ResizeTo(fNumRow, fData->GetSize());
     214    fM.ResizeTo(fNumRow, fData->GetNumEntries());
    203215
    204216    for (int x=0; x<fM.GetNcols(); x++)
    … …  
    287299void MHMatrix::Print(Option_t *) const
    288300{
    289     int n=0;
    290 
    291     TIter Next(fData->GetSize() ? fData : fRules);
    292     MData *data = NULL;
    293     while ((data=(MData*)Next()))
    294     {
    295         *fLog << all << " Column " << setw(3) << n++ << ": " << flush;
    296         data->Print();
    297         *fLog << endl;
    298     }
     301    if (fData)
     302        fData->Print();
     303    else
     304        *fLog << all << "Sorry, no column information available." << endl;
    299305
    300306    fM.Print();
    … …  
    336342}
    337343
     344// --------------------------------------------------------------------------
     345//
     346// Calculated the distance of vector evt from the reference sample
     347// represented by the covariance metrix m.
     348//  - If n<0 the kernel method is applied and
     349//    sum(epx(-d)/n is returned.
     350//  - For n>=0 the n nearest neighbors are summed and
     351//    sqrt(sum(d))/n is returned.
     352//  - if n==0 all distances are summed
     353//
    338354Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const
    339355{
     356    if (num==0) // may later be used for another method
     357    {
     358        TVector d = evt;
     359        d *= m;
     360        return sqrt(d*evt);
     361    }
     362
    340363    const Int_t rows = fM.GetNrows();
    341364    const Int_t cols = fM.GetNcols();
    … …  
    358381
    359382        dists[i] = d2*d; // square of distance
     383
     384        //
     385        // This corrects for numerical uncertanties in cases of very
     386        // small distances...
     387        //
     388        if (dists[i]<0)
     389            dists[i]=0;
    360390    }
    361391
    … …  
    363393    TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);
    364394
    365     const Int_t n = num<rows ? num : rows;
    366 
    367     Double_t res = 0;
    368     for (int i=0; i<n; i++)
    369         res += dists[idx[i]];
    370 
    371     return sqrt(res)/n;
    372 }
    373 
     395    const Int_t n = abs(num)<rows ? abs(num) : rows;
     396
     397    if (num<0)
     398    {
     399        //
     400        // Kernel function sum
     401        //
     402        const Double_t h = 1;
     403
     404        Double_t res = 0;
     405        for (int i=0; i<n; i++)
     406            res += exp(-dists[idx[i]]/h);
     407
     408        return log(res/n);
     409    }
     410    else
     411    {
     412        //
     413        // Nearest Neighbor sum
     414        //
     415        Double_t res = 0;
     416        for (int i=0; i<n; i++)
     417            res += dists[idx[i]];
     418
     419        return sqrt(res/n);
     420    }
     421}
     422
     423// --------------------------------------------------------------------------
     424//
     425// Calls calc dist. In the case of the first call the covariance matrix
     426// fM2 is calculated.
     427//  - If n<0 it is divided by (nrows-1)/h while h is the kernel factor.
     428//
    374429Double_t MHMatrix::CalcDist(const TVector &evt, Int_t num)
    375430{
    … …  
    379434        fM2.ResizeTo(m);
    380435        fM2 = m;
     436        fM2 *= fM.GetNrows()-1;
    381437        delete &m;
    382438    }
    … …  
    391447    fM.ResizeTo(m);
    392448    fM = m;
    393 
    394     TIter Next(fRules);
    395     TObject *obj = NULL;
    396     while ((obj=Next()))
    397         fData->Add(new MDataChain(obj->GetName()));
    398 }
     449}
     450
     451void MHMatrix::StreamPrimitive(ofstream &out) const
     452{
     453    Bool_t data = fData && !TestBit(kIsOwner);
     454
     455    if (data)
     456    {
     457        fData->SavePrimitive(out);
     458        out << endl;
     459    }
     460
     461    out << "   MHMatrix " << GetUniqueName();
     462
     463    if (data || fName!=gsDefName || fTitle!=gsDefTitle)
     464    {
     465        out << "(";
     466        if (data)
     467            out << "&" << fData->GetUniqueName();
     468        if (fName!=gsDefName || fTitle!=gsDefTitle)
     469        {
     470            if (data)
     471                out << ", ";
     472            out << "\"" << fName << "\"";
     473            if (fTitle!=gsDefTitle)
     474                out << ", \"" << fTitle << "\"";
     475        }
     476    }
     477    out << ");" << endl;
     478
     479    if (fData && TestBit(kIsOwner))
     480        for (int i=0; i<fData->GetNumEntries(); i++)
     481            out << "   " << GetUniqueName() << ".AddColumn(\"" << (*fData)[i].GetRule() << "\");" << endl;
     482}

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

    r1355 r1489  
    99#endif
    1010
    11 class MDataChain;
     11class MDataArray;
    1212
    1313class MHMatrix : public MH
    1414{
    1515protected:
    16     Int_t   fNumRow; //! Number of dimensions of histogram
    17     TMatrix fM;      // Matrix to be filled
     16    Int_t   fNumRow;    //! Number of dimensions of histogram
     17    TMatrix fM;         // Matrix to be filled
    1818
    19     TMatrix fM2;     //!
     19    TMatrix fM2;        //!
    2020
    21     TList  *fData;   //! List of data members (columns)
    22     TList  *fRules;  // List of data members as text for storage
     21    MDataArray *fData;  // List of data members (columns)
     22
     23    enum {
     24        kIsOwner = BIT(14)
     25    };
    2326
    2427    void AddRow();
    … …  
    2831    Bool_t Finalize();
    2932
     33    void StreamPrimitive(ofstream &out) const;
     34
    3035public:
    3136    MHMatrix(const char *name=NULL, const char *title=NULL);
     37    MHMatrix(MDataArray *mat, const char *name=NULL, const char *title=NULL);
    3238    ~MHMatrix();
    3339
    3440    void AddColumn(const char *name);
    35     void AddColumns(const MHMatrix *matrix);
     41    void AddColumns(MDataArray *mat);
    3642
    37     //    TMatrix &GetM() { return fM; }
     43    MDataArray *GetColumns() { return fData; }
     44
    3845    const TMatrix &GetM() const { return fM; }
    39     const TList *GetRules() const { return fRules; }
    4046
    4147    //void Draw(Option_t *opt=NULL);
    … …  
    4955    Double_t CalcDist(const TVector &v, Int_t num = 25);
    5056
     57    void SetIOwner(Bool_t b=kTRUE) { b ? SetBit(kIsOwner) : ResetBit(kIsOwner); }
     58
    5159    void Reassign();
    5260
Note: See TracChangeset for help on using the changeset viewer.