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

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