source: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C@ 2475

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