Index: trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1335) +++ trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 1336) @@ -41,4 +41,6 @@ #pragma link C++ class MHMcEnergyMigration+; +#pragma link C++ class MHCompProb+; +#pragma link C++ class MHHadroness+; #endif Index: trunk/MagicSoft/Mars/mhist/MFillH.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 1335) +++ trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 1336) @@ -243,5 +243,9 @@ if (!obj) + { + if (cls==name) + *fLog << inf << "Object '" << fHName << "' not found in parlist... trying to create one." << endl; obj = pList->FindCreateObj(cls, name); + } if (!obj) Index: trunk/MagicSoft/Mars/mhist/MH3.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MH3.cc (revision 1335) +++ trunk/MagicSoft/Mars/mhist/MH3.cc (revision 1336) @@ -319,4 +319,5 @@ p->Draw("same"); p->SetBit(kCanDelete); + p->SetDirectory(NULL); } if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2) @@ -325,4 +326,5 @@ p->Draw("same"); p->SetBit(kCanDelete); + p->SetDirectory(NULL); } @@ -352,4 +354,5 @@ p->Draw("same"); p->SetBit(kCanDelete); + p->SetDirectory(NULL); } if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2) @@ -358,4 +361,5 @@ p->Draw("same"); p->SetBit(kCanDelete); + p->SetDirectory(NULL); } Index: trunk/MagicSoft/Mars/mhist/MHCompProb.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCompProb.cc (revision 1336) +++ trunk/MagicSoft/Mars/mhist/MHCompProb.cc (revision 1336) @@ -0,0 +1,241 @@ +/* ======================================================================== *\ +! +! * +! * 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): Abelardo Moralejo +! Author(s): Thomas Bretz, 5/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////// +// +// MHCompProb +// +// This class contains different histograms of the Hillas parameters +// and composite probabilities based on them. +// +/////////////////////////////////////////////////////////////////////// + +#include "MHCompProb.h" + +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MBinning.h" +#include "MDataChain.h" + +#include "MMcEvt.hxx" + +ClassImp(MHCompProb); + +// -------------------------------------------------------------------------- +// +// Setup histograms +// +MHCompProb::MHCompProb(Int_t nbins, const char *name, const char *title) + : fNumLoop(0) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHCompProb"; + fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms"; + + fData = new TList; + fRules = new TList; + fHists = new TList; + fHistVar = new TList; + + fData->SetOwner(); + fRules->SetOwner(); + fHists->SetOwner(); + fHistVar->SetOwner(); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHCompProb::~MHCompProb() +{ + delete fData; + delete fRules; + delete fHists; + delete fHistVar; +} + +// -------------------------------------------------------------------------- +// +// +// +void MHCompProb::Add(const char *rule, Int_t n, Float_t min, Float_t max) +{ + MDataChain &chain = *new MDataChain(rule); + fData->Add(&chain); + + TNamed &name = *new TNamed(rule, ""); + fRules->Add(&name); + + TH1D &hist = *new TH1D(Form("Hist_%s", rule), rule, n, min, max); + hist.SetXTitle(rule); + hist.SetYTitle("Counts"); + hist.SetDirectory(NULL); + fHists->Add(&hist); + + TH1D &varhist = *new TH1D; + varhist.SetName(Form("Var_%s", rule)); + varhist.SetTitle(rule); + varhist.SetXTitle(rule); + varhist.SetYTitle("Counts"); + varhist.SetDirectory(NULL); + fHistVar->Add(&varhist); +} + +// -------------------------------------------------------------------------- +// +// +// +Bool_t MHCompProb::SetupFill(const MParList *plist) +{ + if (fData->GetSize()==0) + { + *fLog << err << "No data members spcified for usage... aborting." << endl; + return kFALSE; + } + + TIter Next(fData); + MData *data=NULL; + while ((data=(MData*)Next())) + if (!data->PreProcess(plist)) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +// +void MHCompProb::Fill(TList &list) +{ + MData *data = NULL; + + TIter NextD(fData); + TIter NextH(&list); + + while ((data=(MData*)NextD())) + { + TH1D *hist = (TH1D*)NextH(); + hist->Fill(data->GetValue()); + } +} + +// -------------------------------------------------------------------------- +// +// +// +Bool_t MHCompProb::Fill(const MParContainer *par) +{ + const MMcEvt &mcevt = *(MMcEvt*)par; + + switch (fNumLoop) + { + case 0: // First loop : fill the fixed-bin histograms with gammas. + if (mcevt.GetPartId() == kGAMMA) + Fill(*fHists); + return kTRUE; + + case 1: // Second Loop: fill the variable-bin histograms with protons. + if (mcevt.GetPartId() != kGAMMA) + Fill(*fHistVar); + return kTRUE; + default: + *fLog << err << "Error - Invalid Loop Number... aborting." << endl; + return kFALSE; + } +} + +// -------------------------------------------------------------------------- +// +// +// +Bool_t MHCompProb::Finalize() +{ + switch (fNumLoop++) + { + case 0: + *fLog << inf << "Finished filling fixed bin size histograms with gamma-data." << endl; + SetBinningHistVar(); + fHists->Delete(); + return kTRUE; + case 1: + *fLog << inf << "Finished filling variable bin size histogram with proton data." << endl; + return kTRUE; + default: + *fLog << err << "Error - Invalid Loop Number... aborting." << endl; + return kFALSE; + } +} + +// -------------------------------------------------------------------------- +// +// +// +void MHCompProb::SetBinningHistVar() +{ + Int_t nedges = 51; // Number of bins in variable-bin histograms. + + TIter NextH(fHists); + TIter NextV(fHistVar); + TH1D *hist = NULL; + while ((hist=(TH1D*)NextH())) + { + Int_t n = hist->GetNbinsX(); + + TArrayD edges(nedges); + + edges[0] = hist->GetBinLowEdge(1); + edges[nedges-1] = hist->GetBinLowEdge(n+1); + + Float_t newwidth = hist->Integral(1, n)/nedges; + + Int_t jbin = 1; + for (Int_t j=1; jIntegral(1, j) <= jbin*newwidth) + continue; + + edges[jbin++] = hist->GetBinLowEdge(j+1); + } + + MBinning bins; + bins.SetEdges(edges); + + SetBinning((TH1D*)NextV(), &bins); + } +} + Index: trunk/MagicSoft/Mars/mhist/MHCompProb.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCompProb.h (revision 1336) +++ trunk/MagicSoft/Mars/mhist/MHCompProb.h (revision 1336) @@ -0,0 +1,40 @@ +#ifndef MARS_MHCompProb +#define MARS_MHCompProb + +#ifndef MARS_MH +#include "MH.h" +#endif + +class MParList; + +class MHCompProb : public MH +{ +private: + Int_t fNumLoop; //! Counter of the loop (two eventloops needed) + + TList *fRules; // Rules describing the used data sets + + TList *fData; //! MDataChain objects + TList *fHists; //! fixed bin size histograms + TList *fHistVar; // variable bin size histograms + + void Fill(TList &list); + void SetBinningHistVar(); + +public: + MHCompProb(Int_t nbins, const char *name=NULL, const char *title=NULL); + ~MHCompProb(); + + void Add(const char *rule, Int_t n, Float_t min, Float_t max); + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par); + Bool_t Finalize(); + + const TList *GetRules() const { return fRules; } + TList *GetHistVar() const { return fHistVar; } + + ClassDef(MHCompProb, 1) // Histogram to be used for the calculation of the composite probabilities +}; + +#endif Index: trunk/MagicSoft/Mars/mhist/MHHadroness.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHHadroness.cc (revision 1336) +++ trunk/MagicSoft/Mars/mhist/MHHadroness.cc (revision 1336) @@ -0,0 +1,375 @@ +/* ======================================================================== *\ +! +! * +! * 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, 5/2002 +! Author(s): Abelardo Moralejo +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHHadroness +// +// This is histogram is a way to evaluate the quality of a gamma/hadron +// seperation method. It is filled from a MHadroness container, which +// stores a hadroness for the current event. The Value must be in the +// range [0,1]. To fill the histograms correctly the information +// whether it is a gamma or hadron (not a gamma) must be available from +// a MMcEvt container. +// +// In the constructor you can change the number of used bns for the +// evaluation. +// +// The meaning of the histograms (Draw, DrawClone) are the following: +// * Upper Left Corner: +// - black: histogram of all hadronesses for gammas +// - red: histogram of all hadronesses for non gammas +// * Upper Right Corner: +// - black: acceptance of gammas vs. the hadroness +// - red: acceptance of non gammas vs. the hadroness +// - blue: 2D distance of (acceptance_hadrons, acceptances_gammas) +// to optimum (0, 1) +// 1-sqrt(Ag*Ag + (1-Ah)*(1-Ah)) +// * Bottom Left Corner: +// Naive quality factor: Ag/sqrt(Ah) +// * Bottom Right Corner: +// - black: Acceprtance Gammas vs. Acceptance Hadrons +// - blue cross: minimum of distance to (0, 1) +// +//////////////////////////////////////////////////////////////////////////// +#include "MHHadroness.h" + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MBinning.h" +#include "MHadroness.h" + +#include "MMcEvt.hxx" + +ClassImp(MHHadroness); + +// -------------------------------------------------------------------------- +// +// Setup histograms, nbins is the number of bins used for the evaluation. +// The default is 100 bins. +// +MHHadroness::MHHadroness(Int_t nbins, const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHHadroness"; + fTitle = title ? title : "Gamma/Hadron Separation Quality Histograms"; + + fGraph = new TGraph; + fGraph->SetTitle("Acceptance Gammas vs. Hadrons"); + fGraph->SetMaximum(1); + + fGhness = new TH1D("Ghness", "Hadronness", nbins, 0, 1); + fPhness = new TH1D("Phness", "Hadronness", nbins, 0, 1); + fGhness->SetXTitle("Hadroness"); + fPhness->SetXTitle("Hadroness"); + fGhness->SetYTitle("Counts"); + fPhness->SetYTitle("Counts"); + + fIntGhness = new TH1D("AccGammas", "Acceptance", nbins, 0, 1); + fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1); + fIntGhness->SetXTitle("Hadroness"); + fIntPhness->SetXTitle("Hadroness"); + fIntGhness->SetYTitle("Acceptance"); + fIntPhness->SetYTitle("Acceptance"); + + fQfac = new TH1D("Qfac", "Naive Quality factor", nbins, 0, 1); + fQfac->SetXTitle("Hadroness"); + fQfac->SetYTitle("Quality"); + + fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1); + fMinDist->SetXTitle("Hadroness"); + fMinDist->SetYTitle("Distance"); + + fQfac->SetDirectory(NULL); + fGhness->SetDirectory(NULL); + fPhness->SetDirectory(NULL); + fMinDist->SetDirectory(NULL); + fIntGhness->SetDirectory(NULL); + fIntPhness->SetDirectory(NULL); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms. +// +MHHadroness::~MHHadroness() +{ + delete fGhness; + delete fIntGhness; + delete fPhness; + delete fIntPhness; + delete fQfac; + delete fMinDist; + delete fGraph; +} + +// -------------------------------------------------------------------------- +// +// Setup Filling of the histograms. It needs: +// MMcEvt and MHadroness +// +Bool_t MHHadroness::SetupFill(const MParList *plist) +{ + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fHadroness = (MHadroness*)plist->FindObject("MHadroness"); + if (!fHadroness) + { + *fLog << err << dbginf << "MHadroness not found... aborting." << endl; + return kFALSE; + } + + /* + MBinning* bins = (MBinning*)plist->FindObject("BinningHadroness"); + if (!bins) + { + *fLog << err << dbginf << "BinningHadroness [MBinning] not found... aborting." << endl; + return kFALSE; + } + + SetBinning(&fHist, binsalpha, binsenergy, binstheta); + + fHist.Sumw2(); + */ + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the Hadroness from a MHadroness container into the corresponding +// histogram dependant on the particle id. +// +Bool_t MHHadroness::Fill(const MParContainer *par) +{ + if (fMcEvt->GetPartId()==kGAMMA) + fGhness->Fill(fHadroness->GetHadroness()); + else + fPhness->Fill(fHadroness->GetHadroness()); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Finalize the histograms: +// - integrate the hadroness histograms --> acceptance +// - fill the Minimum Distance histogram (formular see class description) +// - fill the Quality histogram (formular see class description) +// +Bool_t MHHadroness::Finalize() +{ + *fLog << inf << "Finished filling hadroness histograms." << endl; + + Int_t n = fGhness->GetNbinsX(); + + fGraph->Set(n); + + const Stat_t sumg = fGhness->Integral(1, n+1); + const Stat_t sump = fPhness->Integral(1, n+1); + + for (Int_t i=1; i<=n; i++) + { + const Stat_t ip = fPhness->Integral(1, i)/sump; + const Stat_t ig = fGhness->Integral(1, i)/sumg; + + fIntPhness->SetBinContent(i, ip); + fIntGhness->SetBinContent(i, ig); + + fGraph->SetPoint(i, ip, ig); + + if (ip<=0) + continue; + + fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1))); + fQfac->SetBinContent(i, ig/sqrt(ip)); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Search the corresponding points for the given hadron acceptance (acchad) +// and interpolate the to pointsd (linear) +// +Double_t MHHadroness::GetGammaAcceptance(Double_t acchad) const +{ + const Int_t n = fGraph->GetN(); + const Double_t *x = fGraph->GetX(); + const Double_t *y = fGraph->GetY(); + + Int_t i = 0; + while (iGetMaximumBin(); + *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " at H=" << fMinDist->GetBinCenter(h) << endl; + *fLog << " Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, "; + *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl; + *fLog << "Maximum Q-Factor: " << GetQfac()->GetMaximum() << endl; + + *fLog << endl; +} + +// -------------------------------------------------------------------------- +// +// Draw clone of all histograms. (For the Meaning see class description) +// +TObject *MHHadroness::DrawClone(Option_t *opt) const +{ + TCanvas &c = *MakeDefCanvas("Hadroness", fTitle); + c.Divide(2, 2); + + gROOT->SetSelectedPad(NULL); + + c.cd(1); + gStyle->SetOptStat(10); + Getghness()->DrawCopy(); + Getphness()->SetLineColor(2); + Getphness()->DrawCopy("same"); + c.cd(4); + TGraph &g = (TGraph&)*fGraph->DrawClone("AC"); + g.SetBit(kCanDelete); + gPad->Modified(); + gPad->Update(); + if (g.GetHistogram()) + { + g.GetXaxis()->SetRangeUser(0, 1); + g.GetXaxis()->SetTitle("Acceptance Hadrons"); + g.GetYaxis()->SetTitle("Acceptance Gammas"); + gPad->Modified(); + gPad->Update(); + } + const Int_t h = fMinDist->GetMaximumBin(); + TMarker *m = new TMarker(fIntPhness->GetBinContent(h), + fIntGhness->GetBinContent(h), kStar); + m->SetMarkerColor(4); + m->SetBit(kCanDelete); + m->Draw(); + c.cd(2); + gStyle->SetOptStat(0); + Getighness()->DrawCopy(); + Getiphness()->SetLineColor(2); + Getiphness()->DrawCopy("same"); + fMinDist->SetLineColor(4); + fMinDist->DrawCopy("same"); + c.cd(3); + GetQfac()->DrawCopy(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw all histograms. (For the Meaning see class description) +// +void MHHadroness::Draw(Option_t *) +{ + if (!gPad) + MakeDefCanvas("Hadroness", fTitle); + + gPad->Divide(2, 2); + + gPad->cd(1); + gStyle->SetOptStat(10); + Getghness()->Draw(); + Getphness()->SetLineColor(2); + Getphness()->Draw("same"); + gPad->cd(4); + fGraph->Draw("AC"); + gPad->Modified(); + gPad->Update(); + if (fGraph->GetHistogram()) + { + fGraph->GetXaxis()->SetRangeUser(0, 1); + fGraph->GetXaxis()->SetTitle("Acceptance Hadrons"); + fGraph->GetYaxis()->SetTitle("Acceptance Gammas"); + gPad->Modified(); + gPad->Update(); + } + const Int_t h = fMinDist->GetMaximumBin(); + TMarker *m = new TMarker(fIntPhness->GetBinContent(h), + fIntGhness->GetBinContent(h), kStar); + m->SetMarkerColor(4); + m->SetBit(kCanDelete); + m->Draw(); + gPad->cd(2); + gStyle->SetOptStat(0); + Getighness()->Draw(); + Getiphness()->SetLineColor(2); + Getiphness()->Draw("same"); + fMinDist->SetLineColor(4); + fMinDist->Draw("same"); + gPad->cd(3); + GetQfac()->Draw(); +} Index: trunk/MagicSoft/Mars/mhist/MHHadroness.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHHadroness.h (revision 1336) +++ trunk/MagicSoft/Mars/mhist/MHHadroness.h (revision 1336) @@ -0,0 +1,54 @@ +#ifndef MARS_MHHadroness +#define MARS_MHHadroness + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TGraph; +class MParList; +class MMcEvt; +class MHadroness; + +class MHHadroness : public MH +{ +private: + const MMcEvt *fMcEvt; //! + const MHadroness *fHadroness; //! + + TH1D* fPhness; // Hadrons Hadroness + TH1D* fGhness; // Gammas Hadroness + TH1D* fIntPhness; // Hadrons Acceptance + TH1D* fIntGhness; // Gammas Acceptance + TH1D* fQfac; // Quality factor + TH1D* fMinDist; // Minimum Distance to optimum acceptance + + TGraph *fGraph; // gamma acceptance vs. hadron acceptance + +public: + MHHadroness(Int_t nbins=100, const char *name=NULL, const char *title=NULL); + ~MHHadroness(); + + Double_t GetGammaAcceptance(Double_t acchad) const; + + TH1D *GetQfac() const { return fQfac; } + TH1D *Getghness() const { return fGhness; } + TH1D *Getphness() const { return fPhness; } + TH1D *Getighness() const { return fIntGhness; } + TH1D *Getiphness() const { return fIntPhness; } + //TH2D *GetHist() const { return fHist; } + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par); + Bool_t Finalize(); + + void Print(Option_t *option="") const; + + void Draw(Option_t *opt=""); + TObject *DrawClone(Option_t *opt="") const; + + ClassDef(MHHadroness, 1) // Gamma/Hadron Separation Quality Histograms +}; + +#endif Index: trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1335) +++ trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1336) @@ -43,8 +43,9 @@ // ///////////////////////////////////////////////////////////////////////////// - #include "MHMatrix.h" #include +#include +#include #include "MLog.h" @@ -61,9 +62,15 @@ // Default Constructor // -MHMatrix::MHMatrix() +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(); } @@ -75,4 +82,5 @@ { delete fData; + delete fRules; } @@ -91,5 +99,5 @@ } - MDataChain &chain = *new MDataChain(name); + MDataChain &chain = *new MDataChain(rule); if (!chain.IsValid()) { @@ -100,4 +108,15 @@ fData->Add(&chain); + + TNamed *name = new TNamed(rule, ""); + fRules->Add(name); +} + +void MHMatrix::AddColumns(const MHMatrix *matrix) +{ + TIter Next(matrix->fRules); + TObject *rule=NULL; + while ((rule=Next())) + AddColumn(rule->GetName()); } @@ -264,9 +283,10 @@ { int n=0; - TIter Next(fData); + + TIter Next(fData->GetSize() ? fData : fRules); MData *data = NULL; while ((data=(MData*)Next())) { - *fLog << all << " Column " << setw(3) << n++ << ": "; + *fLog << all << " Column " << setw(3) << n++ << ": " << flush; data->Print(); *fLog << endl; @@ -275,2 +295,121 @@ fM.Print(); } + +const TMatrix *MHMatrix::InvertPosDef() +{ + /* + ---------------------------------- + Substract Mean of Rows from Rows + ---------------------------------- + + const Int_t rows = fM.GetNrows(); + const Int_t cols = fM.GetNcols(); + + for (int i=0; iInvert(&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; iAdd(new MDataChain(obj->GetName())); +} Index: trunk/MagicSoft/Mars/mhist/MHMatrix.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1335) +++ trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1336) @@ -14,8 +14,11 @@ { protected: - Int_t fNumRow; // Number of dimensions of histogram + Int_t fNumRow; //! Number of dimensions of histogram TMatrix fM; // Matrix to be filled - TList *fData; // List of data members (columns) + TMatrix fM2; //! + + TList *fData; //! List of data members (columns) + TList *fRules; // List of data members as text for storage void AddRow(); @@ -26,11 +29,13 @@ public: - MHMatrix(); + MHMatrix(const char *name=NULL, const char *title=NULL); ~MHMatrix(); void AddColumn(const char *name); + void AddColumns(const MHMatrix *matrix); - TMatrix &GetM() { return fM; } + // TMatrix &GetM() { return fM; } const TMatrix &GetM() const { return fM; } + const TList *GetRules() const { return fRules; } //void Draw(Option_t *opt=NULL); @@ -39,4 +44,11 @@ void Print(Option_t *) const; + const TMatrix *InvertPosDef(); + + Double_t CalcDist(const TMatrix &m, const TVector &v, Int_t num = 25) const; + Double_t CalcDist(const TVector &v, Int_t num = 25); + + void Reassign(); + ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events }; Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 1335) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 1336) @@ -48,4 +48,6 @@ MHTimeDiffTime.cc \ MHTimeDiffTheta.cc \ + MHCompProb.cc \ + MHHadroness.cc \ MHMcEnergy.cc \ MHMcEfficiency.cc \ @@ -53,4 +55,5 @@ MHMcEfficiencyEnergy.cc \ MHMcEnergyImpact.cc \ + MHMcEnergyMigration.cc \ MHThetabarTime.cc \ MHMcRate.cc \