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

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