source: trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C@ 3097

Last change on this file since 3097 was 2798, checked in by wittek, 22 years ago
*** empty log message ***
File size: 100.7 KB
Line 
1
2
3//#include "MagicEgyEst.C"
4
5
6void InitBinnings(MParList *plist)
7{
8 gLog << "InitBinnings" << endl;
9
10 //--------------------------------------------
11 MBinning *binse = new MBinning("BinningE");
12 //binse->SetEdgesLog(30, 1.0e2, 1.0e5);
13
14 //This is Daniel's binning in energy:
15 binse->SetEdgesLog(14, 296.296, 86497.6);
16 plist->AddToList(binse);
17
18 //--------------------------------------------
19
20 MBinning *binssize = new MBinning("BinningSize");
21 binssize->SetEdgesLog(50, 10, 1.0e5);
22 plist->AddToList(binssize);
23
24 MBinning *binsdistc = new MBinning("BinningDist");
25 binsdistc->SetEdges(50, 0, 1.4);
26 plist->AddToList(binsdistc);
27
28 MBinning *binswidth = new MBinning("BinningWidth");
29 binswidth->SetEdges(50, 0, 1.0);
30 plist->AddToList(binswidth);
31
32 MBinning *binslength = new MBinning("BinningLength");
33 binslength->SetEdges(50, 0, 1.0);
34 plist->AddToList(binslength);
35
36 MBinning *binsalpha = new MBinning("BinningAlpha");
37 binsalpha->SetEdges(100, -100, 100);
38 plist->AddToList(binsalpha);
39
40 MBinning *binsasym = new MBinning("BinningAsym");
41 binsasym->SetEdges(50, -1.5, 1.5);
42 plist->AddToList(binsasym);
43
44 MBinning *binsm3l = new MBinning("BinningM3Long");
45 binsm3l->SetEdges(50, -1.5, 1.5);
46 plist->AddToList(binsm3l);
47
48 MBinning *binsm3t = new MBinning("BinningM3Trans");
49 binsm3t->SetEdges(50, -1.5, 1.5);
50 plist->AddToList(binsm3t);
51
52
53 //.....
54 MBinning *binsb = new MBinning("BinningSigmabar");
55 binsb->SetEdges( 100, 0.0, 5.0);
56 plist->AddToList(binsb);
57
58 MBinning *binth = new MBinning("BinningTheta");
59 // this is Daniel's binning in theta
60 //Double_t yedge[8] =
61 // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99};
62 // this is our binning
63 Double_t yedge[9] =
64 {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
65 TArrayD yed;
66 yed.Set(9,yedge);
67 binth->SetEdges(yed);
68 plist->AddToList(binth);
69
70 MBinning *bincosth = new MBinning("BinningCosTheta");
71 Double_t zedge[9];
72 for (Int_t i=0; i<9; i++)
73 {
74 zedge[8-i] = cos(yedge[i]/kRad2Deg);
75 }
76 TArrayD zed;
77 zed.Set(9,zedge);
78 bincosth->SetEdges(zed);
79 plist->AddToList(bincosth);
80
81 MBinning *binsdiff = new MBinning("BinningDiffsigma2");
82 binsdiff->SetEdges(100, -5.0, 20.0);
83 plist->AddToList(binsdiff);
84
85 // robert ----------------------------------------------
86 MBinning *binsalphaf = new MBinning("BinningAlphaFlux");
87 binsalphaf->SetEdges(100, -100, 100);
88 plist->AddToList(binsalphaf);
89
90 MBinning *binsdifftime = new MBinning("BinningTimeDiff");
91 binsdifftime->SetEdges(50, 0, 10);
92 plist->AddToList(binsdifftime);
93
94 MBinning *binstime = new MBinning("BinningTime");
95 binstime->SetEdges(50, 44500, 61000);
96 plist->AddToList(binstime);
97}
98
99
100void DeleteBinnings(MParList *plist)
101{
102 gLog << "DeleteBinnings" << endl;
103
104 TObject *bin;
105
106 //--------------------------------------------
107 bin = plist->FindObject("BinningE");
108 if (bin) delete bin;
109
110 //--------------------------------------------
111
112 bin = plist->FindObject("BinningSize");
113 if (bin) delete bin;
114
115 bin = plist->FindObject("BinningDist");
116 if (bin) delete bin;
117
118 bin = plist->FindObject("BinningWidth");
119 if (bin) delete bin;
120
121 bin = plist->FindObject("BinningLength");
122 if (bin) delete bin;
123
124 bin = plist->FindObject("BinningAlpha");
125 if (bin) delete bin;
126
127 bin = plist->FindObject("BinningAsym");
128 if (bin) delete bin;
129
130 bin = plist->FindObject("BinningM3Long");
131 if (bin) delete bin;
132
133 bin = plist->FindObject("BinningM3Trans");
134 if (bin) delete bin;
135
136 //.....
137 bin = plist->FindObject("BinningSigmabar");
138 if (bin) delete bin;
139
140 bin = plist->FindObject("BinningTheta");
141 if (bin) delete bin;
142
143 bin = plist->FindObject("BinningCosTheta");
144 if (bin) delete bin;
145
146 bin = plist->FindObject("BinningDiffsigma2");
147 if (bin) delete bin;
148
149
150 // robert ----------------------------------------------
151 bin = plist->FindObject("BinningAlphaFlux");
152 if (bin) delete bin;
153
154 bin = plist->FindObject("BinningTimeDiff");
155 if (bin) delete bin;
156
157 bin = plist->FindObject("BinningTime");
158 if (bin) delete bin;
159}
160
161
162
163//************************************************************************
164void ONOFFAnalysis()
165{
166 gLog.SetNoColors();
167
168 if (gRandom)
169 delete gRandom;
170 gRandom = new TRandom3(0);
171
172 //-----------------------------------------------
173 //const char *offfile = "~magican/ct1test/wittek/offdata.preproc";
174 const char *offfile = "/data/MAGIC/rootdata/2003_11_29/20031128_032*_D_OffCrab1_E.root";
175
176 //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc";
177 // const char *onfile = "~magican/ct1test/wittek/mkn421_00-01";
178 const char *onfile = "/data/MAGIC/rootdata/2003_11_29/20031128_031*_D_Crab-Nebula_E.root";
179
180 const char *mcfile = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root";
181 //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_4/gh/0/0/G_M0_00_0_550*.root";
182 //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_5/gh/0/0/G_M0_00_0_550*.root";
183 //-----------------------------------------------
184
185 // path for input for Mars
186 TString inPath = "/data/MAGIC/scratch/wittek/";
187
188 // path for output from Mars
189 TString outPath = "/data/MAGIC/scratch/wittek/";
190
191 //-----------------------------------------------
192
193 //TEnv env("macros/CT1env.rc");
194 //Bool_t printEnv = kFALSE;
195
196 //************************************************************************
197
198 // Job A :
199 // - produce MHSigmaTheta plots for ON, OFF and MC data
200 // - write out (or read in) these MHSigmaTheta plots
201 // - read ON (or OFF or MC) data
202 // - pad the events;
203 // - write root file for ON (or OFF or MC) data (ON1.root, ...);
204
205 Bool_t JobA = kTRUE;
206 Bool_t GPad = kFALSE; // generate padding histograms?
207 Bool_t WPad = kFALSE; // write out padding histograms ?
208 Bool_t RPad = kFALSE; // read in padding histograms ?
209 Bool_t Wout = kTRUE; // write out root file ON1.root
210 // (or OFF1.root or MC1.root)?
211
212
213 // Job B_RF_UP : read ON1.root (OFF1.root or MC1.root) file
214 // - if CTrainRF = TRUE : create matrices of training events
215 // and root files of training and test events
216 // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
217 // - if RTree = TRUE : read in trees, otherwise create trees
218 // - calculate hadroness for method of RANDOM FOREST
219 // - update the input files with the hadronesses (ON2.root, OFF2.root
220 // or MC2.root)
221
222 Bool_t JobB_RF_UP = kFALSE;
223 Bool_t CTrainRF = kFALSE; // create matrices of training events
224 // and root files of training and test events
225 Bool_t RTrainRF = kFALSE; // read in matrices of training events
226 Bool_t RTree = kFALSE; // read in trees (otherwise grow them)
227 Bool_t WRF = kFALSE; // update input root file ?
228
229
230 // Job B_SC_UP : read ON2.root (or MC2.root) file
231 // - depending on WParSC : create (or read in) supercuts parameter values
232 // - calculate hadroness for the SUPERCUTS
233 // - update the input files with the hadroness (==>ON3.root or MC3.root)
234
235 Bool_t JobB_SC_UP = kTRUE;
236 Bool_t CMatrix = kTRUE; // create training and test matrices
237 Bool_t RMatrix = kFALSE; // read training and test matrices from file
238 Bool_t WOptimize = kFALSE; // do optimization using the training sample
239 // and write supercuts parameter values
240 // onto the file parSCfile
241 Bool_t RTest = kFALSE; // test the supercuts using the test matrix
242 Bool_t WSC = kFALSE; // update input root file ?
243
244
245 // Job C:
246 // - read ON3.root and MC3.root files
247 // which should have been updated to contain the hadronnesses
248 // for the method of
249 // RF
250 // SUPERCUTS and
251 // - produce Neyman-Pearson plots
252
253 Bool_t JobC = kFALSE;
254
255
256 // Job D :
257 // - select g/h separation method XX
258 // - read ON3 (or MC3) root file
259 // - apply cuts in hadronness
260 // - make plots
261
262 Bool_t JobD = kFALSE;
263
264
265
266 // Job E_XX : extended version of E_XX (including flux plots)
267 // - select g/h separation method XX
268 // - read MC root file
269   // - calculate eff. collection area
270 // - optimize energy estimator
271 // - read ON root file
272 // - apply final cuts
273 // - calculate flux
274 // - write root file for ON data after final cuts
275
276
277 Bool_t JobE_XX = kFALSE;
278 Bool_t CCollArea= kFALSE; // calculate eff. collection areas
279 Bool_t OEEst = kFALSE; // optimize energy estimator
280 Bool_t WEX = kFALSE; // update root file ?
281 Bool_t WRobert = kFALSE; // write out Robert's file ?
282
283
284
285 //************************************************************************
286
287
288 //---------------------------------------------------------------------
289 // Job A
290 //=========
291
292 // - produce the histograms "sigmabar versus Theta", etc.
293 // for ON, OFF and MC data (to be used for the padding)
294 //
295 // - write root file of padded ON (OFF, MC) events (ON1.root, ...)
296 // (after the standard cuts, before the g/h separation)
297
298
299 if (JobA)
300 {
301 gLog << "=====================================================" << endl;
302 gLog << "Macro ONOFFAnalysis : Start of Job A" << endl;
303 gLog << "" << endl;
304 gLog << "Macro ONOFFAnalysis : JobA, WPad, RPad, Wout = "
305 << (JobA ? "kTRUE" : "kFALSE") << ", "
306 << (WPad ? "kTRUE" : "kFALSE") << ", "
307 << (RPad ? "kTRUE" : "kFALSE") << ", "
308 << (Wout ? "kTRUE" : "kFALSE") << endl;
309
310
311 //--------------------------------------------------
312 // names of ON and OFF files to be read
313 // for generating the histograms to be used in the padding
314 TString fileON = onfile;
315 TString fileOFF = offfile;
316 TString fileMC = mcfile;
317 gLog << "fileON, fileOFF, fileMC = " << fileON << ", "
318 << fileOFF << ", " << fileMC << endl;
319
320 // name of file to conatin the histograms for the padding
321 TString outNameSigTh = outPath;
322 outNameSigTh += "SigmaTheta";
323 outNameSigTh += ".root";
324
325 //--------------------------------------------------
326 // type of data to be padded
327 TString typeInput = "ON";
328 //TString typeInput = "OFF";
329 //TString typeInput = "MC";
330 gLog << "typeInput = " << typeInput << endl;
331
332
333 // name of input root file
334 if (typeInput == "ON")
335 TString filenamein(onfile);
336 else if (typeInput == "OFF")
337 TString filenamein(offfile);
338 else if (typeInput == "MC")
339 TString filenamein(mcfile);
340 gLog << "data to be padded : " << filenamein << endl;
341
342 // name of output root file
343 TString outNameImage = outPath;
344 outNameImage += typeInput;
345 outNameImage += "1.root";
346 gLog << "padded data to be written onto : " << outNameImage << endl;
347
348 //--------------------------------------------------
349
350 //************************************************************
351 // generate histograms to be used in the padding
352 //
353 // read ON, OFF and MC data files
354 // generate (or read in) the padding histograms for ON and OFF data
355 // and merge these histograms
356
357 MPad pad;
358 pad.SetName("MPad");
359 pad.SetDataType(typeInput);
360
361 // generate the padding histograms
362 if (GPad)
363 {
364 gLog << "=====================================================" << endl;
365 gLog << "Start generating the padding histograms" << endl;
366 //-----------------------------------------
367 // ON events
368
369 gLog << "-----------" << endl;
370 gLog << "ON events :" << endl;
371 gLog << "-----------" << endl;
372
373 MTaskList tliston;
374 MParList pliston;
375
376 MReadMarsFile readON("Events", fileON);
377 read.DisableAutoScheme();
378 MCT1ReadPreProc readON(fileON);
379
380 //MFSelBasic selthetaon;
381 //selthetaon.SetCuts(-100.0, 29.5, 35.5);
382 //MContinue contthetaon(&selthetaon);
383
384 MBlindPixelCalc blindon;
385 blindon.SetUseBlindPixels();
386
387 MFSelBasic selbasicon;
388 MContinue contbasicon(&selbasicon);
389
390 MHBlindPixels blindON("BlindPixelsON");
391 MFillH fillblindON("BlindPixelsON[MHBlindPixels]", "MBlindPixels");
392 fillblindON.SetName("FillBlind");
393
394 MSigmabarCalc sigbarcalcon;
395
396 MHSigmaTheta sigthON("SigmaThetaON");
397 MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt");
398 fillsigthetaON.SetName("FillSigTheta");
399
400 //*****************************
401 // entries in MParList
402
403 pliston.AddToList(&tliston);
404 InitBinnings(&pliston);
405 pliston.AddToList(&blindON);
406 pliston.AddToList(&sigthON);
407
408
409 //*****************************
410 // entries in MTaskList
411
412 tliston.AddToList(&readON);
413 //tliston.AddToList(&contthetaon);
414
415 tliston.AddToList(&blindon);
416
417 tliston.AddToList(&contbasicon);
418 tliston.AddToList(&fillblindON);
419 tliston.AddToList(&sigbarcalcon);
420 tliston.AddToList(&fillsigthetaON);
421
422 MProgressBar baron;
423 MEvtLoop evtloopon;
424 evtloopon.SetParList(&pliston);
425 evtloopon.SetProgressBar(&baron);
426
427 Int_t maxeventson = -1;
428 //Int_t maxeventson = 10000;
429 if ( !evtloopon.Eventloop(maxeventson) )
430 return;
431
432 tliston.PrintStatistics(0, kTRUE);
433
434 blindON.DrawClone();
435 sigthON.DrawClone();
436
437 // save the histograms for the padding
438 TH2D *sigthon = sigthON.GetSigmaTheta();
439 TH3D *sigpixthon = sigthON.GetSigmaPixTheta();
440 TH3D *diffpixthon = sigthON.GetDiffPixTheta();
441
442 TH2D *blindidthon = blindON.GetBlindId();
443 TH2D *blindnthon = blindON.GetBlindN();
444
445 //-----------------------------------------
446 // OFF events
447
448 gLog << "------------" << endl;
449 gLog << "OFF events :" << endl;
450 gLog << "------------" << endl;
451
452 MTaskList tlistoff;
453 MParList plistoff;
454
455 MReadMarsFile readOFF("Events", fileOFF);
456 read.DisableAutoScheme();
457 // MCT1ReadPreProc readOFF(fileOFF);
458
459 MFSelBasic selthetaoff;
460 selthetaoff.SetCuts(-100.0, 29.5, 35.5);
461 MContinue contthetaoff(&selthetaoff);
462
463 MBlindPixelCalc blindoff;
464 blindoff.SetUseBlindPixels();
465
466 MFSelBasic selbasicoff;
467 MContinue contbasicoff(&selbasicoff);
468
469 MHBlindPixels blindOFF("BlindPixelsOFF");
470 MFillH fillblindOFF("BlindPixelsOFF[MHBlindPixels]", "MBlindPixels");
471 fillblindOFF.SetName("FillBlindOFF");
472
473 MSigmabarCalc sigbarcalcoff;
474
475 MHSigmaTheta sigthOFF("SigmaThetaOFF");
476 MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt");
477 fillsigthetaOFF.SetName("FillSigThetaOFF");
478
479 //*****************************
480 // entries in MParList
481
482 plistoff.AddToList(&tlistoff);
483 InitBinnings(&plistoff);
484 plistoff.AddToList(&blindOFF);
485 plistoff.AddToList(&sigthOFF);
486
487
488 //*****************************
489 // entries in MTaskList
490
491 tlistoff.AddToList(&readOFF);
492 //tlistoff.AddToList(&contthetaoff);
493
494 tlistoff.AddToList(&blindoff);
495
496 tlistoff.AddToList(&contbasicoff);
497 tlistoff.AddToList(&fillblindOFF);
498 tlistoff.AddToList(&sigbarcalcoff);
499 tlistoff.AddToList(&fillsigthetaOFF);
500
501 MProgressBar baroff;
502 MEvtLoop evtloopoff;
503 evtloopoff.SetParList(&plistoff);
504 evtloopoff.SetProgressBar(&baroff);
505
506 Int_t maxeventsoff = -1;
507 //Int_t maxeventsoff = 20000;
508 if ( !evtloopoff.Eventloop(maxeventsoff) )
509 return;
510
511 tlistoff.PrintStatistics(0, kTRUE);
512
513 blindOFF.DrawClone();
514 sigthOFF.DrawClone();
515
516 // save the histograms for the padding
517 TH2D *sigthoff = sigthOFF.GetSigmaTheta();
518 TH3D *sigpixthoff = sigthOFF.GetSigmaPixTheta();
519 TH3D *diffpixthoff = sigthOFF.GetDiffPixTheta();
520
521 TH2D *blindidthoff = blindOFF.GetBlindId();
522 TH2D *blindnthoff = blindOFF.GetBlindN();
523
524
525 //-----------------------------------------
526 // MC events
527
528 gLog << "------------" << endl;
529 gLog << "MC events :" << endl;
530 gLog << "------------" << endl;
531
532 MTaskList tlistmc;
533 MParList plistmc;
534
535 MReadMarsFile readMC("Events", fileMC);
536 read.DisableAutoScheme();
537 // MCT1ReadPreProc readMC(fileMC);
538
539 MFSelBasic selthetamc;
540 selthetamc.SetCuts(-100.0, 29.5, 35.5);
541 MContinue contthetamc(&selthetamc);
542
543 MBlindPixelCalc blindmc;
544 blindmc.SetUseBlindPixels();
545
546 MFSelBasic selbasicmc;
547 MContinue contbasicmc(&selbasicmc);
548
549 MHBlindPixels blindMC("BlindPixelsMC");
550 MFillH fillblindMC("BlindPixelsMC[MHBlindPixels]", "MBlindPixels");
551 fillblindMC.SetName("FillBlindMC");
552
553 MSigmabarCalc sigbarcalcmc;
554
555 MHSigmaTheta sigthMC("SigmaThetaMC");
556 MFillH fillsigthetaMC ("SigmaThetaMC[MHSigmaTheta]", "MMcEvt");
557 fillsigthetaMC.SetName("FillSigThetaMC");
558
559 //*****************************
560 // entries in MParList
561
562 plistmc.AddToList(&tlistmc);
563 InitBinnings(&plistmc);
564 plistmc.AddToList(&blindMC);
565 plistmc.AddToList(&sigthMC);
566
567
568 //*****************************
569 // entries in MTaskList
570
571 tlistmc.AddToList(&readMC);
572 //tlistmc.AddToList(&contthetamc);
573
574 tlistmc.AddToList(&blindmc);
575
576 tlistmc.AddToList(&contbasicmc);
577 tlistmc.AddToList(&fillblindMC);
578 tlistmc.AddToList(&sigbarcalcmc);
579 tlistmc.AddToList(&fillsigthetaMC);
580
581 MProgressBar barmc;
582 MEvtLoop evtloopmc;
583 evtloopmc.SetParList(&plistmc);
584 evtloopmc.SetProgressBar(&barmc);
585
586 Int_t maxeventsmc = -1;
587 //Int_t maxeventsmc = 20000;
588 if ( !evtloopmc.Eventloop(maxeventsmc) )
589 return;
590
591 tlistmc.PrintStatistics(0, kTRUE);
592
593 blindMC.DrawClone();
594 sigthMC.DrawClone();
595
596 // save the histograms for the padding
597 TH2D *sigthmc = sigthMC.GetSigmaTheta();
598 TH3D *sigpixthmc = sigthMC.GetSigmaPixTheta();
599 TH3D *diffpixthmc = sigthMC.GetDiffPixTheta();
600
601 TH2D *blindidthmc = blindMC.GetBlindId();
602 TH2D *blindnthmc = blindMC.GetBlindN();
603
604
605 //-----------------------------------------
606
607 gLog << "End of generating the padding histograms" << endl;
608 gLog << "=====================================================" << endl;
609
610 pad.MergeONOFFMC(sigthmc, diffpixthmc, sigpixthmc,
611 blindidthmc, blindnthmc,
612 sigthon, diffpixthon, sigpixthon,
613 blindidthon, blindnthon,
614 sigthoff, diffpixthoff, sigpixthoff,
615 blindidthoff, blindnthoff);
616
617
618 if (WPad)
619 {
620 // write the padding histograms onto a file ---------
621 pad.WritePaddingDist(outNameSigTh);
622 }
623 }
624
625 // read the padding histograms ---------------------------
626 if (RPad)
627 {
628 pad.ReadPaddingDist(outNameSigTh);
629 }
630
631
632 //************************************************************
633
634 if (Wout)
635 {
636 gLog << "=====================================================" << endl;
637 gLog << "Start the padding" << endl;
638
639 //-----------------------------------------------------------
640 MTaskList tliston;
641 MParList pliston;
642
643 char *sourceName = "MSrcPosCam";
644 MSrcPosCam source(sourceName);
645
646 // geometry is needed in MHHillas... classes
647 MGeomCam *fGeom =
648 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
649
650 //-------------------------------------------
651 // create the tasks which should be executed
652 //
653
654 MReadMarsFile read("Events", filenamein);
655 read.DisableAutoScheme();
656
657 MGeomApply apply;
658 MMcPedestalCopy pcopy;
659 MMcPedestalNSBAdd pnsb;
660
661 MPedestalWorkaround waround;
662
663 // a way to find out whether one is dealing with MC :
664 MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC
665 f1.SetName("Select MC");
666 MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data
667 f2.SetName("Select Data");
668
669 MCerPhotCalc ncalc;
670 ncalc.SetFilter(&f1);
671 MCerPhotAnal2 nanal;
672 nanal.SetFilter(&f2);
673
674 //if (typeInput == "ON")
675 //{
676 // MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc",
677 // "MCT1PointingCorrCalc");
678 //}
679
680 MBlindPixelCalc blindbeforepad;
681 //blindbeforepad.SetUseBlindPixels();
682 blindbeforepad.SetName("BlindBeforePadding");
683
684 MBlindPixelCalc blind;
685 //blind.SetUseBlindPixels();
686 blind.SetName("BlindAfterPadding");
687
688 MFSelBasic selbasic;
689 MContinue contbasic(&selbasic);
690 contbasic.SetName("SelBasic");
691
692 MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels");
693 fillblind.SetName("HBlind");
694
695 MSigmabarCalc sigbarcalc;
696
697 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
698 fillsigtheta.SetName("HSigmaTheta");
699
700 MImgCleanStd clean;
701
702
703 // calculation of image parameters ---------------------
704 TString fHilName = "MHillas";
705 TString fHilNameExt = "MHillasExt";
706 TString fHilNameSrc = "MHillasSrc";
707 TString fImgParName = "MNewImagePar";
708
709 MHillasCalc hcalc;
710 hcalc.SetNameHillas(fHilName);
711 hcalc.SetNameHillasExt(fHilNameExt);
712 hcalc.SetNameNewImgPar(fImgParName);
713
714 MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc);
715 hsrccalc.SetInput(fHilName);
716
717 MFillH hfill1("MHHillas", fHilName);
718 hfill1.SetName("HHillas");
719
720 MFillH hfill2("MHStarMap", fHilName);
721 hfill2.SetName("HStarMap");
722
723 MFillH hfill3("MHHillasExt", fHilNameSrc);
724 hfill3.SetName("HHillasExt");
725
726 MFillH hfill4("MHHillasSrc", fHilNameSrc);
727 hfill4.SetName("HHillasSrc");
728
729 MFillH hfill5("MHNewImagePar", fImgParName);
730 hfill5.SetName("HNewImagePar");
731 // --------------------------------------------------
732
733 MFSelStandard selstandard(fHilNameSrc);
734 selstandard.SetHillasName(fHilName);
735 selstandard.SetImgParName(fImgParName);
736 selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0);
737 MContinue contstandard(&selstandard);
738 contstandard.SetName("SelStandard");
739
740
741 MWriteRootFile write(outNameImage);
742
743 write.AddContainer("MRawRunHeader", "RunHeaders");
744 write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
745 //write.AddContainer("MTime", "Events");
746 write.AddContainer("MMcEvt", "Events");
747 //write.AddContainer("ThetaOrig", "Events");
748 write.AddContainer("MSrcPosCam", "Events");
749 write.AddContainer("MSigmabar", "Events");
750 write.AddContainer("MHillas", "Events");
751 write.AddContainer("MHillasExt", "Events");
752 write.AddContainer("MHillasSrc", "Events");
753 write.AddContainer("MNewImagePar", "Events");
754
755
756 //*****************************
757 // entries in MParList
758
759 pliston.AddToList(&tliston);
760 InitBinnings(&pliston);
761
762 pliston.AddToList(&source);
763
764 //*****************************
765 // entries in MTaskList
766
767 tliston.AddToList(&read);
768
769 tliston.AddToList(&f1);
770 tliston.AddToList(&f2);
771 tliston.AddToList(&apply);
772 tliston.AddToList(&pcopy);
773 //tliston.AddToList(&waround);
774
775 tliston.AddToList(&pnsb);
776 tliston.AddToList(&ncalc);
777 tliston.AddToList(&nanal);
778
779 tliston.AddToList(&blindbeforepad);
780 // tliston.AddToList(&pad);
781 // if (typeInput == "ON")
782 // tliston.AddToList(&pointcorr);
783 tliston.AddToList(&blind);
784 tliston.AddToList(&contbasic);
785
786 tliston.AddToList(&fillblind);
787 tliston.AddToList(&sigbarcalc);
788 tliston.AddToList(&fillsigtheta);
789 tliston.AddToList(&clean);
790
791 tliston.AddToList(&hcalc);
792 tliston.AddToList(&hsrccalc);
793
794 tliston.AddToList(&hfill1);
795 tliston.AddToList(&hfill2);
796 tliston.AddToList(&hfill3);
797 tliston.AddToList(&hfill4);
798 tliston.AddToList(&hfill5);
799
800 tliston.AddToList(&contstandard);
801 tliston.AddToList(&write);
802
803 //*****************************
804
805 //-------------------------------------------
806 // Execute event loop
807 //
808 MProgressBar bar;
809 MEvtLoop evtloop;
810 evtloop.SetParList(&pliston);
811 //evtloop.ReadEnv(env, "", printEnv);
812 evtloop.SetProgressBar(&bar);
813 // evtloop.Write();
814
815 Int_t maxevents = -1;
816 //Int_t maxevents = 1000;
817 if ( !evtloop.Eventloop(maxevents) )
818 return;
819
820 tliston.PrintStatistics(0, kTRUE);
821
822
823 //-------------------------------------------
824 // Display the histograms
825
826 pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone();
827 pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone();
828
829 pliston.FindObject("MHHillas")->DrawClone();
830 pliston.FindObject("MHHillasExt")->DrawClone();
831 pliston.FindObject("MHHillasSrc")->DrawClone();
832 pliston.FindObject("MHNewImagePar")->DrawClone();
833 pliston.FindObject("MHStarMap")->DrawClone();
834
835 DeleteBinnings(&pliston);
836
837 gLog << "End of padding" << endl;
838 gLog << "=====================================================" << endl;
839 }
840
841
842 gLog << "Macro ONOFFAnalysis : End of Job A" << endl;
843 gLog << "===================================================" << endl;
844 }
845
846
847
848
849
850 //---------------------------------------------------------------------
851 // Job B_RF_UP
852 //============
853
854
855 // - create (or read in) the matrices of training events for gammas
856 // and hadrons
857 // - create (or read in) the trees
858 // - then read ON1.root (or MC1.root) file
859 // - calculate the hadroness for the method of RANDOM FOREST
860 // - update input root file with the hadroness
861
862
863 if (JobB_RF_UP)
864 {
865 gLog << "=====================================================" << endl;
866 gLog << "Macro ONOFFAnalysis : Start of Job B_RF_UP" << endl;
867
868 gLog << "" << endl;
869 gLog << "Macro ONOFFAnalysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
870 << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", "
871 << (RTrainRF ? "kTRUE" : "kFALSE") << ", "
872 << (CTrainRF ? "kTRUE" : "kFALSE") << ", "
873 << (RTree ? "kTRUE" : "kFALSE") << ", "
874 << (WRF ? "kTRUE" : "kFALSE") << endl;
875
876
877 //--------------------------------------------
878 // parameters for the random forest
879 Int_t NumTrees = 100;
880 Int_t NumTry = 3;
881 Int_t NdSize = 1;
882
883
884 TString hadRFName = "HadRF";
885 Float_t maxhadronness = 0.23;
886 Float_t maxalpha = 20.0;
887 Float_t maxdist = 10.0;
888
889 TString fHilName = "MHillas";
890 TString fHilNameExt = "MHillasExt";
891 TString fHilNameSrc = "MHillasSrc";
892 TString fImgParName = "MNewImagePar";
893
894
895 TString extin = "1.root";
896 TString extout = "2.root";
897
898 //--------------------------------------------
899 // for the analysis using ON data only set typeMatrixHadrons = "ON"
900 // ON and OFF data = "OFF"
901 TString typeMatrixHadrons = "OFF";
902 gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
903
904
905 // file to be updated (ON, OFF or MC)
906
907 //TString typeInput = "ON";
908 TString typeInput = "OFF";
909 //TString typeInput = "MC";
910 gLog << "typeInput = " << typeInput << endl;
911
912 // name of input root file
913 TString NameData = outPath;
914 NameData += typeInput;
915 TString inNameData(NameData);
916 inNameData += extin;
917 gLog << "inNameData = " << inNameData << endl;
918
919 // name of output root file
920 TString outNameData(NameData);
921 outNameData += extout;
922 gLog << "outNameData = " << outNameData << endl;
923
924 //--------------------------------------------
925 // files to be read for generating
926 // - the matrices of training events
927 // - and the root files of training and test events
928
929
930 // "hadrons" :
931 TString filenameHad = outPath;
932 filenameHad += typeMatrixHadrons;
933 filenameHad += extin;
934 Int_t howManyHadronsTrain = 12000;
935 Int_t howManyHadronsTest = 12000;
936 gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = "
937 << howManyHadronsTrain << ", howManyHadronsTest = "
938 << howManyHadronsTest << endl;
939
940
941 // "gammas" :
942 TString filenameMC = outPath;
943 filenameMC += "MC";
944 filenameMC += extin;
945 Int_t howManyGammasTrain = 12000;
946 Int_t howManyGammasTest = 12000;
947 gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = "
948 << howManyGammasTrain << ", howManyGammasTest = "
949 << howManyGammasTest << endl;
950
951 //--------------------------------------------
952 // files for the matrices of training events
953
954 TString NameGammas = outPath;
955 NameGammas += "RFmatrix_gammas_Train_";
956 NameGammas += "MC";
957 NameGammas += extin;
958
959 TString NameHadrons = outPath;
960 NameHadrons += "RFmatrix_hadrons_Train_";
961 NameHadrons += typeMatrixHadrons;
962 NameHadrons += extin;
963
964
965 //--------------------------------------------
966 // root files for the training events
967
968 TString NameGammasTrain = outPath;
969 NameGammasTrain += "RF_gammas_Train_";
970 NameGammasTrain += "MC";
971 TString inNameGammasTrain(NameGammasTrain);
972 inNameGammasTrain += extin;
973 TString outNameGammasTrain(NameGammasTrain);
974 outNameGammasTrain += extout;
975
976
977 TString NameHadronsTrain = outPath;
978 NameHadronsTrain += "RF_hadrons_Train_";
979 NameHadronsTrain += typeMatrixHadrons;
980 TString inNameHadronsTrain(NameHadronsTrain);
981 inNameHadronsTrain += extin;
982 TString outNameHadronsTrain(NameHadronsTrain);
983 outNameHadronsTrain += extout;
984
985
986 //--------------------------------------------
987 // root files for the test events
988
989 TString NameGammasTest = outPath;
990 NameGammasTest += "RF_gammas_Test_";
991 NameGammasTest += "MC";
992 TString inNameGammasTest(NameGammasTest);
993 inNameGammasTest += extin;
994 TString outNameGammasTest(NameGammasTest);
995 outNameGammasTest += extout;
996
997 TString NameHadronsTest = outPath;
998 NameHadronsTest += "RF_hadrons_Test_";
999 NameHadronsTest += typeMatrixHadrons;
1000 TString inNameHadronsTest(NameHadronsTest);
1001 inNameHadronsTest += extin;
1002 TString outNameHadronsTest(NameHadronsTest);
1003 outNameHadronsTest += extout;
1004
1005 //--------------------------------------------------------------------
1006
1007
1008 MHMatrix matrixg("MatrixGammas");
1009 matrixg.EnableGraphicalOutput();
1010
1011 matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
1012 matrixg.AddColumn("MSigmabar.fSigmabar");
1013 matrixg.AddColumn("log10(MHillas.fSize)");
1014 matrixg.AddColumn("MHillasSrc.fDist");
1015 matrixg.AddColumn("MHillas.fWidth");
1016 matrixg.AddColumn("MHillas.fLength");
1017 matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
1018 matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
1019 matrixg.AddColumn("MNewImagePar.fConc");
1020 matrixg.AddColumn("MNewImagePar.fLeakage1");
1021
1022 MHMatrix matrixh("MatrixHadrons");
1023 matrixh.EnableGraphicalOutput();
1024
1025 matrixh.AddColumns(matrixg.GetColumns());
1026
1027 //--------------------------------------------
1028 // file of trees of the random forest
1029
1030 TString outRF = outPath;
1031 outRF += "RF.root";
1032
1033
1034 //*************************************************************************
1035 // read in matrices of training events
1036if (RTrainRF)
1037 {
1038 const char* mtxName = "MatrixGammas";
1039
1040 gLog << "" << endl;
1041 gLog << "========================================================" << endl;
1042 gLog << "Get matrix for (gammas)" << endl;
1043 gLog << "matrix name = " << mtxName << endl;
1044 gLog << "name of root file = " << NameGammas << endl;
1045 gLog << "" << endl;
1046
1047
1048 // read in the object with the name 'mtxName' from file 'NameGammas'
1049 //
1050 TFile fileg(NameGammas);
1051
1052 matrixg.Read(mtxName);
1053 matrixg.Print("SizeCols");
1054
1055
1056 //*****************************************************************
1057
1058 const char* mtxName = "MatrixHadrons";
1059
1060 gLog << "" << endl;
1061 gLog << "========================================================" << endl;
1062 gLog << " Get matrix for (hadrons)" << endl;
1063 gLog << "matrix name = " << mtxName << endl;
1064 gLog << "name of root file = " << NameHadrons << endl;
1065 gLog << "" << endl;
1066
1067
1068 // read in the object with the name 'mtxName' from file 'NameHadrons'
1069 //
1070 TFile fileh(NameHadrons);
1071
1072 matrixh.Read(mtxName);
1073 matrixh.Print("SizeCols");
1074 }
1075
1076
1077 //*************************************************************************
1078 // create matrices of training events
1079 // and root files of training and test events
1080
1081if (CTrainRF)
1082 {
1083 gLog << "" << endl;
1084 gLog << "========================================================" << endl;
1085 gLog << " Create matrices of training events and root files of training and test events"
1086 << endl;
1087 gLog << " Gammas :" << endl;
1088 gLog << "---------" << endl;
1089
1090 MParList plistg;
1091 MTaskList tlistg;
1092
1093 MReadMarsFile readg("Events", filenameMC);
1094 readg.DisableAutoScheme();
1095
1096 TString mgname("costhg");
1097 MBinning bing("Binning"+mgname);
1098 bing.SetEdges(10, 0., 1.0);
1099
1100 MH3 gref("cos(MMcEvt.fTelescopeTheta)");
1101 gref.SetName(mgname);
1102 MH::SetBinning(&gref.GetHist(), &bing);
1103 //for (Int_t i=1; i<=gref.GetNbins(); i++)
1104 // gref.GetHist().SetBinContent(i, 1.0);
1105
1106 MFEventSelector2 selectorg(gref);
1107 selectorg.SetNumMax(howManyGammasTrain+howManyGammasTest);
1108 selectorg.SetName("selectGammasTrainTest");
1109 selectorg.SetInverted();
1110 //selectorg.SetUseOrigDistribution(kTRUE);
1111
1112 MContinue contg(&selectorg);
1113 contg.SetName("ContGammas");
1114
1115 Double_t probg = ( (Double_t) howManyGammasTrain )
1116 / ( (Double_t)(howManyGammasTrain+howManyGammasTest) );
1117 MFRandomSplit splitg(probg);
1118
1119 MFillH fillmatg("MatrixGammas");
1120 fillmatg.SetFilter(&splitg);
1121 fillmatg.SetName("fillGammas");
1122
1123 //-----------------------
1124 // for writing the root files of training and test events
1125 // for gammas
1126
1127 MWriteRootFile writetraing(inNameGammasTrain, "RECREATE");
1128 writetraing.SetName("WriteGammasTrain");
1129 writetraing.SetFilter(&splitg);
1130
1131 writetraing.AddContainer("MRawRunHeader", "RunHeaders");
1132 writetraing.AddContainer("MTime", "Events");
1133 writetraing.AddContainer("MMcEvt", "Events");
1134 writetraing.AddContainer("ThetaOrig", "Events");
1135 writetraing.AddContainer("MSrcPosCam", "Events");
1136 writetraing.AddContainer("MSigmabar", "Events");
1137 writetraing.AddContainer("MHillas", "Events");
1138 writetraing.AddContainer("MHillasExt", "Events");
1139 writetraing.AddContainer("MHillasSrc", "Events");
1140 writetraing.AddContainer("MNewImagePar", "Events");
1141
1142 MContinue contgtrain(&splitg);
1143 contgtrain.SetName("ContGammaTrain");
1144
1145 MWriteRootFile writetestg(inNameGammasTest, "RECREATE");
1146 writetestg.SetName("WriteGammasTest");
1147
1148 writetestg.AddContainer("MRawRunHeader", "RunHeaders");
1149 writetestg.AddContainer("MTime", "Events");
1150 writetestg.AddContainer("MMcEvt", "Events");
1151 writetestg.AddContainer("ThetaOrig", "Events");
1152 writetestg.AddContainer("MSrcPosCam", "Events");
1153 writetestg.AddContainer("MSigmabar", "Events");
1154 writetestg.AddContainer("MHillas", "Events");
1155 writetestg.AddContainer("MHillasExt", "Events");
1156 writetestg.AddContainer("MHillasSrc", "Events");
1157 writetestg.AddContainer("MNewImagePar", "Events");
1158
1159 //-----------------------
1160
1161 //***************************** fill gammas ***
1162 // entries in MParList
1163
1164 plistg.AddToList(&tlistg);
1165 InitBinnings(&plistg);
1166
1167 plistg.AddToList(&matrixg);
1168
1169 //*****************************
1170 // entries in MTaskList
1171
1172 tlistg.AddToList(&readg);
1173 tlistg.AddToList(&contg);
1174
1175 tlistg.AddToList(&splitg);
1176 tlistg.AddToList(&fillmatg);
1177 tlistg.AddToList(&writetraing);
1178 tlistg.AddToList(&contgtrain);
1179
1180 tlistg.AddToList(&writetestg);
1181
1182 //*****************************
1183
1184 MProgressBar matrixbar;
1185 MEvtLoop evtloopg;
1186 evtloopg.SetName("FillGammaMatrix");
1187 evtloopg.SetParList(&plistg);
1188 //evtloopg.ReadEnv(env, "", printEnv);
1189 evtloopg.SetProgressBar(&matrixbar);
1190
1191 Int_t maxevents = -1;
1192 if (!evtloopg.Eventloop(maxevents))
1193 return;
1194
1195 tlistg.PrintStatistics(0, kTRUE);
1196
1197 matrixg.Print("SizeCols");
1198 Int_t generatedgTrain = matrixg.GetM().GetNrows();
1199 if ( fabs(generatedgTrain-howManyGammasTrain) >
1200 3.0*sqrt(howManyGammasTrain) )
1201 {
1202 gLog << "ONOFFAnalysis.C : no.of generated gamma training events ("
1203 << generatedgTrain << ") is incompatible with the no.of requested events ("
1204 << howManyGammasTrain << ")" << endl;
1205 }
1206
1207
1208 Int_t generatedgTest = writetestg.GetNumExecutions();
1209 if ( fabs(generatedgTest-howManyGammasTest) >
1210 3.0*sqrt(howManyGammasTest) )
1211 {
1212 gLog << "ONOFFAnalysis.C : no.of generated gamma test events ("
1213 << generatedgTest << ") is incompatible with the no.of requested events ("
1214 << howManyGammasTest << ")" << endl;
1215 }
1216
1217 //***************************** fill hadrons ***
1218 gLog << "---------------------------------------------------------------"
1219 << endl;
1220 gLog << " Hadrons :" << endl;
1221 gLog << "----------" << endl;
1222
1223 MParList plisth;
1224 MTaskList tlisth;
1225
1226 MReadMarsFile readh("Events", filenameHad);
1227 readh.DisableAutoScheme();
1228
1229 TString mhname("costhh");
1230 MBinning binh("Binning"+mhname);
1231 binh.SetEdges(10, 0., 1.0);
1232
1233 //MH3 href("cos(MMcEvt.fTelescopeTheta)");
1234 //href.SetName(mhname);
1235 //MH::SetBinning(&href.GetHist(), &binh);
1236 //for (Int_t i=1; i<=href.GetNbins(); i++)
1237 // href.GetHist().SetBinContent(i, 1.0);
1238
1239 //use the original distribution from the gammas
1240 MH3 &href = *(selectorg.GetHistOrig());
1241
1242 MFEventSelector2 selectorh(href);
1243 selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest);
1244 selectorh.SetName("selectHadronsTrainTest");
1245 selectorh.SetInverted();
1246
1247 MContinue conth(&selectorh);
1248 conth.SetName("ContHadrons");
1249
1250 Double_t probh = ( (Double_t) howManyHadronsTrain )
1251 / ( (Double_t)(howManyHadronsTrain+howManyHadronsTest) );
1252 MFRandomSplit splith(probh);
1253
1254 MFillH fillmath("MatrixHadrons");
1255 fillmath.SetFilter(&splith);
1256 fillmath.SetName("fillHadrons");
1257
1258 //-----------------------
1259 // for writing the root files of training and test events
1260 // for hadrons
1261
1262 MWriteRootFile writetrainh(inNameHadronsTrain, "RECREATE");
1263 writetrainh.SetName("WriteHadronsTrain");
1264 writetrainh.SetFilter(&splith);
1265
1266 writetrainh.AddContainer("MRawRunHeader", "RunHeaders");
1267 writetrainh.AddContainer("MTime", "Events");
1268 writetrainh.AddContainer("MMcEvt", "Events");
1269 writetrainh.AddContainer("ThetaOrig", "Events");
1270 writetrainh.AddContainer("MSrcPosCam", "Events");
1271 writetrainh.AddContainer("MSigmabar", "Events");
1272 writetrainh.AddContainer("MHillas", "Events");
1273 writetrainh.AddContainer("MHillasExt", "Events");
1274 writetrainh.AddContainer("MHillasSrc", "Events");
1275 writetrainh.AddContainer("MNewImagePar", "Events");
1276
1277 MContinue conthtrain(&splith);
1278
1279 MWriteRootFile writetesth(inNameHadronsTest, "RECREATE");
1280 writetesth.SetName("WriteHadronsTest");
1281
1282 writetesth.AddContainer("MRawRunHeader", "RunHeaders");
1283 writetesth.AddContainer("MTime", "Events");
1284 writetesth.AddContainer("MMcEvt", "Events");
1285 writetesth.AddContainer("ThetaOrig", "Events");
1286 writetesth.AddContainer("MSrcPosCam", "Events");
1287 writetesth.AddContainer("MSigmabar", "Events");
1288 writetesth.AddContainer("MHillas", "Events");
1289 writetesth.AddContainer("MHillasExt", "Events");
1290 writetesth.AddContainer("MHillasSrc", "Events");
1291 writetesth.AddContainer("MNewImagePar", "Events");
1292
1293
1294 //*****************************
1295 // entries in MParList
1296
1297 plisth.AddToList(&tlisth);
1298 InitBinnings(&plisth);
1299
1300 plisth.AddToList(&matrixh);
1301
1302 //*****************************
1303 // entries in MTaskList
1304
1305 tlisth.AddToList(&readh);
1306 tlisth.AddToList(&conth);
1307
1308 tlisth.AddToList(&splith);
1309 tlisth.AddToList(&fillmath);
1310 tlisth.AddToList(&writetrainh);
1311 tlisth.AddToList(&conthtrain);
1312
1313 tlisth.AddToList(&writetesth);
1314
1315 //*****************************
1316
1317 MProgressBar matrixbar;
1318 MEvtLoop evtlooph;
1319 evtlooph.SetName("FillHadronMatrix");
1320 evtlooph.SetParList(&plisth);
1321 //evtlooph.ReadEnv(env, "", printEnv);
1322 evtlooph.SetProgressBar(&matrixbar);
1323
1324 Int_t maxevents = -1;
1325 if (!evtlooph.Eventloop(maxevents))
1326 return;
1327
1328 tlisth.PrintStatistics(0, kTRUE);
1329
1330 matrixh.Print("SizeCols");
1331 Int_t generatedhTrain = matrixh.GetM().GetNrows();
1332 if ( fabs(generatedhTrain-howManyHadronsTrain) >
1333 3.0*sqrt(howManyHadronsTrain) )
1334 {
1335 gLog << "ONOFFAnalysis.C : no.of generated hadron training events ("
1336 << generatedhTrain << ") is incompatible with the no.of requested events ("
1337 << howManyHadronsTrain << ")" << endl;
1338 }
1339
1340
1341 Int_t generatedhTest = writetesth.GetNumExecutions();
1342 if ( fabs(generatedhTest-howManyHadronsTest) >
1343 3.0*sqrt(howManyHadronsTest) )
1344 {
1345 gLog << "ONOFFAnalysis.C : no.of generated gamma test events ("
1346 << generatedhTest << ") is incompatible with the no.of requested events ("
1347 << howManyHadronsTest << ")" << endl;
1348 }
1349
1350
1351 //*****************************************************
1352
1353
1354 // write out matrices of training events
1355
1356 gLog << "" << endl;
1357 gLog << "========================================================" << endl;
1358 gLog << "Write out matrices of training events" << endl;
1359
1360
1361 //-------------------------------------------
1362 // "gammas"
1363 gLog << "Gammas :" << endl;
1364 matrixg.Print("SizeCols");
1365
1366 TFile writeg(NameGammas, "RECREATE", "");
1367 matrixg.Write();
1368
1369 gLog << "" << endl;
1370 gLog << "Macro ONOFFAnalysis : matrix of training events for gammas written onto file "
1371 << NameGammas << endl;
1372
1373 //-------------------------------------------
1374 // "hadrons"
1375 gLog << "Hadrons :" << endl;
1376 matrixh.Print("SizeCols");
1377
1378 TFile writeh(NameHadrons, "RECREATE", "");
1379 matrixh.Write();
1380
1381 gLog << "" << endl;
1382 gLog << "Macro ONOFFAnalysis : matrix of training events for hadrons written onto file "
1383 << NameHadrons << endl;
1384
1385 }
1386 //********** end of creating matrices of training events ***********
1387
1388
1389 MRanForest *fRanForest;
1390 MRanTree *fRanTree;
1391 //-----------------------------------------------------------------
1392 // read in the trees of the random forest
1393 if (RTree)
1394 {
1395 MParList plisttr;
1396 MTaskList tlisttr;
1397 plisttr.AddToList(&tlisttr);
1398
1399 MReadTree readtr("TREE", outRF);
1400 readtr.DisableAutoScheme();
1401
1402 MRanForestFill rffill;
1403 rffill.SetNumTrees(NumTrees);
1404
1405 // list of tasks for the loop over the trees
1406
1407 tlisttr.AddToList(&readtr);
1408 tlisttr.AddToList(&rffill);
1409
1410 //-------------------
1411 // Execute tree loop
1412 //
1413 MEvtLoop evtlooptr;
1414 evtlooptr.SetName("ReadRFTrees");
1415 evtlooptr.SetParList(&plisttr);
1416 if (!evtlooptr.Eventloop())
1417 return;
1418
1419 tlisttr.PrintStatistics(0, kTRUE);
1420
1421 gLog << "ONOFFAnalysis : RF trees were read in from file "
1422 << outRF << endl;
1423
1424 // get adresses of objects which are used in the next eventloop
1425 fRanForest = (MRanForest*)plisttr->FindObject("MRanForest");
1426 if (!fRanForest)
1427 {
1428 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1429 return kFALSE;
1430 }
1431
1432 fRanTree = (MRanTree*)plisttr->FindObject("MRanTree");
1433 if (!fRanTree)
1434 {
1435 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1436 return kFALSE;
1437 }
1438
1439 }
1440
1441 //-----------------------------------------------------------------
1442 // grow the trees of the random forest (event loop = tree loop)
1443
1444 if (!RTree)
1445 {
1446
1447 gLog << "" << endl;
1448 gLog << "========================================================" << endl;
1449 gLog << "Macro ONOFFAnalysis : start growing trees" << endl;
1450
1451 MTaskList tlist2;
1452 MParList plist2;
1453 plist2.AddToList(&tlist2);
1454
1455 plist2.AddToList(&matrixg);
1456 plist2.AddToList(&matrixh);
1457
1458 MRanForestGrow rfgrow2;
1459 rfgrow2.SetNumTrees(NumTrees);
1460 rfgrow2.SetNumTry(NumTry);
1461 rfgrow2.SetNdSize(NdSize);
1462
1463 MWriteRootFile rfwrite2(outRF);
1464 rfwrite2.AddContainer("MRanTree", "TREE");
1465
1466 MFillH fillh2("MHRanForestGini");
1467
1468 // list of tasks for the loop over the trees
1469
1470 tlist2.AddToList(&rfgrow2);
1471 tlist2.AddToList(&rfwrite2);
1472 tlist2.AddToList(&fillh2);
1473
1474 //-------------------
1475 // Execute tree loop
1476 //
1477 MEvtLoop treeloop;
1478 treeloop.SetName("GrowRFTrees");
1479 treeloop.SetParList(&plist2);
1480
1481 if ( !treeloop.Eventloop() )
1482 return;
1483
1484 tlist2.PrintStatistics(0, kTRUE);
1485
1486 plist2.FindObject("MHRanForestGini")->DrawClone();
1487
1488
1489 // get adresses of objects which are used in the next eventloop
1490 fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
1491 if (!fRanForest)
1492 {
1493 gLog << err << dbginf << "MRanForest not found... aborting." << endl;
1494 return kFALSE;
1495 }
1496
1497 fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
1498 if (!fRanTree)
1499 {
1500 gLog << err << dbginf << "MRanTree not found... aborting." << endl;
1501 return kFALSE;
1502 }
1503
1504 }
1505 // end of growing the trees of the random forest
1506 //-----------------------------------------------------------------
1507
1508
1509 //-----------------------------------------------------------------
1510 // Update the root files with the RF hadronness
1511 //
1512
1513 if (WRF)
1514 {
1515 //TString fileName(inNameHadronsTrain);
1516 //TString outName(outNameHadronsTrain);
1517
1518 //TString fileName(inNameHadronsTest);
1519 //TString outName(outNameHadronsTest);
1520
1521 //TString fileName(inNameGammasTrain);
1522 //TString outName(outNameGammasTrain);
1523
1524 //TString fileName(inNameGammasTest);
1525 //TString outName(outNameGammasTest);
1526
1527 TString fileName(inNameData);
1528 TString outName(outNameData);
1529
1530
1531
1532 gLog << "" << endl;
1533 gLog << "========================================================" << endl;
1534 gLog << "Update root file '" << fileName
1535 << "' with the RF hadronness; ==> " << outName << endl;
1536
1537
1538 MTaskList tliston;
1539 MParList pliston;
1540
1541
1542 // geometry is needed in MHHillas... classes
1543 MGeomCam *fGeom =
1544 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
1545
1546 //-------------------------------------------
1547 // create the tasks which should be executed
1548 //
1549
1550 MReadMarsFile read("Events", fileName);
1551 read.DisableAutoScheme();
1552
1553
1554 //.......................................................................
1555 // calculate hadronnes for method of RANDOM FOREST
1556
1557
1558 MRanForestCalc rfcalc;
1559 rfcalc.SetHadronnessName(hadRFName);
1560
1561
1562 //.......................................................................
1563
1564 //MWriteRootFile write(outName, "UPDATE");
1565 MWriteRootFile write(outName, "RECREATE");
1566
1567 write.AddContainer("MRawRunHeader", "RunHeaders");
1568 write.AddContainer("MTime", "Events");
1569 write.AddContainer("MMcEvt", "Events");
1570 write.AddContainer("ThetaOrig", "Events");
1571 write.AddContainer("MSrcPosCam", "Events");
1572 write.AddContainer("MSigmabar", "Events");
1573 write.AddContainer("MHillas", "Events");
1574 write.AddContainer("MHillasExt", "Events");
1575 write.AddContainer("MHillasSrc", "Events");
1576 write.AddContainer("MNewImagePar", "Events");
1577
1578 write.AddContainer(hadRFName, "Events");
1579
1580 //-----------------------------------------------------------------
1581
1582
1583 MFSelFinal selfinalgh(fHilNameSrc);
1584 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
1585 selfinalgh.SetHadronnessName(hadRFName);
1586 selfinalgh.SetName("SelFinalgh");
1587 MContinue contfinalgh(&selfinalgh);
1588 contfinalgh.SetName("ContSelFinalgh");
1589
1590 MFillH fillranfor("MHRanForest");
1591 fillranfor.SetName("HRanForest");
1592
1593 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
1594 fillhadrf.SetName("HhadRF");
1595
1596 MFSelFinal selfinal(fHilNameSrc);
1597 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
1598 selfinal.SetHadronnessName(hadRFName);
1599 selfinal.SetName("SelFinal");
1600 MContinue contfinal(&selfinal);
1601 contfinal.SetName("ContSelFinal");
1602
1603 TString mh3name = "abs(Alpha)";
1604 MBinning binsalphaabs("Binning"+mh3name);
1605 binsalphaabs.SetEdges(50, -2.0, 98.0);
1606
1607 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
1608 alphaabs.SetName(mh3name);
1609 MFillH alpha(&alphaabs);
1610 alpha.SetName("FillAlphaAbs");
1611
1612
1613 MFillH hfill1("MHHillas", fHilName);
1614 hfill1.SetName("HHillas");
1615
1616 MFillH hfill2("MHStarMap", fHilName);
1617 hfill2.SetName("HStarMap");
1618
1619 MFillH hfill3("MHHillasExt", fHilNameSrc);
1620 hfill3.SetName("HHillasExt");
1621
1622 MFillH hfill4("MHHillasSrc", fHilNameSrc);
1623 hfill4.SetName("HHillasSrc");
1624
1625 MFillH hfill5("MHNewImagePar", fImgParName);
1626 hfill5.SetName("HNewImagePar");
1627
1628 //*****************************
1629 // entries in MParList
1630
1631 pliston.AddToList(&tliston);
1632 InitBinnings(&pliston);
1633
1634 pliston.AddToList(fRanForest);
1635 pliston.AddToList(fRanTree);
1636
1637 pliston.AddToList(&binsalphaabs);
1638 pliston.AddToList(&alphaabs);
1639
1640
1641 //*****************************
1642 // entries in MTaskList
1643
1644 tliston.AddToList(&read);
1645
1646 tliston.AddToList(&rfcalc);
1647 tliston.AddToList(&fillranfor);
1648 tliston.AddToList(&fillhadrf);
1649
1650 tliston.AddToList(&write);
1651 tliston.AddToList(&contfinalgh);
1652
1653 tliston.AddToList(&alpha);
1654 tliston.AddToList(&hfill1);
1655 tliston.AddToList(&hfill2);
1656 tliston.AddToList(&hfill3);
1657 tliston.AddToList(&hfill4);
1658 tliston.AddToList(&hfill5);
1659
1660 tliston.AddToList(&contfinal);
1661
1662 //*****************************
1663
1664 //-------------------------------------------
1665 // Execute event loop
1666 //
1667 MProgressBar bar;
1668 MEvtLoop evtloop;
1669 evtloop.SetName("UpdateRootFile");
1670 evtloop.SetParList(&pliston);
1671 evtloop.SetProgressBar(&bar);
1672
1673 Int_t maxevents = -1;
1674 if ( !evtloop.Eventloop(maxevents) )
1675 return;
1676
1677 tliston.PrintStatistics(0, kTRUE);
1678
1679
1680 //-------------------------------------------
1681 // Display the histograms
1682 //
1683 pliston.FindObject("MHRanForest")->DrawClone();
1684 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
1685 pliston.FindObject("hadRF", "MHHadronness")->Print();
1686
1687 pliston.FindObject("MHHillas")->DrawClone();
1688 pliston.FindObject("MHHillasExt")->DrawClone();
1689 pliston.FindObject("MHHillasSrc")->DrawClone();
1690 pliston.FindObject("MHNewImagePar")->DrawClone();
1691 pliston.FindObject("MHStarMap")->DrawClone();
1692
1693
1694 //-------------------------------------------
1695 // fit alpha distribution to get the number of excess events and
1696 // calculate significance of gamma signal in the alpha plot
1697
1698 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
1699 TH1 &alphaHist = absalpha->GetHist();
1700 alphaHist.SetXTitle("|alpha| [\\circ]");
1701 alphaHist.SetName("alpha-macro");
1702
1703 Double_t alphasig = 13.1;
1704 Double_t alphamin = 30.0;
1705 Double_t alphamax = 90.0;
1706 Int_t degree = 2;
1707 Double_t significance = -99.0;
1708 Bool_t drawpoly = kTRUE;
1709 Bool_t fitgauss = kTRUE;
1710 Bool_t print = kTRUE;
1711
1712 MHFindSignificance findsig;
1713 findsig.SetRebin(kTRUE);
1714 findsig.SetReduceDegree(kFALSE);
1715
1716 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
1717 alphasig, drawpoly, fitgauss, print);
1718 significance = findsig.GetSignificance();
1719 Float_t alphasi = findsig.GetAlphasi();
1720
1721 gLog << "For file '" << fileName << "' : " << endl;
1722 gLog << "Significance of gamma signal after supercuts : "
1723 << significance << " (for |alpha| < " << alphasi << " degrees)"
1724 << endl;
1725
1726 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
1727
1728 //-------------------------------------------
1729
1730
1731 DeleteBinnings(&pliston);
1732 }
1733
1734 gLog << "Macro ONOFFAnalysis : End of Job B_RF_UP" << endl;
1735 gLog << "=======================================================" << endl;
1736 }
1737 //---------------------------------------------------------------------
1738
1739
1740 //---------------------------------------------------------------------
1741 // Job B_SC_UP
1742 //============
1743
1744 // - create (or read in) optimum supercuts parameter values
1745 //
1746 // - calculate the hadroness for the supercuts
1747 //
1748 // - update input root file, including the hadroness
1749
1750
1751 if (JobB_SC_UP)
1752 {
1753 gLog << "=====================================================" << endl;
1754 gLog << "Macro ONOFFAnalysis : Start of Job B_SC_UP" << endl;
1755
1756 gLog << "" << endl;
1757 gLog << "Macro ONOFFAnalysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = "
1758 << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", "
1759 << (CMatrix ? "kTRUE" : "kFALSE") << ", "
1760 << (RMatrix ? "kTRUE" : "kFALSE") << ", "
1761 << (WOptimize ? "kTRUE" : "kFALSE") << ", "
1762 << (RTest ? "kTRUE" : "kFALSE") << ", "
1763 << (WSC ? "kTRUE" : "kFALSE") << endl;
1764
1765
1766 //--------------------------------------------
1767 // file which contains the initial parameter values for the supercuts
1768 // if parSCinit ="" the initial values are taken from the constructor of
1769 // MSupercuts
1770
1771 TString parSCinit = outPath;
1772 //parSCinit += "parSC_1709d";
1773 parSCinit = "";
1774
1775 gLog << "parSCinit = " << parSCinit << endl;
1776
1777 //---------------
1778 // file onto which the optimal parameter values for the supercuts
1779 // are written
1780
1781 TString parSCfile = outPath;
1782 parSCfile += "parSC_2310a";
1783
1784 gLog << "parSCfile = " << parSCfile << endl;
1785
1786 //--------------------------------------------
1787 // file to be updated (either ON or MC)
1788
1789 //TString typeInput = "ON";
1790 //TString typeInput = "OFF";
1791 TString typeInput = "MC";
1792 gLog << "typeInput = " << typeInput << endl;
1793
1794 // name of input root file
1795 TString filenameData = outPath;
1796 filenameData += typeInput;
1797 filenameData += "2.root";
1798 gLog << "filenameData = " << filenameData << endl;
1799
1800 // name of output root file
1801 TString outNameImage = outPath;
1802 outNameImage += typeInput;
1803 outNameImage += "3.root";
1804
1805
1806 //TString outNameImage = filenameData;
1807
1808 gLog << "outNameImage = " << outNameImage << endl;
1809
1810 //--------------------------------------------
1811 // files to be read for optimizing the supercuts
1812 //
1813 // for the training
1814 TString filenameTrain = outPath;
1815 filenameTrain += "ON";
1816 filenameTrain += "1.root";
1817 Int_t howManyTrain = 800000;
1818 gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = "
1819 << howManyTrain << endl;
1820
1821 // for testing
1822 TString filenameTest = outPath;
1823 filenameTest += "ON";
1824 filenameTest += "1.root";
1825 Int_t howManyTest = 800000;
1826
1827 gLog << "filenameTest = " << filenameTest << ", howManyTest = "
1828 << howManyTest << endl;
1829
1830
1831 //--------------------------------------------
1832 // files to contain the matrices (generated from filenameTrain and
1833 // filenameTest)
1834 //
1835 // for the training
1836 TString fileMatrixTrain = outPath;
1837 fileMatrixTrain += "MatrixTrainSC";
1838 fileMatrixTrain += ".root";
1839 gLog << "fileMatrixTrain = " << fileMatrixTrain << endl;
1840
1841 // for testing
1842 TString fileMatrixTest = outPath;
1843 fileMatrixTest += "MatrixTestSC";
1844 fileMatrixTest += ".root";
1845 gLog << "fileMatrixTest = " << fileMatrixTest << endl;
1846
1847
1848
1849 //---------------------------------------------------------------------
1850 // Training and test matrices :
1851 // - either create them and write them onto a file
1852 // - or read them from a file
1853
1854
1855 MFindSupercuts findsuper;
1856 findsuper.SetFilenameParam(parSCfile);
1857 findsuper.SetHadronnessName("HadSC");
1858 //findsuper.SetUseOrigDistribution(kTRUE);
1859
1860 //--------------------------
1861 // create matrices and write them onto files
1862 if (CMatrix)
1863 {
1864 TString mname("costheta");
1865 MBinning bin("Binning"+mname);
1866 bin.SetEdges(10, 0., 1.0);
1867
1868 MH3 mh3("cos(MMcEvt.fTelescopeTheta)");
1869 mh3.SetName(mname);
1870 MH::SetBinning(&mh3.GetHist(), &bin);
1871 //for (Int_t i=1; i<=mh3.GetNbins(); i++)
1872 // mh3.GetHist().SetBinContent(i, 1.0);
1873
1874
1875 if (filenameTrain == filenameTest)
1876 {
1877 if ( !findsuper.DefineTrainTestMatrix(
1878 filenameTrain, mh3,
1879 howManyTrain, howManyTest,
1880 fileMatrixTrain, fileMatrixTest) )
1881 {
1882 *fLog << "MagicAnalysis.C : DefineTrainTestMatrix failed" << endl;
1883 return;
1884 }
1885
1886 }
1887 else
1888 {
1889 if ( !findsuper.DefineTrainMatrix(filenameTrain, mh3,
1890 howManyTrain, fileMatrixTrain) )
1891 {
1892 *fLog << "MagicAnalysis.C : DefineTrainMatrix failed" << endl;
1893 return;
1894 }
1895
1896 if ( !findsuper.DefineTestMatrix( filenameTest, mh3,
1897 howManyTest, fileMatrixTest) )
1898 {
1899 *fLog << "MagicAnalysis.C : DefineTestMatrix failed" << endl;
1900 return;
1901 }
1902 }
1903 }
1904
1905 //--------------------------
1906 // read matrices from files
1907 //
1908
1909 if (RMatrix)
1910 findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest);
1911 //--------------------------
1912
1913
1914
1915 //---------------------------------------------------------------------
1916 // optimize supercuts using the training sample
1917 //
1918 // the initial values are taken
1919 // - from the file parSCinit (if != "")
1920 // - or from the arrays params and steps (if their sizes are != 0)
1921 // - or from the MSupercuts constructor
1922
1923
1924if (WOptimize)
1925 {
1926 gLog << "MagicAnalysis.C : optimize the supercuts using the training matrix"
1927 << endl;
1928
1929 TArrayD params(0);
1930 TArrayD steps(0);
1931
1932 if (parSCinit == "")
1933 {
1934 Double_t vparams[104] = {
1935 // LengthUp
1936 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362,
1937 0.007388, -0.013463,
1938 // LengthLo
1939 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993,
1940 0.000080, -0.007157,
1941 // WidthUp
1942 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353,
1943 0.020711, -0.016703,
1944 // WidthLo
1945 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101,
1946 0.006126, -0.002849,
1947 // DistUp
1948 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909,
1949 0.189697, 0.0,
1950 // DistLo
1951 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374,
1952 -0.001750, 0.0,
1953 // AsymUp
1954 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1955 0.0, 0.0,
1956 // AsymLo
1957 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1958 0.0, 0.0,
1959 // ConcUp
1960 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1961 0.0, 0.0,
1962 // ConcLo
1963 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1964 0.0, 0.0,
1965 // Leakage1Up
1966 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1967 0.0, 0.0,
1968 // Leakage1Lo
1969 -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0,
1970 0.0, 0.0,
1971 // AlphaUp
1972 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0,
1973 0.0, 0.0 };
1974
1975 Double_t vsteps[104] = {
1976 // LengthUp
1977 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002,
1978 0.0008, 0.002,
1979 // LengthLo
1980 0.02, 0.003, 0.05, 0.006, 0.002, 0.3,
1981 0.0001, 0.0008,
1982 // WidthUp
1983 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006,
1984 0.002, 0.002,
1985 // WidthLo
1986 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002,
1987 0.0007, 0.003,
1988 // DistUp
1989 0.2, 0.0, 0.3, 0.02, 0.0, 0.03,
1990 0.02, 0.0
1991 // DistLo
1992 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005,
1993 0.0002, 0.0
1994 // AsymUp
1995 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1996 0.0, 0.0,
1997 // AsymLo
1998 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1999 0.0, 0.0,
2000 // ConcUp
2001 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
2002 0.0, 0.0,
2003 // ConcLo
2004 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
2005 0.0, 0.0,
2006 // Leakage1Up
2007 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
2008 0.0, 0.0,
2009 // Leakage1Lo
2010 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
2011 0.0, 0.0,
2012 // AlphaUp
2013 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
2014 0.0, 0.0 };
2015
2016 params.Set(104, vparams);
2017 steps.Set (104, vsteps );
2018 }
2019
2020 Bool_t rf;
2021 rf = findsuper.FindParams(parSCinit, params, steps);
2022
2023 if (!rf)
2024 {
2025 gLog << "MagicAnalysis.C : optimization of supercuts failed" << endl;
2026 return;
2027 }
2028 }
2029
2030 //--------------------------------------
2031 // test the supercuts on the test sample
2032 //
2033
2034 if (RTest)
2035 {
2036 gLog << "MagicAnalysis.C : test the supercuts on the test matrix" << endl;
2037 Bool_t rt = findsuper.TestParams();
2038 if (!rt)
2039 {
2040 gLog << "MagicAnalysis.C : test of supercuts on the test matrix failed"
2041 << endl;
2042 }
2043
2044 }
2045
2046
2047 //-----------------------------------------------------------------
2048 // Update the input files with the SC hadronness
2049 //
2050
2051 if (WSC)
2052 {
2053 gLog << "" << endl;
2054 gLog << "========================================================" << endl;
2055 gLog << "Update input file '" << filenameData
2056 << "' with the SC hadronness" << endl;
2057
2058
2059 //----------------------------------------------------
2060 // read in optimum parameter values for the supercuts
2061
2062 TFile inparam(parSCfile);
2063 MSupercuts scin;
2064 scin.Read("MSupercuts");
2065 inparam.Close();
2066
2067 gLog << "Parameter values for supercuts were read in from file '"
2068 << parSCfile << "'" << endl;
2069
2070 TArrayD supercutsPar;
2071 supercutsPar = scin.GetParameters();
2072
2073 TArrayD supercutsStep;
2074 supercutsStep = scin.GetStepsizes();
2075
2076 gLog << "Parameter values for supercuts : " << endl;
2077 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
2078 {
2079 gLog << supercutsPar[i] << ", ";
2080 }
2081 gLog << endl;
2082
2083 gLog << "Step values for supercuts : " << endl;
2084 for (Int_t i=0; i<supercutsStep.GetSize(); i++)
2085 {
2086 gLog << supercutsStep[i] << ", ";
2087 }
2088 gLog << endl;
2089
2090
2091 //----------------------------------------------------
2092 MTaskList tliston;
2093 MParList pliston;
2094
2095 // set the parameters of the supercuts
2096 MSupercuts supercuts;
2097 supercuts.SetParameters(supercutsPar);
2098 gLog << "parameter values for the supercuts used for updating the input file ' "
2099 << filenameData << "'" << endl;
2100 supercutsPar = supercuts.GetParameters();
2101 for (Int_t i=0; i<supercutsPar.GetSize(); i++)
2102 {
2103 gLog << supercutsPar[i] << ", ";
2104 }
2105 gLog << endl;
2106
2107
2108 // geometry is needed in MHHillas... classes
2109 MGeomCam *fGeom =
2110 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
2111
2112 //-------------------------------------------
2113 // create the tasks which should be executed
2114 //
2115
2116 MReadMarsFile read("Events", filenameData);
2117 read.DisableAutoScheme();
2118
2119 TString fHilName = "MHillas";
2120 TString fHilNameExt = "MHillasExt";
2121 TString fHilNameSrc = "MHillasSrc";
2122 TString fImgParName = "MNewImagePar";
2123
2124
2125 //.......................................................................
2126 // calculation of hadroness for the supercuts
2127 // (=0.25 if fullfilled, =0.75 otherwise)
2128
2129 TString hadSCName = "HadSC";
2130 MSupercutsCalc sccalc(fHilName, fHilNameSrc);
2131 sccalc.SetHadronnessName(hadSCName);
2132
2133
2134 //.......................................................................
2135
2136
2137 //MWriteRootFile write(outNameImage, "UPDATE");
2138 //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
2139
2140
2141 MWriteRootFile write(outNameImage, "RECREATE");
2142
2143 write.AddContainer("MRawRunHeader", "RunHeaders");
2144 write.AddContainer("MTime", "Events");
2145 write.AddContainer("MMcEvt", "Events");
2146 write.AddContainer("ThetaOrig", "Events");
2147 write.AddContainer("MSrcPosCam", "Events");
2148 write.AddContainer("MSigmabar", "Events");
2149 write.AddContainer("MHillas", "Events");
2150 write.AddContainer("MHillasExt", "Events");
2151 write.AddContainer("MHillasSrc", "Events");
2152 write.AddContainer("MNewImagePar", "Events");
2153
2154 write.AddContainer("HadRF", "Events");
2155 write.AddContainer(hadSCName, "Events");
2156
2157
2158 //-----------------------------------------------------------------
2159 // geometry is needed in MHHillas... classes
2160 MGeomCam *fGeom =
2161 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
2162
2163 Float_t maxhadronness = 0.40;
2164 Float_t maxalpha = 20.0;
2165 Float_t maxdist = 10.0;
2166
2167 MFSelFinal selfinalgh(fHilNameSrc);
2168 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2169 selfinalgh.SetHadronnessName(hadSCName);
2170 selfinalgh.SetName("SelFinalgh");
2171 MContinue contfinalgh(&selfinalgh);
2172 contfinalgh.SetName("ContSelFinalgh");
2173
2174 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
2175 fillhadsc.SetName("HhadSC");
2176
2177 MFSelFinal selfinal(fHilNameSrc);
2178 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2179 selfinal.SetHadronnessName(hadSCName);
2180 selfinal.SetName("SelFinal");
2181 MContinue contfinal(&selfinal);
2182 contfinal.SetName("ContSelFinal");
2183
2184 TString mh3name = "abs(Alpha)";
2185 MBinning binsalphaabs("Binning"+mh3name);
2186 binsalphaabs.SetEdges(50, -2.0, 98.0);
2187
2188 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
2189 alphaabs.SetName(mh3name);
2190
2191 TH1 &alphahist = alphaabs->GetHist();
2192
2193 MFillH alpha(&alphaabs);
2194 alpha.SetName("FillAlphaAbs");
2195
2196 MFillH hfill1("MHHillas", fHilName);
2197 hfill1.SetName("HHillas");
2198
2199 MFillH hfill2("MHStarMap", fHilName);
2200 hfill2.SetName("HStarMap");
2201
2202 MFillH hfill3("MHHillasExt", fHilNameSrc);
2203 hfill3.SetName("HHillasExt");
2204
2205 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2206 hfill4.SetName("HHillasSrc");
2207
2208 MFillH hfill5("MHNewImagePar", fImgParName);
2209 hfill5.SetName("HNewImagePar");
2210
2211 //*****************************
2212 // entries in MParList
2213
2214 pliston.AddToList(&tliston);
2215 InitBinnings(&pliston);
2216
2217 pliston.AddToList(&supercuts);
2218
2219 pliston.AddToList(&binsalphaabs);
2220 pliston.AddToList(&alphaabs);
2221
2222 //*****************************
2223 // entries in MTaskList
2224
2225 tliston.AddToList(&read);
2226
2227 tliston.AddToList(&sccalc);
2228 tliston.AddToList(&fillhadsc);
2229
2230 tliston.AddToList(&write);
2231 tliston.AddToList(&contfinalgh);
2232
2233 tliston.AddToList(&alpha);
2234 tliston.AddToList(&hfill1);
2235 tliston.AddToList(&hfill2);
2236 tliston.AddToList(&hfill3);
2237 tliston.AddToList(&hfill4);
2238 tliston.AddToList(&hfill5);
2239
2240 tliston.AddToList(&contfinal);
2241
2242 //*****************************
2243
2244 //-------------------------------------------
2245 // Execute event loop
2246 //
2247 MProgressBar bar;
2248 MEvtLoop evtloop;
2249 evtloop.SetParList(&pliston);
2250 evtloop.SetProgressBar(&bar);
2251
2252 Int_t maxevents = -1;
2253 if ( !evtloop.Eventloop(maxevents) )
2254 return;
2255
2256 tliston.PrintStatistics(0, kTRUE);
2257
2258
2259 //-------------------------------------------
2260 // Display the histograms
2261 //
2262 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2263
2264 pliston.FindObject("MHHillas")->DrawClone();
2265 pliston.FindObject("MHHillasExt")->DrawClone();
2266 pliston.FindObject("MHHillasSrc")->DrawClone();
2267 pliston.FindObject("MHNewImagePar")->DrawClone();
2268 pliston.FindObject("MHStarMap")->DrawClone();
2269
2270 //-------------------------------------------
2271 // fit alpha distribution to get the number of excess events and
2272 // calculate significance of gamma signal in the alpha plot
2273
2274 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2275 TH1 &alphaHist = absalpha->GetHist();
2276 alphaHist.SetXTitle("|alpha| [\\circ]");
2277 alphaHist.SetName("alpha-macro");
2278
2279 Double_t alphasig = 13.1;
2280 Double_t alphamin = 30.0;
2281 Double_t alphamax = 90.0;
2282 Int_t degree = 2;
2283 Double_t significance = -99.0;
2284 Bool_t drawpoly = kTRUE;
2285 Bool_t fitgauss = kTRUE;
2286 Bool_t print = kTRUE;
2287
2288 MHFindSignificance findsig;
2289 findsig.SetRebin(kTRUE);
2290 findsig.SetReduceDegree(kFALSE);
2291
2292 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2293 alphasig, drawpoly, fitgauss, print);
2294 significance = findsig.GetSignificance();
2295 Float_t alphasi = findsig.GetAlphasi();
2296
2297 gLog << "For file '" << filenameData << "' : " << endl;
2298 gLog << "Significance of gamma signal after supercuts : "
2299 << significance << " (for |alpha| < " << alphasi << " degrees)"
2300 << endl;
2301
2302 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2303
2304 //-------------------------------------------
2305
2306 DeleteBinnings(&pliston);
2307 }
2308
2309
2310 gLog << "Macro ONOFFAnalysis : End of Job B_SC_UP" << endl;
2311 gLog << "=======================================================" << endl;
2312 }
2313 //---------------------------------------------------------------------
2314
2315
2316
2317 //---------------------------------------------------------------------
2318 // Job C
2319 //======
2320
2321 // - read ON1 and MC1 data files
2322 // which should have been updated to contain the hadronnesses
2323 // for the method of Random Forest and for the SUPERCUTS
2324 // - produce Neyman-Pearson plots
2325
2326 if (JobC)
2327 {
2328 gLog << "=====================================================" << endl;
2329 gLog << "Macro ONOFFAnalysis : Start of Job C" << endl;
2330
2331 gLog << "" << endl;
2332 gLog << "Macro ONOFFAnalysis : JobC = "
2333 << (JobC ? "kTRUE" : "kFALSE") << endl;
2334
2335
2336
2337 TString ext("2.root");
2338 TString extout("2.root");
2339
2340 TString typeHadrons("OFF");
2341 TString typeGammas("MC");
2342
2343 //--------------------------------------------
2344 // name of input data file
2345 TString NameData = outPath;
2346 NameData += typeHadrons;
2347 TString inNameData(NameData);
2348 inNameData += ext;
2349 gLog << "inNameData = " << inNameData << endl;
2350
2351 // name of input MC file
2352 TString NameMC = outPath;
2353 NameMC += typeGammas;
2354 TString inNameMC(NameMC);
2355 inNameMC += ext;
2356 gLog << "inNameMC = " << inNameMC << endl;
2357
2358
2359 //--------------------------------------------
2360 // root files for the training events
2361
2362
2363
2364 TString NameGammasTrain = outPath;
2365 NameGammasTrain += "RF_gammas_Train_";
2366 NameGammasTrain += typeGammas;
2367 TString outNameGammasTrain(NameGammasTrain);
2368 outNameGammasTrain += extout;
2369
2370
2371 TString NameHadronsTrain = outPath;
2372 NameHadronsTrain += "RF_hadrons_Train_";
2373 NameHadronsTrain += typeHadrons;
2374 TString outNameHadronsTrain(NameHadronsTrain);
2375 outNameHadronsTrain += extout;
2376
2377
2378 //--------------------------------------------
2379 // root files for the test events
2380
2381 TString NameGammasTest = outPath;
2382 NameGammasTest += "RF_gammas_Test_";
2383 NameGammasTest += typeGammas;
2384 TString outNameGammasTest(NameGammasTest);
2385 outNameGammasTest += extout;
2386
2387 TString NameHadronsTest = outPath;
2388 NameHadronsTest += "RF_hadrons_Test_";
2389 NameHadronsTest += typeHadrons;
2390 TString outNameHadronsTest(NameHadronsTest);
2391 outNameHadronsTest += extout;
2392
2393 //--------------------------------------------------------------------
2394
2395 //TString filenameData(inNameData);
2396 //TString filenameMC(inNameMC);
2397
2398 //TString filenameData(outNameHadronsTrain);
2399 //TString filenameMC(outNameGammasTrain);
2400
2401 TString filenameData(outNameHadronsTest);
2402 TString filenameMC(outNameGammasTest);
2403
2404 gLog << "filenameData = " << filenameData << endl;
2405 gLog << "filenameMC = " << filenameMC << endl;
2406
2407 //-----------------------------------------------------------------
2408
2409 MTaskList tliston;
2410 MParList pliston;
2411
2412
2413 // geometry is needed in MHHillas... classes
2414 MGeomCam *fGeom =
2415 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
2416
2417 //-------------------------------------------
2418 // create the tasks which should be executed
2419 //
2420
2421 MReadMarsFile read("Events", filenameMC);
2422 read.AddFile(filenameData);
2423 read.DisableAutoScheme();
2424
2425
2426 //.......................................................................
2427 // names of hadronness containers
2428
2429 TString hadSCName = "HadSC";
2430 TString hadRFName = "HadRF";
2431
2432 //.......................................................................
2433
2434
2435 TString fHilName = "MHillas";
2436 TString fHilNameExt = "MHillasExt";
2437 TString fHilNameSrc = "MHillasSrc";
2438 TString fImgParName = "MNewImagePar";
2439
2440 Float_t maxhadronness = 0.40;
2441 Float_t maxalpha = 20.0;
2442 Float_t maxdist = 10.0;
2443
2444 MFSelFinal selfinalgh(fHilNameSrc);
2445 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2446 selfinalgh.SetHadronnessName(hadRFName);
2447 selfinalgh.SetName("SelFinalgh");
2448 MContinue contfinalgh(&selfinalgh);
2449 contfinalgh.SetName("ContSelFinalgh");
2450
2451
2452 //MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
2453 //fillhadsc.SetName("HhadSC");
2454 MFillH fillhadrf("hadRF[MHHadronness]", hadRFName);
2455 fillhadrf.SetName("HhadRF");
2456
2457 MFSelFinal selfinal(fHilNameSrc);
2458 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2459 selfinal.SetHadronnessName(hadRFName);
2460 selfinal.SetName("SelFinal");
2461 MContinue contfinal(&selfinal);
2462 contfinal.SetName("ContSelFinal");
2463
2464
2465 MFillH hfill1("MHHillas", fHilName);
2466 hfill1.SetName("HHillas");
2467
2468 MFillH hfill2("MHStarMap", fHilName);
2469 hfill2.SetName("HStarMap");
2470
2471 MFillH hfill3("MHHillasExt", fHilNameSrc);
2472 hfill3.SetName("HHillasExt");
2473
2474 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2475 hfill4.SetName("HHillasSrc");
2476
2477 MFillH hfill5("MHNewImagePar", fImgParName);
2478 hfill5.SetName("HNewImagePar");
2479
2480
2481 //*****************************
2482 // entries in MParList
2483
2484 pliston.AddToList(&tliston);
2485 InitBinnings(&pliston);
2486
2487
2488 //*****************************
2489 // entries in MTaskList
2490
2491 tliston.AddToList(&read);
2492
2493 //tliston.AddToList(&fillhadsc);
2494 tliston.AddToList(&fillhadrf);
2495
2496 tliston.AddToList(&contfinalgh);
2497 tliston.AddToList(&hfill1);
2498 tliston.AddToList(&hfill2);
2499 tliston.AddToList(&hfill3);
2500 tliston.AddToList(&hfill4);
2501 tliston.AddToList(&hfill5);
2502
2503 tliston.AddToList(&contfinal);
2504
2505 //*****************************
2506
2507 //-------------------------------------------
2508 // Execute event loop
2509 //
2510 MProgressBar bar;
2511 MEvtLoop evtloop;
2512 evtloop.SetParList(&pliston);
2513 evtloop.SetProgressBar(&bar);
2514
2515 Int_t maxevents = -1;
2516 //Int_t maxevents = 35000;
2517 if ( !evtloop.Eventloop(maxevents) )
2518 return;
2519
2520 tliston.PrintStatistics(0, kTRUE);
2521
2522
2523 //-------------------------------------------
2524 // Display the histograms
2525 //
2526
2527 //pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2528 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2529
2530 pliston.FindObject("MHHillas")->DrawClone();
2531 pliston.FindObject("MHHillasExt")->DrawClone();
2532 pliston.FindObject("MHHillasSrc")->DrawClone();
2533 pliston.FindObject("MHNewImagePar")->DrawClone();
2534 pliston.FindObject("MHStarMap")->DrawClone();
2535
2536 DeleteBinnings(&pliston);
2537
2538 gLog << "Macro ONOFFAnalysis : End of Job C" << endl;
2539 gLog << "===================================================" << endl;
2540 }
2541
2542
2543
2544 //---------------------------------------------------------------------
2545 // Job D
2546 //======
2547
2548 // - select g/h separation method XX
2549 // - read ON2 (or MC2) root file
2550 // - apply cuts in hadronness
2551 // - make plots
2552
2553
2554 if (JobD)
2555 {
2556 gLog << "=====================================================" << endl;
2557 gLog << "Macro ONOFFAnalysis : Start of Job D" << endl;
2558
2559 gLog << "" << endl;
2560 gLog << "Macro ONOFFAnalysis : JobD = "
2561 << (JobD ? "kTRUE" : "kFALSE") << endl;
2562
2563
2564 // type of data to be analysed
2565 TString typeData = "ON";
2566 //TString typeData = "OFF";
2567 //TString typeData = "MC";
2568 gLog << "typeData = " << typeData << endl;
2569
2570 TString ext = "3.root";
2571
2572
2573 //------------------------------
2574 // selection of g/h separation method
2575 // and definition of final selections
2576
2577 //TString XX("SC");
2578 TString XX("RF");
2579 TString fhadronnessName("Had");
2580 fhadronnessName += XX;
2581 gLog << "fhadronnessName = " << fhadronnessName << endl;
2582
2583 // maximum values of the hadronness, |ALPHA| and DIST
2584 Float_t maxhadronness = 0.233;
2585 Float_t maxalpha = 20.0;
2586 Float_t maxdist = 10.0;
2587 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2588 << maxhadronness << ", " << maxalpha << ", "
2589 << maxdist << endl;
2590
2591
2592 //------------------------------
2593 // name of data file to be analysed
2594 TString filenameData(outPath);
2595 filenameData += typeData;
2596 filenameData += ext;
2597 gLog << "filenameData = " << filenameData << endl;
2598
2599
2600
2601 //*************************************************************************
2602 //
2603 // Analyse the data
2604 //
2605
2606 MTaskList tliston;
2607 MParList pliston;
2608
2609 // geometry is needed in MHHillas... classes
2610 MGeomCam *fGeom =
2611 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
2612
2613
2614 TString fHilName = "MHillas";
2615 TString fHilNameExt = "MHillasExt";
2616 TString fHilNameSrc = "MHillasSrc";
2617 TString fImgParName = "MNewImagePar";
2618
2619 //-------------------------------------------
2620 // create the tasks which should be executed
2621 //
2622
2623 MReadMarsFile read("Events", filenameData);
2624 read.DisableAutoScheme();
2625
2626
2627 //-----------------------------------------------------------------
2628 // geometry is needed in MHHillas... classes
2629 MGeomCam *fGeom =
2630 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
2631
2632 MFSelFinal selfinalgh(fHilNameSrc);
2633 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
2634 selfinalgh.SetHadronnessName(fhadronnessName);
2635 selfinalgh.SetName("SelFinalgh");
2636 MContinue contfinalgh(&selfinalgh);
2637 contfinalgh.SetName("ContSelFinalgh");
2638
2639 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
2640 fillhadsc.SetName("HhadSC");
2641 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
2642 fillhadrf.SetName("HhadRF");
2643
2644 TString mh3name = "abs(Alpha)";
2645 MBinning binsalphaabs("Binning"+mh3name);
2646 binsalphaabs.SetEdges(50, -2.0, 98.0);
2647
2648 MH3 alphaabs("abs(MHillasSrc.fAlpha)");
2649 alphaabs.SetName(mh3name);
2650
2651 TH1 &alphahist = alphaabs->GetHist();
2652
2653 MFillH alpha(&alphaabs);
2654 alpha.SetName("FillAlphaAbs");
2655
2656 MFillH hfill1("MHHillas", fHilName);
2657 hfill1.SetName("HHillas");
2658
2659 MFillH hfill2("MHStarMap", fHilName);
2660 hfill2.SetName("HStarMap");
2661
2662 MFillH hfill3("MHHillasExt", fHilNameSrc);
2663 hfill3.SetName("HHillasExt");
2664
2665 MFillH hfill4("MHHillasSrc", fHilNameSrc);
2666 hfill4.SetName("HHillasSrc");
2667
2668 MFillH hfill5("MHNewImagePar", fImgParName);
2669 hfill5.SetName("HNewImagePar");
2670
2671 MFSelFinal selfinal(fHilNameSrc);
2672 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
2673 selfinal.SetHadronnessName(fhadronnessName);
2674 selfinal.SetName("SelFinal");
2675 MContinue contfinal(&selfinal);
2676 contfinal.SetName("ContSelFinal");
2677
2678
2679 //*****************************
2680 // entries in MParList
2681
2682 pliston.AddToList(&tliston);
2683 InitBinnings(&pliston);
2684 pliston.AddToList(&binsalphaabs);
2685 pliston.AddToList(&alphaabs);
2686
2687 //*****************************
2688 // entries in MTaskList
2689
2690 tliston.AddToList(&read);
2691
2692 tliston.AddToList(&contfinalgh);
2693
2694 tliston.AddToList(&fillhadsc);
2695 tliston.AddToList(&fillhadrf);
2696
2697 tliston.AddToList(&alpha);
2698 tliston.AddToList(&hfill1);
2699 tliston.AddToList(&hfill2);
2700 tliston.AddToList(&hfill3);
2701 tliston.AddToList(&hfill4);
2702 tliston.AddToList(&hfill5);
2703
2704 tliston.AddToList(&contfinal);
2705
2706 //*****************************
2707
2708 //-------------------------------------------
2709 // Execute event loop
2710 //
2711 MProgressBar bar;
2712 MEvtLoop evtloop;
2713 evtloop.SetParList(&pliston);
2714 evtloop.SetProgressBar(&bar);
2715
2716 Int_t maxevents = -1;
2717 //Int_t maxevents = 10000;
2718 if ( !evtloop.Eventloop(maxevents) )
2719 return;
2720
2721 tliston.PrintStatistics(0, kTRUE);
2722
2723
2724 //-------------------------------------------
2725 // Display the histograms
2726 //
2727
2728 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
2729 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
2730
2731 pliston.FindObject("MHHillas")->DrawClone();
2732 pliston.FindObject("MHHillasExt")->DrawClone();
2733 pliston.FindObject("MHHillasSrc")->DrawClone();
2734 pliston.FindObject("MHNewImagePar")->DrawClone();
2735 pliston.FindObject("MHStarMap")->DrawClone();
2736
2737
2738 //-------------------------------------------
2739
2740 // fit alpha distribution to get the number of excess events and
2741 // calculate significance of gamma signal in the alpha plot
2742
2743 MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
2744 TH1 &alphaHist = absalpha->GetHist();
2745 alphaHist.SetXTitle("|alpha| [\\circ]");
2746 alphaHist.SetName("alpha-JobD");
2747
2748 Double_t alphasig = 13.1;
2749 Double_t alphamin = 30.0;
2750 Double_t alphamax = 90.0;
2751 Int_t degree = 2;
2752 Double_t significance = -99.0;
2753 Bool_t drawpoly = kTRUE;
2754 Bool_t fitgauss = kTRUE;
2755 Bool_t print = kTRUE;
2756
2757 MHFindSignificance findsig;
2758 findsig.SetRebin(kTRUE);
2759 findsig.SetReduceDegree(kFALSE);
2760
2761 findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
2762 alphasig, drawpoly, fitgauss, print);
2763 significance = findsig.GetSignificance();
2764 Float_t alphasi = findsig.GetAlphasi();
2765
2766 gLog << "For file '" << filenameData << "' : " << endl;
2767 gLog << "Significance of gamma signal after supercuts : "
2768 << significance << " (for |alpha| < " << alphasi << " degrees)"
2769 << endl;
2770
2771 findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
2772
2773 //-------------------------------------------
2774
2775
2776 DeleteBinnings(&pliston);
2777
2778 gLog << "Macro ONOFFAnalysis : End of Job D" << endl;
2779 gLog << "=======================================================" << endl;
2780 }
2781 //---------------------------------------------------------------------
2782
2783
2784
2785
2786
2787 //---------------------------------------------------------------------
2788 // Job E_XX
2789 //=========
2790
2791 // - select g/h separation method XX
2792 // - read MC_XX2.root file
2793   // - calculate eff. collection area
2794 // - read ON_XX2.root file
2795 // - apply final cuts
2796 // - calculate flux
2797 // - write root file for ON data after final cuts (ON_XX3.root))
2798
2799
2800 if (JobE_XX)
2801 {
2802 gLog << "=====================================================" << endl;
2803 gLog << "Macro ONOFFAnalysis : Start of Job E_XX" << endl;
2804
2805 gLog << "" << endl;
2806 gLog << "Macro ONOFFAnalysis : JobE_XX, CCollArea, OEEst, WEX = "
2807 << (JobE_XX ? "kTRUE" : "kFALSE") << ", "
2808 << (CCollArea?"kTRUE" : "kFALSE") << ", "
2809 << (OEEst ? "kTRUE" : "kFALSE") << ", "
2810 << (WEX ? "kTRUE" : "kFALSE") << endl;
2811
2812
2813 // type of data to be analysed
2814 //TString typeData = "ON";
2815 //TString typeData = "OFF";
2816 TString typeData = "MC";
2817 gLog << "typeData = " << typeData << endl;
2818
2819 TString typeMC = "MC";
2820 TString ext = "3.root";
2821 TString extout = "4.root";
2822
2823 //------------------------------
2824 // selection of g/h separation method
2825 // and definition of final selections
2826
2827 //TString XX("SC");
2828 TString XX("RF");
2829 TString fhadronnessName("Had");
2830 fhadronnessName += XX;
2831 gLog << "fhadronnessName = " << fhadronnessName << endl;
2832
2833 // maximum values of the hadronness, |ALPHA| and DIST
2834 Float_t maxhadronness = 0.23;
2835 Float_t maxalpha = 20.0;
2836 Float_t maxdist = 10.0;
2837 gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
2838 << maxhadronness << ", " << maxalpha << ", "
2839 << maxdist << endl;
2840
2841 //------------------------------
2842 // name of MC file to be used for optimizing the energy estimator
2843 TString filenameOpt(outPath);
2844 filenameOpt += typeMC;
2845 filenameOpt += ext;
2846 gLog << "filenameOpt = " << filenameOpt << endl;
2847
2848 //------------------------------
2849 // name of file containing the parameters of the energy estimator
2850 TString energyParName(outPath);
2851 energyParName += "energyest_";
2852 energyParName += XX;
2853 energyParName += ".root";
2854 gLog << "energyParName = " << energyParName << endl;
2855
2856 //------------------------------
2857 // name of MC file to be used for calculating the eff. collection areas
2858 TString filenameArea(outPath);
2859 filenameArea += typeMC;
2860 filenameArea += ext;
2861 gLog << "filenameArea = " << filenameArea << endl;
2862
2863 //------------------------------
2864 // name of file containing the eff. collection areas
2865 TString collareaName(outPath);
2866 collareaName += "area_";
2867 collareaName += XX;
2868 collareaName += ".root";
2869 gLog << "collareaName = " << collareaName << endl;
2870
2871 //------------------------------
2872 // name of data file to be analysed
2873 TString filenameData(outPath);
2874 filenameData += typeData;
2875 filenameData += ext;
2876 gLog << "filenameData = " << filenameData << endl;
2877
2878 //------------------------------
2879 // name of output data file (after the final cuts)
2880 TString filenameDataout(outPath);
2881 filenameDataout += typeData;
2882 filenameDataout += "_";
2883 filenameDataout += XX;
2884 filenameDataout += extout;
2885 gLog << "filenameDataout = " << filenameDataout << endl;
2886
2887 //------------------------------
2888 // name of file containing histograms for flux calculastion
2889 TString filenameResults(outPath);
2890 filenameResults += typeData;
2891 filenameResults += "Results_";
2892 filenameResults += XX;
2893 filenameResults += extout;
2894 gLog << "filenameResults = " << filenameResults << endl;
2895
2896
2897 //====================================================================
2898
2899 MHMcCollectionArea collarea;
2900 collarea.SetEaxis(MHMcCollectionArea::kLinear);
2901
2902 MParList parlist;
2903 InitBinnings(&parlist);
2904
2905 if (CCollArea)
2906 {
2907 gLog << "-----------------------------------------------" << endl;
2908 gLog << "Start calculation of effective collection areas" << endl;
2909
2910
2911 MTaskList tasklist;
2912
2913 //---------------------------------------
2914 // Setup the tasks to be executed
2915 //
2916 MReadMarsFile reader("Events", filenameArea);
2917 reader.DisableAutoScheme();
2918
2919 MFSelFinal cuthadrons;
2920 cuthadrons.SetHadronnessName(fhadronnessName);
2921 cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
2922
2923 MContinue conthadrons(&cuthadrons);
2924
2925
2926 MFillH filler("MHMcCollectionArea", "MMcEvt");
2927 filler.SetName("CollectionArea");
2928
2929 //********************************
2930 // entries in MParList
2931
2932 parlist.AddToList(&tasklist);
2933
2934 parlist.AddToList(&collarea);
2935
2936 //********************************
2937 // entries in MTaskList
2938
2939 tasklist.AddToList(&reader);
2940 tasklist.AddToList(&conthadrons);
2941 tasklist.AddToList(&filler);
2942
2943 //********************************
2944
2945 //-----------------------------------------
2946 // Execute event loop
2947 //
2948 MEvtLoop magic;
2949 magic.SetParList(&parlist);
2950
2951 MProgressBar bar;
2952 magic.SetProgressBar(&bar);
2953 if (!magic.Eventloop())
2954 return;
2955
2956 tasklist.PrintStatistics(0, kTRUE);
2957
2958 // Calculate effective collection areas
2959 // and display the histograms
2960 //
2961 //MHMcCollectionArea *collarea =
2962 // (MHMcCollectionArea*)parlist.FindObject("MHMcCollectionArea");
2963 collarea.CalcEfficiency();
2964 collarea.DrawClone();
2965
2966
2967
2968 //---------------------------------------------
2969 // Write histograms to a file
2970 //
2971
2972 TFile f(collareaName, "RECREATE");
2973 //collarea.GetHist()->Write();
2974 //collarea.GetHAll()->Write();
2975 //collarea.GetHSel()->Write();
2976 collarea.Write();
2977
2978 f.Close();
2979
2980 gLog << "Collection area plots written onto file " << collareaName << endl;
2981
2982 gLog << "Calculation of effective collection areas done" << endl;
2983 gLog << "-----------------------------------------------" << endl;
2984 //------------------------------------------------------------------
2985 }
2986
2987 if (!CCollArea)
2988 {
2989 gLog << "-----------------------------------------------" << endl;
2990 gLog << "Read in effective collection areas from file "
2991 << collareaName << endl;
2992
2993 TFile collfile(collareaName);
2994 collfile.ls();
2995 collarea.Read("MHMcCollectionArea");
2996 collarea.DrawClone();
2997
2998 gLog << "Effective collection areas were read in from file "
2999 << collareaName << endl;
3000 gLog << "-----------------------------------------------" << endl;
3001 }
3002
3003
3004 // save binnings for call to MagicEEst
3005 MBinning *binsE = (MBinning*)parlist.FindObject("BinningE");
3006 if (!binsE)
3007 {
3008 gLog << "Object 'BinningE' not found in MParList" << endl;
3009 return;
3010 }
3011 MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
3012 if (!binsTheta)
3013 {
3014 gLog << "Object 'BinningTheta' not found in MParList" << endl;
3015 return;
3016 }
3017
3018 //-------------------------------------
3019 TString fHilName = "MHillas";
3020 TString fHilNameExt = "MHillasExt";
3021 TString fHilNameSrc = "MHillasSrc";
3022 TString fImgParName = "MNewImagePar";
3023
3024
3025 if (OEEst)
3026 {
3027 //===========================================================
3028 //
3029 // Optimization of energy estimator
3030 //
3031 gLog << "Macro ONOFFAnalysis.C : calling MagicEEst" << endl;
3032
3033 TString inpath("");
3034 TString outpath("");
3035 Int_t howMany = 2000;
3036 MagicEEst(inpath, filenameOpt, outpath, energyParName,
3037 fHilName, fHilNameSrc, fhadronnessName,
3038 howMany, maxhadronness, maxalpha, maxdist,
3039 binsE, binsTheta);
3040 gLog << "Macro ONOFFAnalysis.C : returning from MagicEEst" << endl;
3041 }
3042
3043 if (WEX)
3044 {
3045 //-----------------------------------------------------------
3046 //
3047 // Read in parameters of energy estimator ("MMcEnergyEst")
3048 // and migration matrix ("MHMcEnergyMigration")
3049 //
3050 gLog << "================================================================"
3051 << endl;
3052 gLog << "Macro ONOFFAnalysis.C : read parameters of energy estimator and migration matrix from file '"
3053 << energyParName << "'" << endl;
3054 TFile enparam(energyParName);
3055 enparam.ls();
3056 MMcEnergyEst mcest("MMcEnergyEst");
3057 mcest.Read("MMcEnergyEst");
3058
3059 //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
3060 gLog << "Parameters of energy estimator were read in" << endl;
3061
3062
3063 gLog << "Read in Migration matrix" << endl;
3064
3065 MHMcEnergyMigration mighiston("MHMcEnergyMigration");
3066 mighiston.Read("MHMcEnergyMigration");
3067 //MHMcEnergyMigration &mighiston =
3068 // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
3069
3070 gLog << "Migration matrix was read in" << endl;
3071
3072
3073 TArrayD parA(mcest.GetNumCoeffA());
3074 TArrayD parB(mcest.GetNumCoeffB());
3075 for (Int_t i=0; i<parA.GetSize(); i++)
3076 parA[i] = mcest.GetCoeff(i);
3077 for (Int_t i=0; i<parB.GetSize(); i++)
3078 parB[i] = mcest.GetCoeff( i+parA.GetSize() );
3079
3080 //*************************************************************************
3081 //
3082 // Analyse the data
3083 //
3084 gLog << "============================================================"
3085 << endl;
3086 gLog << "Analyse the data" << endl;
3087
3088 MTaskList tliston;
3089 MParList pliston;
3090
3091 // geometry is needed in MHHillas... classes
3092 MGeomCam *fGeom =
3093 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
3094
3095
3096 //-------------------------------------------
3097 // create the tasks which should be executed
3098 //
3099
3100 MReadMarsFile read("Events", filenameData);
3101 read.DisableAutoScheme();
3102
3103 //.......................................................................
3104
3105 gLog << "MagicAnalysis.C : write root file '" << filenameDataout
3106 << "'" << endl;
3107
3108 //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
3109
3110
3111 MWriteRootFile write(filenameDataout, "RECREATE");
3112
3113 write.AddContainer("MRawRunHeader", "RunHeaders");
3114 write.AddContainer("MTime", "Events");
3115 write.AddContainer("MMcEvt", "Events");
3116 write.AddContainer("ThetaOrig", "Events");
3117 write.AddContainer("MSrcPosCam", "Events");
3118 write.AddContainer("MSigmabar", "Events");
3119 write.AddContainer("MHillas", "Events");
3120 write.AddContainer("MHillasExt", "Events");
3121 write.AddContainer("MHillasSrc", "Events");
3122 write.AddContainer("MNewImagePar", "Events");
3123
3124 //write.AddContainer("HadNN", "Events");
3125 write.AddContainer("HadSC", "Events");
3126 write.AddContainer("HadRF", "Events");
3127
3128 write.AddContainer("MEnergyEst", "Events");
3129
3130
3131 //-----------------------------------------------------------------
3132 // geometry is needed in MHHillas... classes
3133 MGeomCam *fGeom =
3134 (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam");
3135
3136 MFSelFinal selfinalgh(fHilNameSrc);
3137 selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
3138 selfinalgh.SetHadronnessName(fhadronnessName);
3139 selfinalgh.SetName("SelFinalgh");
3140 MContinue contfinalgh(&selfinalgh);
3141 contfinalgh.SetName("ContSelFinalgh");
3142
3143 //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
3144 //fillhadnn.SetName("HhadNN");
3145 MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
3146 fillhadsc.SetName("HhadSC");
3147 MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
3148 fillhadrf.SetName("HhadRF");
3149
3150 //---------------------------
3151 // calculate estimated energy
3152
3153 MEnergyEstParam eeston(fHilName);
3154 eeston.Add(fHilNameSrc);
3155
3156 eeston.SetCoeffA(parA);
3157 eeston.SetCoeffB(parB);
3158
3159 //---------------------------
3160 // calculate estimated energy using Daniel's parameters
3161
3162 //MEnergyEstParamDanielMkn421 eeston(fHilName);
3163 //eeston.Add(fHilNameSrc);
3164 //eeston.SetCoeffA(parA);
3165 //eeston.SetCoeffB(parB);
3166
3167
3168 //---------------------------
3169
3170
3171 MFillH hfill1("MHHillas", fHilName);
3172 hfill1.SetName("HHillas");
3173
3174 MFillH hfill2("MHStarMap", fHilName);
3175 hfill2.SetName("HStarMap");
3176
3177 MFillH hfill3("MHHillasExt", fHilNameSrc);
3178 hfill3.SetName("HHillasExt");
3179
3180 MFillH hfill4("MHHillasSrc", fHilNameSrc);
3181 hfill4.SetName("HHillasSrc");
3182
3183 MFillH hfill5("MHNewImagePar", fImgParName);
3184 hfill5.SetName("HNewImagePar");
3185
3186 //---------------------------
3187 // new from Robert
3188
3189 MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
3190 hfill6.SetName("HTimeDiffTheta");
3191
3192 MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
3193 hfill6a.SetName("HTimeDiffTime");
3194
3195 MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
3196 hfill7.SetName("HAlphaEnergyTheta");
3197
3198 MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
3199 hfill7a.SetName("HAlphaEnergyTime");
3200
3201 MFillH hfill7b("MHThetabarTime", fHilNameSrc);
3202 hfill7b.SetName("HThetabarTime");
3203
3204 MFillH hfill7c("MHEnergyTime", "MMcEvt");
3205 hfill7c.SetName("HEnergyTime");
3206
3207
3208 //---------------------------
3209
3210 MFSelFinal selfinal(fHilNameSrc);
3211 selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
3212 selfinal.SetHadronnessName(fhadronnessName);
3213 selfinal.SetName("SelFinal");
3214 MContinue contfinal(&selfinal);
3215 contfinal.SetName("ContSelFinal");
3216
3217
3218 //*****************************
3219 // entries in MParList
3220
3221 pliston.AddToList(&tliston);
3222 InitBinnings(&pliston);
3223
3224
3225 //*****************************
3226 // entries in MTaskList
3227
3228 tliston.AddToList(&read);
3229
3230 // robert
3231 tliston.AddToList(&hfill6); //timediff
3232 tliston.AddToList(&hfill6a); //timediff
3233
3234 tliston.AddToList(&contfinalgh);
3235 tliston.AddToList(&eeston);
3236
3237 tliston.AddToList(&write);
3238
3239 //tliston.AddToList(&fillhadnn);
3240 tliston.AddToList(&fillhadsc);
3241 tliston.AddToList(&fillhadrf);
3242
3243 tliston.AddToList(&hfill1);
3244 tliston.AddToList(&hfill2);
3245 tliston.AddToList(&hfill3);
3246 tliston.AddToList(&hfill4);
3247 tliston.AddToList(&hfill5);
3248
3249 //robert
3250 tliston.AddToList(&hfill7);
3251 tliston.AddToList(&hfill7a);
3252 tliston.AddToList(&hfill7b);
3253 tliston.AddToList(&hfill7c);
3254
3255 tliston.AddToList(&contfinal);
3256
3257 //*****************************
3258
3259 //-------------------------------------------
3260 // Execute event loop
3261 //
3262 MProgressBar bar;
3263 MEvtLoop evtloop;
3264 evtloop.SetParList(&pliston);
3265 evtloop.SetProgressBar(&bar);
3266
3267 Int_t maxevents = -1;
3268 if ( !evtloop.Eventloop(maxevents) )
3269 return;
3270
3271 tliston.PrintStatistics(0, kTRUE);
3272
3273
3274 //-------------------------------------------
3275 // Display the histograms
3276 //
3277
3278 //pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
3279
3280 gLog << "before hadRF" << endl;
3281 pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
3282
3283 gLog << "before hadSC" << endl;
3284 pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
3285
3286 gLog << "before MHHillas" << endl;
3287 pliston.FindObject("MHHillas")->DrawClone();
3288
3289 gLog << "before MHHillasExt" << endl;
3290 pliston.FindObject("MHHillasExt")->DrawClone();
3291
3292 gLog << "before MHHillasSrc" << endl;
3293 pliston.FindObject("MHHillasSrc")->DrawClone();
3294
3295 gLog << "before MHNewImagePar" << endl;
3296 pliston.FindObject("MHNewImagePar")->DrawClone();
3297
3298 gLog << "before MHStarMap" << endl;
3299 pliston.FindObject("MHStarMap")->DrawClone();
3300
3301 gLog << "before DeleteBinnings" << endl;
3302
3303 DeleteBinnings(&pliston);
3304
3305 gLog << "before Robert's code" << endl;
3306
3307
3308//rwagner write all relevant histograms onto a file
3309
3310 if (WRobert)
3311 {
3312 gLog << "=======================================================" << endl;
3313 gLog << "Write results onto file '" << filenameResults << "'" << endl;
3314
3315 TFile outfile(filenameResults,"recreate");
3316
3317 MHHillasSrc* hillasSrc =
3318 (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
3319 TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
3320 alphaHist->Write();
3321 gLog << "Alpha plot has been written out" << endl;
3322
3323
3324 MHAlphaEnergyTheta* aetH =
3325 (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta"));
3326 TH3D* aetHist = (TH3D*)(aetH->GetHist());
3327 aetHist->SetName("aetHist");
3328 aetHist->Write();
3329 gLog << "AlphaEnergyTheta plot has been written out" << endl;
3330
3331 MHAlphaEnergyTime* aetH2 =
3332 (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime"));
3333 TH3D* aetHist2 = (TH3D*)(aetH2->GetHist());
3334 aetHist2->SetName("aetimeHist");
3335// aetHist2->DrawClone();
3336 aetHist2->Write();
3337 gLog << "AlphaEnergyTime plot has been written out" << endl;
3338
3339 MHThetabarTime* tbt =
3340 (MHThetabarTime*)(pliston->FindObject("MHThetabarTime"));
3341 TProfile* tbtHist = (TProfile*)(tbt->GetHist());
3342 tbtHist->SetName("tbtHist");
3343 tbtHist->Write();
3344 gLog << "ThetabarTime plot has been written out" << endl;
3345
3346 MHEnergyTime* ent =
3347 (MHEnergyTime*)(pliston->FindObject("MHEnergyTime"));
3348 TH2D* entHist = (TH2D*)(ent->GetHist());
3349 entHist->SetName("entHist");
3350 entHist->Write();
3351 gLog << "EnergyTime plot has been written out" << endl;
3352
3353 MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta");
3354 TH2D* timeHist = (TH2D*)(time->GetHist());
3355 timeHist->SetName("MHTimeDiffTheta");
3356 timeHist->SetTitle("Time diffs");
3357 timeHist->Write();
3358 gLog << "TimeDiffTheta plot has been written out" << endl;
3359
3360
3361 MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime");
3362 TH2D* timeHist2 = (TH2D*)(time2->GetHist());
3363 timeHist2->SetName("MHTimeDiffTime");
3364 timeHist2->SetTitle("Time diffs");
3365 timeHist2->Write();
3366 gLog << "TimeDiffTime plot has been written out" << endl;
3367
3368//rwagner write also collareas to same file
3369 collarea->GetHist()->Write();
3370 collarea->GetHAll()->Write();
3371 collarea->GetHSel()->Write();
3372 gLog << "Effective collection areas have been written out" << endl;
3373
3374//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data
3375//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file
3376
3377 gLog << "before closing outfile" << endl;
3378
3379 //outfile.Close();
3380 gLog << "Results were written onto file '" << filenameResults
3381 << "'" << endl;
3382 gLog << "=======================================================" << endl;
3383 }
3384
3385 }
3386
3387 gLog << "Macro ONOFFAnalysis : End of Job E_XX" << endl;
3388 gLog << "=======================================================" << endl;
3389 }
3390 //---------------------------------------------------------------------
3391
3392}
3393
3394
3395
3396
3397
3398
3399
3400
3401
3402
3403
3404
3405
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420
3421
Note: See TracBrowser for help on using the repository browser.