source: trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc@ 1895

Last change on this file since 1895 was 1809, checked in by wittek, 22 years ago
*** empty log message ***
File size: 6.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Rudy Bock, 5/2002 <mailto:rkb@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2002
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MMultiDimDistCalc
29//
30// Calculated a multidimensional distance. It calculates the distance to
31// all vectors in a given matrix describing Gammas and another one
32// describing Hadrons (non gammas). The shortest distances are avaraged.
33// How many distances are used for avaraging can be specified in the
34// constructor.
35//
36// * If you want to use the nearest neighbor function for calculation use:
37// MMultiDimDistCalc::SetUseKernelMethod(kFALSE);
38// * If you want to use the kernel function for calculation use:
39// MMultiDimDistCalc::SetUseKernelMethod(kTRUE); <default>
40// * To use only the n next neighbors for your calculation use:
41// MMultiDimDistCalc::SetUseNumRows(n);
42// * To use all reference events set the number to 0 <default>
43//
44////////////////////////////////////////////////////////////////////////////
45#include "MMultiDimDistCalc.h"
46
47#include <fstream.h>
48
49#include "MHMatrix.h" // must be before MLogManip.h
50
51#include "MLog.h"
52#include "MLogManip.h"
53
54#include "MParList.h"
55#include "MDataChain.h"
56#include "MDataArray.h"
57
58#include "MHadronness.h"
59
60ClassImp(MMultiDimDistCalc);
61
62static const TString gsDefName = "MMultiDimDistCalc";
63static const TString gsDefTitle = "Composite Probabilities Loop 1/2";
64// --------------------------------------------------------------------------
65//
66// Setup histograms and the number of distances which are used for
67// avaraging in CalcDist
68//
69MMultiDimDistCalc::MMultiDimDistCalc(const char *name, const char *title)
70 : fNum(0), fUseKernel(kTRUE), fData(NULL)
71{
72 //
73 // set the name and title of this object
74 //
75 fName = name ? name : gsDefName.Data();
76 fTitle = title ? title : gsDefTitle.Data();
77
78 /*
79 fData = new TList;
80 fData->SetOwner();
81 */
82}
83
84// --------------------------------------------------------------------------
85//
86// Delete the data chains
87//
88MMultiDimDistCalc::~MMultiDimDistCalc()
89{
90 // delete fData;
91}
92
93// --------------------------------------------------------------------------
94//
95// Needs:
96// - MatrixGammas [MHMatrix]
97// - MatrixHadrons {MHMatrix]
98// - MHadroness
99// - all data containers used to build the matrixes
100//
101// The matrix object can be filles using MFillH. And must be of the same
102// number of columns (with the same meaning).
103//
104Bool_t MMultiDimDistCalc::PreProcess(MParList *plist)
105{
106 fMGammas = (MHMatrix*)plist->FindObject("MatrixGammas", "MHMatrix");
107 if (!fMGammas)
108 {
109 *fLog << err << dbginf << "MatrixGammas [MHMatrix] not found... aborting." << endl;
110 return kFALSE;
111 }
112
113 fMHadrons = (MHMatrix*)plist->FindObject("MatrixHadrons", "MHMatrix");
114 if (!fMHadrons)
115 {
116 *fLog << err << dbginf << "MatrixHadrons [MHMatrix] not found... aborting." << endl;
117 return kFALSE;
118 }
119
120 if (fMGammas->GetM().GetNcols() != fMHadrons->GetM().GetNcols())
121 {
122 *fLog << err << dbginf << "Error matrices have different numbers of columns... aborting." << endl;
123 return kFALSE;
124 }
125
126 /*
127 TIter Next(fMGammas->GetRules());
128 TObject *data=NULL;
129 while ((data=Next()))
130 {
131 MDataChain *chain = new MDataChain(data->GetName());
132 if (!chain->PreProcess(plist))
133 {
134 delete chain;
135 return kFALSE;
136 }
137 fData->Add(chain);
138 }
139 */
140 fData = fMGammas->GetColumns();
141 if (!fData)
142 {
143 *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl;
144 return kFALSE;
145 }
146
147 if (!fData->PreProcess(plist))
148 {
149 *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl;
150 return kFALSE;
151 }
152
153 fHadronness = (MHadronness*)plist->FindCreateObj("MHadronness");
154 if (!fHadronness)
155 return kFALSE;
156
157 return kTRUE;
158}
159
160// --------------------------------------------------------------------------
161//
162// Evaluates the avarage over the fNum shortest distances in a
163// multidimensional space to a matrix (set of vectors) describing
164// gammas and hadrons.
165// The hadroness of the event is defines as the avarage distance to the
166// set of gammas (dg) divided by the avarage distance to the set of
167// hadrons (dh). Because this value is in teh range [0, inf] it is
168// transformed into [0,1] by:
169// H = exp(-dh/dg);
170//
171Bool_t MMultiDimDistCalc::Process()
172{
173 const Double_t ncols = fMGammas->GetM().GetNcols();
174
175 TVector event(ncols);
176
177 for (int i=0; i<fData->GetNumEntries(); i++)
178 event(i) = (*fData)(i);
179
180 /*
181 Int_t n=0;
182 TIter Next(fData);
183 MData *data = NULL;
184 while ((data=(MData*)Next()))
185 event(n++) = data->GetValue();
186 */
187
188 Double_t numg = fNum;
189 Double_t numh = fNum;
190 if (fNum==0)
191 {
192 numg = fMGammas->GetM().GetNrows();
193 numh = fMHadrons->GetM().GetNrows();
194 }
195 if (fUseKernel)
196 {
197 numg = -numg;
198 numh = -numh;
199 }
200
201 Double_t dg = fMGammas->CalcDist(event, numg);
202 Double_t dh = fMHadrons->CalcDist(event, numh);
203
204 if (dg<0 || dh<0)
205 {
206 *fLog << err << "MHMatrix::CalcDist failed (dg=" << dg << ", dh=" << dh << ")... aborting" << endl;
207 return kERROR;
208 }
209
210 Double_t arg;
211
212 if (dg+dh != 0.0)
213 arg = dg / (dg+dh);
214 else
215 arg = 1.e10;
216 //fHadronness->SetHadronness(arg);
217
218 if (dg != 0.0)
219 arg = exp(-dh/dg);
220 else
221 arg = 0.0;
222 fHadronness->SetHadronness(arg);
223
224
225 return kTRUE;
226}
227
228void MMultiDimDistCalc::StreamPrimitive(ofstream &out) const
229{
230 out << " MMultiDimDist " << GetUniqueName();
231
232 if (fName!=gsDefName || fTitle!=gsDefTitle)
233 {
234 out << "(\"" << fName << "\"";
235 if (fTitle!=gsDefTitle)
236 out << ", \"" << fTitle << "\")";
237 }
238 out << ";" << endl;
239
240 if (fNum!=0)
241 out << " " << GetUniqueName() << ".SetUseNumRows(" << fNum << ");" << endl;
242 if (fUseKernel)
243 out << " " << GetUniqueName() << ".SetUseKernelMethod();" << endl;
244}
Note: See TracBrowser for help on using the repository browser.