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

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