source: trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C@ 2167

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