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

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