source: trunk/MagicSoft/Mars/mtemp/mucm/classes/MMyFindSuperCuts.cc@ 6977

Last change on this file since 6977 was 6317, checked in by marcos, 21 years ago
*** empty log message ***
File size: 38.4 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, 6/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Wolfgang Wittek, 7/2003 <mailto:wittek@mppmu.mpg.de>
20! Author(s): Marcos Lopez, 05/2004 <mailto:marcos@gae.ucm.es>
21!
22! Copyright: MAGIC Software Development, 2000-2003
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28// //
29// MMyFindSuperCuts //
30// //
31// Class for otimizing the parameters of the supercuts //
32// //
33// //
34// //
35/////////////////////////////////////////////////////////////////////////////
36#include "MMyFindSuperCuts.h"
37
38#include <math.h> // fabs
39
40#include <TFile.h>
41#include <TArrayD.h>
42#include <TMinuit.h>
43#include <TCanvas.h>
44#include <TStopwatch.h>
45#include <TVirtualFitter.h>
46
47#include "MBinning.h"
48#include "MContinue.h"
49#include "MMySuperCuts.h"
50#include "MMySuperCutsCalc.h"
51#include "MDataElement.h"
52#include "MDataMember.h"
53
54#include "MEvtLoop.h"
55#include "MFSelFinal.h"
56#include "MF.h"
57#include "MFEventSelector.h"
58#include "MFEventSelector2.h"
59#include "MFillH.h"
60//#include "MGeomCamCT1Daniel.h"
61#include "MFEventSelector.h"
62#include "MGeomCamMagic.h"
63#include "MH3.h"
64#include "MHSupercuts.h"
65#include "MHFindSignificance.h"
66#include "MHMatrix.h"
67#include "MHOnSubtraction.h"
68
69#include "MLog.h"
70#include "MLogManip.h"
71#include "MMatrixLoop.h"
72#include "MMinuitInterface.h"
73#include "MParList.h"
74#include "MProgressBar.h"
75#include "MReadMarsFile.h"
76#include "MReadTree.h"
77#include "MTaskList.h"
78
79#include "TVector2.h"
80#include "MSrcPosCam.h"
81#include "MHillasSrcCalc.h"
82#include "MSrcPosCalc.h"
83#include "MObservatory.h"
84
85ClassImp(MMyFindSuperCuts);
86
87using namespace std;
88
89
90//------------------------------------------------------------------------
91//
92// fcnSupercuts
93//
94// - calculates the quantity to be minimized (using TMinuit)
95//
96// - the quantity to be minimized is (-1)*significance of the gamma signal
97// in the alpha distribution (after cuts)
98//
99// - the parameters to be varied in the minimization are the cut parameters
100// (par)
101//
102// -------------------------
103// - Read the Matrix
104// - Fill a histogram with the alpha (6th column)
105// - Fit the alpha plot
106// - Get result of fit (significance)
107//
108static void fcnSuperCuts(Int_t &npar, Double_t *gin, Double_t &f,
109 Double_t *par, Int_t iflag)
110{
111
112 //
113 // Set by hand the paramters range
114 //
115 if ( par[0]<0 || par[0]>1.8 || par[4]>par[0] || par[4]<0 || par[0]<par[4] ||
116 par[2]>1 || par[3]>1 || par[4]>1 || par[1]>1 || par[6]>1 || par[7]>1 ||
117 par[8]<0 || par[8]>1.8 || par[12]>par[0] || par[12]<0 || par[8]<par[12] ||
118 par[13]<0 || par[13]>1. || par[14]<0 || par[14]>1 || par[15]<0 || par[15]>1 ||
119 par[16]<0 || par[16]>1.8 || par[20]<0 || par[20]>1.8 || par[16]<par[20] ||
120 // par[4]>0.18 ||
121 par[24]>200 || par[28]<-200 ||
122 par[32]>1 || par[36]<0 ||
123 par[40]>1
124 )
125 {
126 gLog << "Parameter out of limit" << endl;
127 f=1e10;
128 cout << f << endl;
129 return;
130 }
131
132 //-------------------------------------------------------------
133 // save pointer to the MINUIT object for optimizing the supercuts
134 // because it will be overwritten
135 // when fitting the alpha distribution in MHFindSignificance
136 TMinuit *savePointer = gMinuit;
137 //-------------------------------------------------------------
138
139
140 MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();
141
142 MParList *plistfcn = (MParList*)evtloopfcn->GetParList();
143
144 MMySuperCuts *super = (MMySuperCuts*)plistfcn->FindObject("MMySuperCuts");
145 if (!super)
146 {
147 gLog << "fcnSupercuts : MMySuperCuts object '" << "MMySuperCuts"
148 << "' not found... aborting" << endl;
149 return;
150 }
151
152 //
153 // transfer current parameter values to MMySuperCuts
154 //
155 // Attention : npar is the number of variable parameters
156 // not the total number of parameters
157 //
158 Double_t fMin, fEdm, fErrdef;
159 Int_t fNpari, fNparx, fIstat;
160 gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
161
162 super->SetParameters(TArrayD(fNparx, par));
163
164
165 //$$$$$$$$$$$$$$$$$$$$$
166 // for testing
167 gLog.Separator("Current params");
168 TArrayD checkparameters = super->GetParameters();
169 gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", "
170 << fNparx << endl;
171 // gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
172 // for (Int_t i=0; i<fNparx; i++)
173// {
174// gLog << i << ", " << par[i] << ", " << checkparameters[i] << endl;
175// }
176 super->Print();
177
178 //$$$$$$$$$$$$$$$$$$$$$
179
180
181 //
182 // EventLoop: Fill the alpha plot with these cuts
183 //
184 gLog.SetNullOutput(kTRUE);
185 evtloopfcn->Eventloop();
186
187 //tasklistfcn->PrintStatistics(0, kTRUE);
188
189
190
191 //
192 // Calculate significance
193 //
194 MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3");
195 if (!alpha)
196 return;
197
198 TH1 &alphaHist = alpha->GetHist();
199 alphaHist.SetName("alpha-fcnSupercuts");
200
201 //-------------------------------------------
202 // set Minuit pointer to zero in order not to destroy the TMinuit
203 // object for optimizing the supercuts
204 gMinuit = NULL;
205
206
207 //=================================================================
208 // fit alpha distribution to get the number of excess events and
209 // calculate significance of gamma signal in the alpha plot
210 //
211 const Double_t alphasig = 10.0;
212 const Double_t alphamin = 30.0;
213 const Double_t alphamax = 90.0;
214 const Int_t degree = 0;
215
216 Bool_t drawpoly;
217 Bool_t fitgauss;
218 TCanvas* c;
219 if (iflag == 3)
220 {
221 c= new TCanvas();
222 drawpoly = kTRUE;
223 fitgauss = kTRUE;
224 }
225 else
226 {
227 drawpoly = kFALSE;
228 fitgauss = kFALSE;
229 }
230 fitgauss = kTRUE;
231 const Bool_t print = kFALSE;
232
233 MHFindSignificance findsig;
234 //findsig.SetRebin(kTRUE);
235 findsig.SetRebin(kFALSE);
236 findsig.SetReduceDegree(kFALSE);
237
238
239
240 const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
241 alphasig, drawpoly, fitgauss, print);
242 gLog.SetNullOutput(kFALSE);
243
244 //=================================================================
245
246 // reset gMinuit to the MINUIT object for optimizing the supercuts
247 gMinuit = savePointer;
248 //-------------------------------------------
249
250 if (!rc)
251 {
252 gLog << "fcnSupercuts : FindSigma() failed" << endl;
253 f = 1.e10;
254 return;
255 }
256
257 // plot some quantities during the optimization
258 MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts");
259 if (plotsuper)
260 plotsuper->Fill(&findsig);
261
262
263 //
264 // Print info
265 //
266 gLog << "(Non, Nbg, Nex, S/N, sigma, Significane) = "
267 << findsig.GetNon() << " "
268 << findsig.GetNbg()<< " "
269 << findsig.GetNex() << " "
270 << findsig.GetNex()/ findsig.GetNbg()<< " "
271 << findsig.GetSigmaGauss() << " "
272 << findsig.GetSignificance() << endl;
273
274 const Double_t sigmagauss = findsig.GetSigmaGauss();
275 // if (sigmagauss>10)
276// {
277// gLog << "fcnSupercuts : Alpha plot too width = " << sigmagauss << endl;
278// f = 1.e10;
279// return;
280// }
281
282
283 //
284 // Get significance
285 //
286 const Double_t significance = findsig.GetSignificance();
287 //const Double_t ratio = findsig.GetNex()/findsig.GetNbg();
288
289 // f = significance>0 ? -significance : 0;
290// f = -TMath::Abs(significance)/sqrt(sigmagauss);
291// f = - significance*sqrt(findsig.GetNex())/sqrt(sigmagauss);
292 f = - significance*significance/sqrt(sigmagauss);
293
294 // f = - significance*findsig.GetProb();
295 cout << f << endl;
296
297// //
298// // optimize signal/background ratio
299// Double_t ratio = findsig.GetNbg()>0.0 ?
300// findsig.GetNex()/findsig.GetNbg() : 0.0;
301
302// //f = -ratio;
303// if(significance < 5)
304// f=1e10;
305
306// f=f*ratio*ratio;
307// cout << f << " " << significance << " " << ratio << endl;
308
309}
310
311
312
313// --------------------------------------------------------------------------
314//
315// Default constructor.
316//
317MMyFindSuperCuts::MMyFindSuperCuts(const char *name, const char *title)
318 : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL),fSizeCutLow(2000.0),fSizeCutUp(10000000), fOptimizationMode(0)
319{
320 fName = name ? name : "MMyFindSuperCuts";
321 fTitle = title ? title : "Optimizer of the supercuts";
322
323 //---------------------------
324 // camera geometry is needed for conversion mm ==> degree
325 //fCam = new MGeomCamCT1Daniel;
326 fCam = new MGeomCamMagic;
327
328 // matrices to contain the training/test samples
329 fMatrixTrain = new MHMatrix("MatrixTrain");
330 fMatrixTest = new MHMatrix("MatrixTest");
331
332 // objects of MMySuperCutsCalc to which these matrices are attached
333 fCalcHadTrain = new MMySuperCutsCalc("SupercutsCalcTrain");
334 fCalcHadTest = new MMySuperCutsCalc("SupercutsCalcTest");
335
336
337 // Define columns of matrices
338 *fLog << inf << "Init Mapping of Train Matrix" << endl;
339 fCalcHadTrain->InitMapping(fMatrixTrain);
340 *fLog << inf << "Init Mapping of Test Matrix" << endl;
341 fCalcHadTest->InitMapping(fMatrixTest);
342}
343
344
345// --------------------------------------------------------------------------
346//
347// Default destructor.
348//
349MMyFindSuperCuts::~MMyFindSuperCuts()
350{
351 delete fCam;
352 delete fMatrixTrain;
353 delete fMatrixTest;
354 delete fCalcHadTrain;
355 delete fCalcHadTest;
356}
357
358
359
360// --------------------------------------------------------------------------
361//
362// Define the matrix 'fMatrixTrain' for the training sample
363//
364// alltogether 'howmanytrain' events are read from file 'nametrain';
365// the events are selected according to a target distribution 'hreftrain'
366//
367//
368Bool_t MMyFindSuperCuts::DefineTrainMatrix(
369 const TString &filenametrain, MH3 &hreftrain,
370 const Int_t howmanytrain, const TString &outfiletrain)
371{
372 if (filenametrain.IsNull() || howmanytrain <= 0)
373 return kFALSE;
374
375 *fLog << endl
376 << "***************************************************************"
377 << endl
378 << " Filling training matrix from file: '" << filenametrain << "'"
379 << endl
380 << " Selecting: " << howmanytrain << " events";
381
382 if (!hreftrain.GetHist().GetEntries()==0)
383 {
384 *fLog << ", according to a distribution given by the MH3 object '"
385 << hreftrain.GetName() << "'" << endl;
386 }
387 else
388 {
389 *fLog << ", Randomly" << endl;
390 }
391
392 *fLog << "***************************************************************"
393 << endl;
394
395
396 //
397 // ParList
398 //
399 MParList plist;
400 MTaskList tlist;
401 plist.AddToList(&tlist);
402
403 // entries in MParList
404 plist.AddToList(fCam);
405 plist.AddToList(fMatrixTrain);
406
407
408 //
409 // TaskList
410 //
411 MReadTree read("Events", filenametrain);
412 //MReadTree read("Parameters", filenametrain);
413 read.DisableAutoScheme();
414
415 // MFEventSelector2 seltrain(hreftrain);
416// seltrain.SetNumMax(howmanytrain);
417// seltrain.SetName("selectTrain");
418
419 // Random selection of events
420 MFEventSelector seltrain;
421 seltrain.SetNumSelectEvts(howmanytrain);
422 seltrain.SetName("selectTrain");
423
424 MFillH filltrain(fMatrixTrain);
425 filltrain.SetFilter(&seltrain);
426 filltrain.SetName("fillMatrixTrain");
427
428// // ReCalculate hillass
429// MObservatory obs;
430// plist.AddToList(&obs);
431
432// MSrcPosCalc srccalc;
433// srccalc.SetPositionXY(3.46945e-17, 0.0487805);
434// MHillasSrcCalc hsrc;
435
436
437
438 // entries in MTaskList
439 tlist.AddToList(&read);
440 tlist.AddToList(&seltrain);
441 // tlist.AddToList(&srccalc);
442// tlist.AddToList(&hsrc);
443 tlist.AddToList(&filltrain);
444
445
446 //
447 // EventLoop
448 //
449 MProgressBar bar;
450 MEvtLoop evtloop;
451 evtloop.SetParList(&plist);
452 evtloop.SetName("EvtLoopMatrixTrain");
453 evtloop.SetProgressBar(&bar);
454
455 if (!evtloop.Eventloop())
456 return kFALSE;
457
458 // if (!evtloop.PreProcess())
459// return kFALSE;
460
461// Int_t i=0;
462// while(tlist.Process())
463// {
464// i++;
465// if (i>howmanytrain)
466// break;
467// }
468
469// evtloop.PostProcess();
470
471 tlist.PrintStatistics(0, kTRUE);
472
473
474 //
475 // Print
476 //
477 fMatrixTrain->Print("SizeCols");
478 Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
479 if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
480 {
481 *fLog << "MMyFindSuperCuts::DefineTrainMatrix; no.of generated events ("
482 << howmanygenerated
483 << ") is incompatible with the no.of requested events ("
484 << howmanytrain << ")" << endl;
485 }
486
487 *fLog << "training matrix was filled" << endl;
488 *fLog << "=============================================" << endl;
489
490
491
492 //
493 // Write out training matrix
494 //
495 if (outfiletrain != "")
496 {
497 TFile filetr(outfiletrain, "RECREATE");
498 fMatrixTrain->Write();
499 filetr.Close();
500
501 *fLog << "MMyFindSuperCuts::DefineTrainMatrix; Training matrix was written onto file '"
502 << outfiletrain << "'" << endl;
503 }
504
505 // fMatrixTrain->Print("data");
506
507 return kTRUE;
508}
509
510// --------------------------------------------------------------------------
511//
512// Define the matrix for the test sample
513//
514// alltogether 'howmanytest' events are read from file 'nametest'
515//
516// the events are selected according to a target distribution 'hreftest'
517//
518//
519Bool_t MMyFindSuperCuts::DefineTestMatrix(
520 const TString &nametest, MH3 &hreftest,
521 const Int_t howmanytest, const TString &filetest)
522{
523 if (nametest.IsNull() || howmanytest<=0)
524 return kFALSE;
525
526 *fLog << "=============================================" << endl;
527 *fLog << "fill test matrix from file '" << nametest
528 << "', select " << howmanytest
529 << " events " << endl;
530 if (!hreftest.GetHist().GetEntries()==0)
531 {
532 *fLog << " according to a distribution given by the MH3 object '"
533 << hreftest.GetName() << "'" << endl;
534 }
535 else
536 {
537 *fLog << " randomly" << endl;
538 }
539
540
541 //
542 // ParList
543 //
544 MParList plist;
545 MTaskList tlist;
546
547 // entries in MParList
548 plist.AddToList(&tlist);
549 plist.AddToList(fCam);
550 plist.AddToList(fMatrixTest);
551
552
553 //
554 // TaskList
555 //
556 MReadTree read("Parameters", nametest);
557 read.DisableAutoScheme();
558
559 // MFEventSelector2 seltest(hreftest);
560// seltest.SetNumMax(howmanytest);
561// seltest.SetName("selectTest");
562
563 MFillH filltest(fMatrixTest);
564 //filltest.SetFilter(&seltest);
565
566 // entries in MTaskList
567 tlist.AddToList(&read);
568 //tlist.AddToList(&seltest);
569 tlist.AddToList(&filltest);
570
571
572 //
573 // EventLoop
574 //
575 MProgressBar bar;
576 MEvtLoop evtloop;
577 evtloop.SetParList(&plist);
578 evtloop.SetName("EvtLoopMatrixTest");
579 evtloop.SetProgressBar(&bar);
580
581
582 if (!evtloop.PreProcess())
583 return kFALSE;
584
585 Int_t i=0;
586 while(tlist.Process())
587 {
588 i++;
589 if (i>100000)
590 break;
591 }
592
593 evtloop.PostProcess();
594
595
596// if (!evtloop.Eventloop())
597// return kFALSE;
598
599 tlist.PrintStatistics(0, kTRUE);
600
601
602 //
603 // Print
604 //
605 fMatrixTest->Print("SizeCols");
606 const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
607 if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
608 {
609 *fLog << "MMyFindSuperCuts::DefineTestMatrix; no.of generated events ("
610 << howmanygenerated
611 << ") is incompatible with the no.of requested events ("
612 << howmanytest << ")" << endl;
613 }
614
615 *fLog << "test matrix was filled" << endl;
616 *fLog << "=============================================" << endl;
617
618
619
620 //
621 // Write out test matrix
622 //
623 if (filetest != "")
624 {
625 TFile filete(filetest, "RECREATE", "");
626 fMatrixTest->Write();
627 filete.Close();
628
629 *fLog << "MMyFindSuperCuts::DefineTestMatrix; Test matrix was written onto file '"
630 << filetest << "'" << endl;
631 }
632
633
634 return kTRUE;
635}
636
637
638
639// --------------------------------------------------------------------------
640//
641// Define the matrices for the training and test sample respectively
642//
643//
644//
645Bool_t MMyFindSuperCuts::DefineTrainTestMatrix(
646 const TString &name, MH3 &href,
647 const Int_t howmanytrain, const Int_t howmanytest,
648 const TString &filetrain, const TString &filetest)
649{
650 *fLog << "=============================================" << endl;
651 *fLog << "fill training and test matrix from file '" << name
652 << "', select " << howmanytrain
653 << " training and " << howmanytest << " test events " << endl;
654 if (!href.GetHist().GetEntries()==0)
655 {
656 *fLog << " according to a distribution given by the MH3 object '"
657 << href.GetName() << "'" << endl;
658 }
659 else
660 {
661 *fLog << " randomly" << endl;
662 }
663
664
665 //
666 // ParList
667 //
668 MParList plist;
669 MTaskList tlist;
670
671 // entries in MParList
672 plist.AddToList(&tlist);
673 plist.AddToList(fCam);
674 plist.AddToList(fMatrixTrain);
675 plist.AddToList(fMatrixTest);
676
677
678
679 //
680 // TaskList
681 //
682 MReadMarsFile read("Events", name);
683 read.DisableAutoScheme();
684
685 MFEventSelector2 selector(href);
686 selector.SetNumMax(howmanytrain+howmanytest);
687 selector.SetName("selectTrainTest");
688 selector.SetInverted();
689
690 MContinue cont(&selector);
691 cont.SetName("ContTrainTest");
692
693 Double_t prob = ( (Double_t) howmanytrain )
694 / ( (Double_t)(howmanytrain+howmanytest) );
695 MFEventSelector split;
696 split.SetSelectionRatio(prob);
697
698 MFillH filltrain(fMatrixTrain);
699 filltrain.SetFilter(&split);
700 filltrain.SetName("fillMatrixTrain");
701
702
703 // consider this event as candidate for a test event
704 // only if event was not accepted as a training event
705
706 MContinue conttrain(&split);
707 conttrain.SetName("ContTrain");
708
709 MFillH filltest(fMatrixTest);
710 filltest.SetName("fillMatrixTest");
711
712
713 // entries in MTaskList
714 tlist.AddToList(&read);
715 tlist.AddToList(&cont);
716 tlist.AddToList(&split);
717 tlist.AddToList(&filltrain);
718 tlist.AddToList(&conttrain);
719 tlist.AddToList(&filltest);
720
721
722 //
723 // EventLoop
724 //
725 MProgressBar bar;
726 MEvtLoop evtloop;
727 evtloop.SetParList(&plist);
728 evtloop.SetName("EvtLoopMatrixTrainTest");
729 evtloop.SetProgressBar(&bar);
730
731 Int_t maxev = -1;
732 if (!evtloop.Eventloop(maxev))
733 return kFALSE;
734
735 tlist.PrintStatistics(0, kTRUE);
736
737 fMatrixTrain->Print("SizeCols");
738 const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
739 if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
740 {
741 *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; no.of generated events ("
742 << generatedtrain
743 << ") is incompatible with the no.of requested events ("
744 << howmanytrain << ")" << endl;
745 }
746
747 fMatrixTest->Print("SizeCols");
748 const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
749 if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
750 {
751 *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; no.of generated events ("
752 << generatedtest
753 << ") is incompatible with the no.of requested events ("
754 << howmanytest << ")" << endl;
755 }
756
757
758 *fLog << "training and test matrix were filled" << endl;
759 *fLog << "=============================================" << endl;
760
761
762
763 //
764 // Write out training matrix
765 //
766 if (filetrain != "")
767 {
768 TFile filetr(filetrain, "RECREATE");
769 fMatrixTrain->Write();
770 filetr.Close();
771
772 *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Training matrix was written onto file '"
773 << filetrain << "'" << endl;
774 }
775
776 //
777 // Write out test matrix
778 //
779 if (filetest != "")
780 {
781 TFile filete(filetest, "RECREATE", "");
782 fMatrixTest->Write();
783 filete.Close();
784
785 *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Test matrix was written onto file '"
786 << filetest << "'" << endl;
787 }
788
789 return kTRUE;
790}
791
792
793
794// --------------------------------------------------------------------------
795//
796// Read training and test matrices from files
797//
798Bool_t MMyFindSuperCuts::ReadMatrix(const TString &filetrain, const TString &filetest)
799{
800 //
801 // Read in training matrix
802 //
803 TFile filetr(filetrain);
804 fMatrixTrain->Read("MatrixTrain");
805 fMatrixTrain->Print("SizeCols");
806
807 *fLog << "MMyFindSuperCuts::ReadMatrix; Training matrix was read in from file '" << filetrain << "'" << endl;
808 filetr.Close();
809
810
811 //
812 // Read in test matrix
813 //
814 TFile filete(filetest);
815 fMatrixTest->Read("MatrixTest");
816 fMatrixTest->Print("SizeCols");
817
818 *fLog <<"MMyFindSuperCuts::ReadMatrix; Test matrix was read in from file '"
819 << filetest << "'" << endl;
820 filete.Close();
821
822
823 return kTRUE;
824}
825
826
827
828//------------------------------------------------------------------------
829//
830// Steering program for optimizing the supercuts
831// ---------------------------------------------
832//
833// the criterion for the 'optimum' is
834//
835// - a maximum significance of the gamma signal in the alpha plot,
836// in which the supercuts have been applied
837//
838// The various steps are :
839//
840// - setup the event loop to be executed for each call to fcnSupercuts
841// - call TMinuit to do the minimization of (-significance) :
842// the fcnSupercuts function calculates the significance
843// for the current cut values
844// for this - the alpha plot is produced in the event loop,
845// in which the cuts have been applied
846// - the significance of the gamma signal in the alpha plot
847// is calculated
848//
849// Needed as input : (to be set by the Set functions)
850//
851// - fFilenameParam name of file to which optimum values of the
852// parameters are written
853//
854// - fHadronnessName name of container where MMySuperCutsCalc will
855// put the supercuts hadronness
856//
857// - for the minimization, the starting values of the parameters are taken
858// - from file parSCinit (if it is != "")
859// - or from the arrays params and/or steps
860//
861//----------------------------------------------------------------------
862Bool_t MMyFindSuperCuts::FindParams(TString parSCinit,
863 TArrayD &params, TArrayD &steps )
864{
865
866 fCalcHadTrain->SetSizeCuts(fSizeCutLow,fSizeCutUp);
867 fCalcHadTest->SetSizeCuts(fSizeCutLow,fSizeCutUp);
868
869
870 // Setup the event loop which will be executed in the
871 // fcnSupercuts function of MINUIT
872 //
873 // parSCinit is the name of the file containing the initial values
874 // of the parameters;
875 // if parSCinit = "" 'params' and 'steps' are taken as initial values
876
877 *fLog << "=============================================" << endl;
878 *fLog << "Setup event loop for fcnSupercuts" << endl;
879
880
881 //
882 // Check that all the necessary containers already exists
883 //
884 if (fHadronnessName.IsNull())
885 {
886 *fLog << "MyFindSuperCuts::FindParams; hadronness name is not defined... aborting"
887 << endl;
888 return kFALSE;
889 }
890
891 if (fMatrixTrain == NULL)
892 {
893 *fLog << "MMyFindSuperCuts::FindParams; training matrix is not defined... aborting"
894 << endl;
895 return kFALSE;
896 }
897
898 if (fMatrixTrain->GetM().GetNrows() <= 0)
899 {
900 *fLog << "MMyFindSuperCuts::FindParams; training matrix has no entries"
901 << endl;
902 return kFALSE;
903 }
904
905
906
907 //
908 // ParList
909 //
910 MParList parlistfcn;
911 MTaskList tasklistfcn;
912 parlistfcn.AddToList(&tasklistfcn);
913
914
915
916 //
917 // create container for the supercut parameters
918 // and set them to their initial values
919 //
920 MMySuperCuts super;
921
922 // take initial values from file parSCinit
923 if (parSCinit != "")
924 {
925 TFile inparam(parSCinit);
926 super.Read("MMySuperCuts");
927 inparam.Close();
928 *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from file " << parSCinit << endl;
929 super.Print();
930 }
931
932 // take initial values from 'params' and/or 'steps'
933 //else if (params.GetSize() != 0 || steps.GetSize() != 0 )
934 if (params.GetSize() != 0 || steps.GetSize() != 0 )
935 {
936 if (params.GetSize() != 0)
937 {
938 *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from 'params'" << endl;
939 super.SetParameters(params);
940 }
941 if (steps.GetSize() != 0)
942 {
943 *fLog << "MMyFindSuperCuts::FindParams; initial step sizes are taken from 'steps'"
944 << endl;
945 super.SetStepsizes(steps);
946 }
947 }
948 else
949 {
950 *fLog << "MMyFindSuperCuts::FindParams; initial values and step sizes are taken from the MMySuperCuts constructor"
951 << endl;
952 }
953
954
955 // Alpha binning
956 const TString mh3Name = "AlphaFcn";
957 MBinning binsalpha("Binning"+mh3Name);
958 // binsalpha.SetEdges(54, -12.0, 96.0);
959 binsalpha.SetEdges(45, 0.0, 90.0);
960
961 // Alpha plot |alpha|
962 MH3 alpha("MatrixTrain[6]"); // |alpha| assumed to be in column 6
963 alpha.SetName(mh3Name);
964
965 *fLog << warn << "WARNING----->ALPHA IS ASSUMED TO BE IN COLUMN 6!!!!!!!!"
966 << endl;
967
968
969 // book histograms to be filled during the optimization
970 // ! not in the event loop !
971 MHSupercuts plotsuper;
972 plotsuper.SetupFill(&parlistfcn);
973
974
975 // Add to ParList
976 parlistfcn.AddToList(&super);
977 parlistfcn.AddToList(fCam);
978 parlistfcn.AddToList(fMatrixTrain);
979 parlistfcn.AddToList(&binsalpha);
980 parlistfcn.AddToList(&alpha);
981 parlistfcn.AddToList(&plotsuper);
982
983
984
985 //
986 // MTaskList
987 //
988 MMatrixLoop loop(fMatrixTrain); // loop over rows of matrix
989
990 // calculate supercuts hadronness
991 //fCalcHadTrain->SetHadronnessName(fHadronnessName);
992
993 // apply the supercuts
994 MF scfilter(fHadronnessName+".fHadronness>0.5");
995 MContinue supercuts(&scfilter);
996
997 // Fill Alpha Plot
998 MFillH fillalpha(&alpha);
999
1000
1001 // Add to TaskList
1002 tasklistfcn.AddToList(&loop);
1003 tasklistfcn.AddToList(fCalcHadTrain);
1004 tasklistfcn.AddToList(&supercuts);
1005 tasklistfcn.AddToList(&fillalpha);
1006
1007
1008
1009 //
1010 // EventLoop
1011 //
1012 MEvtLoop evtloopfcn("EvtLoopFCN");
1013 evtloopfcn.SetParList(&parlistfcn);
1014 *fLog << "Event loop for fcnSupercuts has been setup" << endl;
1015
1016
1017
1018
1019 //-----------------------------------------------------------------------
1020 //
1021 //---------- Start of minimization part --------------------
1022 //
1023 // Do the minimization with MINUIT
1024 //
1025 // Be careful: This is not thread safe
1026 //
1027 *fLog << "========================================================"<< endl;
1028 *fLog << "Start minimization for supercuts" << endl;
1029
1030
1031 // -------------------------------------------
1032 // prepare call to MINUIT
1033 //
1034
1035 // get initial values of parameters
1036 fVinit = super.GetParameters();
1037 fStep = super.GetStepsizes();
1038
1039 TString name[fVinit.GetSize()];
1040 fStep.Set(fVinit.GetSize());
1041 fLimlo.Set(fVinit.GetSize());
1042 fLimup.Set(fVinit.GetSize());
1043 fFix.Set(fVinit.GetSize());
1044
1045 fNpar = fVinit.GetSize();
1046
1047
1048
1049 for (UInt_t i=0; i<fNpar; i++)
1050 {
1051 name[i] = "p";
1052 name[i] += i+1;
1053 //fStep[i] = TMath::Abs(fVinit[i]/10.0);
1054 // fLimlo[i] = -100.0;
1055// fLimup[i] = 100.0;
1056 fLimlo[i] = 0.0;
1057 fLimup[i] = 0.0;
1058 fFix[i] = 0;
1059 }
1060
1061 // these parameters make no sense, fix them at 0.0
1062 // for (UInt_t i=32; i<fNpar; i++)
1063// {
1064// //if( i!=3 && i!=7 && i!=11 && i!=15 && i!=19 && i!=23 )
1065// {
1066
1067// //fVinit[i] = 0.0;
1068// fStep[i] = 0.0;
1069// fFix[i] = 1;
1070// }
1071// // fVinit[3+4*i] = 0.0;
1072// // fStep[3+4*i] = 0.0;
1073// // fFix[3+4*i] = 1;
1074
1075// }
1076 // these parameters make no sense, fix them at 0.0
1077 // for (UInt_t i=0; i<24; i++)
1078// {
1079// //fVinit[i] = 0.0;
1080// fStep[i] = 0.0;
1081// fFix[i] = 1;
1082// }
1083
1084
1085 //-----------------------------------------------------------------
1086 //
1087 // 0=supercuts, 1=dyn from supercuts, 2=supercuts, 3=all
1088 //
1089 gLog << "*********************************" << fOptimizationMode << endl;
1090
1091 switch (fOptimizationMode)
1092 {
1093 case 0:
1094 for (UInt_t i=0; i<fNpar; i++)
1095 {
1096 if(i%4==0) // Only optimaze par0, par4, par8,...
1097 continue;
1098
1099 fStep[i] = 0.0;
1100 fFix[i] = 1;
1101 }
1102 break;
1103
1104 case 1:
1105 for (UInt_t i=0; i<fNpar; i++) // Optimize all but par0, par4, par8
1106 {
1107 if(i%4!=0 && i!=19 && i!=23)
1108 continue;
1109
1110 fStep[i] = 0.0;
1111 fFix[i] = 1;
1112
1113 }
1114
1115 break;
1116 }
1117
1118 //
1119 // DistLo/Up doens't dipend on Dist2
1120 //
1121 fStep[19] = 0.0;
1122 fFix[19] = 1;
1123 fStep[23] = 0.0;
1124 fFix[23] = 1;
1125
1126
1127
1128 //
1129 // Fix all parameters from 24 to the end
1130 //
1131 for (UInt_t i=24; i<fNpar; i++)
1132 {
1133 //fVinit[i] = 0.0;
1134 fStep[i] = 0.0;
1135 fFix[i] = 1;
1136 }
1137
1138
1139 // // Static Cuts
1140// //
1141// for (UInt_t i=0; i<fNpar; i++)
1142// {
1143// if( i!=0 && i!=4 && i!=8 && i!=12 && i!=16 && i!=20 &&
1144// i!=24 && i!=28 && i!=32 && i!=36 && i!=40 )
1145// {
1146// fStep[i] = 0.0;
1147// fFix[i] = 1;
1148// }
1149// }
1150
1151 //
1152 // //
1153 //
1154 // for (UInt_t i=0; i<fNpar; i++)
1155// {
1156// if( i==0 || i==4 || i==8 || i==12 || i==16 || i==20)
1157// {
1158// //fVinit[i] = 0.0;
1159// fStep[i] = 0.0;
1160// fFix[i] = 1;
1161// }
1162// }
1163
1164
1165 // Fix
1166 // for (UInt_t i=0; i<36; i++)
1167// {
1168// fStep[i] = 0.0;
1169// fFix[i] = 1;
1170// }
1171 // for (UInt_t i=44; i<fNpar; i++)
1172// {
1173// fStep[i] = 0.0;
1174// fFix[i] = 1;
1175// }
1176
1177// for (UInt_t i=37; i<=39; i++)
1178// {
1179// //fVinit[i] = 0.0;
1180// fStep[i] = 0.0;
1181// fFix[i] = 1;
1182// }
1183
1184
1185 // fStep[40] = .01;
1186 // fStep[28] = 10;
1187 // fStep[24] = 10;
1188 // fVinit[24] = 100;
1189 // fVinit[28] = -74;
1190 // fVinit[32] = 0.16;
1191// fVinit[36] = 0.03;
1192 // fVinit[40] = 0.21;
1193// fVinit[44] = 0.0;
1194
1195
1196
1197 //fVinit[28] = 0;
1198 // fVinit[40] = 0.5;
1199
1200// fVinit[44] = -1e6;
1201 // fStep[44] = 1;
1202// fVinit[36] = 0.;
1203// //fStep[36] = 1.0;
1204
1205
1206 // vary only first 48 parameters
1207 //for (UInt_t i=0; i<fNpar; i++)
1208 //{
1209 // if (i >= 48)
1210 // {
1211 // fStep[i] = 0.0;
1212 // fFix[i] = 1;
1213 // }
1214 //
1215
1216 // //}
1217// // vary only first 23 parameters
1218// //
1219// for (UInt_t i=0; i<fNpar; i++)
1220// {
1221// if (i >= 24)
1222// {
1223// fStep[i] = 0.0;
1224// fFix[i] = 1;
1225// }
1226// }
1227
1228
1229 // -------------------------------------------
1230 // call MINUIT
1231
1232 TStopwatch clock;
1233 clock.Start();
1234
1235 *fLog << "before calling CallMinuit" << endl;
1236
1237 MMinuitInterface inter;
1238 Bool_t rc = inter.CallMinuit(fcnSuperCuts, name,
1239 fVinit, fStep, fLimlo, fLimup, fFix,
1240 &evtloopfcn, "SIMPLEX", kFALSE);
1241
1242 *fLog << "after calling CallMinuit" << endl;
1243
1244 *fLog << "Time spent for the minimization in MINUIT : " << endl;;
1245 clock.Stop();
1246 clock.Print();
1247
1248 plotsuper.DrawClone();
1249
1250 if (!rc)
1251 return kFALSE;
1252
1253 *fLog << "Minimization for supercuts finished" << endl;
1254 *fLog << "========================================================" << endl;
1255
1256
1257
1258 // -----------------------------------------------------------------
1259 // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put
1260 // into MMySuperCuts
1261
1262 //
1263 // Write optimum parameter values onto file filenameParam
1264 //
1265 TFile outparam(fFilenameParam, "RECREATE");
1266 super.Write();
1267 outparam.Close();
1268
1269 *fLog << "Optimum parameter values for supercuts were written onto file '"
1270 << fFilenameParam << "' :" << endl;
1271
1272 const TArrayD &check = super.GetParameters();
1273 for (Int_t i=0; i<check.GetSize(); i++)
1274 *fLog << check[i] << ", ";
1275 *fLog << endl;
1276
1277
1278
1279 *fLog << "End of supercuts optimization part" << endl;
1280 *fLog << "======================================================" << endl;
1281
1282 return kTRUE;
1283}
1284
1285
1286// -----------------------------------------------------------------------
1287//
1288// Test the supercuts on the test sample
1289//
1290Bool_t MMyFindSuperCuts::TestParams()
1291{
1292 if (fMatrixTest->GetM().GetNrows() <= 0)
1293 {
1294 *fLog << "MMyFindSuperCuts::TestParams; test matrix has no entries"
1295 << endl;
1296 return kFALSE;
1297 }
1298
1299 // ------------- TEST the supercuts ------------------------------
1300 //
1301 *fLog << "Test the supercuts on the test sample" << endl;
1302
1303 // -----------------------------------------------------------------
1304 // read optimum parameter values from file filenameParam
1305 // into array 'supercutsPar'
1306
1307 TFile inparam(fFilenameParam);
1308 MMySuperCuts scin;
1309 scin.Read("MMySuperCuts");
1310 inparam.Close();
1311
1312 *fLog << "Optimum parameter values for supercuts were read from file '";
1313 *fLog << fFilenameParam << "' :" << endl;
1314
1315 const TArrayD &supercutsPar = scin.GetParameters();
1316 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1317 *fLog << supercutsPar[i] << ", ";
1318 *fLog << endl;
1319 //---------------------------
1320
1321
1322 // -----------------------------------------------------------------
1323 if (fHadronnessName.IsNull())
1324 {
1325 *fLog << "MMyFindSuperCuts::TestParams; hadronness name is not defined... aborting";
1326 *fLog << endl;
1327 return kFALSE;
1328 }
1329
1330 MParList parlist2;
1331 MTaskList tasklist2;
1332
1333 MMySuperCuts supercuts;
1334 supercuts.SetParameters(supercutsPar);
1335
1336 fCalcHadTest->SetHadronnessName(fHadronnessName);
1337
1338
1339 //-------------------------------------------
1340
1341 MMatrixLoop loop(fMatrixTest);
1342
1343 // plot alpha before applying the supercuts
1344 const TString mh3NameB = "AlphaBefore";
1345 MBinning binsalphabef("Binning"+mh3NameB);
1346 binsalphabef.SetEdges(54, -12.0, 96.0);
1347
1348 // |alpha| is assumed to be in column 7 of the matrix
1349 MH3 alphabefore("MatrixTest[6]");
1350 alphabefore.SetName(mh3NameB);
1351
1352 TH1 &alphahistb = alphabefore.GetHist();
1353 alphahistb.SetName("alphaBefore-TestParams");
1354
1355 MFillH fillalphabefore(&alphabefore);
1356 fillalphabefore.SetName("FillAlphaBefore");
1357
1358 // apply the supercuts
1359 MF scfilter(fHadronnessName+".fHadronness>0.5");
1360 MContinue applysupercuts(&scfilter);
1361
1362 // plot alpha after applying the supercuts
1363 const TString mh3NameA = "AlphaAfter";
1364 MBinning binsalphaaft("Binning"+mh3NameA);
1365 binsalphaaft.SetEdges(54, -12.0, 96.0);
1366
1367 MH3 alphaafter("MatrixTest[6]");
1368 alphaafter.SetName(mh3NameA);
1369
1370 TH1 &alphahista = alphaafter.GetHist();
1371 alphahista.SetName("alphaAfter-TestParams");
1372
1373 MFillH fillalphaafter(&alphaafter);
1374 fillalphaafter.SetName("FillAlphaAfter");
1375
1376 //******************************
1377 // entries in MParList
1378
1379 parlist2.AddToList(&tasklist2);
1380
1381 parlist2.AddToList(&supercuts);
1382
1383 parlist2.AddToList(fCam);
1384 parlist2.AddToList(fMatrixTest);
1385
1386 parlist2.AddToList(&binsalphabef);
1387 parlist2.AddToList(&alphabefore);
1388
1389 parlist2.AddToList(&binsalphaaft);
1390 parlist2.AddToList(&alphaafter);
1391
1392 //******************************
1393 // entries in MTaskList
1394
1395 tasklist2.AddToList(&loop);
1396 tasklist2.AddToList(&fillalphabefore);
1397
1398 tasklist2.AddToList(fCalcHadTest);
1399 tasklist2.AddToList(&applysupercuts);
1400
1401 tasklist2.AddToList(&fillalphaafter);
1402
1403 //******************************
1404
1405 MProgressBar bar2;
1406 MEvtLoop evtloop2;
1407 evtloop2.SetParList(&parlist2);
1408 evtloop2.SetName("EvtLoopTestParams");
1409 evtloop2.SetProgressBar(&bar2);
1410 Int_t maxevents2 = -1;
1411 if (!evtloop2.Eventloop(maxevents2))
1412 return kFALSE;
1413
1414 tasklist2.PrintStatistics(0, kTRUE);
1415
1416
1417 //-------------------------------------------
1418 // draw the alpha plots
1419
1420 MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
1421 TH1 &alphaHist1 = alphaBefore->GetHist();
1422 alphaHist1.SetXTitle("|\\alpha| [\\circ]");
1423
1424 MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
1425 TH1 &alphaHist2 = alphaAfter->GetHist();
1426 alphaHist2.SetXTitle("|\\alpha| [\\circ]");
1427 alphaHist2.SetName("alpha-TestParams");
1428
1429 TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
1430 c->Divide(2,1);
1431
1432 gROOT->SetSelectedPad(NULL);
1433
1434 c->cd(1);
1435 alphaHist1.DrawCopy();
1436
1437 c->cd(2);
1438 alphaHist2.DrawCopy();
1439
1440
1441
1442 //-------------------------------------------
1443 // fit alpha distribution to get the number of excess events and
1444 // calculate significance of gamma signal in the alpha plot
1445
1446 const Double_t alphasig = 20.0;
1447 const Double_t alphamin = 30.0;
1448 const Double_t alphamax = 90.0;
1449 const Int_t degree = 2;
1450 const Bool_t drawpoly = kTRUE;
1451 const Bool_t fitgauss = kTRUE;
1452 const Bool_t print = kTRUE;
1453
1454 MHFindSignificance findsig;
1455 findsig.SetRebin(kTRUE);
1456 findsig.SetReduceDegree(kFALSE);
1457
1458 findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree,
1459 alphasig, drawpoly, fitgauss, print);
1460
1461 const Double_t significance = findsig.GetSignificance();
1462 const Double_t alphasi = findsig.GetAlphasi();
1463
1464 *fLog << "Significance of gamma signal after supercuts in test sample : "
1465 << significance << " (for |alpha| < " << alphasi << " degrees)"
1466 << endl;
1467 //-------------------------------------------
1468
1469
1470 *fLog << "Test of supercuts part finished" << endl;
1471 *fLog << "======================================================" << endl;
1472
1473 return kTRUE;
1474}
1475
Note: See TracBrowser for help on using the repository browser.