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

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