source: trunk/MagicSoft/Mars/mtools/MTFillMatrix.cc@ 9029

Last change on this file since 9029 was 8907, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 13.5 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, 12/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MTFillMatrix
28//
29// Use this tool to fill eg trainings and test-matrices, while the matrix
30// to be filled can be a real matrix (MHMatrix) or a file (MWriteRootFile)
31// or both.
32//
33// First create a reference histogram (MH3). For more details see
34// MFEventSelector2 which is used to select events according to the
35// reference histogram.
36//
37// If no reference histogram is available the number of events are
38// randomly choosen from the sample with a probability which fits
39// the total destination number of events.
40//
41// Here is an example of how to choose 1000 events somehow distributed in
42// size from a file.
43// -----------------------------------------------------------------------
44//
45// MH3 ref("MHillas.fSize"); // choose a predefined size distribution
46// // Now setup the distribution
47//
48// MHMatrix matrix1; // create matrix to fill
49// matrix.AddColumn("MHillas.fWidth"); // setup columns of your matrix
50//
51// MReadMarsFile read("myfile.root"); // Setup your 'reader'
52// read.DisableAutoScheme(); // make sure everything is read
53//
54// MTFillMatrix fill(&ref); // setup MTFillMatrix
55// fill.SetNumDestEvents1(1000); // setup number of events to select
56// fill.SetDestMatrix1(&matrix1); // setup destination matrix
57// if (!fill.Process()) // check if everything worked
58// return;
59// fill.WriteMatrix1("myoutput.root"); // write matrix to file
60//
61//
62// To get two matrices instead of one (splitted randomly) you can add
63// the following before calling Process():
64// ------------------------------------------------------------------------
65//
66// MHMatrix matrix2;
67// fill.SetNumDestEvents2(500); // setup number of events to select
68// fill.SetDestMatrix2(&matrix2); // setup destination matrix
69// [...]
70// fill.WriteMatrix2("myoutput2.root");
71//
72//
73// To write both matrices into a single file use:
74// ----------------------------------------------
75//
76// fill.WriteMatrices("myfile.root");
77//
78/////////////////////////////////////////////////////////////////////////////
79#include "MTFillMatrix.h"
80
81#include <TFile.h>
82#include <TMath.h>
83
84// environment
85#include "MHMatrix.h"
86
87#include "MLog.h"
88#include "MLogManip.h"
89
90// eventloop
91#include "MParList.h"
92#include "MTaskList.h"
93#include "MEvtLoop.h"
94
95// tasks
96#include "MRead.h"
97#include "MFillH.h"
98#include "MContinue.h"
99
100// filters
101#include "MFDataPhrase.h"
102#include "MFilterList.h"
103#include "MFEventSelector.h"
104#include "MFEventSelector2.h"
105
106ClassImp(MTFillMatrix);
107
108using namespace std;
109
110// --------------------------------------------------------------------------
111//
112// Print Size and contained columns.
113// Check whether the number of generated events is compatible with
114// the number of requested events.
115//
116Bool_t MTFillMatrix::CheckResult(MHMatrix *m, Int_t num) const
117{
118 if (!m || num==0)
119 return kTRUE;
120
121 m->Print("SizeCols");
122
123 const Int_t generated = m->GetM().GetNrows();
124 if (TMath::Abs(generated-num) <= TMath::Sqrt(9.0*num))
125 return kTRUE;
126
127 *fLog << warn << "WARNING - No. of generated events (";
128 *fLog << generated << ") incompatible with requested events (";
129 *fLog << num << ") for " << m->GetDescriptor() << endl;
130
131 return kFALSE;
132}
133
134// --------------------------------------------------------------------------
135//
136// Write the given MHMatrix with its name as default name to a
137// file fname.
138//
139Bool_t MTFillMatrix::WriteMatrix(MHMatrix *m, const TString &fname, Int_t i) const
140{
141 if (!m)
142 {
143 *fLog << "ERROR - Unable to write matrix #" << i << " (=NULL)... ignored." << endl;
144 return kFALSE;
145 }
146 if (fname.IsNull())
147 {
148 *fLog << "ERROR - Unable to write matrix, file name empty." << endl;
149 return kFALSE;
150 }
151
152 TFile file(fname, "RECREATE", m->GetTitle());
153 m->Write();
154 return kTRUE;
155}
156
157void MTFillMatrix::Init(const char *name, const char *title)
158{
159 fName = name ? name : "MFillMatrix";
160 fTitle = title ? title : "Tool to fill MHMatrix from file";
161}
162
163MTFillMatrix::MTFillMatrix(const char *name, const char *title)
164: fReference(0), fReader(0), fDestMatrix1(0), fDestMatrix2(0),
165 fNumDestEvents1(0), fNumDestEvents2(0), fWriteFile1(0), fWriteFile2(0)
166{
167 Init(name, title);
168}
169
170// --------------------------------------------------------------------------
171//
172// Constructor. Takes an MH3 as argument. This MH3 is the reference
173// distribution to fill the matrix. More information can be found
174// at MFEventSelector2 which is used to select the events.
175//
176// If no MH3 *ref is given the events are randomly selected from the
177// total sample - this may result in samples which don't have exactly
178// the predefined size, but it is much faster.
179//
180MTFillMatrix::MTFillMatrix(const MH3 *ref, const char *name, const char *title)
181: fReference(0), fReader(0), fDestMatrix1(0), fDestMatrix2(0),
182 fNumDestEvents1(0), fNumDestEvents2(0), fWriteFile1(0), fWriteFile2(0)
183{
184 Init(name, title);
185 if (ref)
186 fReference = (MH3*)ref->Clone();
187}
188
189MTFillMatrix::~MTFillMatrix()
190{
191 if (fReference)
192 delete fReference;
193}
194
195//------------------------------------------------------------------------
196//
197// Function serving AddPreCuts, AddPreTasks and AddPostTasks
198//
199void MTFillMatrix::Add(const TList &src, const TClass *cls, TList &dest)
200{
201 TIter Next(&src);
202 TObject *obj=0;
203 while ((obj=Next()))
204 if (obj->InheritsFrom(cls))
205 dest.Add(obj);
206}
207
208//------------------------------------------------------------------------
209//
210// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5"
211// (The rule is applied, nit inverted: The matrix is filled with
212// the events fullfilling the condition)
213//
214void MTFillMatrix::AddPreCut(const char *rule)
215{
216 MFilter *f = new MFDataPhrase(rule);
217 f->SetBit(kCanDelete);
218 AddPreCut(f);
219}
220
221//------------------------------------------------------------------------
222//
223// Add a cut which is used to fill the matrix. If kCanDelete is set
224// MJOptimize takes the ownership.
225//
226void MTFillMatrix::AddPreCut(MFilter *f)
227{
228 fPreCuts.Add(f);
229}
230
231//------------------------------------------------------------------------
232//
233// Add all entries deriving from MFilter from list to PreCuts.
234// The ownership is not affected.
235//
236void MTFillMatrix::AddPreCuts(const TList &list)
237{
238 Add(list, MFilter::Class(), fPreCuts);
239}
240
241//------------------------------------------------------------------------
242//
243// Add a task which is executed before the precuts. If kCanDelete is set
244// MJOptimize takes the ownership.
245//
246void MTFillMatrix::AddPreTask(MTask *t)
247{
248 fPreTasks.Add(t);
249}
250
251//------------------------------------------------------------------------
252//
253// Add all entries deriving from MTask from list to PreTasks.
254// The ownership is not affected.
255//
256void MTFillMatrix::AddPreTasks(const TList &list)
257{
258 Add(list, MTask::Class(), fPreTasks);
259}
260
261//------------------------------------------------------------------------
262//
263// Add a task which is executed after the precuts. If kCanDelete is set
264// MJOptimize takes the ownership.
265//
266void MTFillMatrix::AddPostTask(MTask *t)
267{
268 fPostTasks.Add(t);
269}
270
271//------------------------------------------------------------------------
272//
273// Add all entries deriving from MTask from list to PostTasks.
274// The ownership is not affected.
275//
276void MTFillMatrix::AddPostTasks(const TList &list)
277{
278 Add(list, MTask::Class(), fPostTasks);
279}
280
281// --------------------------------------------------------------------------
282//
283// Fill the matrix (FIXME: Flow diagram missing)
284//
285Bool_t MTFillMatrix::Process(const MParList &parlist)
286{
287 if (!fReader)
288 {
289 *fLog << err << "ERROR - No task to read data was given... abort." << endl;
290 return kFALSE;
291 }
292
293 *fLog << inf;
294 fLog->Separator(GetDescriptor());
295 *fLog << "Fill " << fDestMatrix1->GetDescriptor() << " with " << fNumDestEvents1 << endl;
296 if (fDestMatrix2)
297 *fLog << "Fill " << fDestMatrix2->GetDescriptor() << " with " << fNumDestEvents2 << endl;
298 *fLog << "Distribution choosen ";
299 if (fReference && fReference->GetHist().GetEntries()>0)
300 *fLog << "from " << fReference->GetDescriptor();
301 else
302 *fLog << "randomly";
303 *fLog << endl;
304
305 //
306 // Create parameter list and task list, add tasklist to parlist
307 //
308 parlist.Print();
309 MParList plist(parlist);
310 MTaskList tlist;
311 plist.AddToList(&tlist);
312
313 //
314 // A selector to select a given number of events from a sample
315 //
316 // FIXME: Merge MFEventSelector and MFEventSelector2
317 MFilter *selector=0;
318 if (fNumDestEvents1>0 || fNumDestEvents2>0)
319 {
320 if (fReference)
321 {
322 // Case of a reference/nominal distribution
323 // The events must be read before selection
324 MFEventSelector2 *sel = new MFEventSelector2(*fReference);
325 sel->SetNumMax(fNumDestEvents1+fNumDestEvents2);
326 sel->SetInverted();
327
328 selector = sel;
329 }
330 else
331 {
332 // Case of a random distribution
333 // The events can be selected before reading
334 MFEventSelector *sel = new MFEventSelector;
335 sel->SetNumSelectEvts(fNumDestEvents1+fNumDestEvents2);
336 fReader->SetSelector(sel);
337
338 selector = sel;
339 }
340 }
341
342 //
343 // Continue for PreCuts
344 //
345 MFilterList list;
346 list.SetName("PreCuts");
347 if (!list.AddToList(fPreCuts))
348 *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
349
350 MContinue cont0(&list);
351 cont0.SetInverted();
352
353 //
354 // Continue for all events which are not (SetInverted())
355 // selected by the 'selector'
356 //
357 MContinue cont(selector);
358
359 //
360 // Create a filter doing a random split
361 //
362 const Double_t prob = (Double_t)fNumDestEvents1/(fNumDestEvents1+fNumDestEvents2);
363 MFEventSelector split;
364 split.SetSelectionRatio(prob);
365
366 //
367 // Create the logical inverted filter for 'split'
368 //
369 MFilterList invsplit;
370 invsplit.AddToList(&split);
371 invsplit.SetInverted();
372
373 //
374 // The two tasks filling the two matrices
375 //
376 MFillH fill1(fDestMatrix1);
377 MFillH fill2(fDestMatrix2);
378 fill1.SetFilter(&split);
379 fill2.SetFilter(&invsplit);
380
381 // entries in MTaskList
382 tlist.AddToList(fReader); // Read events
383 tlist.AddToList(fPreTasks, 0); // PreTasks
384 if (fPreCuts.GetEntries()>0)
385 tlist.AddToList(&cont0); // PreCuts
386 if (fReference && selector)
387 tlist.AddToList(&cont); // select a sample of events
388 tlist.AddToList(&invsplit); // process invsplit (which implicitly processes split)
389 tlist.AddToList(fPostTasks, 0); // PostTasks
390 if (fDestMatrix1)
391 tlist.AddToList(&fill1); // fill matrix 1
392 if (fDestMatrix2)
393 tlist.AddToList(&fill2); // fill matrix 2
394 if (fWriteFile1)
395 {
396 fWriteFile1->SetFilter(&split);
397 tlist.AddToList(fWriteFile1);
398 }
399 if (fWriteFile2)
400 {
401 fWriteFile2->SetFilter(&invsplit);
402 tlist.AddToList(fWriteFile2);
403 }
404
405 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
406
407 //
408 // Execute the eventloop
409 //
410 MEvtLoop evtloop(fName);
411 evtloop.SetParList(&plist);
412 evtloop.SetDisplay(fDisplay);
413 evtloop.SetLogStream(fLog);
414
415 const Bool_t rc = evtloop.Eventloop();
416
417 if (selector)
418 delete selector;
419
420 if (!rc)
421 {
422 *fLog << err << GetDescriptor() << ": Failed." << endl;
423 return kFALSE;
424 }
425
426 // Check the result of filling...
427 CheckResult(fDestMatrix1, fNumDestEvents1);
428 CheckResult(fDestMatrix2, fNumDestEvents2);
429
430 *fLog << inf << GetDescriptor() << ": Done." << endl;
431 return kTRUE;
432}
433
434// --------------------------------------------------------------------------
435//
436// Write both matrices to a file. Return kFALSE if writing one of them
437// failed or both have the same name.
438//
439Bool_t MTFillMatrix::WriteMatrices(const TString &fname) const
440{
441 if (fDestMatrix1 && fDestMatrix2 &&
442 (TString)fDestMatrix1->GetName()==(TString)fDestMatrix2->GetName())
443 {
444 *fLog << "ERROR - Cannot write matrices (both have the same name)... ignored." << endl;
445 return kFALSE;
446 }
447
448 return WriteMatrix1(fname) && WriteMatrix2(fname);
449}
450
451Int_t MTFillMatrix::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
452{
453 Bool_t rc = kFALSE;
454 if (IsEnvDefined(env, prefix, "NumDestEvents1", print))
455 {
456 rc = kTRUE;
457 SetNumDestEvents1(GetEnvValue(env, prefix, "NumDestEvents1", fNumDestEvents1));
458 }
459 if (IsEnvDefined(env, prefix, "NumDestEvents2", print))
460 {
461 rc = kTRUE;
462 SetNumDestEvents2(GetEnvValue(env, prefix, "NumDestEvents2", fNumDestEvents2));
463 }
464 return rc;
465}
Note: See TracBrowser for help on using the repository browser.