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

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