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

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