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

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