source: trunk/MagicSoft/Mars/mhist/MHMatrix.cc@ 1353

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHMatrix
//
// This is a histogram container which holds a matrix with one column per
// data variable. The data variable can be a complex rule (MDataChain).
// Each event for wich Fill is called (by MFillH) is added as a new
// row to the matrix.
//
// For example:
//   MHMatrix m;
//   m.AddColumn("MHillas.fSize");
//   m.AddColumn("MMcEvt.fImpact/100");
//   m.AddColumn("HillasSource.fDist*MGeomCam.fConvMm2Deg");
//   MFillH fillm(&m);
//   taskliost.AddToList(&fillm);
//   [...]
//   m.Print();
//
/////////////////////////////////////////////////////////////////////////////
#include "MHMatrix.h"

#include <math.h>     // fabs on Alpha

#include <TList.h>
#include <TArrayD.h>
#include <TArrayI.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

#include "MDataChain.h"

ClassImp(MHMatrix);

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHMatrix::MHMatrix(const char *name, const char *title)
    : fNumRow(0)
{
    fName  = name  ? name  : "MHMatrix";
    fTitle = title ? title : "Multidimensional Matrix";

    fData = new TList;
    fData->SetOwner();

    fRules = new TList;
    fRules->SetOwner();
}

// --------------------------------------------------------------------------
//
// Destructor
//
MHMatrix::~MHMatrix()
{
    delete fData;
    delete fRules;
}

// --------------------------------------------------------------------------
//
// Add a new column to the matrix. This can only be done before the first
// event (row) was filled into the matrix. For the syntax of the rule
// see MDataChain.
//
void MHMatrix::AddColumn(const char *rule)
{
    if (fM.IsValid())
    {
        *fLog << warn << "Warning - matrix is already in use. Can't add a new column... skipped." << endl;
        return;
    }

    MDataChain &chain = *new MDataChain(rule);
    if (!chain.IsValid())
    {
        *fLog << err << "Error - Rule cannot be translated... ignored." << endl;
        delete &chain;
        return;
    }
    fData->Add(&chain);

    TNamed *name = new TNamed(rule, rule); // Fimxe, in 3.02/07 the title can't be "", why?
    fRules->Add(name);
}

void MHMatrix::AddColumns(const MHMatrix *matrix)
{
    TIter Next(matrix->fRules);
    TObject *rule=NULL;
    while ((rule=Next()))
        AddColumn(rule->GetName());
}

