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

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