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

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