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

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