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

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