source: trunk/MagicSoft/Mars/macros/CT1Analysis.C@ 2408

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