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

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