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

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