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

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