source: branches/AddingGoogleTestEnvironment/mhbase/MHMatrix.cc@ 18180

Last change on this file since 18180 was 9851, checked in by tbretz, 14 years ago
Changed MH::SetBinning and similar functions to take references instead of pointers as arguments. For convenience wrappers for the old style functions are provided.
File size: 33.7 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 2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Rudy Boeck 2003 <mailto:
20! Wolfgang Wittek2003 <mailto:wittek@mppmu.mpg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2003
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// MHMatrix
30//
31// This is a histogram container which holds a matrix with one column per
32// data variable. The data variable can be a complex rule (MDataPhrase).
33// Each event for wich Fill is called (by MFillH) is added as a new
34// row to the matrix.
35//
36// For example:
37// MHMatrix m;
38// m.AddColumn("MHillas.fSize");
39// m.AddColumn("MMcEvt.fImpact/100");
40// m.AddColumn("HillasSource.fDist*MGeomCam.fConvMm2Deg");
41// MFillH fillm(&m);
42// taskliost.AddToList(&fillm);
43// [...]
44// m.Print();
45//
46/////////////////////////////////////////////////////////////////////////////
47#include "MHMatrix.h"
48
49#include <fstream>
50
51#include <TMath.h>
52
53#include <TList.h>
54#include <TArrayF.h>
55#include <TArrayD.h>
56#include <TArrayI.h>
57
58#include <TH1.h>
59#include <TCanvas.h>
60#include <TRandom3.h>
61
62#include "MLog.h"
63#include "MLogManip.h"
64
65#include "MFillH.h"
66#include "MEvtLoop.h"
67#include "MParList.h"
68#include "MTaskList.h"
69#include "MProgressBar.h"
70
71#include "MData.h"
72#include "MDataArray.h"
73#include "MFilter.h"
74
75ClassImp(MHMatrix);
76
77using namespace std;
78
79const TString MHMatrix::gsDefName = "MHMatrix";
80const TString MHMatrix::gsDefTitle = "Multidimensional Matrix";
81
82// --------------------------------------------------------------------------
83//
84// Default Constructor
85//
86MHMatrix::MHMatrix(const char *name, const char *title)
87 : fNumRows(0), fData(NULL)
88{
89 fName = name ? name : gsDefName.Data();
90 fTitle = title ? title : gsDefTitle.Data();
91}
92
93// --------------------------------------------------------------------------
94//
95// Default Constructor
96//
97MHMatrix::MHMatrix(const TMatrix &m, const char *name, const char *title)
98 : fNumRows(m.GetNrows()), fM(m), fData(NULL)
99{
100 fName = name ? name : gsDefName.Data();
101 fTitle = title ? title : gsDefTitle.Data();
102}
103
104// --------------------------------------------------------------------------
105//
106// Constructor. Initializes the columns of the matrix with the entries
107// from a MDataArray
108//
109MHMatrix::MHMatrix(MDataArray *mat, const char *name, const char *title)
110 : fNumRows(0), fData(mat)
111{
112 fName = name ? name : gsDefName.Data();
113 fTitle = title ? title : gsDefTitle.Data();
114}
115
116// --------------------------------------------------------------------------
117//
118// Constructor. Initializes the columns of the matrix with the entries
119// from a MDataArray
120//
121MHMatrix::MHMatrix(const TMatrix &m, MDataArray *mat, const char *name, const char *title)
122 : fNumRows(m.GetNrows()), fM(m), fData(mat)
123{
124 fName = name ? name : gsDefName.Data();
125 fTitle = title ? title : gsDefTitle.Data();
126}
127
128// --------------------------------------------------------------------------
129//
130// Destructor. Does not deleted a user given MDataArray, except IsOwner
131// was called.
132//
133MHMatrix::~MHMatrix()
134{
135 if (TestBit(kIsOwner) && fData)
136 delete fData;
137}
138
139// --------------------------------------------------------------------------
140//
141Bool_t MHMatrix::SetNumRow(Int_t row)
142{
143 if (row>=fM.GetNrows() || row<0) return kFALSE;
144 fRow = row;
145 return kTRUE;
146}
147
148// --------------------------------------------------------------------------
149//
150// Add a new column to the matrix. This can only be done before the first
151// event (row) was filled into the matrix. For the syntax of the rule
152// see MDataPhrase.
153// Returns the index of the new column, -1 in case of failure.
154// (0, 1, 2, ... for the 1st, 2nd, 3rd, ...)
155//
156Int_t MHMatrix::AddColumn(const char *rule)
157{
158 const Int_t idx = fData ? fData->FindRule(rule) : -1;
159 if (idx>=0)
160 return idx;
161
162 if (IsValid(fM))
163 {
164 *fLog << warn << "Warning - matrix is already in use. Can't add a new column... skipped." << endl;
165 return -1;
166 }
167
168 if (TestBit(kIsLocked))
169 {
170 *fLog << warn << "Warning - matrix is locked. Can't add new column... skipped." << endl;
171 return -1;
172 }
173
174 if (!fData)
175 {
176 fData = new MDataArray;
177 SetBit(kIsOwner);
178 }
179
180 fData->AddEntry(rule);
181 return fData->GetNumEntries()-1;
182}
183
184// --------------------------------------------------------------------------
185//
186void MHMatrix::AddColumns(MDataArray *matrix)
187{
188 if (IsValid(fM))
189 {
190 *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl;
191 return;
192 }
193
194 if (TestBit(kIsLocked))
195 {
196 *fLog << warn << "Warning - matrix is locked. Can't add new columns... skipped." << endl;
197 return;
198 }
199
200 if (fData)
201 *fLog << warn << "Warning - columns already added... replacing." << endl;
202
203 if (fData && TestBit(kIsOwner))
204 {
205 delete fData;
206 ResetBit(kIsOwner);
207 }
208
209 fData = matrix;
210}
211
212void MHMatrix::AddColumns(const TCollection &list)
213{
214 TIter Next(&list);
215 TObject *obj = 0;
216 while ((obj=Next()))
217 AddColumn(obj->GetName());
218}
219
220
221// --------------------------------------------------------------------------
222//
223// Checks whether at least one column is available and PreProcesses all
224// data chains.
225//
226Bool_t MHMatrix::SetupFill(const MParList *plist)
227{
228 if (!fData)
229 {
230 *fLog << err << "Error - No Columns initialized... aborting." << endl;
231 return kFALSE;
232 }
233
234 return fData->PreProcess(plist);
235}
236
237// --------------------------------------------------------------------------
238//
239// If the matrix has not enough rows double the number of available rows.
240//
241void MHMatrix::AddRow()
242{
243 fNumRows++;
244
245 if (fM.GetNrows() > fNumRows)
246 return;
247
248 if (!IsValid(fM))
249 {
250 fM.ResizeTo(1, fData->GetNumEntries());
251 return;
252 }
253
254#if ROOT_VERSION_CODE < ROOT_VERSION(3,05,07)
255 TMatrix m(fM);
256#endif
257
258 // If we exceed ~50MB of memory we start growing slower
259 const Int_t n = fM.GetNrows()>10000000/fM.GetNcols() ? TMath::Nint(fM.GetNrows()*1.2) : fM.GetNrows()*2;
260 fM.ResizeTo(n, fData->GetNumEntries());
261
262#if ROOT_VERSION_CODE < ROOT_VERSION(3,05,07)
263 TVector vold(fM.GetNcols());
264 for (int x=0; x<m.GetNrows(); x++)
265 TMatrixRow(fM, x) = vold = TMatrixRow(m, x);
266#endif
267}
268
269// --------------------------------------------------------------------------
270//
271// Add the values correspoding to the columns to the new row
272//
273Int_t MHMatrix::Fill(const MParContainer *par, const Stat_t w)
274{
275 AddRow();
276
277 for (int col=0; col<fData->GetNumEntries(); col++)
278 fM(fNumRows-1, col) = (*fData)(col);
279
280 return kTRUE;
281}
282
283// --------------------------------------------------------------------------
284//
285// Resize the matrix to a number of rows which corresponds to the number of
286// rows which have really been filled with values.
287//
288Bool_t MHMatrix::Finalize()
289{
290 //
291 // It's not a fatal error so we don't need to stop PostProcessing...
292 //
293 if (fData->GetNumEntries()==0 || fNumRows<1)
294 return kTRUE;
295
296 if (fNumRows != fM.GetNrows())
297 {
298 TMatrix m(fM);
299 CopyCrop(fM, m, fNumRows);
300 }
301
302 return kTRUE;
303}
304
305/*
306// --------------------------------------------------------------------------
307//
308// Draw clone of histogram. So that the object can be deleted
309// and the histogram is still visible in the canvas.
310// The cloned object are deleted together with the canvas if the canvas is
311// destroyed. If you want to handle destroying the canvas you can get a
312// pointer to it from this function
313//
314TObject *MHMatrix::DrawClone(Option_t *opt) const
315{
316 TCanvas &c = *MH::MakeDefCanvas(fHist);
317
318 //
319 // This is necessary to get the expected bahviour of DrawClone
320 //
321 gROOT->SetSelectedPad(NULL);
322
323 fHist->DrawCopy(opt);
324
325 TString str(opt);
326 if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
327 {
328 TProfile *p = ((TH2*)fHist)->ProfileX();
329 p->Draw("same");
330 p->SetBit(kCanDelete);
331 }
332 if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
333 {
334 TProfile *p = ((TH2*)fHist)->ProfileY();
335 p->Draw("same");
336 p->SetBit(kCanDelete);
337 }
338
339 c.Modified();
340 c.Update();
341
342 return &c;
343}
344
345// --------------------------------------------------------------------------
346//
347// Creates a new canvas and draws the histogram into it.
348// Be careful: The histogram belongs to this object and won't get deleted
349// together with the canvas.
350//
351void MHMatrix::Draw(Option_t *opt)
352{
353 if (!gPad)
354 MH::MakeDefCanvas(fHist);
355
356 fHist->Draw(opt);
357
358 TString str(opt);
359 if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
360 {
361 TProfile *p = ((TH2*)fHist)->ProfileX();
362 p->Draw("same");
363 p->SetBit(kCanDelete);
364 }
365 if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
366 {
367 TProfile *p = ((TH2*)fHist)->ProfileY();
368 p->Draw("same");
369 p->SetBit(kCanDelete);
370 }
371
372 gPad->Modified();
373 gPad->Update();
374}
375*/
376
377// --------------------------------------------------------------------------
378//
379// Prints the meaning of the columns and the contents of the matrix.
380// Becareful, this can take a long time for matrices with many rows.
381// Use the option 'size' to print the size of the matrix.
382// Use the option 'cols' to print the culumns
383// Use the option 'data' to print the contents
384//
385void MHMatrix::Print(Option_t *o) const
386{
387 TString str(o);
388
389 *fLog << all << flush;
390
391 if (str.Contains("size", TString::kIgnoreCase))
392 {
393 *fLog << GetDescriptor() << ": NumColumns=" << fM.GetNcols();
394 *fLog << " NumRows=" << fM.GetNrows() << endl;
395 }
396
397 if (!fData && str.Contains("cols", TString::kIgnoreCase))
398 *fLog << "Sorry, no column information available." << endl;
399
400 if (fData && str.Contains("cols", TString::kIgnoreCase))
401 {
402 fData->SetName(fName);
403 fData->Print();
404 }
405
406 if (str.Contains("data", TString::kIgnoreCase))
407 fM.Print();
408}
409
410// --------------------------------------------------------------------------
411//
412const TMatrix *MHMatrix::InvertPosDef()
413{
414 TMatrix m(fM);
415
416 const Int_t rows = m.GetNrows();
417 const Int_t cols = m.GetNcols();
418
419 for (int x=0; x<cols; x++)
420 {
421 Double_t avg = 0;
422 for (int y=0; y<rows; y++)
423 avg += fM(y, x);
424
425 avg /= rows;
426
427#if ROOT_VERSION_CODE > ROOT_VERSION(5,00,00)
428 TMatrixFColumn(m, x) += -avg;
429#else
430 TMatrixColumn(m, x) += -avg;
431#endif
432 }
433
434 TMatrix *m2 = new TMatrix(m, TMatrix::kTransposeMult, m);
435
436 Double_t det;
437 m2->Invert(&det);
438 if (det==0)
439 {
440 *fLog << err << "ERROR - MHMatrix::InvertPosDef failed (Matrix is singular)." << endl;
441 delete m2;
442 return NULL;
443 }
444
445 // m2->Print();
446
447 return m2;
448}
449
450// --------------------------------------------------------------------------
451//
452// Calculated the distance of vector evt from the reference sample
453// represented by the covariance metrix m.
454// - If n<0 the kernel method is applied and
455// -log(sum(epx(-d/h))/n) is returned.
456// - For n>0 the n nearest neighbors are summed and
457// sqrt(sum(d)/n) is returned.
458// - if n==0 all distances are summed
459//
460Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const
461{
462 if (num==0) // may later be used for another method
463 {
464 TVector d = evt;
465 d *= m;
466 return TMath::Sqrt(d*evt);
467 }
468
469 const Int_t rows = fM.GetNrows();
470 const Int_t cols = fM.GetNcols();
471
472 TArrayD dists(rows);
473
474 //
475 // Calculate: v^T * M * v
476 //
477 for (int i=0; i<rows; i++)
478 {
479 TVector col(cols);
480#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
481 col = TMatrixRow(fM, i);
482#else
483 col = TMatrixFRow_const(fM, i);
484#endif
485 TVector d = evt;
486 d -= col;
487
488 TVector d2 = d;
489 d2 *= m;
490
491 dists[i] = d2*d; // square of distance
492
493 //
494 // This corrects for numerical uncertanties in cases of very
495 // small distances...
496 //
497 if (dists[i]<0)
498 dists[i]=0;
499 }
500
501 TArrayI idx(rows);
502 TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);
503
504 Int_t from = 0;
505 Int_t to = TMath::Abs(num)<rows ? TMath::Abs(num) : rows;
506 //
507 // This is a zero-suppression for the case a test- and trainings
508 // sample is identical. This would result in an unwanted leading
509 // zero in the array. To suppress also numerical uncertanties of
510 // zero we cut at 1e-5. Due to Rudy this should be enough. If
511 // you encounter problems we can also use (eg) 1e-25
512 //
513 if (dists[idx[0]]<1e-5)
514 {
515 from++;
516 to ++;
517 if (to>rows)
518 to = rows;
519 }
520
521 if (num<0)
522 {
523 //
524 // Kernel function sum (window size h set according to literature)
525 //
526 const Double_t h = TMath::Power(rows, -1./(cols+4));
527 const Double_t hwin = h*h*2;
528
529 Double_t res = 0;
530 for (int i=from; i<to; i++)
531 res += TMath::Exp(-dists[idx[i]]/hwin);
532
533 return -TMath::Log(res/(to-from));
534 }
535 else
536 {
537 //
538 // Nearest Neighbor sum
539 //
540 Double_t res = 0;
541 for (int i=from; i<to; i++)
542 res += dists[idx[i]];
543
544 return TMath::Sqrt(res/(to-from));
545 }
546}
547
548// --------------------------------------------------------------------------
549//
550// Calls calc dist. In the case of the first call the covariance matrix
551// fM2 is calculated.
552// - If n<0 it is divided by (nrows-1)/h while h is the kernel factor.
553//
554Double_t MHMatrix::CalcDist(const TVector &evt, Int_t num)
555{
556 if (!IsValid(fM2))
557 {
558 if (!IsValid(fM))
559 {
560 *fLog << err << "MHMatrix::CalcDist - ERROR: fM not valid." << endl;
561 return -1;
562 }
563
564 const TMatrix *m = InvertPosDef();
565 if (!m)
566 return -1;
567
568 fM2.ResizeTo(*m);
569 fM2 = *m;
570 fM2 *= fM.GetNrows()-1;
571 delete m;
572 }
573
574 return CalcDist(fM2, evt, num);
575}
576
577// --------------------------------------------------------------------------
578//
579void MHMatrix::Reassign()
580{
581 TMatrix m = fM;
582 fM.ResizeTo(1,1);
583 fM.ResizeTo(m);
584 fM = m;
585}
586
587// --------------------------------------------------------------------------
588//
589// Implementation of SavePrimitive. Used to write the call to a constructor
590// to a macro. In the original root implementation it is used to write
591// gui elements to a macro-file.
592//
593void MHMatrix::StreamPrimitive(ostream &out) const
594{
595 Bool_t data = fData && !TestBit(kIsOwner);
596
597 if (data)
598 {
599 fData->SavePrimitive(out);
600 out << endl;
601 }
602
603 out << " MHMatrix " << GetUniqueName();
604
605 if (data || fName!=gsDefName || fTitle!=gsDefTitle)
606 {
607 out << "(";
608 if (data)
609 out << "&" << fData->GetUniqueName();
610 if (fName!=gsDefName || fTitle!=gsDefTitle)
611 {
612 if (data)
613 out << ", ";
614 out << "\"" << fName << "\"";
615 if (fTitle!=gsDefTitle)
616 out << ", \"" << fTitle << "\"";
617 }
618 }
619 out << ");" << endl;
620
621 if (fData && TestBit(kIsOwner))
622 for (int i=0; i<fData->GetNumEntries(); i++)
623 out << " " << GetUniqueName() << ".AddColumn(\"" << (*fData)[i].GetRule() << "\");" << endl;
624}
625
626// --------------------------------------------------------------------------
627//
628const TArrayI MHMatrix::GetIndexOfSortedColumn(Int_t ncol, Bool_t desc) const
629{
630#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
631 TMatrixColumn col(fM, ncol);
632#else
633 TMatrixFColumn_const col(fM, ncol);
634#endif
635
636 const Int_t n = fM.GetNrows();
637
638 TArrayF array(n);
639
640 for (int i=0; i<n; i++)
641 array[i] = col(i);
642
643 TArrayI idx(n);
644 TMath::Sort(n, array.GetArray(), idx.GetArray(), desc);
645
646 return idx;
647}
648
649// --------------------------------------------------------------------------
650//
651void MHMatrix::SortMatrixByColumn(Int_t ncol, Bool_t desc)
652{
653 TArrayI idx = GetIndexOfSortedColumn(ncol, desc);
654
655 const Int_t n = fM.GetNrows();
656
657#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
658 TVector vold(fM.GetNcols());
659#endif
660
661 TMatrix m(n, fM.GetNcols());
662 for (int i=0; i<n; i++)
663#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
664 TMatrixRow(m, i) = vold = TMatrixRow(fM, idx[i]);
665#else
666 TMatrixFRow(m, i) = TMatrixFRow_const(fM, idx[i]);
667#endif
668 fM = m;
669}
670
671// --------------------------------------------------------------------------
672//
673Bool_t MHMatrix::Fill(MParList *plist, MTask *read, MFilter *filter)
674{
675 //
676 // Read data into Matrix
677 //
678 const Bool_t is = plist->IsOwner();
679 plist->SetOwner(kFALSE);
680
681 MTaskList tlist;
682 plist->Replace(&tlist);
683
684 MFillH fillh(this);
685
686 tlist.AddToList(read);
687
688 if (filter)
689 {
690 tlist.AddToList(filter);
691 fillh.SetFilter(filter);
692 }
693
694 tlist.AddToList(&fillh);
695
696 //MProgressBar bar;
697 MEvtLoop evtloop("MHMatrix::Fill-EvtLoop");
698 evtloop.SetParList(plist);
699 //evtloop.SetProgressBar(&bar);
700
701 if (!evtloop.Eventloop())
702 return kFALSE;
703
704 tlist.PrintStatistics();
705
706 plist->Remove(&tlist);
707 plist->SetOwner(is);
708
709 return kTRUE;
710}
711
712// --------------------------------------------------------------------------
713//
714// Return a comma seperated list of all data members used in the matrix.
715// This is mainly used in MTask::AddToBranchList
716//
717TString MHMatrix::GetDataMember() const
718{
719 return fData ? fData->GetDataMember() : TString("");
720}
721
722// --------------------------------------------------------------------------
723//
724//
725void MHMatrix::ReduceNumberOfRows(UInt_t numrows, const TString opt)
726{
727 UInt_t rows = fM.GetNrows();
728
729 if (rows==numrows)
730 {
731 *fLog << warn << "Matrix has already the correct number of rows..." << endl;
732 return;
733 }
734
735 Float_t ratio = (Float_t)numrows/fM.GetNrows();
736
737 if (ratio>=1)
738 {
739 *fLog << warn << "Matrix cannot be enlarged..." << endl;
740 return;
741 }
742
743 Double_t sum = 0;
744
745 UInt_t oldrow = 0;
746 UInt_t newrow = 0;
747
748#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
749 TVector vold(fM.GetNcols());
750#endif
751 while (oldrow<rows)
752 {
753 sum += ratio;
754
755 if (newrow<=(unsigned int)sum)
756#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
757 TMatrixRow(fM, newrow++) = vold = TMatrixRow(fM, oldrow);
758#else
759 TMatrixFRow(fM, newrow++) = TMatrixFRow_const(fM, oldrow);
760#endif
761
762 oldrow++;
763 }
764}
765
766// ------------------------------------------------------------------------
767//
768// Used in DefRefMatrix to display the result graphically
769//
770void MHMatrix::DrawDefRefInfo(const TH1 &hth, const TH1 &hthd, const TH1 &thsh, Int_t refcolumn)
771{
772 //
773 // Fill a histogram with the distribution after raduction
774 //
775 TH1F hta;
776 hta.SetDirectory(NULL);
777 hta.SetName("hta");
778 hta.SetTitle("Distribution after reduction");
779 CopyBinning(hth, hta);
780
781 for (Int_t i=0; i<fM.GetNrows(); i++)
782 hta.Fill(fM(i, refcolumn));
783
784 TCanvas *th1 = MakeDefCanvas(this);
785 th1->Divide(2,2);
786
787 th1->cd(1);
788 hth.DrawCopy(); // real histogram before
789
790 th1->cd(2);
791 hta.DrawCopy(); // histogram after
792
793 th1->cd(3);
794 hthd.DrawCopy(); // correction factors
795
796 th1->cd(4);
797 thsh.DrawCopy(); // target
798}
799
800// ------------------------------------------------------------------------
801//
802// Resizes th etarget matrix to rows*source.GetNcol() and copies
803// the data from the first (n)rows or the source into the target matrix.
804//
805void MHMatrix::CopyCrop(TMatrix &target, const TMatrix &source, Int_t rows)
806{
807#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
808 TVector v(source.GetNcols());
809#endif
810 target.ResizeTo(rows, source.GetNcols());
811 for (Int_t ir=0; ir<rows; ir++)
812#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
813 TMatrixRow(target, ir) = v = TMatrixRow(source, ir);
814#else
815 TMatrixFRow(target, ir) = TMatrixFRow_const(source, ir);
816#endif
817}
818
819// ------------------------------------------------------------------------
820//
821// Define the reference matrix
822// refcolumn number of the column (starting at 0) containing the variable,
823// for which a target distribution may be given;
824// thsh histogram containing the target distribution of the variable
825// nmaxevts the number of events the reference matrix should have after
826// the renormalization
827// rest a TMatrix conatining the resulting (not choosen)
828// columns of the primary matrix. Maybe NULL if you
829// are not interested in this
830//
831Bool_t MHMatrix::DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh,
832 Int_t nmaxevts, TMatrix *rest)
833{
834 if (!IsValid(fM))
835 {
836 *fLog << err << dbginf << "Matrix not initialized" << endl;
837 return kFALSE;
838 }
839
840 if (thsh.GetMinimum()<0)
841 {
842 *fLog << err << dbginf << "Renormalization not possible: ";
843 *fLog << "Target Distribution has values < 0" << endl;
844 return kFALSE;
845 }
846
847
848 if (nmaxevts>fM.GetNrows())
849 {
850 *fLog << warn << dbginf << "No.requested (" << nmaxevts;
851 *fLog << ") > available events (" << fM.GetNrows() << ")... ";
852 *fLog << "setting equal." << endl;
853 nmaxevts = fM.GetNrows();
854 }
855
856
857 if (nmaxevts<0)
858 {
859 *fLog << err << dbginf << "Number of requested events < 0" << endl;
860 return kFALSE;
861 }
862
863 if (nmaxevts==0)
864 nmaxevts = fM.GetNrows();
865
866 //
867 // refcol is the column number starting at 0; it is >= 0
868 //
869 // number of the column (count from 0) containing
870 // the variable for which the target distribution is given
871 //
872
873 //
874 // Calculate normalization factors
875 //
876 //const int nbins = thsh.GetNbinsX();
877 //const double frombin = thsh.GetBinLowEdge(1);
878 //const double tobin = thsh.GetBinLowEdge(nbins+1);
879 //const double dbin = thsh.GetBinWidth(1);
880
881 const Int_t nrows = fM.GetNrows();
882 const Int_t ncols = fM.GetNcols();
883
884 //
885 // set up the real histogram (distribution before)
886 //
887 //TH1F hth("th", "Distribution before reduction", nbins, frombin, tobin);
888 TH1F hth;
889 hth.SetNameTitle("th", "Distribution before reduction");
890 CopyBinning(thsh, hth);
891 hth.SetDirectory(NULL);
892 for (Int_t j=0; j<nrows; j++)
893 hth.Fill(fM(j, refcolumn));
894
895 //TH1F hthd("thd", "Correction factors", nbins, frombin, tobin);
896 TH1F hthd;
897 hthd.SetNameTitle("thd", "Correction factors");
898 CopyBinning(thsh, hthd);
899 hthd.SetDirectory(NULL);
900 hthd.Divide(&thsh, &hth, 1, 1);
901
902 if (hthd.GetMaximum() <= 0)
903 {
904 *fLog << err << dbginf << "Maximum correction factor <= 0... abort." << endl;
905 return kFALSE;
906 }
907
908 //
909 // ===== obtain correction factors (normalization factors)
910 //
911 hthd.Scale(1/hthd.GetMaximum());
912
913 //
914 // get random access
915 //
916 TArrayI ind(nrows);
917 GetRandomArrayI(ind);
918
919 //
920 // define new matrix
921 //
922 Int_t evtcount1 = -1;
923 Int_t evtcount2 = 0;
924
925 TMatrix mnewtmp(nrows, ncols);
926 TMatrix mrest(nrows, ncols);
927
928 TArrayF cumulweight(nrows); // keep track for each bin how many events
929
930 //
931 // Project values in reference column into [0,1]
932 //
933 TVector v(fM.GetNrows());
934#if ROOT_VERSION_CODE > ROOT_VERSION(5,00,00)
935 v = TMatrixFColumn_const(fM, refcolumn);
936#else
937 v = TMatrixColumn(fM, refcolumn);
938#endif
939 //v += -frombin;
940 //v *= 1/dbin;
941
942 //
943 // select events (distribution after renormalization)
944 //
945 Int_t ir;
946#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
947 TVector vold(fM.GetNcols());
948#endif
949 for (ir=0; ir<nrows; ir++)
950 {
951 // const Int_t indref = (Int_t)v(ind[ir]);
952 const Int_t indref = hthd.FindBin(v(ind[ir])) - 1;
953 cumulweight[indref] += hthd.GetBinContent(indref+1);
954 if (cumulweight[indref]<=0.5)
955 {
956#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
957 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
958#else
959 TMatrixFRow(mrest, evtcount2++) = TMatrixFRow_const(fM, ind[ir]);
960#endif
961 continue;
962 }
963
964 cumulweight[indref] -= 1.;
965 if (++evtcount1 >= nmaxevts)
966 break;
967
968#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
969 TMatrixRow(mnewtmp, evtcount1) = vold = TMatrixRow(fM, ind[ir]);
970#else
971 TMatrixFRow(mnewtmp, evtcount1) = TMatrixFRow_const(fM, ind[ir]);
972#endif
973 }
974
975 for (/*empty*/; ir<nrows; ir++)
976#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
977 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
978#else
979 TMatrixFRow(mrest, evtcount2++) = TMatrixFRow_const(fM, ind[ir]);
980#endif
981
982 //
983 // reduce size
984 //
985 // matrix fM having the requested distribution
986 // and the requested number of rows;
987 // this is the matrix to be used in the g/h separation
988 //
989 CopyCrop(fM, mnewtmp, evtcount1);
990 fNumRows = evtcount1;
991
992 if (evtcount1 < nmaxevts)
993 *fLog << warn << "Reference sample contains less events (" << evtcount1 << ") than requested (" << nmaxevts << ")" << endl;
994
995 if (TestBit(kEnableGraphicalOutput))
996 DrawDefRefInfo(hth, hthd, thsh, refcolumn);
997
998 if (rest)
999 CopyCrop(*rest, mrest, evtcount2);
1000
1001 return kTRUE;
1002}
1003
1004// ------------------------------------------------------------------------
1005//
1006// Returns a array containing randomly sorted indices
1007//
1008void MHMatrix::GetRandomArrayI(TArrayI &ind) const
1009{
1010 const Int_t rows = ind.GetSize();
1011
1012 TArrayF ranx(rows);
1013
1014 TRandom3 rnd(0);
1015 for (Int_t i=0; i<rows; i++)
1016 ranx[i] = rnd.Rndm(i);
1017
1018 TMath::Sort(rows, ranx.GetArray(), ind.GetArray(), kTRUE);
1019}
1020
1021// ------------------------------------------------------------------------
1022//
1023// Define the reference matrix
1024// nmaxevts maximum number of events in the reference matrix
1025// rest a TMatrix conatining the resulting (not choosen)
1026// columns of the primary matrix. Maybe NULL if you
1027// are not interested in this
1028//
1029// the target distribution will be set
1030// equal to the real distribution; the events in the reference
1031// matrix will then be simply a random selection of the events
1032// in the original matrix.
1033//
1034Bool_t MHMatrix::DefRefMatrix(Int_t nmaxevts, TMatrix *rest)
1035{
1036 if (!IsValid(fM))
1037 {
1038 *fLog << err << dbginf << "Matrix not initialized" << endl;
1039 return kFALSE;
1040 }
1041
1042 if (nmaxevts>fM.GetNrows())
1043 {
1044 *fLog << dbginf << "No.of requested events (" << nmaxevts
1045 << ") exceeds no.of available events (" << fM.GetNrows()
1046 << ")" << endl;
1047 *fLog << dbginf
1048 << " set no.of requested events = no.of available events"
1049 << endl;
1050 nmaxevts = fM.GetNrows();
1051 }
1052
1053 if (nmaxevts<0)
1054 {
1055 *fLog << err << dbginf << "Number of requested events < 0" << endl;
1056 return kFALSE;
1057 }
1058
1059 if (nmaxevts==0)
1060 nmaxevts = fM.GetNrows();
1061
1062 const Int_t nrows = fM.GetNrows();
1063 const Int_t ncols = fM.GetNcols();
1064
1065 //
1066 // get random access
1067 //
1068 TArrayI ind(nrows);
1069 GetRandomArrayI(ind);
1070
1071 //
1072 // define new matrix
1073 //
1074 Int_t evtcount1 = 0;
1075 Int_t evtcount2 = 0;
1076
1077 TMatrix mnewtmp(nrows, ncols);
1078 TMatrix mrest(nrows, ncols);
1079
1080 //
1081 // select events (distribution after renormalization)
1082 //
1083#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1084 TVector vold(fM.GetNcols());
1085#endif
1086 for (Int_t ir=0; ir<nmaxevts; ir++)
1087#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1088 TMatrixRow(mnewtmp, evtcount1++) = vold = TMatrixRow(fM, ind[ir]);
1089#else
1090 TMatrixFRow(mnewtmp, evtcount1++) = TMatrixFRow_const(fM, ind[ir]);
1091#endif
1092
1093 for (Int_t ir=nmaxevts; ir<nrows; ir++)
1094#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1095 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
1096#else
1097 TMatrixFRow(mrest, evtcount2++) = TMatrixFRow_const(fM, ind[ir]);
1098#endif
1099
1100 //
1101 // reduce size
1102 //
1103 // matrix fM having the requested distribution
1104 // and the requested number of rows;
1105 // this is the matrix to be used in the g/h separation
1106 //
1107 CopyCrop(fM, mnewtmp, evtcount1);
1108 fNumRows = evtcount1;
1109
1110 if (evtcount1 < nmaxevts)
1111 *fLog << warn << "The reference sample contains less events (" << evtcount1 << ") than requested (" << nmaxevts << ")" << endl;
1112
1113 if (!rest)
1114 return kTRUE;
1115
1116 CopyCrop(*rest, mrest, evtcount2);
1117
1118 return kTRUE;
1119}
1120
1121// --------------------------------------------------------------------------
1122//
1123// overload TOject member function read
1124// in order to reset the name of the object read
1125//
1126Int_t MHMatrix::Read(const char *name)
1127{
1128 Int_t ret = TObject::Read(name);
1129 SetName(name);
1130
1131 return ret;
1132}
1133
1134// --------------------------------------------------------------------------
1135//
1136// Read the setup from a TEnv:
1137// Column0, Column1, Column2, ..., Column10, ..., Column100, ...
1138//
1139// Searching stops if the first key isn't found in the TEnv. Empty
1140// columns are not allowed
1141//
1142// eg.
1143// MHMatrix.Column0: cos(MMcEvt.fTelescopeTheta)
1144// MHMatrix.Column1: MHillasSrc.fAlpha
1145//
1146Int_t MHMatrix::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
1147{
1148 if (IsValid(fM))
1149 {
1150 *fLog << err << "ERROR - matrix is already in use. Can't add a new column from TEnv... skipped." << endl;
1151 return kERROR;
1152 }
1153
1154 if (TestBit(kIsLocked))
1155 {
1156 *fLog << err << "ERROR - matrix is locked. Can't add new column from TEnv... skipped." << endl;
1157 return kERROR;
1158 }
1159
1160 //
1161 // Search (beginning with 0) all keys
1162 //
1163 int i=0;
1164 while (1)
1165 {
1166 TString idx = "Column";
1167 idx += i;
1168
1169 // Output if print set to kTRUE
1170 if (!IsEnvDefined(env, prefix, idx, print))
1171 break;
1172
1173 // Try to get the file name
1174 TString name = GetEnvValue(env, prefix, idx, "");
1175 if (name.IsNull())
1176 {
1177 *fLog << warn << prefix+"."+idx << " empty." << endl;
1178 continue;
1179 }
1180
1181 if (i==0)
1182 {
1183 if (fData)
1184 {
1185 *fLog << inf << "Removing all existing columns in " << GetDescriptor() << endl;
1186 fData->Delete();
1187 }
1188 else
1189 {
1190 fData = new MDataArray;
1191 SetBit(kIsOwner);
1192 }
1193 }
1194
1195 fData->AddEntry(name);
1196 i++;
1197 }
1198
1199 return i!=0;
1200}
1201
1202// --------------------------------------------------------------------------
1203//
1204// ShuffleEvents. Shuffles the order of the matrix rows.
1205//
1206//
1207void MHMatrix::ShuffleRows(UInt_t seed)
1208{
1209 TRandom rnd(seed);
1210
1211 TVector v(fM.GetNcols());
1212#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1213 TVector tmp(fM.GetNcols());
1214#endif
1215 for (Int_t irow = 0; irow<fNumRows; irow++)
1216 {
1217 const Int_t jrow = rnd.Integer(fNumRows);
1218
1219#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1220 v = TMatrixRow(fM, irow);
1221 TMatrixRow(fM, irow) = tmp = TMatrixRow(fM, jrow);
1222 TMatrixRow(fM, jrow) = v;
1223#else
1224 v = TMatrixFRow_const(fM, irow);
1225 TMatrixFRow(fM, irow) = TMatrixFRow_const(fM, jrow);
1226 TMatrixFRow(fM, jrow) = v;
1227#endif
1228 }
1229
1230 *fLog << warn << GetDescriptor() << ": Attention! Matrix rows have been shuffled." << endl;
1231}
1232
1233// --------------------------------------------------------------------------
1234//
1235// Reduces the number of rows to the given number num by cutting out the
1236// last rows.
1237//
1238void MHMatrix::ReduceRows(UInt_t num)
1239{
1240 if ((Int_t)num>=fM.GetNrows())
1241 {
1242 *fLog << warn << GetDescriptor() << ": Warning - " << num << " >= rows=" << fM.GetNrows() << endl;
1243 return;
1244 }
1245
1246#if ROOT_VERSION_CODE < ROOT_VERSION(3,05,07)
1247 const TMatrix m(fM);
1248#endif
1249 fM.ResizeTo(num, fM.GetNcols());
1250
1251#if ROOT_VERSION_CODE < ROOT_VERSION(3,05,07)
1252 TVector tmp(fM.GetNcols());
1253 for (UInt_t irow=0; irow<num; irow++)
1254 TMatrixRow(fM, irow) = tmp = TMatrixRow(m, irow);
1255#endif
1256}
1257
1258// --------------------------------------------------------------------------
1259//
1260// Remove rows which contains numbers not fullfilling TMath::Finite
1261//
1262Bool_t MHMatrix::RemoveInvalidRows()
1263{
1264 TMatrix m(fM);
1265
1266 const Int_t ncol=fM.GetNcols();
1267
1268#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1269 TVector vold(ncol);
1270#endif
1271
1272 int irow=0;
1273
1274 for (int i=0; i<m.GetNrows(); i++)
1275 {
1276#if ROOT_VERSION_CODE > ROOT_VERSION(5,00,00)
1277 const TMatrixFRow_const &row = TMatrixFRow_const(m, i);
1278#else
1279 const TMatrixRow &row = TMatrixRow(m, i);
1280#endif
1281 // finite (-> math.h) checks for NaN as well as inf
1282 int jcol;
1283 for (jcol=0; jcol<ncol; jcol++)
1284 if (!TMath::Finite(row(jcol)))
1285 break;
1286
1287 if (jcol==ncol)
1288#if ROOT_VERSION_CODE < ROOT_VERSION(4,00,8)
1289 TMatrixRow(fM, irow++) = vold = row;
1290#else
1291 TMatrixFRow(fM, irow++) = row;
1292#endif
1293 else
1294 *fLog << warn << "Warning - MHMatrix::RemoveInvalidRows: row #" << i<< " removed." << endl;
1295 }
1296
1297 // Do not use ResizeTo (in older root versions this doesn't save the contents
1298 ReduceRows(irow);
1299
1300 return kTRUE;
1301}
1302
1303// --------------------------------------------------------------------------
1304//
1305// Returns the row pointed to by fNumRow into TVector v
1306//
1307void MHMatrix::GetRow(TVector &v) const
1308{
1309 Int_t ncols = fM.GetNcols();
1310
1311 v.ResizeTo(ncols);
1312
1313 while (ncols--)
1314 v(ncols) = fM(fRow, ncols);
1315}
Note: See TracBrowser for help on using the repository browser.