Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1488) +++ /trunk/MagicSoft/Mars/Changelog (revision 1489) @@ -1,3 +1,56 @@ -*-*- END -*-*- + + 2002/08/08: Thomas Bretz + + * manalysis/MHillasSrc.cc: + - use double dist instead of single fDist for calculation + + * manalysis/MMultiDimDistCalc.[h,cc]: + - added support for the kernel method + - added stream primitive + - changed version number to 1 + - adapted to new MHMatrix (using MDataArray) + + * mdata/MDataArray.[h,cc]: + - added + + * mdata/DataLinkDef.h, madata/Makefile: + - added MDataArray + + * mfileio/MWriteRootFile.cc: + - fixed some bugs in StreamPrimitive + - StreamPrimtive doesn't write the default name/title anymore + + * mhist/MHMatrix.[h,cc]: + - replaced the Arrays for the rules by a MDataArray + - implemented StreamPrimitive + - implement the use of the kernel function for num<0 + - multiply fM2 by nevts-1 + - added sanity check in case of dists[i]<0 + + * mhist/MHHillas.[h,cc]: + - added fUsedPix, fCorePix + - added fUsedPix, fCorePix to plots + - changed layout of plots + - changed name and title of MakeDefCanvas + + * mhist/MHHillasSrc.[h,cc]: + - changed plot of Alpha from fabs(fAlpha) to fAlpha + - changed name and title of MakeDefCanvas + + * mhist/MHillasExt.[h,cc]: + - changed layout of plots + - changed name and title of MakeDefCanvas + + * mbase/MTask.cc: + - fixed wrong streaming of filter name + + * macros/starplot.C: + - added + + * macros/dohtml.C: + - added starplot.C + + 2002/08/07: Thomas Bretz Index: /trunk/MagicSoft/Mars/macros/dohtml.C =================================================================== --- /trunk/MagicSoft/Mars/macros/dohtml.C (revision 1488) +++ /trunk/MagicSoft/Mars/macros/dohtml.C (revision 1489) @@ -52,4 +52,5 @@ html.Convert("trigrate.C", "MARS - Calculate Trigger Rate from a MC root file"); html.Convert("star.C", "MARS - (St)andard (A)nalysis and (R)econstruction"); + html.Convert("starplot.C", "MARS - Plot parameters from file created with star.C"); html.Convert("comprob.C", "MARS - Calculation of composite probabilities for G/H-Seperation"); html.Convert("multidimdist.C", "MARS - Calculation of multidimensional distances for G/H-Seperation"); Index: /trunk/MagicSoft/Mars/macros/starplot.C =================================================================== --- /trunk/MagicSoft/Mars/macros/starplot.C (revision 1489) +++ /trunk/MagicSoft/Mars/macros/starplot.C (revision 1489) @@ -0,0 +1,114 @@ +/* ======================================================================== *\ +! +! * +! * 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 et al, 08/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + + +void starplot(const char *filename="Gamma_*.root") +{ + // + // This is a demonstration program which plots the Hillas + // parameter from a file created with star.C + // + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + + MTaskList tlist; + plist.AddToList(&tlist); + + // + // The geometry container must be created by yourself to make sure + // that you don't choos a wrong geometry by chance + // + MGeomCamMagic geomcam; + plist.AddToList(&geomcam); + + // + // Use this if you want to change the binning of one of + // the histograms. You can use: + // BinningConc, BinningConc1, BinningAsym, BinningM3Long, + // BinningM3Trans, BinningWidth, BinningLength, BinningDist, + // BinningHeadTail, BinningAlpha, BinningSize, BinningDelta, + // BinningPixels and BinningCamera + // + // For more information see MBinning and the corresponding + // histograms + // + // MBinning binsalpha("BinningAlpha"); + // binsalpha.SetEdges(90, 0, 90); // 90 bins from 0 to 90 deg + // plist.AddToList(&binsalpha); + + // MBinning binssize("BinningSize"); + // binssize.SetEdgesLog(50, 1, 1e7); + // plist.AddToList(&binssize); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + // The first argument is the tree you want to read. + // Events: Cosmic ray events + // PedEvents: Pedestal Events + // CalEvents: Calibration Events + // + MReadMarsFile read("Events", filename); + read.DisableAutoScheme(); + + MFillH hfill("MHHillasExt", "MHillas"); + MFillH sfill("MHStarMap", "MHillas"); + MFillH hfill2s("HistSource [MHHillasSrc]", "HillasSource"); + MFillH hfill2a("HistAntiSrc [MHHillasSrc]", "HillasAntiSrc"); + + tlist.AddToList(&read); + tlist.AddToList(&hfill); + tlist.AddToList(&sfill); + tlist.AddToList(&hfill2s); + tlist.AddToList(&hfill2a); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // + // After the analysis is finished we can display the histograms + // + plist.FindObject("HistSource")->DrawClone(); + plist.FindObject("MHHillasExt")->DrawClone(); + plist.FindObject("HistAntiSrc")->DrawClone(); + plist.FindObject("MHStarMap")->DrawClone(); +} + Index: /trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc (revision 1488) +++ /trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc (revision 1489) @@ -69,5 +69,5 @@ { fName = name ? name : "MHillasSrc"; - fTitle = title ? title : "parameters depending in source position"; + fTitle = title ? title : "Parameters depending in source position"; } Index: /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc (revision 1488) +++ /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc (revision 1489) @@ -52,4 +52,5 @@ #include "MParList.h" #include "MDataChain.h" +#include "MDataArray.h" #include "MHadroness.h" @@ -65,5 +66,5 @@ // MMultiDimDistCalc::MMultiDimDistCalc(const char *name, const char *title) - : fNum(0), fUseKernel(kFALSE) + : fNum(0), fUseKernel(kFALSE), fData(NULL) { // @@ -73,6 +74,8 @@ fTitle = title ? title : gsDefTitle.Data(); - fData = new TList; - fData->SetOwner(); + /* + fData = new TList; + fData->SetOwner(); + */ } @@ -83,5 +86,5 @@ MMultiDimDistCalc::~MMultiDimDistCalc() { - delete fData; + // delete fData; } @@ -119,4 +122,5 @@ } + /* TIter Next(fMGammas->GetRules()); TObject *data=NULL; @@ -130,4 +134,17 @@ } fData->Add(chain); + } + */ + fData = fMGammas->GetColumns(); + if (!fData) + { + *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl; + return kFALSE; + } + + if (!fData->PreProcess(plist)) + { + *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl; + return kFALSE; } @@ -156,4 +173,8 @@ TVector event(ncols); + for (int i=0; iGetNumEntries(); i++) + event(i) = (*fData)(i); + + /* Int_t n=0; TIter Next(fData); @@ -161,4 +182,5 @@ while ((data=(MData*)Next())) event(n++) = data->GetValue(); + */ Double_t numg = fNum; Index: /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h (revision 1488) +++ /trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.h (revision 1489) @@ -9,4 +9,5 @@ class MParList; class MHadroness; +class MDataArray; class MMultiDimDistCalc : public MTask @@ -21,5 +22,5 @@ MHadroness *fHadroness; //! Output container for calculated hadroness - TList *fData; //! Used to store the MDataChains to get the event values + MDataArray *fData; //! Used to store the MDataChains to get the event values void StreamPrimitive(ofstream &out) const; Index: /trunk/MagicSoft/Mars/mbase/MTask.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTask.cc (revision 1488) +++ /trunk/MagicSoft/Mars/mbase/MTask.cc (revision 1489) @@ -281,5 +281,5 @@ fFilter->SavePrimitive(out); - out << " " << ToLower(fName) << ".SetFilter(&" << ToLower(fFilter->GetName()) <<");" << endl; */ -} + out << " " << GetUniqueName() << ".SetFilter(&" << fFilter->GetUniqueName() <<");" << endl; +} Index: /trunk/MagicSoft/Mars/mhist/MHHillas.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillas.cc (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHHillas.cc (revision 1489) @@ -67,8 +67,10 @@ fTitle = title ? title : "Source independent image parameters"; - fLength = new TH1F("Length", "Length of Ellipse", 100, 0, 296.7); - fWidth = new TH1F("Width", "Width of Ellipse", 100, 0, 296.7); - fDistC = new TH1F("DistC", "Distance from center of camera", 100, 0, 445); - fDelta = new TH1F("Delta", "Angle (Main axis - x-axis)", 101, -90, 90); + fLength = new TH1F("Length", "Length of Ellipse", 100, 0, 296.7); + fWidth = new TH1F("Width", "Width of Ellipse", 100, 0, 296.7); + fDistC = new TH1F("DistC", "Distance from center of camera", 100, 0, 445); + fDelta = new TH1F("Delta", "Angle (Main axis - x-axis)", 101, -90, 90); + fUsedPix = new TH1F("UsedPix", "Number of used pixels", 150, 0, 150); + fCorePix = new TH1F("CorePix", "Number of core pixels", 150, 0, 150); fLength->SetDirectory(NULL); @@ -76,4 +78,6 @@ fDistC->SetDirectory(NULL); fDelta->SetDirectory(NULL); + fUsedPix->SetDirectory(NULL); + fCorePix->SetDirectory(NULL); fLength->SetXTitle("Length [mm]"); @@ -81,4 +85,6 @@ fDistC->SetXTitle("Distance [mm]"); fDelta->SetXTitle("Delta [\\circ]"); + fUsedPix->SetXTitle("Number of Pixels"); + fCorePix->SetXTitle("Number of Pixels"); fLength->SetYTitle("Counts"); @@ -86,4 +92,6 @@ fDistC->SetYTitle("Counts"); fDelta->SetYTitle("Counts"); + fUsedPix->SetYTitle("Counts"); + fCorePix->SetYTitle("Counts"); MBinning bins; @@ -122,4 +130,7 @@ delete fSize; delete fCenter; + + delete fUsedPix; + delete fCorePix; } @@ -150,4 +161,6 @@ ApplyBinning(*plist, "Delta", fDelta); ApplyBinning(*plist, "Size", fSize); + ApplyBinning(*plist, "Pixels", fUsedPix); + ApplyBinning(*plist, "Pixels", fCorePix); const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera"); @@ -235,10 +248,12 @@ const Double_t scale = fUseMmScale ? 1 : fMm2Deg; - fLength->Fill(scale*h.GetLength()); - fWidth ->Fill(scale*h.GetWidth()); - fDistC ->Fill(scale*d); - fCenter->Fill(scale*h.GetMeanX(), scale*h.GetMeanY()); - fDelta ->Fill(kRad2Deg*h.GetDelta()); - fSize ->Fill(h.GetSize()); + fLength ->Fill(scale*h.GetLength()); + fWidth ->Fill(scale*h.GetWidth()); + fDistC ->Fill(scale*d); + fCenter ->Fill(scale*h.GetMeanX(), scale*h.GetMeanY()); + fDelta ->Fill(kRad2Deg*h.GetDelta()); + fSize ->Fill(h.GetSize()); + fUsedPix->Fill(h.GetNumUsedPixels()); + fCorePix->Fill(h.GetNumCorePixels()); return kTRUE; @@ -269,5 +284,5 @@ TObject *MHHillas::DrawClone(Option_t *opt) const { - TCanvas *c = MakeDefCanvas("Hillas", fTitle, 720, 810); + TCanvas *c = MakeDefCanvas(this, 720, 810); c->Divide(2,3); @@ -275,12 +290,17 @@ c->cd(1); - fLength->DrawCopy(); + fWidth->DrawCopy(); + fLength->SetLineColor(kBlue); + fLength->DrawCopy("same"); c->cd(2); - fWidth->DrawCopy(); - - c->cd(3); gPad->SetLogx(); fSize->DrawCopy(); + + c->cd(3); + fCorePix->SetLineColor(kRed); + fCorePix->DrawCopy(); + fUsedPix->SetLineColor(kGreen); + fUsedPix->DrawCopy("same"); c->cd(4); @@ -309,17 +329,22 @@ { if (!gPad) - MakeDefCanvas("Hillas", fTitle, 720, 810); + MakeDefCanvas(this, 720, 810); gPad->Divide(2,3); gPad->cd(1); - fLength->Draw(); + fWidth->Draw(); + fLength->SetLineColor(kBlue); + fLength->Draw("same"); gPad->cd(2); - fWidth->Draw(); - - gPad->cd(3); gPad->SetLogx(); fSize->Draw(); + + gPad->cd(3); + fCorePix->SetLineColor(kRed); + fCorePix->Draw(); + fUsedPix->SetLineColor(kGreen); + fUsedPix->Draw("same"); gPad->cd(4); Index: /trunk/MagicSoft/Mars/mhist/MHHillas.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillas.h (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHHillas.h (revision 1489) @@ -22,4 +22,7 @@ TH1F *fSize; //-> TH2F *fCenter; //-> + + TH1F *fUsedPix; //-> + TH1F *fCorePix; //-> void SetColors() const; Index: /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc (revision 1489) @@ -62,5 +62,5 @@ // fName = name ? name : "MHHillasExt"; - fTitle = title ? title : "Container for Hillas (ext) histograms"; + fTitle = title ? title : "Container for extended Hillas histograms"; // @@ -69,4 +69,7 @@ // connect all the histogram with the container fHist // + fHConc1.SetLineColor(kBlue); + fHConc.SetFillStyle(0); + fHConc.SetDirectory(NULL); fHConc1.SetDirectory(NULL); @@ -220,27 +223,22 @@ TObject *MHHillasExt::DrawClone(Option_t *opt) const { - TCanvas &c = *MakeDefCanvas("Hillas", "Histograms of Hillas Parameters", - 720, 810); - c.Divide(2, 3); + TCanvas &c = *MakeDefCanvas(this, 720, 540); + c.Divide(2, 2); gROOT->SetSelectedPad(NULL); - // - // This is necessary to get the expected bahviour of DrawClone - // + c.cd(1); - ((TH1F&)fHConc).DrawCopy(); + ((TH1F)fHConc1).DrawCopy(); + ((TH1F)fHConc).DrawCopy("same"); c.cd(2); - ((TH1F&)fHConc1).DrawCopy(); - - c.cd(5); - ((TH1F&)fHAsym).DrawCopy(); + ((TH1F)fHAsym).DrawCopy(); c.cd(3); - ((TH1F&)fHM3Long).DrawCopy(); + ((TH1F)fHM3Long).DrawCopy(); c.cd(4); - ((TH1F&)fHM3Trans).DrawCopy(); + ((TH1F)fHM3Trans).DrawCopy(); c.Modified(); @@ -261,23 +259,21 @@ { if (!gPad) - MakeDefCanvas("Hillas", "Histograms of Hillas Parameters", - 720, 810); - - gPad->Divide(2, 3); + MakeDefCanvas(this, 720, 540); + + gPad->Divide(2, 2); gPad->cd(1); - fHConc.DrawCopy(); + fHConc1.SetLineColor(kBlue); + fHConc1.Draw(); + fHConc.Draw("same"); gPad->cd(2); - fHConc1.DrawCopy(); - - gPad->cd(5); - fHAsym.DrawCopy(); + fHAsym.Draw(); gPad->cd(3); - fHM3Long.DrawCopy(); + fHM3Long.Draw(); gPad->cd(4); - fHM3Trans.DrawCopy(); + fHM3Trans.Draw(); gPad->Modified(); Index: /trunk/MagicSoft/Mars/mhist/MHHillasSrc.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillasSrc.cc (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHHillasSrc.cc (revision 1489) @@ -68,5 +68,5 @@ // connect all the histogram with the container fHist // - fAlpha = new TH1F("Alpha", "Alpha of Ellipse", 90, 0, 90); + fAlpha = new TH1F("Alpha", "Alpha of Ellipse", 181, -90, 90); fDist = new TH1F("Dist", "Dist of Ellipse", 100, 0, 445); fHeadTail = new TH1F("HeadTail", "HeadTail of Ellipse", 101, -445, 445); @@ -81,5 +81,5 @@ fDist->SetXTitle("Dist [mm]"); fHeadTail->SetXTitle("Head-Tail [mm]"); - fCosDA->SetXTitle("cos(\\delta,\\alpha) [mm]"); + fCosDA->SetXTitle("cos(\\delta,\\alpha)"); fAlpha->SetYTitle("Counts"); @@ -138,5 +138,5 @@ const MHillasSrc &h = *(MHillasSrc*)par; - fAlpha ->Fill(fabs(h.GetAlpha())); + fAlpha ->Fill(h.GetAlpha()); fDist ->Fill(fUseMmScale ? h.GetDist() : fMm2Deg*h.GetDist()); fHeadTail->Fill(fUseMmScale ? h.GetHeadTail() : fMm2Deg*h.GetHeadTail()); @@ -211,6 +211,5 @@ TObject *MHHillasSrc::DrawClone(Option_t *opt) const { - TCanvas *c = MakeDefCanvas("HillasSrc", "Histograms of Source dependant Parameters", - 700, 500); + TCanvas *c = MakeDefCanvas(this, 700, 500); c->Divide(2, 2); @@ -250,6 +249,5 @@ { if (!gPad) - MakeDefCanvas("HillasSrc", "Histograms of Src dependant Parameters", - 700, 500); + MakeDefCanvas(this, 700, 500); // FIXME: Display Source position Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1489) @@ -45,4 +45,6 @@ #include "MHMatrix.h" +#include + #include #include @@ -54,33 +56,44 @@ #include "MParList.h" -#include "MDataChain.h" +#include "MData.h" +#include "MDataArray.h" ClassImp(MHMatrix); -// -------------------------------------------------------------------------- -// -// Default Constructor +static const TString gsDefName = "MHMatrix"; +static const TString gsDefTitle = "Multidimensional Matrix"; + +// -------------------------------------------------------------------------- +// +// 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 + : fNumRow(0), fData(NULL) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +// -------------------------------------------------------------------------- +// +// Constructor. Initializes the columns of the matrix with the entries +// from a MDataArray +// +MHMatrix::MHMatrix(MDataArray *mat, const char *name, const char *title) + : fNumRow(0), fData(mat) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +// -------------------------------------------------------------------------- +// +// Destructor. Does not deleted a user given MDataArray, except IsOwner +// was called. // MHMatrix::~MHMatrix() { - delete fData; - delete fRules; + if (TestBit(kIsOwner) && fData) + delete fData; } @@ -99,23 +112,31 @@ } - MDataChain &chain = *new MDataChain(rule); - if (!chain.IsValid()) - { - *fLog << err << "Error - Rule cannot be translated... ignored." << endl; - delete &chain; + if (!fData) + { + fData = new MDataArray; + SetBit(kIsOwner); + } + + fData->AddEntry(rule); +} + +void MHMatrix::AddColumns(MDataArray *matrix) +{ + if (fM.IsValid()) + { + *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl; 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()); + + if (fData) + *fLog << warn << "Warning - columns already added... replacing." << endl; + + if (fData && TestBit(kIsOwner)) + { + delete fData; + ResetBit(kIsOwner); + } + + fData = matrix; } @@ -127,17 +148,11 @@ Bool_t MHMatrix::SetupFill(const MParList *plist) { - if (fData->GetSize()==0) - { - *fLog << err << "Error - No Column specified... aborting." << endl; + if (!fData) + { + *fLog << err << "Error - No Columns initialized... aborting." << endl; return kFALSE; } - TIter Next(fData); - MData *data = NULL; - while ((data=(MData*)Next())) - if (!data->PreProcess(plist)) - return kFALSE; - - return kTRUE; + return fData->PreProcess(plist); } @@ -155,5 +170,5 @@ if (!fM.IsValid()) { - fM.ResizeTo(1, fData->GetSize()); + fM.ResizeTo(1, fData->GetNumEntries()); return; } @@ -161,5 +176,5 @@ TMatrix m(fM); - fM.ResizeTo(fM.GetNrows()*2, fData->GetSize()); + fM.ResizeTo(fM.GetNrows()*2, fData->GetNumEntries()); for (int x=0; xGetValue(); + for (int col=0; colGetNumEntries(); col++) + fM(fNumRow-1, col) = (*fData)(col); return kTRUE; @@ -195,10 +207,10 @@ // It's not a fatal error so we don't need to stop PostProcessing... // - if (fData->GetSize()<1 || fNumRow<1) + if (fData->GetNumEntries()==0 || fNumRow<1) return kTRUE; TMatrix m(fM); - fM.ResizeTo(fNumRow, fData->GetSize()); + fM.ResizeTo(fNumRow, fData->GetNumEntries()); for (int x=0; xGetSize() ? fData : fRules); - MData *data = NULL; - while ((data=(MData*)Next())) - { - *fLog << all << " Column " << setw(3) << n++ << ": " << flush; - data->Print(); - *fLog << endl; - } + if (fData) + fData->Print(); + else + *fLog << all << "Sorry, no column information available." << endl; fM.Print(); @@ -336,6 +342,23 @@ } +// -------------------------------------------------------------------------- +// +// Calculated the distance of vector evt from the reference sample +// represented by the covariance metrix m. +// - If n<0 the kernel method is applied and +// sum(epx(-d)/n is returned. +// - For n>=0 the n nearest neighbors are summed and +// sqrt(sum(d))/n is returned. +// - if n==0 all distances are summed +// Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const { + if (num==0) // may later be used for another method + { + TVector d = evt; + d *= m; + return sqrt(d*evt); + } + const Int_t rows = fM.GetNrows(); const Int_t cols = fM.GetNcols(); @@ -358,4 +381,11 @@ dists[i] = d2*d; // square of distance + + // + // This corrects for numerical uncertanties in cases of very + // small distances... + // + if (dists[i]<0) + dists[i]=0; } @@ -363,13 +393,38 @@ TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE); - const Int_t n = numAdd(new MDataChain(obj->GetName())); -} +} + +void MHMatrix::StreamPrimitive(ofstream &out) const +{ + Bool_t data = fData && !TestBit(kIsOwner); + + if (data) + { + fData->SavePrimitive(out); + out << endl; + } + + out << " MHMatrix " << GetUniqueName(); + + if (data || fName!=gsDefName || fTitle!=gsDefTitle) + { + out << "("; + if (data) + out << "&" << fData->GetUniqueName(); + if (fName!=gsDefName || fTitle!=gsDefTitle) + { + if (data) + out << ", "; + out << "\"" << fName << "\""; + if (fTitle!=gsDefTitle) + out << ", \"" << fTitle << "\""; + } + } + out << ");" << endl; + + if (fData && TestBit(kIsOwner)) + for (int i=0; iGetNumEntries(); i++) + out << " " << GetUniqueName() << ".AddColumn(\"" << (*fData)[i].GetRule() << "\");" << endl; +} Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1488) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1489) @@ -9,16 +9,19 @@ #endif -class MDataChain; +class MDataArray; class MHMatrix : public MH { protected: - Int_t fNumRow; //! Number of dimensions of histogram - TMatrix fM; // Matrix to be filled + Int_t fNumRow; //! Number of dimensions of histogram + TMatrix fM; // Matrix to be filled - TMatrix fM2; //! + TMatrix fM2; //! - TList *fData; //! List of data members (columns) - TList *fRules; // List of data members as text for storage + MDataArray *fData; // List of data members (columns) + + enum { + kIsOwner = BIT(14) + }; void AddRow(); @@ -28,14 +31,17 @@ Bool_t Finalize(); + void StreamPrimitive(ofstream &out) const; + public: MHMatrix(const char *name=NULL, const char *title=NULL); + MHMatrix(MDataArray *mat, const char *name=NULL, const char *title=NULL); ~MHMatrix(); void AddColumn(const char *name); - void AddColumns(const MHMatrix *matrix); + void AddColumns(MDataArray *mat); - // TMatrix &GetM() { return fM; } + MDataArray *GetColumns() { return fData; } + const TMatrix &GetM() const { return fM; } - const TList *GetRules() const { return fRules; } //void Draw(Option_t *opt=NULL); @@ -49,4 +55,6 @@ Double_t CalcDist(const TVector &v, Int_t num = 25); + void SetIOwner(Bool_t b=kTRUE) { b ? SetBit(kIsOwner) : ResetBit(kIsOwner); } + void Reassign();