source: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.cc@ 6723

Last change on this file since 6723 was 6355, checked in by mazin, 20 years ago
*** empty log message ***
File size: 118.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Wolfgang Wittek 9/2003 <mailto:wittek@mppmu.mpg.de>
19! David Paneque 11/2003 <mailto:dpaneque@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MFindSupercutsONOFF //
29// //
30// Class for optimizing the parameters of the supercuts //
31// Using ON and OFF data //
32// //
33// //
34/////////////////////////////////////////////////////////////////////////////
35#include "MFindSupercutsONOFF.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 "MSupercuts.h"
49#include "MSupercutsCalcONOFF.h"
50#include "MDataElement.h"
51#include "MDataMember.h"
52
53
54#include <TPostScript.h>
55
56#include "MEvtLoop.h"
57#include "MFCT1SelFinal.h"
58#include "MF.h"
59#include "MFEventSelector.h"
60#include "MFEventSelector2.h"
61#include "MFillH.h"
62//#include "MGeomCamCT1Daniel.h"
63//#include "MGeomCamCT1.h"
64#include "MGeomCamMagic.h"
65#include "MFRandomSplit.h"
66#include "MH3.h"
67#include "MHCT1Supercuts.h"
68#include "MHFindSignificance.h" // To be removed at some point...
69#include "MHFindSignificanceONOFF.h"
70#include "MTSupercutsApplied.h"
71#include "MHMatrix.h"
72#include "MHOnSubtraction.h"
73#include "MDataValue.h"
74// #include "MDataString.h"
75
76#include "MLog.h"
77#include "MLogManip.h"
78#include "MMatrixLoop.h"
79#include "MMinuitInterface.h"
80#include "MParList.h"
81#include "MProgressBar.h"
82#include "MReadMarsFile.h"
83#include "MReadTree.h"
84#include "MTaskList.h"
85
86
87ClassImp(MFindSupercutsONOFF);
88
89using namespace std;
90
91
92
93// Function that computes the normalization factor using COUNTED events in alpha histogram
94// for ON data and alpha histogram for OFF data;
95// in alpha region defined by AlphaBkgMin and AlphaBkgMax
96
97// It is defined outside the class MFindSupercutsONOFF so that it can be used
98// in function fcnSupercuts
99
100
101static Double_t ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF,
102 Double_t AlphaBkgMin,
103 Double_t AlphaBkgMax)
104{
105
106 Double_t NormFactor = 0.0;
107 Double_t ONEvents = 0.0;
108 Double_t OFFEvents = 0.0;
109
110 const Double_t SmallQuantity = 0.01;
111
112 Double_t xlo = 0.0;
113 Double_t xup = 0.0;
114 Double_t width = 0.0;
115
116 Int_t BinCounterOFF = 0;
117 Int_t BinCounterON = 0;
118
119
120 // I make a copy of the histograms so that nothing happens to the
121 // histograms used in argument
122
123 TH1* HistON = (TH1*) histON->Clone();
124 TH1* HistOFF = (TH1*) histOFF->Clone();
125
126 if ( !HistON )
127 {
128 gLog << "ComputeNormFactorFromAlphaBkg; " << endl
129 << "Clone of ON histogram could not be generated"
130 << endl;
131 return 0.0;
132 }
133
134
135
136 if ( !HistOFF )
137 {
138 gLog << "ComputeNormFactorFromAlphaBkg; " << endl
139 << " Clone of OFF histogram could not be generated"
140 << endl;
141 return 0.0;
142 }
143
144 // Calculate the number of OFF events in the Bkg region
145 // defined by AlphaBkgMin and AlphaBkgMax
146 // ___________________________________________________
147
148
149 Int_t nbinsOFF = HistOFF -> GetNbinsX();
150 Double_t binwidthOFF = HistOFF -> GetBinWidth(1);
151
152 for (Int_t i=1; i<=nbinsOFF; i++)
153 {
154 xlo = HistOFF->GetBinLowEdge(i);
155 xup = HistOFF->GetBinLowEdge(i+1);
156
157 // bin must be completely contained in the bkg region
158 if ( xlo >= (AlphaBkgMin-SmallQuantity) && xup <= (AlphaBkgMax+SmallQuantity) )
159 {
160 width = fabs(xup-xlo);
161 if (fabs(width-binwidthOFF) > SmallQuantity)
162 {
163 gLog << "ComputeNormFactorFromAlphaBkg; "
164 << endl << " HistOFF has variable binning, which is not allowed"
165 << endl;
166 return 0.0;
167 }
168
169 BinCounterOFF++;
170
171 OFFEvents += HistOFF->GetBinContent(i);
172
173 }
174 }
175
176
177 // Calculate the number of ON events in the Bkg region
178 // defined by AlphaBkgMin and AlphaBkgMax
179 // ___________________________________________________
180
181
182 Int_t nbinsON = HistON -> GetNbinsX();
183 Double_t binwidthON = HistON -> GetBinWidth(1);
184
185 for (Int_t i=1; i<=nbinsON; i++)
186 {
187 xlo = HistON->GetBinLowEdge(i);
188 xup = HistON->GetBinLowEdge(i+1);
189
190 // bin must be completely contained in the bkg region
191 if ( xlo >= (AlphaBkgMin-SmallQuantity) && xup <= (AlphaBkgMax+SmallQuantity) )
192 {
193 width = fabs(xup-xlo);
194 if (fabs(width-binwidthON) > SmallQuantity)
195 {
196 gLog << "ComputeNormFactorFromAlphaBkg; "
197 << endl << " HistON has variable binning, which is not allowed"
198 << endl;
199 return 0.0;
200 }
201
202 BinCounterON++;
203 ONEvents += HistON->GetBinContent(i);
204
205 }
206 }
207
208
209
210
211 // NormFactor is computed
212
213 if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity)
214 {
215 gLog << "ComputeNormFactorFromAlphaBkg; "
216 << endl
217 << "ONEvents or OFFEvents computed in bkg region are < "
218 << SmallQuantity << endl;
219 return 0.0;
220
221 }
222
223
224
225
226 NormFactor = ONEvents/OFFEvents;
227
228 Double_t error = 1/ONEvents + 1/OFFEvents;
229 error = TMath::Sqrt(error);
230 error = error * NormFactor;
231
232 // tmp info
233 gLog << "ComputeNormFactorFromAlphaBkg;" << endl
234 << "ON Events in bkg region = " << ONEvents
235 << " (" << BinCounterON << " bins)"
236 << " , OFF Events in bkg region = " << OFFEvents
237 << " (" << BinCounterOFF << " bins)" << endl
238 <<"NormFactor computed from bkg region = " << NormFactor
239 << " +/- " << error << endl;
240 // end temp
241
242 return NormFactor;
243
244}
245
246
247
248
249
250
251
252
253//------------------------------------------------------------------------
254//
255// fcnSupercuts
256//
257// - calculates the quantity to be minimized (using TMinuit)
258//
259// - the quantity to be minimized is (-1)*significance of the gamma signal
260// in the alpha distribution (after cuts)
261//
262// - the parameters to be varied in the minimization are the cut parameters
263// (par)
264//
265static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f,
266 Double_t *par, Int_t iflag)
267{
268 //cout << "entry fcnSupercuts" << endl;
269
270 //-------------------------------------------------------------
271 // save pointer to the MINUIT object for optimizing the supercuts
272 // because it will be overwritten
273 // when fitting the alpha distribution in MHFindSignificance
274 TMinuit *savePointer = gMinuit;
275 //-------------------------------------------------------------
276
277
278 MEvtLoop* evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();
279
280 // Event loops for ON and OFF data are recovered.
281
282 MEvtLoop* ONDataevtloopfcn = &evtloopfcn[0];
283 MEvtLoop* OFFDataevtloopfcn = &evtloopfcn[1];
284
285
286 // Parameter list from event loops for ON and OFF data are recovered
287
288 MParList *ONDataplistfcn = (MParList*) ONDataevtloopfcn->GetParList();
289 MParList *OFFDataplistfcn = (MParList*) OFFDataevtloopfcn->GetParList();
290
291 // MTaskList *ONDataTasklistfcn = (MTaskList*) ONDataevtloopfcn->GetTaskList();
292 // MTaskList *OFFDataTasklistfcn = (MTaskList*) OFFDataevtloopfcn->GetTaskList();
293
294
295
296 // Container for supercuts is retrieved from ONDataplistfcn.
297 // NOTE: The same supercuts parameter container is used in OFFDataplistfcn
298
299 MSupercuts *super = (MSupercuts*) ONDataplistfcn->FindObject("MSupercuts");
300 if (!super)
301 {
302 gLog << "fcnSupercuts : MSupercuts object '" << "MSupercuts"
303 << "' not found... aborting" << endl;
304 return;
305 }
306
307
308 // Normalization factor for train sample (Train ON before cuts/ train OFF before cuts)
309 // is retrieved from the ONDataplistfcn
310
311 Double_t NormalizationFactorTrain = 0.0;
312
313 MDataValue* NormFactorContainer = (MDataValue*) ONDataplistfcn->FindObject("NormFactorTrain");
314 NormalizationFactorTrain = NormFactorContainer -> GetValue();
315
316
317 gLog << "fcnSupercuts : Normalization factor retrieved from ONDataplistfcn: " << endl
318 << "NormalizationFactorTrain = " << NormalizationFactorTrain << endl;
319
320
321 // Degree of polynomials used to fit the ON and the OFF data are retrieved from
322 // the event loop
323
324 Int_t degreeON;
325 Int_t degreeOFF;
326
327 MDataValue* DegreeONContainer = (MDataValue*) ONDataplistfcn->FindObject("DegreeON");
328 degreeON = int(DegreeONContainer -> GetValue());
329
330 MDataValue* DegreeOFFContainer = (MDataValue*) OFFDataplistfcn->FindObject("DegreeOFF");
331 degreeOFF = int(DegreeOFFContainer -> GetValue());
332
333 gLog << "fcnSupercuts : fDegreeON and fDegreeOFF retrieved from ONDataplistfcn and OFFDataplistfcn: " << endl
334 << "fDegreeON, fDegreeOFF: = " << degreeON << " ," << degreeOFF << endl;
335
336
337 Double_t AlphaSig = 0.0;
338 MDataValue* AlphaSigContainer = (MDataValue*) ONDataplistfcn->FindObject("AlphaSigValue");
339 AlphaSig = AlphaSigContainer -> GetValue();
340
341 Double_t AlphaBkgMin = 0.0;
342 MDataValue* AlphaBkgMinContainer =
343 (MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMinValue");
344 AlphaBkgMin = AlphaBkgMinContainer -> GetValue();
345
346 Double_t AlphaBkgMax = 0.0;
347 MDataValue* AlphaBkgMaxContainer =
348 (MDataValue*) ONDataplistfcn->FindObject("AlphaBkgMaxValue");
349 AlphaBkgMax = AlphaBkgMaxContainer -> GetValue();
350
351
352 gLog << "fcnSupercuts : AlphaSig and AlphaBkgMin-AlphaBkgMax retrieved from ONDataplistfcn: "
353 << endl
354 << "AlphaSig = " << AlphaSig << "; AlphaBkgMin-AlphaBkgMax = "
355 << AlphaBkgMin << "-" << AlphaBkgMax << endl;
356
357
358
359 // Variable fUseFittedQuantities is retrieved from ONDataplistfcn
360
361 Bool_t UseFittedQuantities;
362 MDataValue* UseFittedQuantitiesContainer =
363 (MDataValue*) ONDataplistfcn->FindObject("UseFittedQuantitiesValue");
364 UseFittedQuantities = bool(UseFittedQuantitiesContainer -> GetValue());
365
366 if (UseFittedQuantities)
367 {
368 gLog << "fcnSupercuts : UseFittedQuantities variable set to kTRUE" << endl;
369 }
370 else
371 {
372 gLog << "fcnSupercuts : UseFittedQuantities variable set to kFALSE" << endl;
373 }
374
375
376 Bool_t UseNormFactorFromAlphaBkg;
377 MDataValue* UseNormFactorFromAlphaBkgContainer =
378 (MDataValue*) ONDataplistfcn -> FindObject("UseNormFactorFromAlphaBkgValue");
379 UseNormFactorFromAlphaBkg = bool(UseNormFactorFromAlphaBkgContainer -> GetValue());
380
381
382 if (UseNormFactorFromAlphaBkg)
383 {
384 gLog << "fcnSupercuts : UseNormFactorFromAlphaBkg variable set to kTRUE" << endl;
385 }
386 else
387 {
388 gLog << "fcnSupercuts : UseNormFactorFromAlphaBkg variable set to kFALSE" << endl;
389 }
390
391
392
393
394 //
395 // transfer current parameter values to MSupercuts
396 //
397 // Attention : npar is the number of variable parameters
398 // not the total number of parameters
399 //
400 Double_t fMin, fEdm, fErrdef;
401 Int_t fNpari, fNparx, fIstat;
402 gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
403
404 super->SetParameters(TArrayD(fNparx, par));
405
406 //$$$$$$$$$$$$$$$$$$$$$
407 // for testing
408 // TArrayD checkparameters = super->GetParameters();
409 //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", "
410 // << fNparx << endl;
411 //gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
412 //for (Int_t i=0; i<fNparx; i++)
413 //{
414 //gLog << i << ", " << par[i] << ", " << checkparameters[i] << endl;
415 // }
416 //$$$$$$$$$$$$$$$$$$$$$
417
418 //
419 // plot alpha for ON data with the current cuts
420 //
421 if (!ONDataevtloopfcn->Eventloop())
422 {
423 gLog << "fcnsupercuts : ONDataevtloopfcn->Eventloop() failed" << endl;
424 }
425
426 // Somehow (??) I can not use the function MTaskList::PrintStatistics...
427 // it produces a segmentation fault...
428
429 //ONDataTasklistfcn->PrintStatistics(0, kFALSE);
430
431 MH3* alpha = (MH3*)ONDataplistfcn->FindObject("AlphaFcn", "MH3");
432 if (!alpha)
433 return;
434
435 TH1 &alphaHist = alpha->GetHist();
436 alphaHist.SetName("alpha-fcnSupercuts");
437
438
439 //
440 // plot alpha for OFF data with the current cuts
441 //
442 if(!OFFDataevtloopfcn->Eventloop())
443 {
444 gLog << "fcnsupercuts : OFFDataevtloopfcn->Eventloop() failed" << endl;
445 }
446
447 // Somehow (??) I can not use the function MTaskList::PrintStatistics...
448 // it produces a segmentation fault...
449
450 //OFFDataTasklistfcn->PrintStatistics(0, kFALSE);
451
452 MH3* alphaOFF = (MH3*) OFFDataplistfcn->FindObject("AlphaOFFFcn", "MH3");
453 if (!alphaOFF)
454 return;
455
456 TH1 &alphaHistOFF = alphaOFF->GetHist();
457 alphaHistOFF.SetName("alphaOFF-fcnSupercuts");
458
459
460
461
462
463 if (UseNormFactorFromAlphaBkg)
464 {
465 // Normalization factor computed using alpha bkg region
466 Double_t NewNormFactor =
467 ComputeNormFactorFromAlphaBkg(&alphaHist, &alphaHistOFF,
468 AlphaBkgMin, AlphaBkgMax);
469
470 gLog << "Normalization factor computed from alpha plot (after cuts) " << endl
471 << "using counted number of ON and OFF events in alpha region " << endl
472 << "defined by range " << AlphaBkgMin << "-" << AlphaBkgMax << endl
473 << "Normalization factor = " << NewNormFactor << endl;
474
475 gLog << "Normalization factor used is the one computed in bkg region; " << endl
476 << "i.e. " << NewNormFactor << " instead of " << NormalizationFactorTrain << endl;
477
478 NormalizationFactorTrain = NewNormFactor;
479
480
481 }
482
483
484
485
486
487
488
489 //-------------------------------------------
490 // set Minuit pointer to zero in order not to destroy the TMinuit
491 // object for optimizing the supercuts
492 gMinuit = NULL;
493
494 //=================================================================
495 // fit alpha distribution to get the number of excess events and
496 // calculate significance of gamma signal in the alpha plot
497
498 const Double_t alphasig = AlphaSig;
499 const Double_t alphamin = AlphaBkgMin;
500 // alpha min for bkg region in ON data
501 const Double_t alphamax = AlphaBkgMax; // alpha max for bkg region in ON data
502
503
504 Bool_t drawpoly;
505 Bool_t fitgauss;
506 Bool_t saveplots;
507
508
509 if (iflag == 3)
510 {// Even though minimization finished successfully, I will NOT produce
511 // final plots now... i'll do it later via the function
512 // "MFindSupercutsONOFF::OutputOptimizationOnTrainSample()"
513
514 drawpoly = kFALSE;
515 fitgauss = kFALSE;
516 saveplots = kFALSE;
517 }
518 else
519 {
520 drawpoly = kFALSE;
521 fitgauss = kFALSE;
522 saveplots = kFALSE;
523
524
525 }
526
527
528 const Bool_t print = kTRUE;
529
530 MHFindSignificanceONOFF findsig;
531 findsig.SetRebin(kTRUE);
532 findsig.SetReduceDegree(kFALSE);
533 findsig.SetUseFittedQuantities(UseFittedQuantities);
534
535 // TPostScript* DummyPs = new TPostScript("dummy.ps");
536
537 TString DummyPs = ("Dummy");
538
539
540 const Bool_t rc = findsig.FindSigmaONOFF(&alphaHist,&alphaHistOFF,
541 NormalizationFactorTrain,
542 alphamin, alphamax,
543 degreeON, degreeOFF,
544 alphasig, drawpoly, fitgauss,
545 print, saveplots,
546 DummyPs);
547 //DummyPs -> Close();
548 //delete DummyPs;
549
550
551 //=================================================================
552
553 // reset gMinuit to the MINUIT object for optimizing the supercuts
554 gMinuit = savePointer;
555 //-------------------------------------------
556
557 if (!rc)
558 {
559 gLog << "fcnSupercuts : FindSigmaONOFF() failed" << endl;
560 f = 1.e10;
561 return;
562 }
563
564 /*
565
566 // plot some quantities during the optimization
567 MHCT1Supercuts *plotsuper = (MHCT1Supercuts*)ONDataplistfcn->FindObject("MHCT1Supercuts");
568 if (plotsuper)
569 plotsuper->Fill(&findsig);
570
571
572
573 */
574
575 //------------------------
576 // get significance
577 const Double_t significance = findsig.GetSignificance();
578 f = significance>0 ? -significance : 0;
579
580
581 gLog << " ///****************************************** ///" << endl;
582 gLog << "Significance (Li&Ma)is now: " << f << endl;
583 gLog << " ///****************************************** ///" << endl;
584
585 //------------------------
586 // optimize signal/background ratio
587 //Double_t ratio = findsig.GetNbg()>0.0 ?
588 // findsig.GetNex()/findsig.GetNbg() : 0.0;
589 //f = -ratio;
590
591 //-------------------------------------------
592 // final calculations
593 //if (iflag == 3)
594 //{
595 //
596 //}
597
598 //-------------------------------------------------------------
599}
600
601
602// --------------------------------------------------------------------------
603//
604// Default constructor.
605//
606MFindSupercutsONOFF::MFindSupercutsONOFF(const char *name, const char *title)
607 : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL)
608{
609 fName = name ? name : "MFindSupercutsONOFF";
610 fTitle = title ? title : "Optimizer of the supercuts";
611
612 //---------------------------
613 // camera geometry is needed for conversion mm ==> degree
614 //fCam = new MGeomCamCT1Daniel;
615 // fCam = new MGeomCamCT1;
616 fCam = new MGeomCamMagic;
617
618 // matrices to contain the training/test samples
619 fMatrixTrain = new MHMatrix("MatrixTrain");
620 fMatrixTest = new MHMatrix("MatrixTest");
621 fMatrixTrainOFF = new MHMatrix("MatrixTrainOFF");
622 fMatrixTestOFF = new MHMatrix("MatrixTestOFF");
623
624
625 // objects of MSupercutsCalcONOFF to which these matrices are attached
626 fCalcHadTrain = new MSupercutsCalcONOFF("SupercutsCalcTrain");
627 fCalcHadTest = new MSupercutsCalcONOFF("SupercutsCalcTest");
628 fCalcHadTrainOFF = new MSupercutsCalcONOFF("SupercutsCalcTrainOFF");
629 fCalcHadTestOFF = new MSupercutsCalcONOFF("SupercutsCalcTestOFF");
630
631 // Define columns of matrices
632 fCalcHadTrain->InitMapping(fMatrixTrain);
633 fCalcHadTest->InitMapping(fMatrixTest);
634 fCalcHadTrainOFF->InitMapping(fMatrixTrainOFF);
635 fCalcHadTestOFF->InitMapping(fMatrixTestOFF);
636
637 // For the time being weights method is not implemented
638 //UseWeights = kFALSE;
639
640 // Normalization factors are initialized to -1
641 // Functions determining Number of excess events will not work
642 // with negative norm factors;
643 // This will ensure that Norm factors are computed before running
644 // DetExcessONOFF() function.
645
646 fNormFactorTrain = -1.0;
647 fNormFactorTest = -1.0;
648
649 // SigmaLiMa and Nex member variables are initialized to 0
650
651 fSigmaLiMaTrain = 0.0;
652 fSigmaLiMaTest = 0.0;
653 fNexTrain = 0.0;
654 fNexTest = 0.0;
655
656
657 // Cuts (low and up) in variable MPointingPos.fZd
658 // (at some point real theta)
659 // The default is not cut, i.e. all values (0-1) are taken
660
661 //fThetaMin = 0; // in miliradians // FIXME: change name
662 //fThetaMax = 1570; // in miliradians // FIXME: change name
663
664 fThetaMin = 0; // degrees (new in magic)
665 fThetaMax = 90; // degrees (new in magic)
666
667
668
669 fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1
670
671 // By default, bkg region is set to alpha range 30-90 for CT1
672 fAlphaBkgMin = 30;
673 fAlphaBkgMax = 90;
674
675
676 // By default, bining for alpha plots is set to
677
678 fNAlphaBins = 54;
679 fAlphaBinLow = -12;
680 fAlphaBinUp = 96;
681
682
683
684 fNormFactorFromAlphaBkg = kTRUE; // By default normalization factor is
685 // computed using bkg region in alpha histograms (after cuts)
686
687// fActualCosThetaBinCenter = 0;
688
689 fTuneNormFactor = kTRUE; // Norm factors will be corrected using
690 //the total amount of OFF events before cuts and the estimated excess events
691 // fNormFactorTrain = fNormFactorTrain - Ngammas/EventsInTrainMatrixOFF
692
693
694 // use quantities computed from the fits
695 // The variable allows the user to NOT use these quantities when there is not
696 // enough statistics and fit not always is possible.
697 // Default value is kTRUE
698 fUseFittedQuantities = kTRUE;
699
700
701
702
703 // Boolean variable that allows the user to set some limits to the
704 // some of the Minuit parameters. For the time being, only the limits
705 // for the parameters which do NOT depend in size, dist and theta are set,
706 // i.e. static limits. The value of this boolean variable is set in the
707 // constructor of the class.
708
709 fSetLimitsToSomeMinuitParams = kTRUE;
710
711 // Limits for the Minuit parameters. For the time being the values are set in the constructor
712 // of the class. One MUST be very careful to set limits such that the expected final values
713 // (optimized values) are far away from the limits.
714
715 // Values in degrees
716
717
718 fMinuitDistUPUpperLimit = 2.0;
719 fMinuitDistUPLowerLimit = 0.5;
720 fMinuitLengthUPUpperLimit = 0.8;
721 fMinuitLengthUPLowerLimit = 0.0;
722 fMinuitWidthUPUpperLimit = 0.5;
723 fMinuitWidthUPLowerLimit = 0.0;
724 fMinuitLeakage1UPUpperLimit = 1.5;
725 fMinuitLeakage1UPLowerLimit = -0.5;
726
727 fMinuitDistLOWUpperLimit = 1.0;
728 fMinuitDistLOWLowerLimit = 0.0;
729 fMinuitLengthLOWUpperLimit = 0.5;
730 fMinuitLengthLOWLowerLimit = -0.3;
731 fMinuitWidthLOWUpperLimit = 0.4;
732 fMinuitWidthLOWLowerLimit = -0.3;
733
734
735
736 // Boolean variable that controls wether the optimization of the
737 // parameters (MMinuitInterface::CallMinuit(..) in function FindParams(..))
738 // takes place or not. kTRUE will skip such optimization.
739 // This variable is useful to test the optmized parameters (previously found
740 // and stored in root file) on the TRAIN sample.
741
742 fSkipOptimization = kFALSE;
743
744 // Boolean variable that allows the user to write the initial parameters
745 // into the root file that will be used to store the optimum cuts.
746 // If fUseInitialSCParams = kTRUE , parameters are written.
747 // In this way, the initial SC parameters can be applied on the data (train/test)
748
749 // The initial parameters are ONLY written to the root file if
750 // there is NO SC params optimization, i.e., if variable
751 // fSkipOptimization = kTRUE;
752
753 // The default value is obviously kFALSE.
754
755 fUseInitialSCParams = kTRUE;
756
757
758 // Set wether to use or not hillas dist
759 fUseDist = kTRUE;
760
761
762
763
764 fGammaEfficiency = 0.5; // Fraction of gammas that remain after cuts
765 // Quantity that will have to be determined with MC, yet for the
766 // time being I set it to 0.5 (standard value)
767
768 fPsFilename = NULL;
769 fPsFilename2 = NULL;
770
771
772 ////////////////////////////////////////////////////
773 // TMP
774
775 // There are quite some problems during the data preprocessing.
776 // For the time being, I will add some cuts to the functions
777 // DefineTrainTestMatrixThetaRange and for OFF, so that I can
778 // make a kind of preprocess on my own. This allows me
779 // to make a very silly preprocess with wolfgangs macro, which
780 // might be free of corrupted data, and then I can do on my own.
781
782 fSizeCutLow = 0.1; // To prevent empty events
783 fSizeCutUp = 10000000;
784
785 // Angular cuts are converted to mm, which
786 // are the units of the preprocessed data....
787
788
789 // Angular cuts not yet implemented ...
790 fConvMMToDeg = 0.00337034;
791
792 fDistCutLow = 0.4/fConvMMToDeg;
793 fDistCutUp = 1.5/fConvMMToDeg;
794
795 fLengthCutLow = 0.1/fConvMMToDeg;
796 fLengthCutUp = 1/fConvMMToDeg;
797
798 fWidthCutLow = 0.07/fConvMMToDeg;
799 fWidthCutUp = 1/fConvMMToDeg;
800
801 // ENDTMP
802 ////////////////////////////////////////////////////
803
804
805 // Degree of polynomials used to fit the ON OFF data is
806 // set initially to 2.
807
808 fDegreeON = 2;
809 fDegreeOFF = 2;
810
811
812}
813
814// --------------------------------------------------------------------------
815//
816// Default destructor.
817//
818MFindSupercutsONOFF::~MFindSupercutsONOFF()
819{
820
821 *fLog << "destructor of MFindSupercutsONOFF is called" << endl;
822
823 fPsFilename = NULL;
824 fPsFilename2 = NULL;
825
826
827 delete fCam;
828 delete fMatrixTrain;
829 delete fMatrixTest;
830 delete fCalcHadTrain;
831 delete fCalcHadTest;
832 delete fMatrixTrainOFF;
833 delete fMatrixTestOFF;
834 delete fCalcHadTrainOFF;
835 delete fCalcHadTestOFF;
836
837 *fLog << "destructor of MFindSupercutsONOFF finished successfully" << endl;
838
839}
840
841// --------------------------------------------------------------------------
842//
843// Function that sets the name of the PostScript file where alpha distributions
844// for the different Theta bins will be stored.
845// It also initializes
846
847
848
849void MFindSupercutsONOFF::SetPostScriptFile (TPostScript* PsFile)
850{
851 fPsFilename = PsFile;
852
853 *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file "
854 << fPsFilename -> GetName() << endl;
855
856}
857
858void MFindSupercutsONOFF::SetPostScriptFile2 (TPostScript &PsFile)
859{
860 fPsFilename2 = new TPostScript (PsFile);
861
862 *fLog << "MFindSupercutsONOFF : Results (alpha distributions with excess and significances) will be stored in PostScript file "
863 << fPsFilename2 -> GetName() << endl;
864
865}
866
867void MFindSupercutsONOFF::SetPsFilenameString (const TString filename)
868{
869 fPsFilenameString = filename;
870}
871
872void MFindSupercutsONOFF::SetSkipOptimization(Bool_t b)
873{
874 fSkipOptimization = b;
875 if (fSkipOptimization)
876 {
877 *fLog << "MFindSupercutsONOFF :: SetSkipOptimization " << endl
878 << "Variable fSkipOptimization is kTRUE, and therefore "
879 << "the optimization of supercuts is skipped. Hope that's "
880 << "what you want... " << endl;
881 }
882
883}
884
885
886void MFindSupercutsONOFF::SetUseInitialSCParams(Bool_t b)
887{
888 fUseInitialSCParams = b;
889 if (fUseInitialSCParams)
890 {
891 *fLog << "MFindSupercutsONOFF :: SetUseInitialSCParams " << endl
892 << "Variable fUseInitialSCParams is kTRUE. " << endl;
893
894 if (fSkipOptimization)
895 {
896 *fLog << "The Initial SC Parameters will be applied on the selected data."
897 << endl;
898 }
899
900 else
901 {
902 *fLog << "However, fSkipOptimization = kFALSE, and therefore, the "
903 << "the supercuts will be optimized. The final cuts that "
904 << "will be applied to the data will NOT be the initial SC parameters."
905 << endl;
906
907 }
908
909
910
911 }
912
913}
914
915
916
917Bool_t MFindSupercutsONOFF::SetAlphaSig(Double_t alphasig)
918{
919 // check that alpha is within the limits 0-90
920 if (alphasig <= 0 || alphasig > 90)
921 {
922 *fLog << "MFindSupercutsONOFF ::SetAlphaSig; "
923 << "value " << alphasig << " is not within the the "
924 << "logical limits of alpha; 0-90" << endl;
925 return kFALSE;
926 }
927
928
929 fAlphaSig = alphasig;
930
931 return kTRUE;
932}
933
934Bool_t MFindSupercutsONOFF::SetAlphaBkgMin(Double_t alpha)
935{
936 // check that alpha is within the limits 0-90
937 if (alpha <= 0 || alpha >= 90)
938 {
939 *fLog << "MFindSupercutsONOFF ::SetAlphaBkgMin; "
940 << "value " << alpha << " is not within the the "
941 << "logical limits of alpha; 0-90" << endl;
942 return kFALSE;
943 }
944
945
946 fAlphaBkgMin = alpha;
947
948 return kTRUE;
949}
950
951
952Bool_t MFindSupercutsONOFF::SetAlphaBkgMax(Double_t alpha)
953{
954 // check that alpha is within the limits 0-90
955 if (alpha <= 0 || alpha > 90.001)
956 {
957 *fLog << "MFindSupercutsONOFF ::SetAlphaBkgMax; "
958 << "value " << alpha << " is not within the the "
959 << "logical limits of alpha; 0-90" << endl;
960 return kFALSE;
961 }
962
963
964 fAlphaBkgMax = alpha;
965
966 return kTRUE;
967}
968
969
970// Function that checks that the values of the member data
971// fAlphaSig, fAlphaBkgMin and fAlphaBkgMax make sense
972// (ie, fAlphaSig < fAlphaBkgMin < fAlphaBkgMax)
973
974Bool_t MFindSupercutsONOFF::CheckAlphaSigBkg()
975{
976
977 if (fAlphaSig > fAlphaBkgMin)
978 {
979 *fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl
980 << "fAlphaSig > fAlphaBkgMin, which should not occur..." << endl
981 << "fAlphaSig = " << fAlphaSig << ", fAlphaBkgMin = " << fAlphaBkgMin
982 << endl;
983
984 return kFALSE;
985 }
986
987 if (fAlphaBkgMax < fAlphaBkgMin)
988 {
989 *fLog << "MFindSupercutsONOFF ::CheckAlphaSigBkg(); " << endl
990 << "fAlphaBkgMin > fAlphaBkgMax, which should not occur..." << endl
991 << "fAlphaBkgMin = " << fAlphaBkgMin << ", fAlphaBkgMax = " << fAlphaBkgMax
992 << endl;
993
994 return kFALSE;
995 }
996
997 return kTRUE;
998
999}
1000
1001
1002/*
1003// Function that computes the normalization factor using COUNTED events ON and OFF
1004// in alpha region defined by fAlphaBkgMin and fAlphaBkgMax
1005
1006
1007Double_t MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg(TH1 *histON, TH1 *histOFF)
1008{
1009
1010 Double_t NormFactor = 0.0;
1011 Double_t ONEvents = 0.0;
1012 Double_t OFFEvents = 0.0;
1013
1014 const Double_t SmallQuantity = 0.01;
1015
1016 Double_t xlo = 0.0;
1017 Double_t xup = 0.0;
1018 Double_t width = 0.0;
1019
1020 Int_t BinCounterOFF = 0;
1021 Int_t BinCounterON = 0;
1022
1023
1024 // I make a copy of the histograms so that nothing happens to the
1025 // histograms used in argument
1026
1027 TH1* HistON = (TH1*) histON->Clone();
1028 TH1* HistOFF = (TH1*) histOFF->Clone();
1029
1030 if ( !HistON )
1031 {
1032 *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl
1033 << "Clone of ON histogram could not be generated"
1034 << endl;
1035 return 0.0;
1036 }
1037
1038
1039
1040 if ( !HistOFF )
1041 {
1042 *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; " << endl
1043 << " Clone of OFF histogram could not be generated"
1044 << endl;
1045 return 0.0;
1046 }
1047
1048 // Calculate the number of OFF events in the Bkg region
1049 // defined by fAlphaBkgMin and fAlphaBkgMax
1050 // ___________________________________________________
1051
1052
1053 Int_t nbinsOFF = HistOFF -> GetNbinsX();
1054 Double_t binwidthOFF = HistOFF -> GetBinWidth(1);
1055
1056 for (Int_t i=1; i<=nbinsOFF; i++)
1057 {
1058 xlo = HistOFF->GetBinLowEdge(i);
1059 xup = HistOFF->GetBinLowEdge(i+1);
1060
1061 // bin must be completely contained in the bkg region
1062 if ( xlo >= (fAlphaBkgMin-SmallQuantity) && xup <= (fAlphaBkgMax+SmallQuantity) )
1063 {
1064 width = fabs(xup-xlo);
1065 if (fabs(width-binwidthOFF) > SmallQuantity)
1066 {
1067 *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; "
1068 << endl << " HistOFF has variable binning, which is not allowed"
1069 << endl;
1070 return 0.0;
1071 }
1072
1073 BinCounterOFF++;
1074
1075 OFFEvents += HistOFF->GetBinContent(i);
1076
1077 }
1078 }
1079
1080
1081 // Calculate the number of ON events in the Bkg region
1082 // defined by fAlphaBkgMin and fAlphaBkgMax
1083 // ___________________________________________________
1084
1085
1086 Int_t nbinsON = HistON -> GetNbinsX();
1087 Double_t binwidthON = HistON -> GetBinWidth(1);
1088
1089 for (Int_t i=1; i<=nbinsON; i++)
1090 {
1091 xlo = HistON->GetBinLowEdge(i);
1092 xup = HistON->GetBinLowEdge(i+1);
1093
1094 // bin must be completely contained in the bkg region
1095 if ( xlo >= (fAlphaBkgMin-SmallQuantity) && xup <= (fAlphaBkgMax+SmallQuantity) )
1096 {
1097 width = fabs(xup-xlo);
1098 if (fabs(width-binwidthON) > SmallQuantity)
1099 {
1100 *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; "
1101 << endl << " HistON has variable binning, which is not allowed"
1102 << endl;
1103 return 0.0;
1104 }
1105
1106 BinCounterON++;
1107 ONEvents += HistON->GetBinContent(i);
1108
1109 }
1110 }
1111
1112
1113
1114
1115 // NormFactor is computed
1116
1117 if (ONEvents < SmallQuantity || OFFEvents < SmallQuantity)
1118 {
1119 *fLog << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg; "
1120 << endl
1121 << "ONEvents or OFFEvents computed in bkg region are < "
1122 << SmallQuantity << endl;
1123 return 0.0;
1124
1125 }
1126
1127
1128
1129
1130 NormFactor = ONEvents/OFFEvents;
1131
1132 Double_t error = 1/ONEvents + 1/OFFEvents;
1133 error = TMath::Sqrt(error);
1134 error = error * NormFactor;
1135
1136 // tmp info
1137 cout << "MFindSupercutsONOFF::ComputeNormFactorFromAlphaBkg;" << endl
1138 << "ON Events in bkg region = " << ONEvents
1139 << " (" << BinCounterON << " bins)"
1140 << " , OFF Events in bkg region = " << OFFEvents
1141 << " (" << BinCounterOFF << " bins)" << endl
1142 <<"NormFactor computed from bkg region = " << NormFactor
1143 << " +/- " << error << endl;
1144 // end temp
1145
1146 return NormFactor;
1147
1148}
1149
1150
1151*/
1152
1153
1154// Function that set the values of fThetaMin and fThetaMax and also
1155// fThetaRangeString (with value in miliradians);
1156// data members that are used basically to print/plot
1157// information.
1158
1159Bool_t MFindSupercutsONOFF::SetThetaRange(Double_t ThetaMin, Double_t ThetaMax)
1160{
1161
1162// Check that values are reasonable... well ... i guess this was done
1163// in previous functions...
1164
1165 // fThetaMin = int(ThetaMin*1000.0);
1166 // fThetaMax = int(ThetaMax*1000.0);
1167
1168
1169 // in new magic theta is given in deg
1170 // 0.5 is added so that the rounding to integer is correct
1171 fThetaMin = int(ThetaMin + 0.5);
1172 fThetaMax = int(ThetaMax + 0.5);
1173
1174 fThetaRangeString = ("ThetaRange");
1175 fThetaRangeString += (fThetaMin);
1176 fThetaRangeString += ("_");
1177 fThetaRangeString += (fThetaMax);
1178 // fThetaRangeString += ("mRad");
1179 fThetaRangeString += ("_Degrees"); // new in magic
1180
1181
1182 return kTRUE;
1183
1184}
1185
1186// Function that sets Size range
1187Bool_t MFindSupercutsONOFF::SetSizeRange(Double_t SizeMin, Double_t SizeMax)
1188{
1189
1190
1191 fSizeCutLow = SizeMin;
1192 fSizeCutUp = SizeMax;
1193
1194 *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl
1195 << "Data matrices will be filled with events whose MHillas.fSize " << endl
1196 << "is in the range "
1197 << fSizeCutLow <<"-"<<fSizeCutUp << endl;
1198
1199
1200
1201 return kTRUE;
1202}
1203
1204
1205Bool_t MFindSupercutsONOFF::SetFilters(Double_t LeakageMax, Double_t DistMax, Double_t DistMin)
1206{
1207
1208 fDistCutLow = DistMin/fConvMMToDeg;
1209 fDistCutUp = DistMax/fConvMMToDeg;
1210 fLeakageMax = LeakageMax;
1211
1212 *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl
1213 << "Data matrices will be filled with events whose MHillasSrc.fDist " << endl
1214 << "is in the range "
1215 << fDistCutLow <<"-"<< fDistCutUp << " degrees" << endl
1216 << "and fLeakage " << "< " << fLeakageMax << endl;
1217
1218 return kTRUE;
1219}
1220
1221
1222
1223// Function that sets the names of all parameter containers
1224// used to store the supercuts applied to ON/OFF Train/Test samples
1225
1226void MFindSupercutsONOFF::SetSupercutsAppliedTreeNames()
1227{
1228
1229 char* sc = {"SupercutsApplied"};
1230
1231 fTrainONSupercutsAppliedName = (sc);
1232 fTrainONSupercutsAppliedName += ("TrainON");
1233
1234 fTrainOFFSupercutsAppliedName = (sc);
1235 fTrainOFFSupercutsAppliedName += ("TrainOFF");
1236
1237 fTestONSupercutsAppliedName = (sc);
1238 fTestONSupercutsAppliedName += ("TestON");
1239
1240 fTestOFFSupercutsAppliedName = (sc);
1241 fTestOFFSupercutsAppliedName += ("TestOFF");
1242
1243
1244 if (fThetaRangeString.IsNull())
1245 {
1246 *fLog << "MFindSupercutsONOFF::SetSupercutsAppliedTreeNames; "
1247 << " fThetaRangeString is NOT defined..." << endl;
1248
1249 fTrainONSupercutsAppliedName += ("ThetaRangeStringUndefined");
1250 fTrainOFFSupercutsAppliedName += ("ThetaRangeStringUndefined");
1251 fTestONSupercutsAppliedName += ("ThetaRangeStringUndefined");
1252 fTestOFFSupercutsAppliedName += ("ThetaRangeStringUndefined");
1253 }
1254 else
1255 {
1256 fTrainONSupercutsAppliedName += (fThetaRangeString);
1257 fTrainOFFSupercutsAppliedName += (fThetaRangeString);
1258 fTestONSupercutsAppliedName += (fThetaRangeString);
1259 fTestOFFSupercutsAppliedName += (fThetaRangeString);
1260 }
1261
1262 // Info about names
1263
1264 *fLog << "MFindSupercutsONOFF::SetSupercutsAppliedTreeNames; " << endl
1265 << "Names of the MTSupercutsApplied Trees for Train (ON/OFF) "
1266 << " and Test (ON/OFF) samples are the following ones: " << endl
1267 << fTrainONSupercutsAppliedName << ", "
1268 << fTrainOFFSupercutsAppliedName << ", "
1269 << fTestONSupercutsAppliedName << ", "
1270 << fTestOFFSupercutsAppliedName << endl;
1271
1272}
1273
1274
1275void MFindSupercutsONOFF::SetUseOrigDistribution(Bool_t b)
1276{
1277 fUseOrigDistribution = b;
1278
1279 if (fUseOrigDistribution == kTRUE)
1280 {
1281
1282 *fLog << "MFindSupercutsONOFF : when defining training and test matrices use the original distribution"
1283 << endl;
1284 }
1285
1286 else
1287 {
1288 *fLog << "MFindSupercutsONOFF : when defining training and test matrices, events will be selected according to distribution given by a MH3 object"
1289 << endl;
1290 }
1291}
1292
1293
1294
1295// --------------------------------------------------------------------------
1296//
1297// Define the matrix 'fMatrixTrain' for the training sample
1298//
1299// alltogether 'howmanytrain' events are read from file 'nametrain';
1300// the events are selected according to a target distribution 'hreftrain'
1301//
1302//
1303Bool_t MFindSupercutsONOFF::DefineTrainMatrix(const TString &nametrain, MH3 &hreftrain,
1304 const Int_t howmanytrain,
1305 const TString &filetrain)
1306{
1307 if (nametrain.IsNull() || howmanytrain <= 0)
1308 return kFALSE;
1309
1310 *fLog << "=============================================" << endl;
1311 *fLog << "fill training matrix from file '" << nametrain
1312 << "', select " << howmanytrain
1313 << " events " << endl;
1314
1315 if (!fUseOrigDistribution)
1316 {
1317 *fLog << " according to a distribution given by the MH3 object '"
1318 << hreftrain.GetName() << "'" << endl;
1319 }
1320 else
1321 {
1322 *fLog << " randomly" << endl;
1323 }
1324
1325
1326
1327 MParList plist;
1328 MTaskList tlist;
1329
1330 MReadMarsFile read("Events", nametrain);
1331 read.DisableAutoScheme();
1332
1333
1334 MFEventSelector2 seltrain(hreftrain);
1335 seltrain.SetNumMax(howmanytrain);
1336 seltrain.SetName("selectTrain");
1337 if (fUseOrigDistribution)
1338 {
1339 *fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
1340 << " DOES NOT EXIST NOW..." << endl;
1341 // seltrain.SetUseOrigDistribution(kTRUE);
1342 }
1343
1344
1345 MFillH filltrain(fMatrixTrain);
1346 filltrain.SetFilter(&seltrain);
1347 filltrain.SetName("fillMatrixTrain");
1348
1349 //******************************
1350 // entries in MParList
1351
1352 plist.AddToList(&tlist);
1353 plist.AddToList(fCam);
1354 plist.AddToList(fMatrixTrain);
1355
1356 //******************************
1357 // entries in MTaskList
1358
1359 tlist.AddToList(&read);
1360 tlist.AddToList(&seltrain);
1361 tlist.AddToList(&filltrain);
1362
1363 //******************************
1364
1365 MProgressBar bar;
1366 MEvtLoop evtloop;
1367 evtloop.SetParList(&plist);
1368 evtloop.SetName("EvtLoopMatrixTrain");
1369 evtloop.SetProgressBar(&bar);
1370
1371 if (!evtloop.Eventloop())
1372 return kFALSE;
1373
1374 tlist.PrintStatistics(0, kTRUE);
1375
1376 fMatrixTrain->Print("SizeCols");
1377 Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
1378 if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
1379 {
1380 *fLog << "MFindSupercutsONOFF::DefineTrainMatrix; no.of generated events ("
1381 << howmanygenerated
1382 << ") is incompatible with the no.of requested events ("
1383 << howmanytrain << ")" << endl;
1384 }
1385
1386 *fLog << "training matrix was filled" << endl;
1387 *fLog << "=============================================" << endl;
1388
1389 //--------------------------
1390 // write out training matrix
1391
1392 if (filetrain != "")
1393 {
1394 TFile filetr(filetrain, "RECREATE");
1395 fMatrixTrain->Write();
1396 filetr.Close();
1397
1398 *fLog << "MFindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
1399 << filetrain << "'" << endl;
1400 }
1401
1402
1403 return kTRUE;
1404}
1405
1406
1407// --------------------------------------------------------------------------
1408
1409
1410
1411// --------------------------------------------------------------------------
1412//
1413// Define the matrix for the test sample
1414//
1415// alltogether 'howmanytest' events are read from file 'nametest'
1416//
1417// the events are selected according to a target distribution 'hreftest'
1418//
1419//
1420Bool_t MFindSupercutsONOFF::DefineTestMatrix(const TString &nametest, MH3 &hreftest,
1421 const Int_t howmanytest, const TString &filetest)
1422{
1423 if (nametest.IsNull() || howmanytest<=0)
1424 return kFALSE;
1425
1426 *fLog << "=============================================" << endl;
1427 *fLog << "fill test matrix from file '" << nametest
1428 << "', select " << howmanytest << " events " << endl;
1429
1430
1431 if (!fUseOrigDistribution)
1432 {
1433 *fLog << " according to a distribution given by the MH3 object '"
1434 << hreftest.GetName() << "'" << endl;
1435 }
1436 else
1437 {
1438 *fLog << " randomly" << endl;
1439 }
1440
1441
1442 MParList plist;
1443 MTaskList tlist;
1444
1445 MReadMarsFile read("Events", nametest);
1446 read.DisableAutoScheme();
1447
1448 MFEventSelector2 seltest(hreftest);
1449 seltest.SetNumMax(howmanytest);
1450 seltest.SetName("selectTest");
1451 if (fUseOrigDistribution)
1452 {
1453 *fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
1454 << " DOES NOT EXIST NOW..." << endl;
1455 // seltest.SetUseOrigDistribution(kTRUE);
1456 }
1457
1458
1459 MFillH filltest(fMatrixTest);
1460 filltest.SetFilter(&seltest);
1461 filltest.SetName("fillMatrixTest");
1462
1463
1464
1465 //******************************
1466 // entries in MParList
1467
1468 plist.AddToList(&tlist);
1469 plist.AddToList(fCam);
1470 plist.AddToList(fMatrixTest);
1471
1472 //******************************
1473 // entries in MTaskList
1474
1475 tlist.AddToList(&read);
1476 tlist.AddToList(&seltest);
1477 tlist.AddToList(&filltest);
1478
1479 //******************************
1480
1481 MProgressBar bar;
1482 MEvtLoop evtloop;
1483 evtloop.SetParList(&plist);
1484 evtloop.SetName("EvtLoopMatrixTest");
1485 evtloop.SetProgressBar(&bar);
1486
1487 if (!evtloop.Eventloop())
1488 return kFALSE;
1489
1490 tlist.PrintStatistics(0, kTRUE);
1491
1492 fMatrixTest->Print("SizeCols");
1493 const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
1494 if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
1495 {
1496 *fLog << "MFindSupercutsONOFF::DefineTestMatrix; no.of generated events ("
1497 << howmanygenerated
1498 << ") is incompatible with the no.of requested events ("
1499 << howmanytest << ")" << endl;
1500 }
1501
1502 *fLog << "test matrix was filled" << endl;
1503 *fLog << "=============================================" << endl;
1504
1505 //--------------------------
1506 // write out test matrix
1507
1508 if (filetest != "")
1509 {
1510 TFile filete(filetest, "RECREATE", "");
1511 fMatrixTest->Write();
1512 filete.Close();
1513
1514 *fLog << "MFindSupercutsONOFF::DefineTestMatrix; Test matrix was written onto file '"
1515 << filetest << "'" << endl;
1516 }
1517
1518
1519
1520 return kTRUE;
1521}
1522
1523
1524
1525// --------------------------------------------------------------------------
1526//
1527// Define the matrices for the training and test sample respectively
1528//
1529//
1530//
1531Bool_t MFindSupercutsONOFF::DefineTrainTestMatrix(
1532 const TString &name, MH3 &href,
1533 const Int_t howmanytrain, const Int_t howmanytest,
1534 const TString &filetrain, const TString &filetest)
1535{
1536 *fLog << "=============================================" << endl;
1537 *fLog << "fill training and test matrix from file '" << name
1538 << "', select " << howmanytrain
1539 << " training and " << howmanytest << " test events " << endl;
1540 if (!fUseOrigDistribution)
1541 {
1542 *fLog << " according to a distribution given by the MH3 object '"
1543 << href.GetName() << "'" << endl;
1544 }
1545 else
1546 {
1547 *fLog << " randomly" << endl;
1548 }
1549
1550
1551 MParList plist;
1552 MTaskList tlist;
1553
1554 MReadMarsFile read("Events", name);
1555 read.DisableAutoScheme();
1556
1557 MFEventSelector2 selector(href);
1558 selector.SetNumMax(howmanytrain+howmanytest);
1559 selector.SetName("selectTrainTest");
1560 selector.SetInverted();
1561 if (fUseOrigDistribution)
1562 {
1563 *fLog << "MFindSupercutsONFF; MFEventSelector2::SetUseOrigDistribution(Bool)"
1564 << " DOES NOT EXIST NOW..." << endl;
1565 // selector.SetUseOrigDistribution(kTRUE);
1566 }
1567
1568 MContinue cont(&selector);
1569 cont.SetName("ContTrainTest");
1570
1571 Double_t prob = ( (Double_t) howmanytrain )
1572 / ( (Double_t)(howmanytrain+howmanytest) );
1573 MFRandomSplit split(prob);
1574
1575 MFillH filltrain(fMatrixTrain);
1576 filltrain.SetFilter(&split);
1577 filltrain.SetName("fillMatrixTrain");
1578
1579
1580 // consider this event as candidate for a test event
1581 // only if event was not accepted as a training event
1582
1583 MContinue conttrain(&split);
1584 conttrain.SetName("ContTrain");
1585
1586 MFillH filltest(fMatrixTest);
1587 filltest.SetName("fillMatrixTest");
1588
1589
1590 //******************************
1591 // entries in MParList
1592
1593 plist.AddToList(&tlist);
1594 plist.AddToList(fCam);
1595 plist.AddToList(fMatrixTrain);
1596 plist.AddToList(fMatrixTest);
1597
1598 //******************************
1599
1600
1601 //******************************
1602 // entries in MTaskList
1603
1604 tlist.AddToList(&read);
1605 tlist.AddToList(&cont);
1606
1607 tlist.AddToList(&split);
1608 tlist.AddToList(&filltrain);
1609 tlist.AddToList(&conttrain);
1610
1611 tlist.AddToList(&filltest);
1612
1613 //******************************
1614
1615 MProgressBar bar;
1616 MEvtLoop evtloop;
1617 evtloop.SetParList(&plist);
1618 evtloop.SetName("EvtLoopMatrixTrainTest");
1619 evtloop.SetProgressBar(&bar);
1620
1621 Int_t maxev = -1;
1622 if (!evtloop.Eventloop(maxev))
1623 return kFALSE;
1624
1625 tlist.PrintStatistics(0, kTRUE);
1626
1627 fMatrixTrain->Print("SizeCols");
1628
1629 const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
1630 if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
1631 {
1632 *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events ("
1633 << generatedtrain
1634 << ") is incompatible with the no.of requested events ("
1635 << howmanytrain << ")" << endl;
1636 }
1637
1638 fMatrixTest->Print("SizeCols");
1639 const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
1640 if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
1641 {
1642 *fLog << "MFindSupercuts::DefineTrainTestMatrix; no.of generated events ("
1643 << generatedtest
1644 << ") is incompatible with the no.of requested events ("
1645 << howmanytest << ")" << endl;
1646 }
1647
1648
1649 *fLog << "training and test matrix were filled" << endl;
1650 *fLog << "=============================================" << endl;
1651
1652
1653 //--------------------------
1654 // write out training matrix
1655
1656 if (filetrain != "")
1657 {
1658 TFile filetr(filetrain, "RECREATE");
1659 fMatrixTrain->Write();
1660 filetr.Close();
1661
1662 *fLog << "MFindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
1663 << filetrain << "'" << endl;
1664 }
1665
1666 //--------------------------
1667 // write out test matrix
1668
1669 if (filetest != "")
1670 {
1671 TFile filete(filetest, "RECREATE", "");
1672 fMatrixTest->Write();
1673 filete.Close();
1674
1675 *fLog << "MFindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
1676 << filetest << "'" << endl;
1677 }
1678
1679 return kTRUE;
1680}
1681
1682// --------------------------------------------------------------------------
1683//
1684// Define the matrices for the training and test OFF sample respectively
1685//
1686//
1687//
1688
1689Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange(
1690 const TString &name,
1691 const Double_t whichfractiontrain,
1692 const Double_t whichfractiontest,
1693 Double_t ThetaMin, Double_t ThetaMax,
1694 const TString &filetrain, const TString &filetest)
1695
1696
1697
1698{
1699 *fLog << "=============================================" << endl;
1700 *fLog << "Fill training and testing OFF matrices from file '" << name
1701 << "', select a fraction of " << whichfractiontrain
1702 << " events for the training and a fraction of " << endl
1703 << whichfractiontest << " for the testing" << endl;
1704
1705
1706 MParList plist;
1707 MTaskList tlist;
1708
1709 MReadMarsFile read("Events", name);
1710 read.DisableAutoScheme();
1711
1712
1713
1714
1715
1716 // Cuts in Theta
1717
1718 // TString ThetaCutMinString ("ThetaOrig.fVal>");
1719 TString ThetaCutMinString ("MPointingPos.fZd>"); // new! // for magic
1720 ThetaCutMinString += ThetaMin;
1721 MContinue ThetaCutMin(ThetaCutMinString);
1722 ThetaCutMin.SetInverted();
1723
1724 //TString ThetaCutMaxString ("ThetaOrig.fVal<");
1725 TString ThetaCutMaxString ("MPointingPos.fZd<"); // new! // for magic
1726 ThetaCutMaxString += ThetaMax;
1727 MContinue ThetaCutMax(ThetaCutMaxString);
1728 ThetaCutMax.SetInverted();
1729
1730
1731
1732 // Cuts in Size,
1733
1734 TString SizeCutMinString ("MHillas.fSize>");
1735 SizeCutMinString += fSizeCutLow;
1736 MContinue SizeCutMin(SizeCutMinString);
1737 SizeCutMin.SetInverted();
1738
1739
1740 TString SizeCutMaxString ("MHillas.fSize<");
1741 SizeCutMaxString += fSizeCutUp;
1742 MContinue SizeCutMax(SizeCutMaxString);
1743 SizeCutMax.SetInverted();
1744
1745 // Cuts in Dist
1746 TString DistCutMinString ("MHillasSrc.fDist>");
1747 DistCutMinString += fDistCutLow;
1748 MContinue DistCutMin(DistCutMinString);
1749 DistCutMin.SetInverted();
1750
1751 TString DistCutMaxString ("MHillasSrc.fDist<");
1752 DistCutMaxString += fDistCutUp;
1753 MContinue DistCutMax(DistCutMaxString);
1754 DistCutMax.SetInverted();
1755
1756
1757 // Cuts in leakage
1758 TString LeakCutMaxString ("MNewImagePar.fLeakage1<");
1759 LeakCutMaxString += fLeakageMax;
1760 MContinue LeakCutMax(LeakCutMaxString);
1761 LeakCutMax.SetInverted();
1762
1763
1764
1765
1766 Double_t prob = whichfractiontrain /(whichfractiontrain+whichfractiontest);
1767
1768
1769
1770
1771 MFRandomSplit split(prob);
1772
1773 MFillH filltrain(fMatrixTrainOFF);
1774 filltrain.SetName("fillMatrixTrainOFF");
1775 filltrain.SetFilter(&split);
1776
1777
1778 // consider this event as candidate for a test event
1779 // only if event was not accepted as a training event
1780
1781 MContinue conttrain(&split);
1782 conttrain.SetName("ContTrain");
1783
1784 MFillH filltest(fMatrixTestOFF);
1785 filltest.SetName("fillMatrixTestOFF");
1786
1787
1788 //******************************
1789 // entries in MParList
1790
1791 plist.AddToList(&tlist);
1792 plist.AddToList(fCam);
1793 plist.AddToList(fMatrixTrainOFF);
1794 plist.AddToList(fMatrixTestOFF);
1795
1796 //******************************
1797 // entries in MTaskList
1798
1799 tlist.AddToList(&read);
1800 tlist.AddToList(&ThetaCutMin);
1801 tlist.AddToList(&ThetaCutMax);
1802
1803
1804 tlist.AddToList(&SizeCutMin);
1805 tlist.AddToList(&SizeCutMax);
1806
1807 tlist.AddToList(&DistCutMin);
1808 tlist.AddToList(&DistCutMax);
1809 tlist.AddToList(&LeakCutMax);
1810
1811
1812 tlist.AddToList(&split);
1813 tlist.AddToList(&filltrain);
1814
1815 tlist.AddToList(&conttrain);
1816 tlist.AddToList(&filltest);
1817
1818 //******************************
1819
1820 MProgressBar bar;
1821 MEvtLoop evtloop;
1822 evtloop.SetParList(&plist);
1823 evtloop.SetName("EvtLoopMatrixTrainTestOFF");
1824 evtloop.SetProgressBar(&bar);
1825
1826 Int_t maxev = -1;
1827 if (!evtloop.Eventloop(maxev))
1828 return kFALSE;
1829
1830 tlist.PrintStatistics(0, kTRUE);
1831
1832 fMatrixTrainOFF->Print("SizeCols");
1833
1834
1835 fMatrixTestOFF->Print("SizeCols");
1836
1837
1838 *fLog << "train and test matrix OFF were filled" << endl;
1839 *fLog << "=============================================" << endl;
1840
1841
1842 //--------------------------
1843 // write out training matrix
1844
1845 if (filetrain != "")
1846 {
1847 TFile filetr(filetrain, "RECREATE");
1848 fMatrixTrainOFF->Write();
1849 filetr.Close();
1850
1851 *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Training matrix was written onto file '"
1852 << filetrain << "'" << endl;
1853 }
1854
1855 //--------------------------
1856 // write out test matrix
1857
1858 if (filetest != "")
1859 {
1860 TFile filete(filetest, "RECREATE", "");
1861 fMatrixTestOFF->Write();
1862 filete.Close();
1863
1864 *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixOFFThetaRange; Test matrix was written onto file '"
1865 << filetest << "'" << endl;
1866 }
1867
1868 return kTRUE;
1869}
1870
1871
1872// --------------------------------------------------------------------------
1873//
1874// Define the matrices for the training and test sample respectively
1875//
1876//
1877//
1878Bool_t MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange(
1879 const TString &name,
1880 const Double_t whichfractiontrain,
1881 const Double_t whichfractiontest,
1882 Double_t ThetaMin, Double_t ThetaMax,
1883 const TString &filetrain, const TString &filetest)
1884
1885
1886{
1887 *fLog << "=============================================" << endl;
1888 *fLog << "Fill training and testing ON matrices from file '" << name
1889 << "', select a fraction of " << whichfractiontrain
1890 << " events for the training and a fraction of " << endl
1891 << whichfractiontest << " for the testing" << endl;
1892
1893
1894 MParList plist;
1895 MTaskList tlist;
1896
1897
1898 MReadMarsFile read("Events", name);
1899 read.DisableAutoScheme();
1900
1901
1902 // Cuts in Theta
1903
1904 //TString ThetaCutMinString ("ThetaOrig.fVal>");
1905 TString ThetaCutMinString ("MPointingPos.fZd>"); // for magic
1906 ThetaCutMinString += ThetaMin;
1907 MContinue ThetaCutMin(ThetaCutMinString);
1908 ThetaCutMin.SetInverted();
1909
1910 //TString ThetaCutMaxString ("ThetaOrig.fVal<");
1911 TString ThetaCutMaxString ("MPointingPos.fZd<"); // for magic
1912 ThetaCutMaxString += ThetaMax;
1913 MContinue ThetaCutMax(ThetaCutMaxString);
1914 ThetaCutMax.SetInverted();
1915
1916
1917
1918 // Cuts in Size,
1919
1920 TString SizeCutMinString ("MHillas.fSize>");
1921 SizeCutMinString += fSizeCutLow;
1922 MContinue SizeCutMin(SizeCutMinString);
1923 SizeCutMin.SetInverted();
1924
1925
1926 TString SizeCutMaxString ("MHillas.fSize<");
1927 SizeCutMaxString += fSizeCutUp;
1928 MContinue SizeCutMax(SizeCutMaxString);
1929 SizeCutMax.SetInverted();
1930
1931
1932 // Cuts in Dist
1933 TString DistCutMinString ("MHillasSrc.fDist>");
1934 DistCutMinString += fDistCutLow;
1935 MContinue DistCutMin(DistCutMinString);
1936 DistCutMin.SetInverted();
1937
1938 TString DistCutMaxString ("MHillasSrc.fDist<");
1939 DistCutMaxString += fDistCutUp;
1940 MContinue DistCutMax(DistCutMaxString);
1941 DistCutMax.SetInverted();
1942
1943
1944 // Cuts in leakage
1945 TString LeakCutMaxString ("MNewImagePar.fLeakage1<");
1946 LeakCutMaxString += fLeakageMax;
1947 MContinue LeakCutMax(LeakCutMaxString);
1948 LeakCutMax.SetInverted();
1949
1950
1951 Double_t prob = whichfractiontrain/(whichfractiontrain + whichfractiontest);
1952
1953
1954
1955 MFRandomSplit split(prob);
1956
1957 MFillH filltrain(fMatrixTrain);
1958 filltrain.SetName("fillMatrixTrain");
1959 filltrain.SetFilter(&split);
1960
1961
1962 // consider this event as candidate for a test event
1963 // only if event was not accepted as a training event
1964
1965 MContinue conttrain(&split);
1966 conttrain.SetName("ContTrain");
1967
1968 MFillH filltest(fMatrixTest);
1969 filltest.SetName("fillMatrixTest");
1970
1971
1972
1973
1974 //******************************
1975 // entries in MParList
1976
1977 plist.AddToList(&tlist);
1978 plist.AddToList(fCam);
1979 plist.AddToList(fMatrixTrain);
1980 plist.AddToList(fMatrixTest);
1981
1982 //******************************
1983 // entries in MTaskList
1984
1985 tlist.AddToList(&read);
1986 tlist.AddToList(&ThetaCutMin);
1987 tlist.AddToList(&ThetaCutMax);
1988
1989
1990
1991 tlist.AddToList(&SizeCutMin);
1992 tlist.AddToList(&SizeCutMax);
1993
1994 tlist.AddToList(&DistCutMin);
1995 tlist.AddToList(&DistCutMax);
1996 tlist.AddToList(&LeakCutMax);
1997
1998
1999
2000 tlist.AddToList(&split);
2001 tlist.AddToList(&filltrain);
2002
2003 tlist.AddToList(&conttrain);
2004 tlist.AddToList(&filltest);
2005
2006 //******************************
2007
2008 MProgressBar bar;
2009 MEvtLoop evtloop;
2010 evtloop.SetParList(&plist);
2011 evtloop.SetName("EvtLoopMatrixTrainTest");
2012 evtloop.SetProgressBar(&bar);
2013
2014 Int_t maxev = -1;
2015 if (!evtloop.Eventloop(maxev))
2016 return kFALSE;
2017
2018 tlist.PrintStatistics(0, kTRUE);
2019
2020 fMatrixTrain->Print("SizeCols");
2021
2022 fMatrixTest->Print("SizeCols");
2023
2024
2025 *fLog << "training and test matrix were filled" << endl;
2026 *fLog << "=============================================" << endl;
2027
2028
2029 //--------------------------
2030 // write out training matrix
2031
2032 if (filetrain != "")
2033 {
2034 TFile filetr(filetrain, "RECREATE");
2035 fMatrixTrain->Write();
2036 filetr.Close();
2037
2038 *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Train matrix was written onto file '"
2039 << filetrain << "'" << endl;
2040 }
2041
2042
2043 //--------------------------
2044 // write out test matrix
2045
2046 if (filetest != "")
2047 {
2048 TFile filete(filetest, "RECREATE", "");
2049 fMatrixTest->Write();
2050 filete.Close();
2051
2052 *fLog << "MFindSupercutsONOFF::DefineTrainTestMatrixThetaRange; Test matrix was written onto file '"
2053 << filetest << "'" << endl;
2054 }
2055
2056 return kTRUE;
2057}
2058
2059
2060
2061
2062
2063/// **********///
2064
2065// --------------------------------------------------------------------------
2066//
2067// Read only training matrices ON and OFF
2068//
2069//
2070
2071Bool_t MFindSupercutsONOFF::ReadMatrixTrain(const TString &filetrainON, const TString &filetrainOFF)
2072{
2073 //--------------------------
2074 // read in training matrix ON sample
2075
2076 TFile filetr(filetrainON);
2077 fMatrixTrain->Read("MatrixTrain");
2078 fMatrixTrain->Print("SizeCols");
2079
2080 *fLog << "MFindSupercuts::ReadMatrixTrain; Training matrix for ON sample was read in from file '"
2081 << filetrainON << "'" << endl;
2082 filetr.Close();
2083
2084
2085 // read in training matrix OFF sample
2086
2087 TFile filetrOFF(filetrainOFF);
2088 fMatrixTrainOFF->Read("MatrixTrainOFF");
2089 fMatrixTrainOFF->Print("SizeCols");
2090
2091 *fLog << "MFindSupercutsONOFF::ReadMatrixTrain; Training matrix for OFF sample was read in from file '"
2092 << filetrainOFF << "'" << endl;
2093 filetrOFF.Close();
2094
2095 return kTRUE;
2096
2097}
2098
2099// --------------------------------------------------------------------------
2100//
2101// Read only test matrices ON and OFF
2102//
2103//
2104
2105Bool_t MFindSupercutsONOFF::ReadMatrixTest(const TString &filetestON, const TString &filetestOFF)
2106{
2107 //--------------------------
2108 // read in testing matrix ON sample
2109
2110//--------------------------
2111 // read in test matrix for ON sample
2112
2113 TFile filete(filetestON);
2114 fMatrixTest->Read("MatrixTest");
2115 fMatrixTest->Print("SizeCols");
2116
2117 *fLog << "MFindSupercuts::ReadMatrixTest; Test matrix for ON sample was read in from file '"
2118 << filetestON << "'" << endl;
2119 filete.Close();
2120
2121
2122 //--------------------------
2123 // read in test matrix for OFF sample
2124
2125 TFile fileteOFF(filetestOFF);
2126 fMatrixTestOFF->Read("MatrixTestOFF");
2127 fMatrixTestOFF->Print("SizeCols");
2128
2129 *fLog << "MFindSupercutsONOFF::ReadMatrixTest; Test matrix for OFF sample was read in from file '"
2130 << filetestOFF << "'" << endl;
2131 filete.Close();
2132
2133
2134
2135 return kTRUE;
2136
2137}
2138
2139
2140/// **********///
2141
2142// --------------------------------------------------------------------------
2143//
2144// Read training and test matrices from files
2145//
2146//
2147
2148Bool_t MFindSupercutsONOFF::ReadMatrix(const TString &filetrain, const TString &filetest)
2149{
2150 //--------------------------
2151 // read in training matrix
2152
2153 TFile filetr(filetrain);
2154 fMatrixTrain->Read("MatrixTrain");
2155 fMatrixTrain->Print("SizeCols");
2156
2157 *fLog << "MFindSupercuts::ReadMatrix; Training matrix was read in from file '"
2158 << filetrain << "'" << endl;
2159 filetr.Close();
2160
2161
2162 //--------------------------
2163 // read in test matrix
2164
2165 TFile filete(filetest);
2166 fMatrixTest->Read("MatrixTest");
2167 fMatrixTest->Print("SizeCols");
2168
2169 *fLog << "MFindSupercuts::ReadMatrix; Test matrix was read in from file '"
2170 << filetest << "'" << endl;
2171 filete.Close();
2172
2173 return kTRUE;
2174}
2175
2176
2177// Read training and test matrices OFF from files
2178//
2179//
2180
2181Bool_t MFindSupercutsONOFF::ReadMatrixOFF(const TString &filetrain, const TString &filetest)
2182{
2183 //--------------------------
2184 // read in training matrix
2185
2186 TFile filetr(filetrain);
2187 fMatrixTrainOFF->Read("MatrixTrainOFF");
2188 fMatrixTrainOFF->Print("SizeCols");
2189
2190 *fLog << "MFindSupercutsONOFF::ReadMatrixOFF; Training matrix OFF was read in from file '"
2191 << filetrain << "'" << endl;
2192 filetr.Close();
2193
2194
2195 //--------------------------
2196 // read in test matrix
2197
2198 TFile filete(filetest);
2199 fMatrixTestOFF->Read("MatrixTestOFF");
2200 fMatrixTestOFF->Print("SizeCols");
2201
2202 *fLog << "MFindSupercutsONOFF::ReadMatrixONOFF; Test matrix OFF was read in from file '"
2203 << filetest << "'" << endl;
2204 filete.Close();
2205
2206 return kTRUE;
2207}
2208
2209
2210
2211// Function to compute the normalization factor for Train sample.
2212// The normalization factor is defined as the ratio of OFF/ON events.
2213
2214Bool_t MFindSupercutsONOFF::ComputeNormFactorTrain()
2215{
2216 Int_t EventsInTrainMatrixON = fMatrixTrain->GetM().GetNrows();
2217 Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows();
2218
2219 // Info...
2220 *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; " << endl
2221 << "EventsInTrainMatrixON = " << EventsInTrainMatrixON
2222 << " , EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF << endl;
2223
2224
2225 fNormFactorTrain = double(EventsInTrainMatrixON)/double(EventsInTrainMatrixOFF);
2226
2227
2228 *fLog << "MFindSupercutsONOFF::ComputeNormFactorTrain; fNormFactorTrain is : "
2229 << fNormFactorTrain << endl;
2230
2231 return kTRUE;
2232}
2233
2234
2235
2236
2237// Function to compute the normalization factor for Test sample.
2238// The normalization factor is defined as the ratio of OFF/ON events.
2239
2240Bool_t MFindSupercutsONOFF::ComputeNormFactorTest()
2241{
2242 Int_t EventsInTestMatrixON = fMatrixTest->GetM().GetNrows();
2243 Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows();
2244
2245 // Info...
2246 *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; " << endl
2247 << "EventsInTestMatrixON = " << EventsInTestMatrixON
2248 << " , EventsInTestMatrixOFF = " << EventsInTestMatrixOFF << endl;
2249
2250 fNormFactorTest = double(EventsInTestMatrixON)/double(EventsInTestMatrixOFF);
2251
2252 *fLog << "MFindSupercutsONOFF::ComputeNormFactorTest; fNormFactorTest is : "
2253 << fNormFactorTest << endl;
2254
2255 return kTRUE;
2256}
2257
2258
2259Bool_t MFindSupercutsONOFF::SetGammaEfficiency(Double_t gammaeff)
2260{
2261
2262 if (gammaeff < 0.0 || gammaeff >1.0)
2263 {
2264 *fLog << "MFindSupercutsONOFF::SetGammaEfficiency; " << endl
2265 << "The argument of SetGammaEfficiency ("
2266 << gammaeff << ") is outside the range 0-1" << endl
2267 << "This is not allowed..." << endl;
2268 return kFALSE;
2269
2270 }
2271
2272 fGammaEfficiency = gammaeff;
2273
2274 // TEMP INFO
2275 cout << "fGammaEfficiency has been set to : " << fGammaEfficiency << endl;
2276 // END TEMP
2277 return kTRUE;
2278
2279
2280}
2281
2282//------------------------------------------------------------------------
2283//
2284// Steering program for optimizing the supercuts
2285// ---------------------------------------------
2286//
2287// the criterion for the 'optimum' is
2288//
2289// - a maximum significance of the gamma signal in the alpha plot,
2290// in which the supercuts have been applied
2291//
2292// The various steps are :
2293//
2294// - setup the 2 event loops (ON and OFF) to be executed for each call to fcnSupercuts
2295// - call TMinuit to do the minimization of (-significance) :
2296// the fcnSupercuts function calculates the significance
2297// for the current cut values
2298// for this - the 2 alpha plots (ON and OFF) are produced in the 2 event loops,
2299// in which the SAME cuts have been applied
2300// - the significance of the gamma signal is extracted from the
2301// the 2 alpha plots (ON OFF) using MHFindSignificanceONOFF::FindSigmaONOFF.
2302//
2303//
2304// Needed as input : (to be set by the Set functions)
2305//
2306// - fFilenameParam name of file to which optimum values of the
2307// parameters are written
2308//
2309// - fHadronnessName name of container where MSupercutsCalcONOFF will
2310// put the supercuts hadronness for ON data
2311
2312// - fHadronnessNameOFF name of container where MSupercutsCalcONOFF will
2313// put the supercuts hadronness for OFF data
2314//
2315// - fAlphaSig, fAlphaBkgMin, fAlphaBkgMax ; Signal and Bkg region in alpha histo.
2316
2317// - for the minimization, the starting values of the parameters are taken
2318// - from file parSCinit (if it is != "")
2319// - or from the arrays params and/or steps
2320//
2321//----------------------------------------------------------------------
2322Bool_t MFindSupercutsONOFF::FindParams(TString parSCinit,
2323 TArrayD &params, TArrayD &steps)
2324{
2325 // Setup the event loops which will be executed in the
2326 // fcnSupercuts function of MINUIT
2327 //
2328 // parSCinit is the name of the file containing the initial values
2329 // of the parameters;
2330 // if parSCinit = "" 'params' and 'steps' are taken as initial values
2331
2332 *fLog << "=============================================" << endl;
2333 *fLog << "Setup event loop for fcnSupercuts" << endl;
2334
2335 if (fHadronnessName.IsNull())
2336 {
2337 *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for ON data is not defined... aborting"
2338 << endl;
2339 return kFALSE;
2340 }
2341
2342
2343 if (fHadronnessNameOFF.IsNull())
2344 {
2345 *fLog << "MFindSupercutsONOFF::FindParams; hadronness name for OFF data is not defined... aborting"
2346 << endl;
2347 return kFALSE;
2348 }
2349
2350
2351 if (fMatrixTrain == NULL)
2352 {
2353 *fLog << "MFindSupercuts::FindParams; training matrix is not defined... aborting"
2354 << endl;
2355 return kFALSE;
2356 }
2357
2358 if (fMatrixTrain->GetM().GetNrows() <= 0)
2359 {
2360 *fLog << "MFindSupercuts::FindParams; training matrix has no entries"
2361 << endl;
2362 return kFALSE;
2363 }
2364
2365
2366
2367 if (fMatrixTrainOFF == NULL)
2368 {
2369 *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF is not defined... aborting"
2370 << endl;
2371 return kFALSE;
2372 }
2373
2374 if (fMatrixTrainOFF->GetM().GetNrows() <= 0)
2375 {
2376 *fLog << "MFindSupercutsONOFF::FindParams; training matrix OFF has no entries"
2377 << endl;
2378 return kFALSE;
2379 }
2380
2381 if (fAlphaDistributionsRootFilename.IsNull())
2382 {
2383 *fLog << "MFindSupercutsONOFF::FindParams; fAlphaDistributionsRootFilename is not defined... program aborting..." << endl;
2384
2385 return kFALSE;
2386 }
2387
2388
2389
2390 //--------------------------------
2391 // create container for the supercut parameters
2392 // and set them to their initial values
2393 MSupercuts super;
2394
2395 // take initial values from file parSCinit
2396 if (parSCinit != "")
2397 {
2398 TFile inparam(parSCinit);
2399 super.Read("MSupercuts");
2400 inparam.Close();
2401 *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
2402 << parSCinit << endl;
2403 }
2404
2405 // take initial values from 'params' and/or 'steps'
2406 else if (params.GetSize() != 0 || steps.GetSize() != 0 )
2407 {
2408 if (params.GetSize() != 0)
2409 {
2410 *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from 'params'"
2411 << endl;
2412 super.SetParameters(params);
2413 }
2414 if (steps.GetSize() != 0)
2415 {
2416 *fLog << "MFindSupercutsONOFF::FindParams; initial step sizes are taken from 'steps'"
2417 << endl;
2418 super.SetStepsizes(steps);
2419 }
2420
2421 /*
2422 // TMP
2423 // Print parameters
2424 if (params.GetSize() == steps.GetSize())
2425 {
2426 *fLog << "MFindSupercutsONOFF; SC parameters and Setps are: " << endl;
2427 for (Int_t z = 0; z < params.GetSize(); z++)
2428 {
2429 cout << params[z] << setw(20) << steps[z] << endl;
2430 }
2431 }
2432
2433 // ENDTMP
2434
2435 */
2436
2437 /*
2438 // TMP2
2439
2440 TArrayD paramsBis = super.GetParameters();
2441 TArrayD stepsBis = super.GetStepsizes();
2442 if (paramsBis.GetSize() == stepsBis.GetSize())
2443 {
2444 *fLog << "MFindSupercutsONOFF; SC parametersBis and SetpsBis are: " << endl;
2445 for (Int_t z = 0; z < paramsBis.GetSize(); z++)
2446 {
2447 cout << paramsBis[z] << setw(20) << stepsBis[z] << endl;
2448 }
2449 }
2450
2451 // ENDTMP2
2452 */
2453
2454 else
2455 {
2456 *fLog << "MFindSupercutsONOFF; ERROR: params.GetSize() != steps.GetSize() "
2457 << endl;
2458
2459 }
2460
2461
2462
2463 }
2464 else
2465 {
2466 *fLog << "MFindSupercutsONOFF::FindParams; initial values and step sizes are taken from the MSupercuts constructor"
2467 << endl;
2468 }
2469
2470
2471 // Computation of fNormFactorTrain and creation of a container for
2472 // storing this value and include it in the MParList for being
2473 // used by function fcnSupercuts
2474
2475 if (!ComputeNormFactorTrain())
2476 {
2477 *fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl;
2478 return kFALSE;
2479 }
2480
2481
2482 *fLog << "MFindSupercutsONOFF::FindParams; fNormFactorTrain = " << fNormFactorTrain << endl;
2483
2484
2485
2486 MDataValue NormFactorCont(fNormFactorTrain);
2487 NormFactorCont.SetName("NormFactorTrain");
2488
2489
2490 // Value of fAlphaSig, fAlphaBkgMin and fAlphaBkgMax
2491 // are stored in MDataValue containers,
2492 // then they will be passed to the MParList and will be
2493 // accessible to fcnsupercuts
2494
2495 MDataValue AlphaSigCont(fAlphaSig);
2496 AlphaSigCont.SetName("AlphaSigValue");
2497
2498 MDataValue AlphaBkgMinCont(fAlphaBkgMin);
2499 AlphaBkgMinCont.SetName("AlphaBkgMinValue");
2500
2501 MDataValue AlphaBkgMaxCont(fAlphaBkgMax);
2502 AlphaBkgMaxCont.SetName("AlphaBkgMaxValue");
2503
2504
2505 MDataValue DegreeONCont (fDegreeON);
2506 DegreeONCont.SetName("DegreeON");
2507
2508 MDataValue DegreeOFFCont (fDegreeOFF);
2509 DegreeOFFCont.SetName("DegreeOFF");
2510
2511
2512
2513
2514
2515 // Value of fUseFittedQuantities is stored in container
2516 // and passed to the MParList and will be
2517 // accessible to fcnsupercuts
2518
2519
2520
2521 MDataValue UseFittedQuantitiesCont(fUseFittedQuantities);
2522 UseFittedQuantitiesCont.SetName("UseFittedQuantitiesValue");
2523
2524
2525 // Value of fNormFactorFromAlphaBkg is stored in container
2526 // and passed to the MParList and will be
2527 // accessible to fcnsupercuts
2528
2529 MDataValue UseNormFactorFromAlphaBkgCont(fNormFactorFromAlphaBkg);
2530 UseNormFactorFromAlphaBkgCont.SetName("UseNormFactorFromAlphaBkgValue");
2531
2532
2533
2534 // A vector of pointers to objects of the class MEvtLoop is defined.
2535 // One of the pointed loops will be used to compute ON data alpha plot
2536 // and the other for computing OFF data alpha plot.
2537
2538 MEvtLoop evtloopfcn[2] = {MEvtLoop("ONDataEvtLoopFCN"),
2539 MEvtLoop("OFFDataEvtLoopFCN")};
2540
2541
2542
2543
2544 // ******************************************************************
2545 // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR ON DATA IS SET
2546 // ******************************************************************
2547
2548 // -----------------------------------------------------------------
2549
2550 MParList parlistfcn;
2551 MTaskList tasklistfcn;
2552
2553 // loop over rows of matrix
2554 MMatrixLoop loop(fMatrixTrain);
2555
2556
2557 //--------------------------------
2558 // calculate supercuts hadronness
2559 fCalcHadTrain->SetHadronnessName(fHadronnessName);
2560
2561 // Boolean variable that controls in class MSupercutsCalcONOFF
2562 // wether the supercuts are stored or not
2563 // is set to the default value, kFALSE. (no storage of supercuts)
2564
2565 fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE);
2566
2567
2568
2569 // Set boolean variable that controls wether cuts are
2570 // dynamic or static
2571
2572 fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts);
2573
2574 if (!fUseStaticCuts)
2575 {
2576 // Set boolean variable that controls wether the theta variable
2577 // is used or not in the computation of the dynamical cuts
2578
2579 fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta);
2580 }
2581
2582 // apply the supercuts
2583 MF scfilter(fHadronnessName+".fHadronness>0.5");
2584 MContinue supercuts(&scfilter);
2585
2586 // plot |alpha|
2587 const TString mh3Name = "AlphaFcn";
2588 MBinning binsalpha("Binning"+mh3Name);
2589 //binsalpha.SetEdges(54, -12.0, 96.0);
2590 binsalpha.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
2591
2592
2593 *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;
2594
2595 // |alpha| is assumed to be in column 7 of the matrix
2596 MH3 alpha("MatrixTrain[7]");
2597 alpha.SetName(mh3Name);
2598
2599 MFillH fillalpha(&alpha);
2600
2601 // book histograms to be filled during the optimization
2602 // ! not in the event loop !
2603 MHCT1Supercuts plotsuper;
2604 plotsuper.SetupFill(&parlistfcn);
2605
2606 //******************************
2607 // entries in MParList (extension of old MParList)
2608
2609 parlistfcn.AddToList(&tasklistfcn);
2610 parlistfcn.AddToList(&super);
2611 parlistfcn.AddToList(&NormFactorCont);
2612 parlistfcn.AddToList(&AlphaSigCont);
2613 parlistfcn.AddToList(&AlphaBkgMinCont);
2614 parlistfcn.AddToList(&AlphaBkgMaxCont);
2615 parlistfcn.AddToList(&DegreeONCont);
2616 parlistfcn.AddToList(&UseFittedQuantitiesCont);
2617 parlistfcn.AddToList(&UseNormFactorFromAlphaBkgCont);
2618
2619 parlistfcn.AddToList(fCam);
2620 parlistfcn.AddToList(fMatrixTrain);
2621
2622 parlistfcn.AddToList(&binsalpha);
2623 parlistfcn.AddToList(&alpha);
2624
2625 parlistfcn.AddToList(&plotsuper);
2626
2627 //******************************
2628 // entries in MTaskList
2629
2630 tasklistfcn.AddToList(&loop);
2631 tasklistfcn.AddToList(fCalcHadTrain);
2632 tasklistfcn.AddToList(&supercuts);
2633 tasklistfcn.AddToList(&fillalpha);
2634
2635 //******************************
2636
2637
2638 // &evtloopfcn[0] = new MEvtLoop("ONDataEvtLoopFCN");
2639 evtloopfcn[0].SetParList(&parlistfcn);
2640 // MEvtLoop evtloopfcn("EvtLoopFCN");
2641 // evtloopfcn.SetParList(&parlistfcn);
2642 *fLog << "Event loop for computing ON data alpha histo in fcnSupercuts has been setup" << endl;
2643
2644
2645//---- End of setting of loop for compuing ON data alpha plot ---------------
2646
2647
2648
2649 // ******************************************************************
2650 // EVENT LOOP FOR COMPUTING ALPHA HISTOGRAM FOR OFF DATA IS SET
2651 // ******************************************************************
2652
2653 // -----------------------------------------------------------------
2654
2655 MParList parlistfcnOFF;
2656 MTaskList tasklistfcnOFF;
2657
2658 // loop over rows of matrix
2659 MMatrixLoop loopOFF(fMatrixTrainOFF);
2660
2661
2662 //--------------------------------
2663 // calculate supercuts hadronness
2664 fCalcHadTrainOFF->SetHadronnessName(fHadronnessNameOFF);
2665
2666
2667 // Boolean variable that controls in class MSupercutsCalcONOFF
2668 // wether the supercuts are stored or not
2669 // is set to the default value, kFALSE. (no storage of supercuts)
2670
2671 fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE);
2672
2673
2674
2675 // Set boolean variable that controls wether cuts are
2676 // dynamic or static
2677
2678 fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts);
2679
2680
2681 if (!fUseStaticCuts)
2682 {
2683 // Set boolean variable that controls wether the theta variable
2684 // is used or not in the computation of the dynamica cuts
2685
2686 fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta);
2687 }
2688
2689 // apply the supercuts
2690 MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
2691 MContinue supercutsOFF(&scfilterOFF);
2692
2693 // plot |alpha|
2694 const TString mh3NameOFF = "AlphaOFFFcn";
2695 MBinning binsalphaOFF("Binning"+mh3NameOFF);
2696 //binsalphaOFF.SetEdges(54, -12.0, 96.0);
2697 binsalphaOFF.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
2698
2699
2700 *fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;
2701
2702 // |alpha| is assumed to be in column 7 of the matrix
2703 MH3 alphaOFF("MatrixTrainOFF[7]");
2704 alphaOFF.SetName(mh3NameOFF);
2705
2706 MFillH fillalphaOFF(&alphaOFF);
2707
2708
2709
2710 //******************************
2711 // entries in MParList (extension of old MParList)
2712
2713 parlistfcnOFF.AddToList(&tasklistfcnOFF);
2714 parlistfcnOFF.AddToList(&super);
2715 parlistfcnOFF.AddToList(fCam);
2716 parlistfcnOFF.AddToList(fMatrixTrainOFF);
2717
2718 parlistfcnOFF.AddToList(&binsalphaOFF);
2719 parlistfcnOFF.AddToList(&alphaOFF);
2720
2721 parlistfcnOFF.AddToList(&DegreeOFFCont);
2722
2723
2724 //******************************
2725 // entries in MTaskList
2726
2727 tasklistfcnOFF.AddToList(&loopOFF);
2728 tasklistfcnOFF.AddToList(fCalcHadTrainOFF);
2729 tasklistfcnOFF.AddToList(&supercutsOFF);
2730 tasklistfcnOFF.AddToList(&fillalphaOFF);
2731
2732 //******************************
2733
2734
2735 // &evtloopfcn[1] = new MEvtLoop("OFFDataEvtLoopFCN");
2736 evtloopfcn[1].SetParList(&parlistfcnOFF);
2737 // MEvtLoop evtloopfcn("EvtLoopFCN");
2738 // evtloopfcn.SetParList(&parlistfcn);
2739 *fLog << "Event loop for computing OFF data alpha histo in fcnSupercuts has been setup" << endl;
2740
2741
2742//---- End of setting of loop for compuing OFF data alpha plot ---------------
2743
2744
2745 // address of evtloopfcn is used in CallMinuit()
2746
2747
2748 //-----------------------------------------------------------------------
2749 //
2750 //---------- Start of minimization part --------------------
2751 //
2752 // Do the minimization with MINUIT
2753 //
2754 // Be careful: This is not thread safe
2755 //
2756 *fLog << "========================================================" << endl;
2757 *fLog << "Start minimization for supercuts" << endl;
2758
2759
2760 // -------------------------------------------
2761 // prepare call to MINUIT
2762 //
2763
2764 // get initial values of parameters
2765 fVinit = super.GetParameters();
2766 fStep = super.GetStepsizes();
2767
2768 TString name[fVinit.GetSize()];
2769 fStep.Set(fVinit.GetSize());
2770 fLimlo.Set(fVinit.GetSize());
2771 fLimup.Set(fVinit.GetSize());
2772 fFix.Set(fVinit.GetSize());
2773
2774 fNpar = fVinit.GetSize();
2775
2776 for (UInt_t i=0; i<fNpar; i++)
2777 {
2778 name[i] = "p";
2779 name[i] += i+1;
2780 //fStep[i] = TMath::Abs(fVinit[i]/10.0);
2781 fLimlo[i] = -100.0;
2782 fLimup[i] = 100.0;
2783 fFix[i] = 0;
2784 }
2785
2786
2787
2788
2789
2790 if(fSetLimitsToSomeMinuitParams)
2791 {
2792
2793 // Limits for some Minuit parameters are set. For the time being the values are set in the constructor
2794 // of the class. One MUST be very careful to set limits such that the expected final values
2795 // (optimized values) are far away from the limits.
2796 // Limits are set to those Minuit parameters not depending in Size,dist and theta (static values).
2797 // Limits are set only to Dist, Length and Width.
2798
2799
2800 *fLog << endl
2801 <<"Limits set to hillas parameters Length, Width, Dist and Leakage1 for the optimization with Minuit."
2802 << endl;
2803
2804
2805
2806 fLimup[0] = fMinuitLengthUPUpperLimit;
2807 fLimlo[0] = fMinuitLengthUPLowerLimit;
2808 fLimup[8] = fMinuitLengthLOWUpperLimit;
2809 fLimlo[8] = fMinuitLengthLOWLowerLimit;
2810 fLimup[16] = fMinuitWidthUPUpperLimit;
2811 fLimlo[16] = fMinuitWidthUPLowerLimit;
2812 fLimup[24] = fMinuitWidthLOWUpperLimit;
2813 fLimlo[24] = fMinuitWidthLOWLowerLimit;
2814 fLimup[32] = fMinuitDistUPUpperLimit;
2815 fLimlo[32] = fMinuitDistUPLowerLimit;
2816 fLimup[40] = fMinuitDistLOWUpperLimit;
2817 fLimlo[40] = fMinuitDistLOWLowerLimit;
2818
2819 fLimup[80] = fMinuitLeakage1UPUpperLimit;
2820 fLimlo[80] = fMinuitLeakage1UPLowerLimit;
2821
2822
2823
2824 }
2825
2826
2827
2828 // Wolfgang sets manyally some parameters.
2829 // That might give an error if parameters are changed...
2830 // To be changed in future...
2831
2832 // these parameters make no sense, fix them at 0.0
2833 fVinit[33] = 0.0;
2834 fStep[33] = 0.0;
2835 fFix[33] = 1;
2836
2837 fVinit[36] = 0.0;
2838 fStep[36] = 0.0;
2839 fFix[36] = 1;
2840
2841 fVinit[39] = 0.0;
2842 fStep[39] = 0.0;
2843 fFix[39] = 1;
2844
2845 fVinit[41] = 0.0;
2846 fStep[41] = 0.0;
2847 fFix[41] = 1;
2848
2849 fVinit[44] = 0.0;
2850 fStep[44] = 0.0;
2851 fFix[44] = 1;
2852
2853 fVinit[47] = 0.0;
2854 fStep[47] = 0.0;
2855 fFix[47] = 1;
2856
2857
2858
2859 // ADDITIONAL PARAMETERS ARE FIXED ACCORDING TO THE
2860 // VALUES OF THE BOOLEAN VARIABLES fUseStaticCuts, AND
2861 // fNotUseTheta
2862
2863
2864 if (fUseStaticCuts)
2865 {
2866 // Static cuts will be used; all parameters with index
2867 // 1-7 are fixed for ALL variables (width, length, ...)
2868
2869 for (UInt_t i=0; i<fNpar; i = i+8)
2870 {
2871 // At some point, the number of paraemters in the
2872 // dynamical cuts parameterezattion (currently 7)
2873 // will become a variable that can be set from
2874 // outside the class
2875
2876 for (UInt_t j=i+1; j<=i+7; j++)
2877 {
2878 fVinit[j] = 0.0;
2879 fStep[j] = 0.0;
2880 fFix[j] = 1;
2881 }
2882 }
2883 }
2884 else
2885 {
2886 if(fNotUseTheta)
2887 {
2888 // Theta is not used in the parameterization of the dynamical cut
2889 // Parameters with index 2 and 5 are fixed for all variables
2890 // (width, length, ...)
2891
2892 for (UInt_t i=0; i<fNpar; i = i+8)
2893 {
2894 // At some point, the number of paraemters in the
2895 // dynamical cuts parameterezattion (currently 7)
2896 // will become a variable that can be set from
2897 // outside the class
2898
2899 for (UInt_t j=i+2; j<=i+7; j = j+3)
2900 {
2901 fVinit[j] = 0.0;
2902 fStep[j] = 0.0;
2903 fFix[j] = 1;
2904 }
2905 }
2906
2907
2908 }
2909
2910 if(!fUseDist)
2911 {
2912 // parameterization using dist is removed
2913 // Parameters with index 1 and 4 and 7 are fixed for all variables
2914 // (width, length, ...)
2915 for (UInt_t i=0; i<fNpar; i = i+8)
2916 {
2917 // At some point, the number of paraemters in the
2918 // dynamical cuts parameterezattion (currently 7)
2919 // will become a variable that can be set from
2920 // outside the class
2921
2922 for (UInt_t j=i+1; j<=i+7; j = j+3)
2923 {
2924 fVinit[j] = 0.0;
2925 fStep[j] = 0.0;
2926 fFix[j] = 1;
2927 }
2928 }
2929
2930
2931
2932 }
2933
2934
2935
2936
2937 }
2938
2939
2940
2941
2942
2943
2944 // vary only first 48 parameters
2945 //for (UInt_t i=0; i<fNpar; i++)
2946 //{
2947 // if (i >= 48)
2948 // {
2949 // fStep[i] = 0.0;
2950 // fFix[i] = 1;
2951 // }
2952 //}
2953
2954 // -------------------------------------------
2955 // call MINUIT
2956
2957 Bool_t rc;
2958
2959 if(fSkipOptimization)
2960 {
2961 *fLog << "Parameter optimization is skipped. Using previously optimized parameters."
2962 << endl;
2963
2964
2965 rc = kTRUE;
2966
2967
2968 if(fUseInitialSCParams)
2969 {
2970
2971 // write Initial parameter values onto root file filenameParam
2972 // so that later they can be applied on the data
2973
2974 *fLog << "Initial SC parameter values are written onto file '"
2975 << fFilenameParam << "' :" << endl;
2976
2977
2978 TFile outparam(fFilenameParam, "RECREATE");
2979 super.Write();
2980 outparam.Close();
2981
2982
2983
2984 const TArrayD &check = super.GetParameters();
2985 for (Int_t i=0; i<check.GetSize(); i++)
2986 *fLog << check[i] << ", ";
2987 *fLog << endl;
2988
2989
2990 }
2991
2992
2993 }
2994 else
2995 {
2996 TStopwatch clock;
2997 clock.Start();
2998
2999 *fLog << "before calling CallMinuit" << endl;
3000
3001 MMinuitInterface inter;
3002 rc = inter.CallMinuit(fcnSupercuts, name,
3003 fVinit, fStep, fLimlo, fLimup, fFix,
3004 evtloopfcn, "SIMPLEX", kFALSE);
3005
3006 *fLog << "after calling CallMinuit" << endl;
3007
3008 *fLog << "Time spent for the minimization in MINUIT : " << endl;;
3009 clock.Stop();
3010 clock.Print();
3011
3012 // plotsuper.DrawClone();
3013
3014 if (!rc)
3015 {
3016 *fLog << "************* ERROR !!! *********************" << endl
3017 << "Minimization could not finish properly..." << endl
3018 << "MMinuitInterface.CallMinuit() returned kFALSE" << endl;
3019 return kFALSE;
3020 }
3021
3022
3023
3024
3025 *fLog << "Minimization for supercuts finished" << endl;
3026 *fLog << "========================================================" << endl;
3027
3028
3029 // -----------------------------------------------------------------
3030 // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put
3031 // into MSupercuts
3032
3033 // write optimum parameter values onto file filenameParam
3034
3035 TFile outparam(fFilenameParam, "RECREATE");
3036 super.Write();
3037 outparam.Close();
3038
3039 *fLog << "Optimum parameter values for supercuts were written onto file '"
3040 << fFilenameParam << "' :" << endl;
3041
3042 const TArrayD &check = super.GetParameters();
3043 for (Int_t i=0; i<check.GetSize(); i++)
3044 *fLog << check[i] << ", ";
3045 *fLog << endl;
3046
3047
3048
3049
3050 }
3051
3052
3053
3054
3055
3056 *fLog << "End of supercuts optimization part" << endl;
3057 *fLog << "======================================================" << endl;
3058
3059
3060 *fLog << "Applying the optimized supercuts on the TRAIN sample and "
3061 << "producing plots related to the output results" << endl;
3062
3063 *fLog << "======================================================" << endl;
3064
3065 if (!TestParamsOnTrainSample())
3066 {
3067 *fLog << "MFindSupercutsONOFF::FindParams;"
3068 << "TestParamsOnTrainSample failed"
3069 << endl;
3070
3071 return kFALSE;
3072 }
3073
3074
3075
3076 return kTRUE;
3077}
3078
3079
3080
3081
3082
3083
3084
3085// -----------------------------------------------------------------------
3086//
3087// Test the supercuts on the test sample using ON and OFF data
3088// 2 loops are used to fill the 2 alpha distributions
3089// Supercuts applied (LengthUp, LengthLow...)
3090// to all the individual events, as well as
3091// the features of the shower images (Length,Width...) are
3092// stored in 2 TTree objects (for ON and OFF data) in the
3093// root file specified by variable "fAlphaDistributionsRootFilename"
3094
3095
3096
3097Bool_t MFindSupercutsONOFF::TestParamsOnTestSample()
3098{
3099 if (fMatrixTest->GetM().GetNrows() <= 0)
3100 {
3101 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix ON has no entries"
3102 << endl;
3103 return kFALSE;
3104 }
3105
3106
3107 if (fMatrixTestOFF->GetM().GetNrows() <= 0)
3108 {
3109 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; test matrix OFF has no entries"
3110 << endl;
3111 return kFALSE;
3112 }
3113
3114 // ------------- TEST the supercuts ------------------------------
3115 //
3116 *fLog << "Test the supercuts on the test sample ON and OFF" << endl;
3117
3118 // -----------------------------------------------------------------
3119 // read optimum parameter values from file filenameParam
3120 // into array 'supercutsPar'
3121
3122 TFile inparam(fFilenameParam);
3123 MSupercuts scin;
3124 scin.Read("MSupercuts");
3125 inparam.Close();
3126
3127 *fLog << "Optimum parameter values for supercuts were read from file '";
3128 *fLog << fFilenameParam << "' :" << endl;
3129
3130 const TArrayD &supercutsPar = scin.GetParameters();
3131 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
3132 *fLog << supercutsPar[i] << ", ";
3133 *fLog << endl;
3134 //---------------------------
3135
3136
3137 // -----------------------------------------------------------------
3138 if (fHadronnessName.IsNull())
3139 {
3140 *fLog << "MFindSupercutsONOFF::TestParamsONOFF; hadronness name for ON data is not defined... aborting";
3141 *fLog << endl;
3142 return kFALSE;
3143 }
3144
3145 if (fHadronnessNameOFF.IsNull())
3146 {
3147 *fLog << "MFindSupercutsONOFF::TestParamsONOFF; hadronness name for OFF data is not defined... aborting";
3148 *fLog << endl;
3149 return kFALSE;
3150 }
3151
3152 if(fTestONSupercutsAppliedName.IsNull())
3153 {
3154 *fLog << "MFindSupercutsONOFF::TestParamsONOFF; "
3155 << " MTSupercutsApplied tree name for ON Test data is not defined... aborting"
3156 << endl;
3157 return kFALSE;
3158 }
3159
3160 if(fTestOFFSupercutsAppliedName.IsNull())
3161 {
3162 *fLog << "MFindSupercutsONOFF::TestParamsONOFF; "
3163 << " MTSupercutsApplied container name for OFF Test data is not defined... aborting"
3164 << endl;
3165 return kFALSE;
3166 }
3167
3168
3169 MSupercuts supercuts;
3170 supercuts.SetParameters(supercutsPar);
3171
3172 // Objects containing the trees used to store supercuts applied
3173 // TTree Branches are also created
3174 MTSupercutsApplied supapptestON(fTestONSupercutsAppliedName.Data(), "");
3175 supapptestON.CreateTreeBranches();
3176
3177 MTSupercutsApplied supapptestOFF(fTestOFFSupercutsAppliedName.Data(), "");
3178 supapptestOFF.CreateTreeBranches();
3179
3180
3181
3182
3183 fCalcHadTestOFF -> SetHadronnessName(fHadronnessNameOFF);
3184 fCalcHadTest -> SetHadronnessName(fHadronnessName);
3185
3186
3187 // Settings related to the storage of teh supercuts applied
3188 fCalcHadTest -> SetSupercutsAppliedName(fTestONSupercutsAppliedName);
3189 fCalcHadTest -> SetStoreAppliedSupercuts(kTRUE);
3190
3191 fCalcHadTestOFF -> SetSupercutsAppliedName(fTestOFFSupercutsAppliedName);
3192 fCalcHadTestOFF -> SetStoreAppliedSupercuts(kTRUE);
3193
3194
3195
3196 // Set boolean variable that controls wether cuts are
3197 // dynamic or static
3198
3199 fCalcHadTest -> SetVariableUseStaticCuts(fUseStaticCuts);
3200 fCalcHadTestOFF -> SetVariableUseStaticCuts(fUseStaticCuts);
3201
3202 if (!fUseStaticCuts)
3203 {
3204 // Set boolean variable that controls wether the theta variable
3205 // is used or not in the computation of the dynamical cuts
3206
3207 fCalcHadTest -> SetVariableNotUseTheta(fNotUseTheta);
3208 fCalcHadTestOFF -> SetVariableNotUseTheta(fNotUseTheta);
3209 }
3210
3211
3212 // apply the supercuts to OFF data
3213
3214 MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
3215 MContinue applysupercutsOFF(&scfilterOFF);
3216
3217
3218 // apply the supercuts to ON data
3219 MF scfilter(fHadronnessName+".fHadronness>0.5");
3220 MContinue applysupercuts(&scfilter);
3221
3222
3223
3224 // ****************************************************
3225 // Filling OFF alpha distribution
3226
3227 MParList parlistOFF;
3228 MTaskList tasklistOFF;
3229
3230 MMatrixLoop loopOFF(fMatrixTestOFF);
3231
3232 // plot alpha before applying the supercuts
3233 const TString mh3NameOFFB = "AlphaOFFBefore";
3234 MBinning binsalphaOFFbef("Binning"+mh3NameOFFB);
3235 // binsalphaOFFbef.SetEdges(54, -12.0, 96.0);
3236 binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3237
3238 // |alpha| is assumed to be in column 7 of the matrix
3239 MH3 alphaOFFbefore("MatrixTestOFF[7]");
3240 alphaOFFbefore.SetName(mh3NameOFFB);
3241
3242
3243 MFillH fillalphaOFFbefore(&alphaOFFbefore);
3244 fillalphaOFFbefore.SetName("FillAlphaOFFBefore");
3245
3246
3247
3248 // plot alpha OFF after applying the supercuts
3249 const TString mh3NameOFFA = "AlphaOFFAfter";
3250 MBinning binsalphaOFFaft("Binning"+mh3NameOFFA);
3251 // binsalphaOFFaft.SetEdges(54, -12.0, 96.0);
3252 binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3253
3254 MH3 alphaOFFafter("MatrixTestOFF[7]");
3255 alphaOFFafter.SetName(mh3NameOFFA);
3256
3257
3258 MFillH fillalphaOFFafter(&alphaOFFafter);
3259 fillalphaOFFafter.SetName("FillAlphaOFFAfter");
3260
3261 //******************************
3262 // entries in MParList
3263
3264 parlistOFF.AddToList(&tasklistOFF);
3265
3266 parlistOFF.AddToList(&supercuts);
3267
3268 parlistOFF.AddToList(&supapptestOFF);
3269
3270 parlistOFF.AddToList(fCam);
3271 parlistOFF.AddToList(fMatrixTestOFF);
3272
3273 parlistOFF.AddToList(&binsalphaOFFbef);
3274 parlistOFF.AddToList(&alphaOFFbefore);
3275
3276 parlistOFF.AddToList(&binsalphaOFFaft);
3277 parlistOFF.AddToList(&alphaOFFafter);
3278
3279
3280 //******************************
3281 // Entries in MtaskList
3282 tasklistOFF.AddToList(&loopOFF);
3283 tasklistOFF.AddToList(&fillalphaOFFbefore);
3284
3285 tasklistOFF.AddToList(fCalcHadTestOFF);
3286 tasklistOFF.AddToList(&applysupercutsOFF);
3287
3288 tasklistOFF.AddToList(&fillalphaOFFafter);
3289
3290 //******************************
3291
3292 MProgressBar barOFF;
3293 MEvtLoop evtloopOFF;
3294 evtloopOFF.SetParList(&parlistOFF);
3295 evtloopOFF.SetName("EvtLoopTestParamsOFF");
3296 evtloopOFF.SetProgressBar(&barOFF);
3297 Int_t maxeventsOFF = -1;
3298 if (!evtloopOFF.Eventloop(maxeventsOFF))
3299 return kFALSE;
3300
3301 tasklistOFF.PrintStatistics(0, kTRUE);
3302
3303
3304 //-------------------------------------------
3305 // draw the alpha plots
3306
3307 MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3");
3308 TH1 &alphaOFFHistb = alphaOFFBefore->GetHist();
3309 TString alphaOFFHistbName ("TestAlphaOFFBeforeCuts");
3310
3311 if (fThetaRangeString.IsNull())
3312 {
3313 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
3314 << " fThetaRangeString is NOT defined..." << endl;
3315 alphaOFFHistbName += ("ThetaRangeStringUndefined");
3316 }
3317 else
3318 {
3319 alphaOFFHistbName += (fThetaRangeString);
3320 }
3321
3322
3323 alphaOFFHistb.SetXTitle("|\\alpha| [\\circ]");
3324 alphaOFFHistb.SetName(alphaOFFHistbName);
3325
3326
3327 MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3");
3328 TH1 &alphaOFFHista = alphaOFFAfter->GetHist();
3329 TString alphaOFFHistaName ("TestAlphaOFFAfterCuts");
3330
3331 if (fThetaRangeString.IsNull())
3332 {
3333 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
3334 << " fThetaRangeString is NOT defined..." << endl;
3335 alphaOFFHistaName += ("ThetaRangeStringUndefined");
3336 }
3337 else
3338 {
3339 alphaOFFHistaName += (fThetaRangeString);
3340 }
3341
3342
3343 alphaOFFHista.SetXTitle("|\\alpha| [\\circ]");
3344 alphaOFFHista.SetName(alphaOFFHistaName);
3345
3346
3347
3348
3349
3350 // *****************************************************
3351
3352
3353
3354 //-------------------------------------------
3355
3356
3357 MParList parlist2;
3358 MTaskList tasklist2;
3359
3360 MMatrixLoop loop(fMatrixTest);
3361
3362 // plot alpha before applying the supercuts
3363 const TString mh3NameB = "AlphaBefore";
3364 MBinning binsalphabef("Binning"+mh3NameB);
3365 //binsalphabef.SetEdges(54, -12.0, 96.0);
3366 binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3367
3368 // |alpha| is assumed to be in column 7 of the matrix
3369 MH3 alphabefore("MatrixTest[7]");
3370 alphabefore.SetName(mh3NameB);
3371
3372
3373 MFillH fillalphabefore(&alphabefore);
3374 fillalphabefore.SetName("FillAlphaBefore");
3375
3376
3377 // plot alpha after applying the supercuts
3378 const TString mh3NameA = "AlphaAfter";
3379 MBinning binsalphaaft("Binning"+mh3NameA);
3380 // binsalphaaft.SetEdges(54, -12.0, 96.0);
3381 binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3382
3383 MH3 alphaafter("MatrixTest[7]");
3384 alphaafter.SetName(mh3NameA);
3385
3386
3387 MFillH fillalphaafter(&alphaafter);
3388 fillalphaafter.SetName("FillAlphaAfter");
3389
3390 //******************************
3391 // entries in MParList
3392
3393 parlist2.AddToList(&tasklist2);
3394
3395 parlist2.AddToList(&supercuts);
3396
3397 parlist2.AddToList(&supapptestON);
3398
3399 parlist2.AddToList(fCam);
3400 parlist2.AddToList(fMatrixTest);
3401
3402 parlist2.AddToList(&binsalphabef);
3403 parlist2.AddToList(&alphabefore);
3404
3405 parlist2.AddToList(&binsalphaaft);
3406 parlist2.AddToList(&alphaafter);
3407
3408 //******************************
3409 // entries in MTaskList
3410
3411 tasklist2.AddToList(&loop);
3412 tasklist2.AddToList(&fillalphabefore);
3413
3414 tasklist2.AddToList(fCalcHadTest);
3415 tasklist2.AddToList(&applysupercuts);
3416
3417 tasklist2.AddToList(&fillalphaafter);
3418
3419 //******************************
3420
3421 MProgressBar bar2;
3422 MEvtLoop evtloop2;
3423 evtloop2.SetParList(&parlist2);
3424 evtloop2.SetName("EvtLoopTestParams");
3425 evtloop2.SetProgressBar(&bar2);
3426 Int_t maxevents2 = -1;
3427 if (!evtloop2.Eventloop(maxevents2))
3428 return kFALSE;
3429
3430 tasklist2.PrintStatistics(0, kTRUE);
3431
3432
3433 //-------------------------------------------
3434 // draw the alpha plots
3435
3436 MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
3437 TH1 &alphaHist1 = alphaBefore->GetHist();
3438 TString alphaHist1Name ("TestAlphaBeforeCuts");
3439
3440 if (fThetaRangeString.IsNull())
3441 {
3442 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
3443 << " fThetaRangeString is NOT defined..." << endl;
3444 alphaHist1Name += ("ThetaRangeStringUndefined");
3445 }
3446 else
3447 {
3448 alphaHist1Name += (fThetaRangeString);
3449 }
3450
3451
3452
3453
3454 alphaHist1.SetName(alphaHist1Name);
3455 alphaHist1.SetXTitle("|\\alpha| [\\circ]");
3456
3457
3458
3459 MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
3460 TH1 &alphaHist2 = alphaAfter->GetHist();
3461 TString alphaHist2Name ("TestAlphaAfterCuts");
3462
3463
3464 if (fThetaRangeString.IsNull())
3465 {
3466 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
3467 << " fThetaRangeString is NOT defined..." << endl;
3468 alphaHist2Name += ("ThetaRangeStringUndefined");
3469 }
3470 else
3471 {
3472 alphaHist2Name += (fThetaRangeString);
3473 }
3474
3475 alphaHist2.SetXTitle("|\\alpha| [\\circ]");
3476 alphaHist2.SetName(alphaHist2Name);
3477
3478
3479
3480 //fPsFilename -> NewPage();
3481
3482
3483 // Canvas is deleted after saving the plot in ps file
3484 // This is to allow for GRID analysis
3485
3486 TCanvas *c = new TCanvas("TestAlphaBeforeAfterSC",
3487 "TestAlphaBeforeAfterSC",
3488 600, 600);
3489 c->Divide(2,2);
3490
3491
3492
3493 c->cd(1);
3494 alphaOFFHistb.DrawCopy();
3495
3496 c->cd(2);
3497 alphaOFFHista.DrawCopy();
3498
3499
3500
3501 c->cd(3);
3502 alphaHist1.DrawCopy();
3503
3504 c->cd(4);
3505 alphaHist2.DrawCopy();
3506
3507 c->Modified();
3508 c-> Update();
3509
3510
3511 // ********************************************
3512 // TMP solutions while the TPostScript thing is not working
3513 // PsFileName for storing these plots is derived
3514 // from PsFilenameString.
3515
3516
3517 cout << "Alpha distributions will be saved in PostScript file " ;
3518
3519
3520 if (!fPsFilenameString.IsNull())
3521 {
3522 TString filename = (fPsFilenameString);
3523 filename += ("AlphaTESTSampleONOFFBeforeAfterCuts.ps");
3524 cout << filename << endl;
3525 c -> SaveAs(filename);
3526 }
3527
3528 // END OF TEMPORAL SOLUTION
3529 // ********************************************
3530
3531
3532 // Canvas is deleted after saving the plot in ps file
3533 // This is to allow for GRID analysis
3534 delete c;
3535
3536 //-------------------------------------------
3537 // fit alpha distribution to get the number of excess events and
3538 // calculate significance of gamma signal in the alpha plot
3539
3540 const Double_t alphasig = fAlphaSig;
3541 const Double_t alphamin = fAlphaBkgMin;
3542 // alpha min for bkg region in ON data
3543
3544 const Double_t alphamax = fAlphaBkgMax;
3545// const Int_t degree = 2;
3546// const Int_t degreeOFF = 2;
3547 const Bool_t drawpoly = kTRUE;
3548 const Bool_t fitgauss = kTRUE;
3549 const Bool_t print = kTRUE;
3550
3551 const Bool_t saveplots = kTRUE; // Save plots in Psfile
3552
3553
3554
3555 if (!ComputeNormFactorTest())
3556 {
3557 *fLog << "Normalization factor for test sample (ON-OFF) couold not be computed. Aborting..." << endl;
3558 return kFALSE;
3559 }
3560
3561
3562
3563 // Normalization factor factor is corrected using the estimated
3564 // number of excess events and the total number of ON events
3565
3566
3567
3568 if (fTuneNormFactor)
3569 {
3570
3571
3572 // Run MHFindSignificanceONOFF::FindSigmaONOFF
3573 // to get estimation of number of excess gamma events.
3574 // This number will be used in the correction of the
3575 // normalization factor
3576
3577 // No plots will be produced this time...
3578
3579 const Bool_t drawpolyTest = kFALSE;
3580 const Bool_t fitgaussTest = kFALSE;
3581 const Bool_t printTest = kFALSE;
3582
3583 const Bool_t saveplotsTest = kFALSE; // Save plots in Psfile
3584
3585 // TPostScript* DummyPs = new TPostScript("dummy.ps");
3586
3587 TString DummyPs = ("Dummy");
3588
3589 MHFindSignificanceONOFF findsigTest;
3590 findsigTest.SetRebin(kTRUE);
3591 findsigTest.SetReduceDegree(kFALSE);
3592 findsigTest.SetUseFittedQuantities(fUseFittedQuantities);
3593
3594
3595 if (findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista,
3596 fNormFactorTest,
3597 alphamin, alphamax,
3598 fDegreeON, fDegreeOFF,
3599 alphasig, drawpolyTest,
3600 fitgaussTest, printTest,
3601 saveplotsTest, DummyPs))
3602 {
3603
3604 // Update values of Nex and SigmaLiMa for Test sammple
3605
3606 fSigmaLiMaTest = double(findsigTest.GetSignificance());
3607
3608 if(fUseFittedQuantities)
3609 {
3610 fNexTest = double(findsigTest.GetNexONOFFFitted());
3611 }
3612 else
3613 {
3614 fNexTest = double(findsigTest.GetNexONOFF());
3615 }
3616
3617
3618 }
3619 else
3620 {
3621
3622 *fLog << "Normalization Factor tuning for Train sample is not possible"
3623 << endl;
3624 fSigmaLiMaTest = 0.0;
3625 fNexTest = 0.0;
3626 }
3627
3628
3629// DummyPs -> Close();
3630 //DummyPs = NULL;
3631 // delete DummyPs;
3632
3633
3634
3635
3636
3637 Int_t EventsInTestMatrixOFF = fMatrixTestOFF->GetM().GetNrows();
3638 Double_t Ngammas = fNexTest/fGammaEfficiency;
3639 Double_t GammaFraction = Ngammas/EventsInTestMatrixOFF;
3640
3641
3642 *fLog << "MFindSupercutsONOFF::TestParamsOnTestSample; "
3643 << "fNormFactorTest is corrected with fraction of gammas in ON "
3644 << "Data Test sample and the number of OFF events before cuts."
3645 << endl
3646 << "EventsInTestMatrixOFF = " << EventsInTestMatrixOFF
3647 << " Excess events in Test ON sample = " << fNexTest
3648 << " fGammaEfficiency = " << fGammaEfficiency << endl
3649 << "fNormFactorTest Value before correction is "
3650 << fNormFactorTest << endl;
3651
3652 fNormFactorTest = fNormFactorTest - GammaFraction;
3653
3654 *fLog << "fNormFactorTest Value after correction is "
3655 << fNormFactorTest << endl;
3656
3657
3658 }
3659
3660
3661
3662
3663 if (fNormFactorFromAlphaBkg)
3664 {
3665 // Normalization factor computed using alpha bkg region
3666 Double_t NewNormFactor =
3667 ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista,
3668 fAlphaBkgMin, fAlphaBkgMax);
3669
3670 *fLog << "Normalization factor computed from alpha plot (after cuts) " << endl
3671 << "using counted number of ON and OFF events in alpha region " << endl
3672 << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl
3673 << "Normalization factor = " << NewNormFactor << endl;
3674
3675 *fLog << "Normalization factor used is the one computed in bkg region; " << endl
3676 << "i.e. " << NewNormFactor << " instead of " << fNormFactorTest << endl;
3677
3678 fNormFactorTest = NewNormFactor;
3679
3680
3681 }
3682
3683
3684
3685 // Significance is found using MHFindSignificanceONOFF::FindSigmaONOFF
3686
3687
3688 MHFindSignificanceONOFF findsig;
3689 findsig.SetRebin(kTRUE);
3690 findsig.SetReduceDegree(kFALSE);
3691 findsig.SetUseFittedQuantities(fUseFittedQuantities);
3692
3693 TString psfilename; // Name of psfile where Alpha plot will be stored
3694 if (!fPsFilenameString.IsNull())
3695 {
3696 psfilename += (fPsFilenameString);
3697 psfilename += ("TESTSample");
3698 }
3699 else
3700 {
3701 psfilename += ("TESTSample");
3702 }
3703
3704
3705 findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista,
3706 fNormFactorTest,
3707 alphamin, alphamax,
3708 fDegreeON, fDegreeOFF,
3709 alphasig, drawpoly, fitgauss, print,
3710 saveplots, psfilename);
3711
3712 const Double_t significance = findsig.GetSignificance();
3713 const Double_t alphasi = findsig.GetAlphasi();
3714
3715
3716 // Update values of Nex and SigmaLiMa for Test sammple
3717
3718 fSigmaLiMaTest = double(findsig.GetSignificance());
3719
3720 if(fUseFittedQuantities)
3721 {
3722 fNexTest = double(findsig.GetNexONOFFFitted());
3723 }
3724 else
3725 {
3726 fNexTest = double(findsig.GetNexONOFF());
3727 }
3728
3729
3730 // Alpha distributions and containers with teh supercuts applied
3731 // are written into root file
3732
3733 TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
3734 "Alpha Distributions for several Theta bins");
3735
3736 alphaHist2.Write();
3737 alphaOFFHista.Write();
3738
3739 (supapptestON.GetTreePointer())->Write();
3740 (supapptestOFF.GetTreePointer())->Write();
3741
3742 *fLog << "TTree objects containing supercuts applied to TEST sample are "
3743 << "written into root file " << fAlphaDistributionsRootFilename << endl;
3744
3745 (supapptestON.GetTreePointer())->Print();
3746 (supapptestOFF.GetTreePointer())->Print();
3747
3748 //rootfile.Close();
3749
3750
3751 // Boolean variable that controls in class MSupercutsCalcONOFF
3752 // wether the supercuts are stored or not
3753 // is set to the default value, kFALSE.
3754
3755 fCalcHadTest -> SetStoreAppliedSupercuts(kFALSE);
3756 fCalcHadTestOFF -> SetStoreAppliedSupercuts(kFALSE);
3757
3758
3759
3760 *fLog << "Significance of gamma signal after supercuts in Test sample : "
3761 << significance << " (for |alpha| < " << alphasi << " degrees)"
3762 << endl;
3763 //-------------------------------------------
3764
3765
3766 *fLog << "Test of supercuts on TEST sample finished" << endl;
3767 *fLog << "======================================================" << endl;
3768
3769 return kTRUE;
3770}
3771
3772
3773
3774
3775
3776// Function that applies the optimized supercuts on the TEST sample.
3777
3778// Supercuts applied (LengthUp, LengthLow...)
3779// to all the individual events, as well as
3780// the features of the shower images (Length,Width...) are
3781// stored in 2 TTree objects (for ON and OFF data) in the
3782// root file specified by variable "fAlphaDistributionsRootFilename"
3783
3784
3785Bool_t MFindSupercutsONOFF::TestParamsOnTrainSample()
3786{
3787 if (fMatrixTrain->GetM().GetNrows() <= 0)
3788 {
3789 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix ON has no entries"
3790 << endl;
3791 return kFALSE;
3792 }
3793
3794
3795 if (fMatrixTrainOFF->GetM().GetNrows() <= 0)
3796 {
3797 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; train matrix OFF has no entries"
3798 << endl;
3799 return kFALSE;
3800 }
3801
3802 // ------------- TEST the supercuts on TRAIN sample --------------------------
3803 //
3804 *fLog << "Producing output of the SUPERCUTS ON-OFF optimization " << endl;
3805
3806 // -----------------------------------------------------------------
3807 // read optimum parameter values from file filenameParam
3808 // into array 'supercutsPar'
3809
3810 TFile inparam(fFilenameParam);
3811 MSupercuts scin;
3812 scin.Read("MSupercuts");
3813 inparam.Close();
3814
3815 *fLog << "Optimum parameter values for supercuts were read from file '";
3816 *fLog << fFilenameParam << "' :" << endl;
3817
3818 const TArrayD &supercutsPar = scin.GetParameters();
3819 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
3820 *fLog << supercutsPar[i] << ", ";
3821 *fLog << endl;
3822 //---------------------------
3823
3824
3825 // -----------------------------------------------------------------
3826 if (fHadronnessName.IsNull())
3827 {
3828 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; hadronness name for ON data is not defined... aborting";
3829 *fLog << endl;
3830 return kFALSE;
3831 }
3832
3833 if (fHadronnessNameOFF.IsNull())
3834 {
3835 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; hadronness name for OFF data is not defined... aborting";
3836 *fLog << endl;
3837 return kFALSE;
3838 }
3839
3840
3841 if(fTrainONSupercutsAppliedName.IsNull())
3842 {
3843 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
3844 << " MTSupercutsApplied Tree name for ON Train data is not defined... aborting"
3845 << endl;
3846 return kFALSE;
3847 }
3848
3849 if(fTrainOFFSupercutsAppliedName.IsNull())
3850 {
3851 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
3852 << " MTSupercutsApplied Tree name for OFF Train data is not defined... aborting"
3853 << endl;
3854 return kFALSE;
3855 }
3856
3857
3858 MSupercuts supercuts;
3859 supercuts.SetParameters(supercutsPar);
3860
3861 // Objects containing the trees used to store supercuts applied
3862 // TTree Branches are also created
3863 MTSupercutsApplied supapptrainON(fTrainONSupercutsAppliedName,"");
3864 supapptrainON.CreateTreeBranches();
3865
3866 MTSupercutsApplied supapptrainOFF(fTrainOFFSupercutsAppliedName,"");
3867 supapptrainOFF.CreateTreeBranches();
3868
3869 fCalcHadTrainOFF -> SetHadronnessName(fHadronnessNameOFF);
3870 fCalcHadTrain -> SetHadronnessName(fHadronnessName);
3871
3872
3873 // Settings related to the storage of teh supercuts applied
3874 fCalcHadTrain -> SetSupercutsAppliedName(fTrainONSupercutsAppliedName);
3875 fCalcHadTrain -> SetStoreAppliedSupercuts(kTRUE);
3876
3877 fCalcHadTrainOFF -> SetSupercutsAppliedName(fTrainOFFSupercutsAppliedName);
3878 fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kTRUE);
3879
3880
3881 // Set boolean variable that controls wether cuts are
3882 // dynamic or static
3883
3884 fCalcHadTrain -> SetVariableUseStaticCuts(fUseStaticCuts);
3885 fCalcHadTrainOFF -> SetVariableUseStaticCuts(fUseStaticCuts);
3886
3887
3888 if (!fUseStaticCuts)
3889 {
3890 // Set boolean variable that controls wether the theta variable
3891 // is used or not in the computation of the dynamical cuts
3892
3893 fCalcHadTrain -> SetVariableNotUseTheta(fNotUseTheta);
3894 fCalcHadTrainOFF -> SetVariableNotUseTheta(fNotUseTheta);
3895 }
3896
3897
3898 // apply the supercuts to OFF data
3899
3900 MF scfilterOFF(fHadronnessNameOFF+".fHadronness>0.5");
3901 MContinue applysupercutsOFF(&scfilterOFF);
3902
3903 // apply the supercuts to ON data
3904 MF scfilter(fHadronnessName+".fHadronness>0.5");
3905 MContinue applysupercuts(&scfilter);
3906
3907
3908
3909
3910 // ****************************************************
3911 // Filling OFF alpha distribution
3912
3913 MParList parlistOFF;
3914 MTaskList tasklistOFF;
3915
3916 MMatrixLoop loopOFF(fMatrixTrainOFF);
3917
3918 // plot alpha before applying the supercuts
3919 const TString mh3NameOFFB = "AlphaOFFBefore";
3920 MBinning binsalphaOFFbef("Binning"+mh3NameOFFB);
3921 //binsalphaOFFbef.SetEdges(54, -12.0, 96.0);
3922 binsalphaOFFbef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3923
3924
3925 // |alpha| is assumed to be in column 7 of the matrix
3926 MH3 alphaOFFbefore("MatrixTrainOFF[7]");
3927 alphaOFFbefore.SetName(mh3NameOFFB);
3928
3929
3930 MFillH fillalphaOFFbefore(&alphaOFFbefore);
3931 fillalphaOFFbefore.SetName("FillAlphaOFFBefore");
3932
3933
3934
3935 // fill OFF alpha after applying the supercuts
3936 const TString mh3NameOFFA = "AlphaOFFAfter";
3937 MBinning binsalphaOFFaft("Binning"+mh3NameOFFA);
3938 //binsalphaOFFaft.SetEdges(54, -12.0, 96.0);
3939 binsalphaOFFaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
3940
3941 MH3 alphaOFFafter("MatrixTrainOFF[7]");
3942 alphaOFFafter.SetName(mh3NameOFFA);
3943
3944 //TH1 &alphaOFFhista = alphaOFFafter.GetHist();
3945 //alphaOFFhista.SetName("alphaOFFAfter-OptimizationParamOutput");
3946
3947 MFillH fillalphaOFFafter(&alphaOFFafter);
3948 fillalphaOFFafter.SetName("FillAlphaOFFAfter");
3949
3950 //******************************
3951 // entries in MParList
3952
3953 parlistOFF.AddToList(&tasklistOFF);
3954
3955 parlistOFF.AddToList(&supercuts);
3956
3957 parlistOFF.AddToList(&supapptrainOFF);
3958
3959 parlistOFF.AddToList(fCam);
3960 parlistOFF.AddToList(fMatrixTrainOFF);
3961
3962 parlistOFF.AddToList(&binsalphaOFFbef);
3963 parlistOFF.AddToList(&alphaOFFbefore);
3964
3965 parlistOFF.AddToList(&binsalphaOFFaft);
3966 parlistOFF.AddToList(&alphaOFFafter);
3967
3968
3969 //******************************
3970 // Entries in MtaskList
3971 tasklistOFF.AddToList(&loopOFF);
3972 tasklistOFF.AddToList(&fillalphaOFFbefore);
3973
3974 tasklistOFF.AddToList(fCalcHadTrainOFF);
3975 tasklistOFF.AddToList(&applysupercutsOFF);
3976
3977 tasklistOFF.AddToList(&fillalphaOFFafter);
3978
3979 //******************************
3980
3981 MProgressBar barOFF;
3982 MEvtLoop evtloopOFF;
3983 evtloopOFF.SetParList(&parlistOFF);
3984 evtloopOFF.SetName("EvtLoopOptimizationParamOutputOFF");
3985 evtloopOFF.SetProgressBar(&barOFF);
3986 Int_t maxeventsOFF = -1;
3987 if (!evtloopOFF.Eventloop(maxeventsOFF))
3988 return kFALSE;
3989
3990 tasklistOFF.PrintStatistics(0, kTRUE);
3991
3992
3993 //-------------------------------------------
3994 // draw the alpha plots
3995
3996 MH3* alphaOFFBefore = (MH3*)parlistOFF.FindObject(mh3NameOFFB, "MH3");
3997 TH1 &alphaOFFHistb = alphaOFFBefore->GetHist();
3998 alphaOFFHistb.SetXTitle("|\\alpha| [\\circ]");
3999
4000 TString alphaOFFHistbName ("TrainAlphaOFFBeforeCuts");
4001 if (fThetaRangeString.IsNull())
4002 {
4003 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
4004 << " fThetaRangeString is NOT defined..." << endl;
4005 alphaOFFHistbName += ("ThetaRangeStringUndefined");
4006 }
4007 else
4008 {
4009 alphaOFFHistbName += (fThetaRangeString);
4010 }
4011
4012 alphaOFFHistb.SetName(alphaOFFHistbName);
4013
4014 MH3* alphaOFFAfter = (MH3*)parlistOFF.FindObject(mh3NameOFFA, "MH3");
4015 TH1 &alphaOFFHista = alphaOFFAfter->GetHist();
4016 TString alphaOFFHistaName ("TrainAlphaOFFAfterCuts");
4017 if (fThetaRangeString.IsNull())
4018 {
4019 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
4020 << " fThetaRangeString is NOT defined..." << endl;
4021 alphaOFFHistaName += ("ThetaRangeStringUndefined");
4022 }
4023 else
4024 {
4025 alphaOFFHistaName += (fThetaRangeString);
4026 }
4027
4028
4029 alphaOFFHista.SetXTitle("|\\alpha| [\\circ]");
4030 alphaOFFHista.SetName(alphaOFFHistaName);
4031
4032
4033 // *****************************************************
4034
4035
4036
4037 //-------------------------------------------
4038
4039
4040 MParList parlist2;
4041 MTaskList tasklist2;
4042
4043 MMatrixLoop loop(fMatrixTrain);
4044
4045 // plot alpha before applying the supercuts
4046 const TString mh3NameB = "AlphaBefore";
4047 MBinning binsalphabef("Binning"+mh3NameB);
4048 //binsalphabef.SetEdges(54, -12.0, 96.0);
4049 binsalphabef.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
4050
4051 // |alpha| is assumed to be in column 7 of the matrix
4052 MH3 alphabefore("MatrixTrain[7]");
4053 alphabefore.SetName(mh3NameB);
4054
4055 MFillH fillalphabefore(&alphabefore);
4056 fillalphabefore.SetName("FillAlphaBefore");
4057
4058
4059 // plot alpha after applying the supercuts
4060 const TString mh3NameA = "AlphaAfter";
4061 MBinning binsalphaaft("Binning"+mh3NameA);
4062 //binsalphaaft.SetEdges(54, -12.0, 96.0);
4063 binsalphaaft.SetEdges(fNAlphaBins, fAlphaBinLow, fAlphaBinUp);
4064
4065 MH3 alphaafter("MatrixTrain[7]");
4066 alphaafter.SetName(mh3NameA);
4067
4068
4069 MFillH fillalphaafter(&alphaafter);
4070 fillalphaafter.SetName("FillAlphaAfter");
4071
4072 //******************************
4073 // entries in MParList
4074
4075 parlist2.AddToList(&tasklist2);
4076
4077 parlist2.AddToList(&supercuts);
4078
4079 parlist2.AddToList(&supapptrainON);
4080
4081 parlist2.AddToList(fCam);
4082 parlist2.AddToList(fMatrixTrain);
4083
4084 parlist2.AddToList(&binsalphabef);
4085 parlist2.AddToList(&alphabefore);
4086
4087 parlist2.AddToList(&binsalphaaft);
4088 parlist2.AddToList(&alphaafter);
4089
4090 //******************************
4091 // entries in MTaskList
4092
4093 tasklist2.AddToList(&loop);
4094 tasklist2.AddToList(&fillalphabefore);
4095
4096 tasklist2.AddToList(fCalcHadTrain);
4097 tasklist2.AddToList(&applysupercuts);
4098
4099 tasklist2.AddToList(&fillalphaafter);
4100
4101 //******************************
4102
4103 MProgressBar bar2;
4104 MEvtLoop evtloop2;
4105 evtloop2.SetParList(&parlist2);
4106 evtloop2.SetName("EvtLoopOptimizationParamOutput");
4107 evtloop2.SetProgressBar(&bar2);
4108 Int_t maxevents2 = -1;
4109 if (!evtloop2.Eventloop(maxevents2))
4110 return kFALSE;
4111
4112 tasklist2.PrintStatistics(0, kTRUE);
4113
4114
4115 //-------------------------------------------
4116 // draw the alpha plots
4117
4118 MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
4119 TH1 &alphaHist1 = alphaBefore->GetHist();
4120 TString alphaHist1Name ("TrainAlphaBeforeCuts");
4121 if (fThetaRangeString.IsNull())
4122 {
4123 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
4124 << " fThetaRangeString is NOT defined..." << endl;
4125 alphaHist1Name += ("ThetaRangeStringUndefined");
4126 }
4127 else
4128 {
4129 alphaHist1Name += (fThetaRangeString);
4130 }
4131
4132
4133 alphaHist1.SetXTitle("|\\alpha| [\\circ]");
4134 alphaHist1.SetName(alphaHist1Name);
4135
4136 MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
4137 TH1 &alphaHist2 = alphaAfter->GetHist();
4138 TString alphaHist2Name ("TrainAlphaAfterCuts");
4139 if (fThetaRangeString.IsNull())
4140 {
4141 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
4142 << " fThetaRangeString is NOT defined..." << endl;
4143 alphaHist2Name += ("ThetaRangeStringUndefined");
4144 }
4145 else
4146 {
4147 alphaHist2Name += (fThetaRangeString);
4148 }
4149
4150 alphaHist2.SetXTitle("|\\alpha| [\\circ]");
4151 alphaHist2.SetName(alphaHist2Name);
4152
4153
4154 // fPsFilename -> NewPage();
4155
4156
4157 // Canvas is deleted after saving the plot in ps file
4158 // This is to allow for GRID analysis
4159
4160 TCanvas *c = new TCanvas("TrainAlphaBeforeAfterSC",
4161 "TrainAlphaBeforeAfterSC",
4162 600, 600);
4163
4164 c->Divide(2,2);
4165
4166
4167
4168 c->cd(1);
4169 alphaOFFHistb.DrawCopy();
4170
4171 c->cd(2);
4172 alphaOFFHista.DrawCopy();
4173
4174
4175 //c->Modified();
4176 //c-> Update();
4177
4178
4179 // fPsFilename -> NewPage();
4180
4181 c->cd(3);
4182 alphaHist1.DrawCopy();
4183
4184 c->cd(4);
4185 alphaHist2.DrawCopy();
4186
4187
4188 c -> Modified();
4189 c -> Update();
4190
4191 // ********************************************
4192 // TMP solutions while the TPostScript thing is not working
4193 // PsFileName for storing these plots is derived
4194 // from PsFilenameString.
4195
4196 cout << "Alpha distributions will be saved in PostScript file " ;
4197
4198 if (!fPsFilenameString.IsNull())
4199 {
4200 TString filename = (fPsFilenameString);
4201 filename += ("AlphaTRAINSampleONOFFBeforeAfterCuts.ps");
4202 cout << filename << endl;
4203 c -> SaveAs(filename);
4204 }
4205
4206 // END OF TEMPORAL SOLUTION
4207 // ********************************************
4208
4209
4210
4211
4212 // fPsFilename -> NewPage();
4213
4214
4215 // Canvas is deleted after saving the plot in ps file
4216 // This is to allow for GRID analysis
4217
4218 delete c;
4219
4220
4221
4222
4223 //-------------------------------------------
4224 // fit alpha distribution to get the number of excess events and
4225 // calculate significance of gamma signal in the alpha plot
4226
4227 const Double_t alphasig = fAlphaSig;
4228 const Double_t alphamin = fAlphaBkgMin;
4229 // alpha min for bkg region in ON data
4230
4231 const Double_t alphamax = fAlphaBkgMax;
4232// const Int_t degree = 2;
4233// const Int_t degreeOFF = 2;
4234 const Bool_t drawpoly = kTRUE;
4235 const Bool_t fitgauss = kTRUE;
4236 const Bool_t print = kTRUE;
4237
4238 Bool_t saveplots = kTRUE; // Save plots in Psfile
4239
4240
4241
4242
4243
4244
4245 if (!ComputeNormFactorTrain())
4246 {
4247 *fLog << "Normalization factor for train sample (ON-OFF) could not be computed. Aborting..." << endl;
4248 return kFALSE;
4249 }
4250
4251
4252 // Normalization factor factor is corrected using the estimated
4253 // number of excess events and the total number of OFF events
4254 // fNormFactor = fNormFactor - Ngammas/EventsInMatrixOFF
4255
4256
4257 if (fTuneNormFactor)
4258 {
4259
4260
4261 // Run MHFindSignificanceONOFF::FindSigmaONOFF
4262 // to get estimation of number of excess gamma events.
4263 // This number will be used in the correction of the
4264 // normalization factor
4265
4266 // No plots will be produced this time...
4267
4268 const Bool_t drawpolyTest = kFALSE;
4269 const Bool_t fitgaussTest = kFALSE;
4270 const Bool_t printTest = kFALSE;
4271
4272 const Bool_t saveplotsTest = kFALSE; // Save plots in Psfile
4273
4274 //TPostScript* DummyPs = new TPostScript("dummy.ps");
4275
4276 TString DummyPs = ("Dummy");
4277
4278 MHFindSignificanceONOFF findsigTest;
4279 findsigTest.SetRebin(kTRUE);
4280 findsigTest.SetReduceDegree(kFALSE);
4281 findsigTest.SetUseFittedQuantities(fUseFittedQuantities);
4282
4283 if(findsigTest.FindSigmaONOFF(&alphaHist2, &alphaOFFHista,
4284 fNormFactorTrain,
4285 alphamin, alphamax,
4286 fDegreeON, fDegreeOFF,
4287 alphasig, drawpolyTest,
4288 fitgaussTest, printTest,
4289 saveplotsTest, DummyPs))
4290
4291 {
4292
4293 // Update values of Nex and SigmaLiMa for Test sammple
4294
4295 fSigmaLiMaTrain = double(findsigTest.GetSignificance());
4296
4297 if(fUseFittedQuantities)
4298 {
4299 fNexTrain = double(findsigTest.GetNexONOFFFitted());
4300 }
4301 else
4302 {
4303 fNexTrain = double(findsigTest.GetNexONOFF());
4304 }
4305
4306 }
4307
4308 else
4309 {
4310 *fLog << "Normalization Factor tuning for Train sample is not possible"
4311 << endl;
4312
4313 fSigmaLiMaTrain = 0.0;
4314 fNexTrain = 0.0;
4315
4316 }
4317
4318 //DummyPs -> Close();
4319 //DummyPs = NULL;
4320 //delete DummyPs;
4321
4322
4323
4324
4325 Int_t EventsInTrainMatrixOFF = fMatrixTrainOFF->GetM().GetNrows();
4326 Double_t Ngammas = fNexTrain/fGammaEfficiency;
4327 Double_t GammaFraction = Ngammas/EventsInTrainMatrixOFF;
4328
4329
4330 *fLog << "MFindSupercutsONOFF::TestParamsOnTrainSample; "
4331 << "fNormFactorTrain is corrected with fraction of gammas in ON "
4332 << "Data Train sample and the number of OFF events before cuts."
4333 << endl
4334 << "EventsInTrainMatrixOFF = " << EventsInTrainMatrixOFF
4335 << " Excess events in Train ON sample = " << fNexTrain
4336 << " fGammaEfficiency = " << fGammaEfficiency << endl
4337 << "fNormFactorTrain Value before correction is "
4338 << fNormFactorTrain << endl;
4339
4340 fNormFactorTrain = fNormFactorTrain - GammaFraction;
4341
4342 *fLog << "fNormFactorTrain Value after correction is "
4343 << fNormFactorTrain << endl;
4344
4345
4346
4347 }
4348
4349 if (fNormFactorFromAlphaBkg)
4350 {
4351 Double_t NewNormFactor =
4352 ComputeNormFactorFromAlphaBkg(&alphaHist2, &alphaOFFHista,
4353 fAlphaBkgMin, fAlphaBkgMax);
4354
4355 *fLog << "Normalization factor computed from alpha plot (after cuts) " << endl
4356 << "using counted number of ON and OFF events in alpha region " << endl
4357 << "defined by range " << fAlphaBkgMin << "-" << fAlphaBkgMax << endl
4358 << "Normalization factor = " << NewNormFactor << endl;
4359
4360 *fLog << "Normalization factor used is the one computed in bkg region; " << endl
4361 << "i.e. " << NewNormFactor << " instead of " << fNormFactorTrain << endl;
4362
4363 fNormFactorTrain = NewNormFactor;
4364
4365
4366 }
4367
4368
4369 MHFindSignificanceONOFF findsig;
4370 findsig.SetRebin(kTRUE);
4371 findsig.SetReduceDegree(kFALSE);
4372 findsig.SetUseFittedQuantities(fUseFittedQuantities);
4373
4374 TString psfilename; // Name of psfile where Alpha plot will be stored
4375 if (!fPsFilenameString.IsNull())
4376 {
4377 psfilename += (fPsFilenameString);
4378 psfilename += ("TRAINSample");
4379 }
4380 else
4381 {
4382 psfilename += ("TRAINSample");
4383 }
4384
4385
4386 findsig.FindSigmaONOFF(&alphaHist2, &alphaOFFHista,
4387 fNormFactorTrain,
4388 alphamin, alphamax,
4389 fDegreeON, fDegreeOFF,
4390 alphasig, drawpoly, fitgauss, print,
4391 saveplots, psfilename);
4392
4393 const Double_t significance = findsig.GetSignificance();
4394 const Double_t alphasi = findsig.GetAlphasi();
4395
4396 // Update Train values Nex, SigmaLiMa
4397
4398 fSigmaLiMaTrain = double(findsig.GetSignificance());
4399
4400 if(fUseFittedQuantities)
4401 {
4402 fNexTrain = double(findsig.GetNexONOFFFitted());
4403 }
4404 else
4405 {
4406 fNexTrain = double(findsig.GetNexONOFF());
4407 }
4408
4409
4410 // Alpha distributions and containers with the supercuts applied
4411 // are written into root file
4412
4413 TFile rootfile(fAlphaDistributionsRootFilename, "UPDATE",
4414 "Alpha Distributions for several Theta bins");
4415
4416 alphaHist2.Write();
4417 alphaOFFHista.Write();
4418
4419 (supapptrainON.GetTreePointer())->Write();
4420 (supapptrainOFF.GetTreePointer())->Write();
4421
4422 *fLog << "TTree objects containing supercuts applied to TEST sample are "
4423 << "written into root file " << fAlphaDistributionsRootFilename << endl;
4424
4425 (supapptrainON.GetTreePointer())->Print();
4426 (supapptrainOFF.GetTreePointer())->Print();
4427
4428
4429 //rootfile.Close();
4430
4431
4432 // Boolean variable that controls in class MSupercutsCalcONOFF
4433 // wether the supercuts are stored or not
4434 // is set to the default value, kFALSE.
4435
4436 fCalcHadTrain -> SetStoreAppliedSupercuts(kFALSE);
4437 fCalcHadTrainOFF -> SetStoreAppliedSupercuts(kFALSE);
4438
4439
4440
4441
4442
4443 *fLog << "Significance of gamma signal after supercuts in train sample : "
4444 << significance << " (for |alpha| < " << alphasi << " degrees)"
4445 << endl;
4446 //-------------------------------------------
4447
4448
4449 *fLog << "Test of supercuts on TRAIN sample finished" << endl;
4450 *fLog << "======================================================" << endl;
4451
4452 return kTRUE;
4453}
4454
4455
4456
4457
4458
4459
4460
4461
4462
4463
4464
4465
4466
4467
4468
4469
4470
4471
4472
4473
4474
Note: See TracBrowser for help on using the repository browser.