source: trunk/MagicSoft/Mars/manalysisct1/MCT1FindSupercuts.cc@ 8123

Last change on this file since 8123 was 4457, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 32.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, 6/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Wolfgang Wittek, 7/2003 <mailto:wittek@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MCT1FindSupercuts //
29// //
30// Class for otimizing the parameters of the supercuts //
31// //
32// //
33// //
34/////////////////////////////////////////////////////////////////////////////
35#include "MCT1FindSupercuts.h"
36
37#include <math.h> // fabs
38
39#include <TFile.h>
40#include <TArrayD.h>
41#include <TMinuit.h>
42#include <TCanvas.h>
43#include <TStopwatch.h>
44#include <TVirtualFitter.h>
45
46#include "MBinning.h"
47#include "MContinue.h"
48#include "MCT1Supercuts.h"
49#include "MCT1SupercutsCalc.h"
50#include "MDataElement.h"
51#include "MDataMember.h"
52
53#include "MEvtLoop.h"
54#include "MFCT1SelFinal.h"
55#include "MF.h"
56#include "MFEventSelector.h"
57#include "MFEventSelector2.h"
58#include "MFillH.h"
59//#include "MGeomCamCT1Daniel.h"
60#include "MFEventSelector.h"
61#include "MGeomCamCT1.h"
62#include "MH3.h"
63#include "MHCT1Supercuts.h"
64#include "MHFindSignificance.h"
65#include "MHMatrix.h"
66#include "MHOnSubtraction.h"
67
68#include "MLog.h"
69#include "MLogManip.h"
70#include "MMatrixLoop.h"
71#include "MMinuitInterface.h"
72#include "MParList.h"
73#include "MProgressBar.h"
74#include "MReadMarsFile.h"
75#include "MReadTree.h"
76#include "MTaskList.h"
77
78
79ClassImp(MCT1FindSupercuts);
80
81using namespace std;
82
83
84//------------------------------------------------------------------------
85//
86// fcnSupercuts
87//
88// - calculates the quantity to be minimized (using TMinuit)
89//
90// - the quantity to be minimized is (-1)*significance of the gamma signal
91// in the alpha distribution (after cuts)
92//
93// - the parameters to be varied in the minimization are the cut parameters
94// (par)
95//
96static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f,
97 Double_t *par, Int_t iflag)
98{
99 //cout << "entry fcnSupercuts" << endl;
100
101 //-------------------------------------------------------------
102 // save pointer to the MINUIT object for optimizing the supercuts
103 // because it will be overwritten
104 // when fitting the alpha distribution in MHFindSignificance
105 TMinuit *savePointer = gMinuit;
106 //-------------------------------------------------------------
107
108
109 MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();
110
111 MParList *plistfcn = (MParList*)evtloopfcn->GetParList();
112 //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList");
113
114 MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts");
115 if (!super)
116 {
117 gLog << "fcnSupercuts : MCT1Supercuts object '" << "MCT1Supercuts"
118 << "' not found... aborting" << endl;
119 return;
120 }
121
122 //
123 // transfer current parameter values to MCT1Supercuts
124 //
125 // Attention : npar is the number of variable parameters
126 // not the total number of parameters
127 //
128 Double_t fMin, fEdm, fErrdef;
129 Int_t fNpari, fNparx, fIstat;
130 gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
131
132 super->SetParameters(TArrayD(fNparx, par));
133
134 //$$$$$$$$$$$$$$$$$$$$$
135 // for testing
136 //TArrayD checkparameters = super->GetParameters();
137 //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", "
138 // << fNparx << endl;
139 //gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
140 //for (Int_t i=0; i<fNparx; i++)
141 //{
142 // gLog << i << ", " << par[i] << ", " << checkparameters[i] << endl;
143 //}
144 //$$$$$$$$$$$$$$$$$$$$$
145
146 //
147 // plot alpha with the current cuts
148 //
149 evtloopfcn->Eventloop();
150
151 //tasklistfcn->PrintStatistics(0, kTRUE);
152
153 MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3");
154 if (!alpha)
155 return;
156
157 TH1 &alphaHist = alpha->GetHist();
158 alphaHist.SetName("alpha-fcnSupercuts");
159
160 //-------------------------------------------
161 // set Minuit pointer to zero in order not to destroy the TMinuit
162 // object for optimizing the supercuts
163 gMinuit = NULL;
164
165 //=================================================================
166 // fit alpha distribution to get the number of excess events and
167 // calculate significance of gamma signal in the alpha plot
168
169 const Double_t alphasig = 20.0;
170 const Double_t alphamin = 30.0;
171 const Double_t alphamax = 90.0;
172 const Int_t degree = 2;
173
174 Bool_t drawpoly;
175 Bool_t fitgauss;
176 if (iflag == 3)
177 {
178 drawpoly = kTRUE;
179 fitgauss = kTRUE;
180 }
181 else
182 {
183 drawpoly = kFALSE;
184 fitgauss = kFALSE;
185 }
186 //drawpoly = kFALSE;
187 //fitgauss = kFALSE;
188
189 const Bool_t print = kTRUE;
190
191 MHFindSignificance findsig;
192 findsig.SetRebin(kTRUE);
193 findsig.SetReduceDegree(kFALSE);
194
195 const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
196 alphasig, drawpoly, fitgauss, print);
197
198 //=================================================================
199
200 // reset gMinuit to the MINUIT object for optimizing the supercuts
201 gMinuit = savePointer;
202 //-------------------------------------------
203
204 if (!rc)
205 {
206 gLog << "fcnSupercuts : FindSigma() failed" << endl;
207 f = 1.e10;
208 return;
209 }
210
211 // plot some quantities during the optimization
212 MHCT1Supercuts *plotsuper = (MHCT1Supercuts*)plistfcn->FindObject("MHCT1Supercuts");
213 if (plotsuper)
214 plotsuper->Fill(&findsig);
215
216 //------------------------
217 // get significance
218 const Double_t significance = findsig.GetSignificance();
219 f = significance>0 ? -significance : 0;
220
221
222 //------------------------
223 // optimize signal/background ratio
224 //Double_t ratio = findsig.GetNbg()>0.0 ?
225 // findsig.GetNex()/findsig.GetNbg() : 0.0;
226 //f = -ratio;
227
228 //-------------------------------------------
229 // final calculations
230 //if (iflag == 3)
231 //{
232 //
233 //}
234
235 //-------------------------------------------------------------
236}
237
238
239// --------------------------------------------------------------------------
240//
241// Default constructor.
242//
243MCT1FindSupercuts::MCT1FindSupercuts(const char *name, const char *title)
244: fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL)
245{
246 fName = name ? name : "MCT1FindSupercuts";
247 fTitle = title ? title : "Optimizer of the supercuts";
248
249 //---------------------------
250 // camera geometry is needed for conversion mm ==> degree
251 //fCam = new MGeomCamCT1Daniel;
252 fCam = new MGeomCamCT1;
253
254 // matrices to contain the training/test samples
255 fMatrixTrain = new MHMatrix("MatrixTrain");
256 fMatrixTest = new MHMatrix("MatrixTest");
257
258 // objects of MCT1SupercutsCalc to which these matrices are attached
259 fCalcHadTrain = new MCT1SupercutsCalc("SupercutsCalcTrain");
260 fCalcHadTest = new MCT1SupercutsCalc("SupercutsCalcTest");
261
262 // Define columns of matrices
263 fCalcHadTrain->InitMapping(fMatrixTrain);
264 fCalcHadTest->InitMapping(fMatrixTest);
265}
266
267// --------------------------------------------------------------------------
268//
269// Default destructor.
270//
271MCT1FindSupercuts::~MCT1FindSupercuts()
272{
273 delete fCam;
274 delete fMatrixTrain;
275 delete fMatrixTest;
276 delete fCalcHadTrain;
277 delete fCalcHadTest;
278}
279
280// --------------------------------------------------------------------------
281//
282// Define the matrix 'fMatrixTrain' for the training sample
283//
284// alltogether 'howmanytrain' events are read from file 'nametrain';
285// the events are selected according to a target distribution 'hreftrain'
286//
287//
288Bool_t MCT1FindSupercuts::DefineTrainMatrix(
289 const TString &nametrain, MH3 &hreftrain,
290 const Int_t howmanytrain, const TString &filetrain)
291{
292 if (nametrain.IsNull() || howmanytrain <= 0)
293 return kFALSE;
294
295 *fLog << "=============================================" << endl;
296 *fLog << "fill training matrix from file '" << nametrain
297 << "', select " << howmanytrain
298 << " events " << endl;
299 if (!hreftrain.GetHist().GetEntries()==0)
300 {
301 *fLog << " according to a distribution given by the MH3 object '"
302 << hreftrain.GetName() << "'" << endl;
303 }
304 else
305 {
306 *fLog << " randomly" << endl;
307 }
308
309
310 MParList plist;
311 MTaskList tlist;
312
313 MReadMarsFile read("Events", nametrain);
314 read.DisableAutoScheme();
315
316 MFEventSelector2 seltrain(hreftrain);
317 seltrain.SetNumMax(howmanytrain);
318 seltrain.SetName("selectTrain");
319
320 MFillH filltrain(fMatrixTrain);
321 filltrain.SetFilter(&seltrain);
322 filltrain.SetName("fillMatrixTrain");
323
324 //******************************
325 // entries in MParList
326
327 plist.AddToList(&tlist);
328 plist.AddToList(fCam);
329 plist.AddToList(fMatrixTrain);
330
331 //******************************
332 // entries in MTaskList
333
334 tlist.AddToList(&read);
335 tlist.AddToList(&seltrain);
336 tlist.AddToList(&filltrain);
337
338 //******************************
339
340 MProgressBar bar;
341 MEvtLoop evtloop;
342 evtloop.SetParList(&plist);
343 evtloop.SetName("EvtLoopMatrixTrain");
344 evtloop.SetProgressBar(&bar);
345
346 if (!evtloop.Eventloop())
347 return kFALSE;
348
349 tlist.PrintStatistics(0, kTRUE);
350
351 fMatrixTrain->Print("SizeCols");
352 Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
353 if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
354 {
355 *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
356 << howmanygenerated
357 << ") is incompatible with the no.of requested events ("
358 << howmanytrain << ")" << endl;
359 }
360
361 *fLog << "training matrix was filled" << endl;
362 *fLog << "=============================================" << endl;
363
364 //--------------------------
365 // write out training matrix
366
367 if (filetrain != "")
368 {
369 TFile filetr(filetrain, "RECREATE");
370 fMatrixTrain->Write();
371 filetr.Close();
372
373 *fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
374 << filetrain << "'" << endl;
375 }
376
377
378 return kTRUE;
379}
380
381// --------------------------------------------------------------------------
382//
383// Define the matrix for the test sample
384//
385// alltogether 'howmanytest' events are read from file 'nametest'
386//
387// the events are selected according to a target distribution 'hreftest'
388//
389//
390Bool_t MCT1FindSupercuts::DefineTestMatrix(
391 const TString &nametest, MH3 &hreftest,
392 const Int_t howmanytest, const TString &filetest)
393{
394 if (nametest.IsNull() || howmanytest<=0)
395 return kFALSE;
396
397 *fLog << "=============================================" << endl;
398 *fLog << "fill test matrix from file '" << nametest
399 << "', select " << howmanytest
400 << " events " << endl;
401 if (!hreftest.GetHist().GetEntries()==0)
402 {
403 *fLog << " according to a distribution given by the MH3 object '"
404 << hreftest.GetName() << "'" << endl;
405 }
406 else
407 {
408 *fLog << " randomly" << endl;
409 }
410
411
412 MParList plist;
413 MTaskList tlist;
414
415 MReadMarsFile read("Events", nametest);
416 read.DisableAutoScheme();
417
418 MFEventSelector2 seltest(hreftest);
419 seltest.SetNumMax(howmanytest);
420 seltest.SetName("selectTest");
421
422 MFillH filltest(fMatrixTest);
423 filltest.SetFilter(&seltest);
424
425 //******************************
426 // entries in MParList
427
428 plist.AddToList(&tlist);
429 plist.AddToList(fCam);
430 plist.AddToList(fMatrixTest);
431
432 //******************************
433 // entries in MTaskList
434
435 tlist.AddToList(&read);
436 tlist.AddToList(&seltest);
437 tlist.AddToList(&filltest);
438
439 //******************************
440
441 MProgressBar bar;
442 MEvtLoop evtloop;
443 evtloop.SetParList(&plist);
444 evtloop.SetName("EvtLoopMatrixTest");
445 evtloop.SetProgressBar(&bar);
446
447 if (!evtloop.Eventloop())
448 return kFALSE;
449
450 tlist.PrintStatistics(0, kTRUE);
451
452 fMatrixTest->Print("SizeCols");
453 const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
454 if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
455 {
456 *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events ("
457 << howmanygenerated
458 << ") is incompatible with the no.of requested events ("
459 << howmanytest << ")" << endl;
460 }
461
462 *fLog << "test matrix was filled" << endl;
463 *fLog << "=============================================" << endl;
464
465 //--------------------------
466 // write out test matrix
467
468 if (filetest != "")
469 {
470 TFile filete(filetest, "RECREATE", "");
471 fMatrixTest->Write();
472 filete.Close();
473
474 *fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '"
475 << filetest << "'" << endl;
476 }
477
478
479 return kTRUE;
480}
481
482// --------------------------------------------------------------------------
483//
484// Define the matrices for the training and test sample respectively
485//
486//
487//
488Bool_t MCT1FindSupercuts::DefineTrainTestMatrix(
489 const TString &name, MH3 &href,
490 const Int_t howmanytrain, const Int_t howmanytest,
491 const TString &filetrain, const TString &filetest)
492{
493 *fLog << "=============================================" << endl;
494 *fLog << "fill training and test matrix from file '" << name
495 << "', select " << howmanytrain
496 << " training and " << howmanytest << " test events " << endl;
497 if (!href.GetHist().GetEntries()==0)
498 {
499 *fLog << " according to a distribution given by the MH3 object '"
500 << href.GetName() << "'" << endl;
501 }
502 else
503 {
504 *fLog << " randomly" << endl;
505 }
506
507
508 MParList plist;
509 MTaskList tlist;
510
511 MReadMarsFile read("Events", name);
512 read.DisableAutoScheme();
513
514 MFEventSelector2 selector(href);
515 selector.SetNumMax(howmanytrain+howmanytest);
516 selector.SetName("selectTrainTest");
517 selector.SetInverted();
518
519 MContinue cont(&selector);
520 cont.SetName("ContTrainTest");
521
522 Double_t prob = ( (Double_t) howmanytrain )
523 / ( (Double_t)(howmanytrain+howmanytest) );
524 MFEventSelector split;
525 split.SetSelectionRatio(prob);
526
527 MFillH filltrain(fMatrixTrain);
528 filltrain.SetFilter(&split);
529 filltrain.SetName("fillMatrixTrain");
530
531
532 // consider this event as candidate for a test event
533 // only if event was not accepted as a training event
534
535 MContinue conttrain(&split);
536 conttrain.SetName("ContTrain");
537
538 MFillH filltest(fMatrixTest);
539 filltest.SetName("fillMatrixTest");
540
541
542 //******************************
543 // entries in MParList
544
545 plist.AddToList(&tlist);
546 plist.AddToList(fCam);
547 plist.AddToList(fMatrixTrain);
548 plist.AddToList(fMatrixTest);
549
550 //******************************
551 // entries in MTaskList
552
553 tlist.AddToList(&read);
554 tlist.AddToList(&cont);
555
556 tlist.AddToList(&split);
557 tlist.AddToList(&filltrain);
558 tlist.AddToList(&conttrain);
559
560 tlist.AddToList(&filltest);
561
562 //******************************
563
564 MProgressBar bar;
565 MEvtLoop evtloop;
566 evtloop.SetParList(&plist);
567 evtloop.SetName("EvtLoopMatrixTrainTest");
568 evtloop.SetProgressBar(&bar);
569
570 Int_t maxev = -1;
571 if (!evtloop.Eventloop(maxev))
572 return kFALSE;
573
574 tlist.PrintStatistics(0, kTRUE);
575
576 fMatrixTrain->Print("SizeCols");
577 const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
578 if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
579 {
580 *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
581 << generatedtrain
582 << ") is incompatible with the no.of requested events ("
583 << howmanytrain << ")" << endl;
584 }
585
586 fMatrixTest->Print("SizeCols");
587 const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
588 if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
589 {
590 *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
591 << generatedtest
592 << ") is incompatible with the no.of requested events ("
593 << howmanytest << ")" << endl;
594 }
595
596
597 *fLog << "training and test matrix were filled" << endl;
598 *fLog << "=============================================" << endl;
599
600
601 //--------------------------
602 // write out training matrix
603
604 if (filetrain != "")
605 {
606 TFile filetr(filetrain, "RECREATE");
607 fMatrixTrain->Write();
608 filetr.Close();
609
610 *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
611 << filetrain << "'" << endl;
612 }
613
614 //--------------------------
615 // write out test matrix
616
617 if (filetest != "")
618 {
619 TFile filete(filetest, "RECREATE", "");
620 fMatrixTest->Write();
621 filete.Close();
622
623 *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
624 << filetest << "'" << endl;
625 }
626
627 return kTRUE;
628}
629
630// --------------------------------------------------------------------------
631//
632// Read training and test matrices from files
633//
634//
635
636Bool_t MCT1FindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest)
637{
638 //--------------------------
639 // read in training matrix
640
641 TFile filetr(filetrain);
642 fMatrixTrain->Read("MatrixTrain");
643 fMatrixTrain->Print("SizeCols");
644
645 *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
646 << filetrain << "'" << endl;
647 filetr.Close();
648
649
650 //--------------------------
651 // read in test matrix
652
653 TFile filete(filetest);
654 fMatrixTest->Read("MatrixTest");
655 fMatrixTest->Print("SizeCols");
656
657 *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
658 << filetest << "'" << endl;
659 filete.Close();
660
661 return kTRUE;
662}
663
664
665//------------------------------------------------------------------------
666//
667// Steering program for optimizing the supercuts
668// ---------------------------------------------
669//
670// the criterion for the 'optimum' is
671//
672// - a maximum significance of the gamma signal in the alpha plot,
673// in which the supercuts have been applied
674//
675// The various steps are :
676//
677// - setup the event loop to be executed for each call to fcnSupercuts
678// - call TMinuit to do the minimization of (-significance) :
679// the fcnSupercuts function calculates the significance
680// for the current cut values
681// for this - the alpha plot is produced in the event loop,
682// in which the cuts have been applied
683// - the significance of the gamma signal in the alpha plot
684// is calculated
685//
686// Needed as input : (to be set by the Set functions)
687//
688// - fFilenameParam name of file to which optimum values of the
689// parameters are written
690//
691// - fHadronnessName name of container where MCT1SupercutsCalc will
692// put the supercuts hadronness
693//
694// - for the minimization, the starting values of the parameters are taken
695// - from file parSCinit (if it is != "")
696// - or from the arrays params and/or steps
697//
698//----------------------------------------------------------------------
699Bool_t MCT1FindSupercuts::FindParams(TString parSCinit,
700 TArrayD &params, TArrayD &steps)
701{
702 // Setup the event loop which will be executed in the
703 // fcnSupercuts function of MINUIT
704 //
705 // parSCinit is the name of the file containing the initial values
706 // of the parameters;
707 // if parSCinit = "" 'params' and 'steps' are taken as initial values
708
709 *fLog << "=============================================" << endl;
710 *fLog << "Setup event loop for fcnSupercuts" << endl;
711
712 if (fHadronnessName.IsNull())
713 {
714 *fLog << "MCT1FindSupercuts::FindParams; hadronness name is not defined... aborting"
715 << endl;
716 return kFALSE;
717 }
718
719 if (fMatrixTrain == NULL)
720 {
721 *fLog << "MCT1FindSupercuts::FindParams; training matrix is not defined... aborting"
722 << endl;
723 return kFALSE;
724 }
725
726 if (fMatrixTrain->GetM().GetNrows() <= 0)
727 {
728 *fLog << "MCT1FindSupercuts::FindParams; training matrix has no entries"
729 << endl;
730 return kFALSE;
731 }
732
733 MParList parlistfcn;
734 MTaskList tasklistfcn;
735
736 // loop over rows of matrix
737 MMatrixLoop loop(fMatrixTrain);
738
739 //--------------------------------
740 // create container for the supercut parameters
741 // and set them to their initial values
742 MCT1Supercuts super;
743
744 // take initial values from file parSCinit
745 if (parSCinit != "")
746 {
747 TFile inparam(parSCinit);
748 super.Read("MCT1Supercuts");
749 inparam.Close();
750 *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file "
751 << parSCinit << endl;
752 }
753
754 // take initial values from 'params' and/or 'steps'
755 else if (params.GetSize() != 0 || steps.GetSize() != 0 )
756 {
757 if (params.GetSize() != 0)
758 {
759 *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'"
760 << endl;
761 super.SetParameters(params);
762 }
763 if (steps.GetSize() != 0)
764 {
765 *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'"
766 << endl;
767 super.SetStepsizes(steps);
768 }
769 }
770 else
771 {
772 *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor"
773 << endl;
774 }
775
776
777 //--------------------------------
778 // calculate supercuts hadronness
779 fCalcHadTrain->SetHadronnessName(fHadronnessName);
780
781 // apply the supercuts
782 MF scfilter(fHadronnessName+".fHadronness>0.5");
783 MContinue supercuts(&scfilter);
784
785 // plot |alpha|
786 const TString mh3Name = "AlphaFcn";
787 MBinning binsalpha("Binning"+mh3Name);
788 binsalpha.SetEdges(54, -12.0, 96.0);
789
790 *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;
791
792 // |alpha| is assumed to be in column 7 of the matrix
793 MH3 alpha("MatrixTrain[7]");
794 alpha.SetName(mh3Name);
795
796 MFillH fillalpha(&alpha);
797
798 // book histograms to be filled during the optimization
799 // ! not in the event loop !
800 MHCT1Supercuts plotsuper;
801 plotsuper.SetupFill(&parlistfcn);
802
803 //******************************
804 // entries in MParList (extension of old MParList)
805
806 parlistfcn.AddToList(&tasklistfcn);
807 parlistfcn.AddToList(&super);
808 parlistfcn.AddToList(fCam);
809 parlistfcn.AddToList(fMatrixTrain);
810
811 parlistfcn.AddToList(&binsalpha);
812 parlistfcn.AddToList(&alpha);
813
814 parlistfcn.AddToList(&plotsuper);
815
816 //******************************
817 // entries in MTaskList
818
819 tasklistfcn.AddToList(&loop);
820 tasklistfcn.AddToList(fCalcHadTrain);
821 tasklistfcn.AddToList(&supercuts);
822 tasklistfcn.AddToList(&fillalpha);
823
824
825 //******************************
826
827 MEvtLoop evtloopfcn("EvtLoopFCN");
828 evtloopfcn.SetParList(&parlistfcn);
829 *fLog << "Event loop for fcnSupercuts has been setup" << endl;
830
831 // address of evtloopfcn is used in CallMinuit()
832
833
834 //-----------------------------------------------------------------------
835 //
836 //---------- Start of minimization part --------------------
837 //
838 // Do the minimization with MINUIT
839 //
840 // Be careful: This is not thread safe
841 //
842 *fLog << "========================================================" << endl;
843 *fLog << "Start minimization for supercuts" << endl;
844
845
846 // -------------------------------------------
847 // prepare call to MINUIT
848 //
849
850 // get initial values of parameters
851 fVinit = super.GetParameters();
852 fStep = super.GetStepsizes();
853
854 TString name[fVinit.GetSize()];
855 fStep.Set(fVinit.GetSize());
856 fLimlo.Set(fVinit.GetSize());
857 fLimup.Set(fVinit.GetSize());
858 fFix.Set(fVinit.GetSize());
859
860 fNpar = fVinit.GetSize();
861
862 for (UInt_t i=0; i<fNpar; i++)
863 {
864 name[i] = "p";
865 name[i] += i+1;
866 //fStep[i] = TMath::Abs(fVinit[i]/10.0);
867 fLimlo[i] = -100.0;
868 fLimup[i] = 100.0;
869 fFix[i] = 0;
870 }
871
872 // these parameters make no sense, fix them at 0.0
873 fVinit[33] = 0.0;
874 fStep[33] = 0.0;
875 fFix[33] = 1;
876
877 fVinit[36] = 0.0;
878 fStep[36] = 0.0;
879 fFix[36] = 1;
880
881 fVinit[39] = 0.0;
882 fStep[39] = 0.0;
883 fFix[39] = 1;
884
885 fVinit[41] = 0.0;
886 fStep[41] = 0.0;
887 fFix[41] = 1;
888
889 fVinit[44] = 0.0;
890 fStep[44] = 0.0;
891 fFix[44] = 1;
892
893 fVinit[47] = 0.0;
894 fStep[47] = 0.0;
895 fFix[47] = 1;
896
897 // vary only first 48 parameters
898 //for (UInt_t i=0; i<fNpar; i++)
899 //{
900 // if (i >= 48)
901 // {
902 // fStep[i] = 0.0;
903 // fFix[i] = 1;
904 // }
905 //}
906
907
908 // -------------------------------------------
909 // call MINUIT
910
911 TStopwatch clock;
912 clock.Start();
913
914 *fLog << "before calling CallMinuit" << endl;
915
916 MMinuitInterface inter;
917 Bool_t rc = inter.CallMinuit(fcnSupercuts, name,
918 fVinit, fStep, fLimlo, fLimup, fFix,
919 &evtloopfcn, "SIMPLEX", kFALSE);
920
921 *fLog << "after calling CallMinuit" << endl;
922
923 *fLog << "Time spent for the minimization in MINUIT : " << endl;;
924 clock.Stop();
925 clock.Print();
926
927 plotsuper.DrawClone();
928
929 if (!rc)
930 return kFALSE;
931
932 *fLog << "Minimization for supercuts finished" << endl;
933 *fLog << "========================================================" << endl;
934
935
936 // -----------------------------------------------------------------
937 // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put
938 // into MCT1Supercuts
939
940 // write optimum parameter values onto file filenameParam
941
942 TFile outparam(fFilenameParam, "RECREATE");
943 super.Write();
944 outparam.Close();
945
946 *fLog << "Optimum parameter values for supercuts were written onto file '"
947 << fFilenameParam << "' :" << endl;
948
949 const TArrayD &check = super.GetParameters();
950 for (Int_t i=0; i<check.GetSize(); i++)
951 *fLog << check[i] << ", ";
952 *fLog << endl;
953
954
955
956 *fLog << "End of supercuts optimization part" << endl;
957 *fLog << "======================================================" << endl;
958
959 return kTRUE;
960}
961
962
963// -----------------------------------------------------------------------
964//
965// Test the supercuts on the test sample
966//
967
968Bool_t MCT1FindSupercuts::TestParams()
969{
970 if (fMatrixTest->GetM().GetNrows() <= 0)
971 {
972 *fLog << "MCT1FindSupercuts::TestParams; test matrix has no entries"
973 << endl;
974 return kFALSE;
975 }
976
977 // ------------- TEST the supercuts ------------------------------
978 //
979 *fLog << "Test the supercuts on the test sample" << endl;
980
981 // -----------------------------------------------------------------
982 // read optimum parameter values from file filenameParam
983 // into array 'supercutsPar'
984
985 TFile inparam(fFilenameParam);
986 MCT1Supercuts scin;
987 scin.Read("MCT1Supercuts");
988 inparam.Close();
989
990 *fLog << "Optimum parameter values for supercuts were read from file '";
991 *fLog << fFilenameParam << "' :" << endl;
992
993 const TArrayD &supercutsPar = scin.GetParameters();
994 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
995 *fLog << supercutsPar[i] << ", ";
996 *fLog << endl;
997 //---------------------------
998
999
1000 // -----------------------------------------------------------------
1001 if (fHadronnessName.IsNull())
1002 {
1003 *fLog << "MCT1FindSupercuts::TestParams; hadronness name is not defined... aborting";
1004 *fLog << endl;
1005 return kFALSE;
1006 }
1007
1008 MParList parlist2;
1009 MTaskList tasklist2;
1010
1011 MCT1Supercuts supercuts;
1012 supercuts.SetParameters(supercutsPar);
1013
1014 fCalcHadTest->SetHadronnessName(fHadronnessName);
1015
1016
1017 //-------------------------------------------
1018
1019 MMatrixLoop loop(fMatrixTest);
1020
1021 // plot alpha before applying the supercuts
1022 const TString mh3NameB = "AlphaBefore";
1023 MBinning binsalphabef("Binning"+mh3NameB);
1024 binsalphabef.SetEdges(54, -12.0, 96.0);
1025
1026 // |alpha| is assumed to be in column 7 of the matrix
1027 MH3 alphabefore("MatrixTest[7]");
1028 alphabefore.SetName(mh3NameB);
1029
1030 TH1 &alphahistb = alphabefore.GetHist();
1031 alphahistb.SetName("alphaBefore-TestParams");
1032
1033 MFillH fillalphabefore(&alphabefore);
1034 fillalphabefore.SetName("FillAlphaBefore");
1035
1036 // apply the supercuts
1037 MF scfilter(fHadronnessName+".fHadronness>0.5");
1038 MContinue applysupercuts(&scfilter);
1039
1040 // plot alpha after applying the supercuts
1041 const TString mh3NameA = "AlphaAfter";
1042 MBinning binsalphaaft("Binning"+mh3NameA);
1043 binsalphaaft.SetEdges(54, -12.0, 96.0);
1044
1045 MH3 alphaafter("MatrixTest[7]");
1046 alphaafter.SetName(mh3NameA);
1047
1048 TH1 &alphahista = alphaafter.GetHist();
1049 alphahista.SetName("alphaAfter-TestParams");
1050
1051 MFillH fillalphaafter(&alphaafter);
1052 fillalphaafter.SetName("FillAlphaAfter");
1053
1054 //******************************
1055 // entries in MParList
1056
1057 parlist2.AddToList(&tasklist2);
1058
1059 parlist2.AddToList(&supercuts);
1060
1061 parlist2.AddToList(fCam);
1062 parlist2.AddToList(fMatrixTest);
1063
1064 parlist2.AddToList(&binsalphabef);
1065 parlist2.AddToList(&alphabefore);
1066
1067 parlist2.AddToList(&binsalphaaft);
1068 parlist2.AddToList(&alphaafter);
1069
1070 //******************************
1071 // entries in MTaskList
1072
1073 tasklist2.AddToList(&loop);
1074 tasklist2.AddToList(&fillalphabefore);
1075
1076 tasklist2.AddToList(fCalcHadTest);
1077 tasklist2.AddToList(&applysupercuts);
1078
1079 tasklist2.AddToList(&fillalphaafter);
1080
1081 //******************************
1082
1083 MProgressBar bar2;
1084 MEvtLoop evtloop2;
1085 evtloop2.SetParList(&parlist2);
1086 evtloop2.SetName("EvtLoopTestParams");
1087 evtloop2.SetProgressBar(&bar2);
1088 Int_t maxevents2 = -1;
1089 if (!evtloop2.Eventloop(maxevents2))
1090 return kFALSE;
1091
1092 tasklist2.PrintStatistics(0, kTRUE);
1093
1094
1095 //-------------------------------------------
1096 // draw the alpha plots
1097
1098 MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
1099 TH1 &alphaHist1 = alphaBefore->GetHist();
1100 alphaHist1.SetXTitle("|\\alpha| [\\circ]");
1101
1102 MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
1103 TH1 &alphaHist2 = alphaAfter->GetHist();
1104 alphaHist2.SetXTitle("|\\alpha| [\\circ]");
1105 alphaHist2.SetName("alpha-TestParams");
1106
1107 TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
1108 c->Divide(2,1);
1109
1110 gROOT->SetSelectedPad(NULL);
1111
1112 c->cd(1);
1113 alphaHist1.DrawCopy();
1114
1115 c->cd(2);
1116 alphaHist2.DrawCopy();
1117
1118
1119
1120 //-------------------------------------------
1121 // fit alpha distribution to get the number of excess events and
1122 // calculate significance of gamma signal in the alpha plot
1123
1124 const Double_t alphasig = 20.0;
1125 const Double_t alphamin = 30.0;
1126 const Double_t alphamax = 90.0;
1127 const Int_t degree = 2;
1128 const Bool_t drawpoly = kTRUE;
1129 const Bool_t fitgauss = kTRUE;
1130 const Bool_t print = kTRUE;
1131
1132 MHFindSignificance findsig;
1133 findsig.SetRebin(kTRUE);
1134 findsig.SetReduceDegree(kFALSE);
1135
1136 findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree,
1137 alphasig, drawpoly, fitgauss, print);
1138
1139 const Double_t significance = findsig.GetSignificance();
1140 const Double_t alphasi = findsig.GetAlphasi();
1141
1142 *fLog << "Significance of gamma signal after supercuts in test sample : "
1143 << significance << " (for |alpha| < " << alphasi << " degrees)"
1144 << endl;
1145 //-------------------------------------------
1146
1147
1148 *fLog << "Test of supercuts part finished" << endl;
1149 *fLog << "======================================================" << endl;
1150
1151 return kTRUE;
1152}
1153
Note: See TracBrowser for help on using the repository browser.