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

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