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

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