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

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