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

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