// --------------------------------------------------------------------------
//
// Checks whether at least one column is available and PreProcesses all
// data chains.
//
Bool_t MHMatrix::SetupFill(const MParList *plist)
{
    if (fData->GetSize()==0)
    {
        *fLog << err << "Error - No Column specified... aborting." << endl;
        return kFALSE;
    }

    TIter Next(fData);
    MData *data = NULL;
    while ((data=(MData*)Next()))
        if (!data->PreProcess(plist))
            return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// If the matrix has not enough rows double the number of available rows.
//
void MHMatrix::AddRow()
{
    fNumRow++;

    if (fM.GetNrows() > fNumRow)
        return;

    if (!fM.IsValid())
    {
        fM.ResizeTo(1, fData->GetSize());
        return;
    }

    TMatrix m(fM);

    fM.ResizeTo(fM.GetNrows()*2, fData->GetSize());

    for (int x=0; x<m.GetNcols(); x++)
        for (int y=0; y<m.GetNrows(); y++)
            fM(y, x) = m(y, x);
}

// --------------------------------------------------------------------------
//
// Add the values correspoding to the columns to the new row
//
Bool_t MHMatrix::Fill(const MParContainer *par)
{
    AddRow();

    int col=0;
    TIter Next(fData);
    MData *data = NULL;
    while ((data=(MData*)Next()))
        fM(fNumRow-1, col++) = data->GetValue();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Resize the matrix to a number of rows which corresponds to the number of
// rows which have really been filled with values.
//
Bool_t MHMatrix::Finalize()
{
    TMatrix m(fM);

    fM.ResizeTo(fNumRow, fData->GetSize());

    for (int x=0; x<fM.GetNcols(); x++)
        for (int y=0; y<fM.GetNrows(); y++)
            fM(y, x) = m(y, x);

    return kTRUE;
}
/*
// --------------------------------------------------------------------------
//
// Draw clone of histogram. So that the object can be deleted
// and the histogram is still visible in the canvas.
// The cloned object are deleted together with the canvas if the canvas is
// destroyed. If you want to handle destroying the canvas you can get a
// pointer to it from this function
//
TObject *MHMatrix::DrawClone(Option_t *opt) const
{
    TCanvas &c = *MH::MakeDefCanvas(fHist);

    //
    // This is necessary to get the expected bahviour of DrawClone
    //
    gROOT->SetSelectedPad(NULL);

    fHist->DrawCopy(opt);

    TString str(opt);
    if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileX();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }
    if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileY();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }

    c.Modified();
    c.Update();

    return &c;
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the histogram into it.
// Be careful: The histogram belongs to this object and won't get deleted
// together with the canvas.
//
void MHMatrix::Draw(Option_t *opt)
{
    if (!gPad)
        MH::MakeDefCanvas(fHist);

    fHist->Draw(opt);

    TString str(opt);
    if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileX();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }
    if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileY();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }

    gPad->Modified();
    gPad->Update();
}
*/

// --------------------------------------------------------------------------
//
// Prints the meaning of the columns and the contents of the matrix.
// Becareful, this can take a long time for matrices with many rows.
//
void MHMatrix::Print(Option_t *) const
{
    int n=0;

    TIter Next(fData->GetSize() ? fData : fRules);
    MData *data = NULL;
    while ((data=(MData*)Next()))
    {
        *fLog << all << " Column " << setw(3) << n++ << ": " << flush;
        data->Print();
        *fLog << endl;
    }

    fM.Print();
}

const TMatrix *MHMatrix::InvertPosDef()
{
    TMatrix m(fM);

    const Int_t rows = m.GetNrows();
    const Int_t cols = m.GetNcols();

    for (int i=0; i<cols; i++)
    {
        Double_t avg = 0;
        for (int j=0; j<rows; j++)
            avg += fM(j, i);

        avg /= rows;

        for (int j=0; j<rows; j++)
            m(j, i) -= avg;
    }

    TMatrix *m2 = new TMatrix(m, TMatrix::kTransposeMult, m);

    Double_t det;
    m2->Invert(&det);
    if (det==0)
    {
        *fLog << err << "ERROR - MHMatrix::InvertPosDef failed (Matrix is sigular)." << endl;
        delete m2;
        return NULL;
    }

    // m2->Print();

    return m2;
}

Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const
{
    const Int_t rows = fM.GetNrows();
    const Int_t cols = fM.GetNcols();

    TArrayD dists(rows);

    //
    // Calculate:  v^T * M * v
    //
    for (int i=0; i<rows; i++)
    {
        TVector col(cols);
        col = TMatrixRow(fM, i);

        TVector d = evt;
        d -= col;

        TVector d2 = d;
        d2 *= m;

        dists[i] = fabs(d2*d);
    }

    TArrayI idx(rows);
    TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);

    const Int_t n = num<rows ? num : rows;

    Double_t res = 0;
    for (int i=0; i<n; i++)
        res += dists[idx[i]];

    return res/n;
}

Double_t MHMatrix::CalcDist(const TVector &evt, Int_t num)
{
    if (!fM2.IsValid())
    {
        const TMatrix &m = *InvertPosDef();
        fM2.ResizeTo(m);
        fM2 = m;
        delete &m;
    }

    return CalcDist(fM2, evt, num);
}

void MHMatrix::Reassign()
{
    TMatrix m = fM;
    fM.ResizeTo(1,1);
    fM.ResizeTo(m);
    fM = m;

    TIter Next(fRules);
    TObject *obj = NULL;
    while ((obj=Next()))
        fData->Add(new MDataChain(obj->GetName()));
}