source: trunk/MagicSoft/Mars/macros/ONAnalysis.C@ 5111

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