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

Last change on this file since 2029 was 1963, checked in by wittek, 23 years ago
*** empty log message ***
File size: 66.9 KB
Line 
1
2#include "CT1EEst.C"
3
4void InitBinnings(MParList *plist)
5 {
6 cout << "InitBinnings" << endl;
7
8 MBinning *binssize = new MBinning("BinningSize");
9 binssize->SetEdgesLog(50, 10, 1.0e5);
10 plist->AddToList(binssize);
11
12 MBinning *binsdistc = new MBinning("BinningDist");
13 binsdistc->SetEdges(50, 0, 1.4);
14 plist->AddToList(binsdistc);
15
16 MBinning *binswidth = new MBinning("BinningWidth");
17 binswidth->SetEdges(50, 0, 1.0);
18 plist->AddToList(binswidth);
19
20 MBinning *binslength = new MBinning("BinningLength");
21 binslength->SetEdges(50, 0, 1.0);
22 plist->AddToList(binslength);
23
24 MBinning *binsalpha = new MBinning("BinningAlpha");
25 binsalpha->SetEdges(100, -100, 100);
26 plist->AddToList(binsalpha);
27
28 MBinning *binsasym = new MBinning("BinningAsym");
29 binsasym->SetEdges(50, -1.5, 1.5);
30 plist->AddToList(binsasym);
31
32 MBinning *binsht = new MBinning("BinningHeadTail");
33 binsht->SetEdges(50, -1.5, 1.5);
34 plist->AddToList(binsht);
35
36 MBinning *binsm3l = new MBinning("BinningM3Long");
37 binsm3l->SetEdges(50, -1.5, 1.5);
38 plist->AddToList(binsm3l);
39
40 MBinning *binsm3t = new MBinning("BinningM3Trans");
41 binsm3t->SetEdges(50, -1.5, 1.5);
42 plist->AddToList(binsm3t);
43
44
45 //.....
46 MBinning *binsb = new MBinning("BinningSigmabar");
47 binsb->SetEdges( 100, 0.0, 5.0);
48 plist->AddToList(binsb);
49
50 Int_t nbins = 8;
51 TArrayD xbins(nbins+1);
52 xbins[0] = 0.0;
53 xbins[1] = 7.5;
54 xbins[2] = 17.5;
55 xbins[3] = 23.5;
56 xbins[4] = 29.5;
57 xbins[5] = 35.5;
58 xbins[6] = 42.0;
59 xbins[7] = 55.0;
60 xbins[8] = 68.0;
61
62 MBinning *binth = new MBinning("BinningTheta");
63 binth->SetEdges(xbins);
64 plist->AddToList(binth);
65
66 MBinning *binsdiff = new MBinning("BinningDiffsigma2");
67 binsdiff->SetEdges(100, -5.0, 20.0);
68 plist->AddToList(binsdiff);
69 }
70
71void DeleteBinnings(MParList *plist)
72 {
73 cout << "DeleteBinnings" << endl;
74
75 TObject *bin;
76
77 bin = plist->FindObject("BinningSize");
78 if (!bin) delete bin;
79 bin = plist->FindObject("BinningDist");
80 if (!bin) delete bin;
81 bin = plist->FindObject("BinningWidth");
82 if (!bin) delete bin;
83 bin = plist->FindObject("BinningLength");
84 if (!bin) delete bin;
85
86 bin = plist->FindObject("BinningAlpha");
87 if (!bin) delete bin;
88 bin = plist->FindObject("BinningAsym");
89 if (!bin) delete bin;
90 bin = plist->FindObject("BinningHeadTail");
91 if (!bin) delete bin;
92 bin = plist->FindObject("BinningM3Long");
93 if (!bin) delete bin;
94 bin = plist->FindObject("BinningM3Trans");
95 if (!bin) delete bin;
96
97 //.....
98 bin = plist->FindObject("BinningSigmabar");
99 if (!bin) delete bin;
100 bin = plist->FindObject("BinningTheta");
101 if (!bin) delete bin;
102 bin = plist->FindObject("BinningDiffsigma2");
103 if (!bin) delete bin;
104 }
105
106void CT1Analysis()
107{
108 //-----------------------------------------------
109 const char *offfile = "~magican/ct1test/wittek/offdata.preproc";
110
111 //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc";
112 const char *onfile = "~magican/ct1test/wittek/mkn421_00-01";
113
114 const char *mcfile = "~magican/ct1test/wittek/mkn421_mc_pedrms_0.001.preproc";
115 //-----------------------------------------------
116
117 // path for input for Mars
118 TString inPath = "~magican/ct1test/wittek/marsoutput/optionB/";
119
120 // path for output from Mars
121 TString outPath = "~magican/ct1test/wittek/marsoutput/optionB/";
122
123
124
125 //************************************************************************
126
127 // Job A_ON : read ON data
128 // - generate sigmabar vs. Theta plot;
129 // - write root file for ON data (ON1.root);
130 // - generate hadron matrix for g/h separation
131
132 Bool_t JobA_ON = kFALSE;
133 Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ?
134 Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ?
135 Bool_t WON1 = kFALSE; // write out root file ON1.root ?
136
137
138 // Job A_MC : read MC gamma data,
139 // - read sigmabar vs. Theta plot from ON data
140 // - do padding;
141 // - write root file for MC gammas (MC1.root);
142 // - generate gamma matrix for g/h separation)
143
144 Bool_t JobA_MC = kFALSE;
145 Bool_t WLookMC = kFALSE; // write out look-alike events for gammas ?
146 Bool_t WMC1 = kFALSE; // write out root file MC1.root ?
147
148
149 // Job B_NN_UP : read ON1.root (or MC1.root) file
150 // - read hadron and gamma matrix
151 // - calculate hadroness for method of Nearest Neighbors
152 // and for the supercuts
153 // - update the input files with the hadronesses (ON2.root or MC2.root)
154
155 Bool_t JobB_NN_UP = kFALSE;
156 Bool_t WNN = kFALSE; // write root file ?
157
158
159 // Job B_EST_UP :
160 // - read MC2.root file
161 // - select g/h separation method XX
162 // - optimize energy estimation for events passing the final cuts
163 // - write parameters of energy estimator onto file
164 // - update ON2.root and MC2.root files with estimated energy
165 // (ON_XX2.root and MC_XX2.root)
166
167 Bool_t JobB_EST_UP = kFALSE;
168 Bool_t WESTUP = kFALSE; // update root file ?
169
170
171 // Job C_NN : read ON1.root and MC1.root files
172 // - read look-alike events for hadrons and gammas
173 // - calculate hadronness for method of NEAREST NEIGHBOR
174 // - produce Neyman-Pearson plot
175
176 Bool_t JobC_NN = kFALSE;
177
178
179 // Job C_RF : read ON1.root and MC1.root files
180 // - read look-alike events for hadrons and gammas
181 // - calculate hadronness for method of RANDOM FOREST
182 // - produce Neyman-Pearson plot
183
184 Bool_t JobC_RF = kFALSE;
185
186
187 // Job D_XX :
188 // - select g/h separation method XX
189 // - read MC_XX2.root file
190   // - calculate eff. collection area
191 // - read ON_XX2.root file
192 // - apply final cuts
193 // - calculate flux
194 // - write root file for ON data after final cuts (ON3.root))
195 Bool_t JobD_XX = kTRUE;
196 Bool_t WXXD = kFALSE; // write out root file ON3.root ?
197
198
199 //************************************************************************
200
201
202 //---------------------------------------------------------------------
203 // Job A_ON
204 //=========
205 // read ON data file
206
207 // - produce the 2D-histogram "sigmabar versus Theta"
208 // (SigmaTheta_ON.root) for ON data
209 // (to be used for the padding of the MC gamma data)
210
211 // - write a file of look-alike events (matrix_hadrons_ON.root)
212 // (to be used later together with the corresponding MC gamma file
213 // for the g/h separation in the analysis of the ON data)
214
215 // - write a file of ON events (ON1.root)
216 // (after the standard cuts, before the g/h separation)
217 // (to be used together with the corresponding MC gamma file (MC1.root)
218 // for the optimization of the g/h separation)
219
220
221 if (JobA_ON)
222 {
223 cout << "=====================================================" << endl;
224 cout << "Macro CT1Analysis : Start of Job A_ON" << endl;
225 cout << "" << endl;
226 cout << "Macro CT1Analysis : JobA_ON, WHistON, WLookON, WON1 = "
227 << JobA_ON << ", " << WHistON << ", " << WLookON << ", "
228 << WON1 << endl;
229
230 TString typeInput = "ON";
231 cout << "typeInput = " << typeInput << endl;
232
233 // name of input root file
234 TString filenamein(onfile);
235
236 // name of output root file
237 TString outNameImage = outPath;
238 outNameImage += typeInput;
239 outNameImage += "1.root";
240
241
242 //-----------------------------------------------------------
243 MTaskList tliston;
244 MParList pliston;
245 MFilterList flist;
246
247 char *sourceName = "MSrcPosCam";
248 MSrcPosCam source(sourceName);
249
250
251 //-------------------------------------------
252 // create the tasks which should be executed
253 //
254
255 MCT1ReadPreProc read(filenamein);
256
257 MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc",
258 "MCT1PointingCorrCalc");
259
260 MBlindPixelCalc blind;
261
262 MFCT1SelBasic selbasic;
263 MContinue contbasic(&selbasic);
264
265 MSigmabarCalc sigbarcalc;
266
267 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
268 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
269
270 MImgCleanStd clean;
271
272 // calculate also extended image parameters
273 TString fHilName = "MHillas";
274 MHillasExt hext(fHilName);
275
276 MHillasCalc hcalc(fHilName);
277
278 TString fHilSrcName = "MHillasSrc";
279 MHillasSrc hilsrc(fHilSrcName);
280
281 MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
282 csrc1.SetInput(fHilName);
283
284 TString fnewparName = "MNewImagePar";
285 MNewImagePar newpar(fnewparName);
286
287 MNewImageParCalc newimage(sourceName, fnewparName);
288 newimage.SetInput(fHilName);
289
290 MFCT1SelStandard selstandard(fHilName, fHilSrcName);
291 MContinue contstandard(&selstandard);
292
293 MFillH hfill1("MHHillas", fHilName);
294 hfill1.SetTitle("Task to plot the source independent Hillas parameters");
295
296 MFillH hfill2("MHStarMap", fHilName);
297 hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");
298
299 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
300
301 MFillH hfill3("MHHillasExt", fHilSrcName);
302 hfill3.SetTitle("Task to plot the extended Hillas parameters");
303
304 MFillH hfill4("MHHillasSrc", fHilSrcName);
305 hfill4.SetTitle("Task to plot the source dependent Hillas parameters");
306
307 MFillH hfill5("MHNewImagePar", fnewparName);
308 hfill5.SetTitle("Task to plot the new image parameters");
309
310 //+++++++++++++++++++++++++++++++++++++++++++++++++++
311 // prepare writing of look-alike events (for hadrons)
312
313 const char* mtxName = "MatrixHadrons";
314 Int_t howMany = 500000;
315 TString outName = outPath;
316 outName += "matrix_hadrons_";
317 outName += typeInput;
318 outName += ".root";
319
320
321 cout << "" << endl;
322 cout << "========================================================" << endl;
323 cout << "Matrix for (hadrons)" << endl;
324 cout << "input file = " << filenamein<< endl;
325 cout << "matrix name = " << mtxName << endl;
326 cout << "max. no.of look-alike events = " << howMany << endl;
327 cout << "name of output root file = " << outName << endl;
328 cout << "" << endl;
329
330
331 // matrix limitation for look-alike events (approximate number)
332 MFEventSelector selector("", "");
333 selector.SetNumSelectEvts(howMany);
334
335 //
336 // --- setup the matrix and add it to the parlist
337 //
338 MHMatrix *pmatrix = new MHMatrix(mtxName);
339 MHMatrix& matrix = *pmatrix;
340
341 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
342 matrix.AddColumn("MSigmabar.fSigmabar");
343 matrix.AddColumn("log10(MHillas.fSize)");
344 matrix.AddColumn("MHillasSrc.fDist");
345 matrix.AddColumn("MHillas.fWidth");
346 matrix.AddColumn("MHillas.fLength");
347 matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
348 matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
349 matrix.AddColumn("MHillas.fConc");
350 matrix.AddColumn("MNewImagePar.fLeakage1");
351
352 matrix.Print("cols");
353
354 MFillH fillmatg(mtxName);
355 fillmatg.SetFilter(&flist);
356 //+++++++++++++++++++++++++++++++++++++++++++++++++++
357
358 if (WON1)
359 {
360 MWriteRootFile write(outNameImage);
361
362 write.AddContainer("MRawRunHeader", "RunHeaders");
363 write.AddContainer("MTime", "Events");
364 write.AddContainer("MMcEvt", "Events");
365 write.AddContainer("MSrcPosCam", "Events");
366 write.AddContainer("MSigmabar", "Events");
367 write.AddContainer("MHillas", "Events");
368 write.AddContainer("MHillasSrc", "Events");
369 write.AddContainer("MNewImagePar", "Events");
370 }
371
372 //*****************************
373 // entries in MFilterList
374
375 flist.AddToList(&selector);
376
377
378 //*****************************
379 // entries in MParList
380
381 pliston.AddToList(&tliston);
382 InitBinnings(&pliston);
383
384 pliston.AddToList(&source);
385 pliston.AddToList(&sigbarcalc);
386
387 pliston.AddToList(&hext);
388 pliston.AddToList(&hilsrc);
389 pliston.AddToList(&newpar);
390 pliston.AddToList(&hhilext);
391
392 pliston.AddToList(pmatrix);
393
394
395 //*****************************
396 // entries in MTaskList
397
398 tliston.AddToList(&read);
399 tliston.AddToList(&pointcorr);
400 tliston.AddToList(&blind);
401
402 tliston.AddToList(&contbasic);
403 tliston.AddToList(&sigbarcalc);
404 tliston.AddToList(&fillsigtheta);
405 tliston.AddToList(&clean);
406
407 tliston.AddToList(&hcalc);
408 tliston.AddToList(&csrc1);
409 tliston.AddToList(&newimage);
410
411 tliston.AddToList(&flist);
412 tliston.AddToList(&fillmatg);
413
414 tliston.AddToList(&hfill1);
415 tliston.AddToList(&hfill2);
416 tliston.AddToList(&hfill3);
417 tliston.AddToList(&hfill4);
418 tliston.AddToList(&hfill5);
419
420 tliston.AddToList(&contstandard);
421 if (WON1)
422 tliston.AddToList(&write);
423
424 //*****************************
425
426 //-------------------------------------------
427 // Execute event loop
428 //
429 MProgressBar bar;
430 MEvtLoop evtloop;
431 evtloop.SetParList(&pliston);
432 evtloop.SetProgressBar(&bar);
433 //evtloop.Write();
434
435 //Int_t maxevents = 1000000000;
436 Int_t maxevents = 20000;
437 if ( !evtloop.Eventloop(maxevents) )
438 return;
439
440 tliston.PrintStatistics(0, kTRUE);
441 DeleteBinnings(&pliston);
442
443 //-------------------------------------------
444 // Display the histograms
445 //
446 TObject *c;
447 c = ((MHSigmaTheta*)pliston.FindObject("SigmaTheta", "MHSigmaTheta"))
448 ->DrawClone();
449
450 c = pliston.FindObject("MHHillas")->DrawClone();
451 c = pliston.FindObject("MHHillasExt")->DrawClone();
452 c = pliston.FindObject("MHHillasSrc")->DrawClone();
453 c = pliston.FindObject("MHNewImagePar")->DrawClone();
454 c = pliston.FindObject("MHStarMap")->DrawClone();
455
456
457 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
458 //
459 // ----------------------------------------------------------
460 // Definition of the reference sample of look-alike events
461 // (this is a subsample of the original sample)
462 // ----------------------------------------------------------
463 //
464 cout << "" << endl;
465 cout << "========================================================" << endl;
466 // Select a maximum of nmaxevts events from the sample of look-alike
467 // events. They will form the reference sample.
468 Int_t nmaxevts = 2000;
469
470 // target distribution for the variable in column refcolumn (start at 0);
471 Int_t refcolumn = 0;
472 Int_t nbins = 5;
473 Float_t frombin = 0.5;
474 Float_t tobin = 1.0;
475 TH1F *thsh = new TH1F("thsh","target distribution",
476 nbins, frombin, tobin);
477 thsh->SetXTitle("cos( \\Theta)");
478 thsh->SetYTitle("Counts");
479 Float_t dbin = (tobin-frombin)/nbins;
480 Float_t lbin = frombin +dbin*0.5;
481 for (Int_t j=1; j<=nbins; j++)
482 {
483 thsh->Fill(lbin,1.0);
484 lbin += dbin;
485 }
486
487 cout << "" << endl;
488 cout << "========================================================" << endl;
489 cout << "Macro CT1Analysis : define reference sample for hadrons" << endl;
490 cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
491 << refcolumn << ", " << nmaxevts << endl;
492
493 matrix.EnableGraphicalOutput();
494 Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);
495
496 if ( !matrixok )
497 {
498 cout << "Macro CT1Analysis : Reference matrix for hadrons cannot be defined" << endl;
499 return;
500 }
501 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
502
503 //-------------------------------------------
504 // write out look-alike events (for hadrons)
505 //
506 if (WLookON)
507 {
508 TFile *writejens = new TFile(outName, "RECREATE", "");
509 matrix.Write();
510
511 cout << "" << endl;
512 cout << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
513 << outName << endl;
514
515 writejens->Close();
516 delete writejens;
517 }
518
519
520 //-------------------------------------------
521 // Write histograms onto a file
522 if (WHistON)
523 {
524 MHSigmaTheta *sigtheta =
525 (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
526 if (!sigtheta)
527 {
528 cout << "Object 'SigmaTheta' not found" << endl;
529 return;
530 }
531 TH2D *fHSigTh = sigtheta->GetSigmaTheta();
532 TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
533 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
534
535 TString outNameSigTh = outPath;
536 outNameSigTh += "SigmaTheta_";
537 outNameSigTh += typeInput;
538 outNameSigTh += ".root";
539
540 TFile outfile(outNameSigTh, "RECREATE");
541 fHSigTh->Write();
542 fHSigPixTh->Write();
543 fHDifPixTh->Write();
544 outfile.Close();
545
546 cout << "" << endl;
547 cout << "File " << outNameSigTh << " was written out" << endl;
548 }
549
550
551 cout << "Macro CT1Analysis : End of Job A_ON" << endl;
552 cout << "===================================================" << endl;
553 }
554
555
556 //---------------------------------------------------------------------
557 // Job A_MC
558 //=========
559
560 // read MC gamma data
561
562 // - to pad them
563 // (using the 2D-histogram "sigmabar versus Theta"
564 // (SigmaTheta_ON.root) of the ON data)
565
566 // - to write a file of look-alike events (matrix_gammas_MC.root)
567 // (to be used later together with the corresponding hadron file
568 // for the g/h separation in the analysis of the ON data)
569
570 // - to write a file of padded MC gamma events (MC1.root)
571 // (after the standard cuts, before the g/h separation)
572 // (to be used together with the corresponding hadron file
573 // for the optimization of the g/h separation)
574
575
576 if (JobA_MC)
577 {
578 cout << "=====================================================" << endl;
579 cout << "Macro CT1Analysis : Start of Job A_MC" << endl;
580
581 cout << "" << endl;
582 cout << "Macro CT1Analysis : JobA_MC, WLookMC, WMC1 = "
583 << JobA_MC << ", " << WLookMC << ", " << WMC1 << endl;
584
585 TString typeInput = "MC";
586 cout << "typeInput = " << typeInput << endl;
587
588
589 // name of input root file
590 TString filenamein(mcfile);
591
592 // name of output root file
593 TString outNameImage = outPath;
594 outNameImage += typeInput;
595 outNameImage += "1.root";
596
597
598 // use for padding sigmabar vs. Theta from ON data
599 TString typeHist = "ON";
600 cout << "typeHist = " << typeHist << endl;
601
602 //------------------------------------
603 // Get the histograms "2D-ThetaSigmabar"
604 // and "3D-ThetaPixSigma"
605 // and "3D-ThetaPixDiff"
606
607
608 TString outNameSigTh = outPath;
609 outNameSigTh += "SigmaTheta_";
610 outNameSigTh += typeHist;
611 outNameSigTh += ".root";
612
613
614 cout << "Reading in file " << outNameSigTh << endl;
615
616 TFile *infile = new TFile(outNameSigTh);
617 infile->ls();
618
619 TH2D *fHSigmaTheta =
620 (TH2D*) gROOT.FindObject("2D-ThetaSigmabar");
621 if (!fHSigmaTheta)
622 {
623 cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
624 return;
625 }
626 cout << "Object '2D-ThetaSigmabar' was read in" << endl;
627
628 TH3D *fHSigmaPixTheta =
629 (TH3D*) gROOT.FindObject("3D-ThetaPixSigma");
630 if (!fHSigmaPixTheta)
631 {
632 cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
633 return;
634 }
635 cout << "Object '3D-ThetaPixSigma' was read in" << endl;
636
637 TH3D *fHDiffPixTheta =
638 (TH3D*) gROOT.FindObject("3D-ThetaPixDiff");
639 if (!fHSigmaPixTheta)
640 {
641 cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
642 return;
643 }
644 cout << "Object '3D-ThetaPixDiff' was read in" << endl;
645
646 //------------------------------------
647
648 MTaskList tlist;
649 MParList plist;
650 MFilterList flist;
651
652 char *sourceName = "MSrcPosCam";
653 MSrcPosCam source(sourceName);
654
655 //-------------------------------------------
656 // create the tasks which should be executed
657 //
658
659 MCT1ReadPreProc read(filenamein);
660
661 MBlindPixelCalc blind;
662
663 MFCT1SelBasic selbasic;
664 MContinue contbasic(&selbasic);
665
666
667 // There are 2 options for Thomas Schweizer's padding
668 // fPadFlag = 1 get Sigmabar from fHSigmaTheta
669 // and Sigma from fHDiffPixTheta
670 // fPadFlag = 2 get Sigma from fHSigmaPixTheta
671
672 MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)",
673 fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta);
674 padthomas.SetPadFlag(1);
675
676 //...........................................
677
678 MSigmabarCalc sigbarcalc;
679
680 MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt");
681 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
682
683 MImgCleanStd clean;
684
685 // calculate also extended image parameters
686 TString fHilName = "MHillas";
687 MHillasExt hext(fHilName);
688
689 MHillasCalc hcalc(fHilName);
690
691 TString fHilSrcName = "MHillasSrc";
692 MHillasSrc hilsrc(fHilSrcName);
693
694 MHillasSrcCalc csrc1(sourceName, fHilSrcName, "MHillas", "");
695 csrc1.SetInput(fHilName);
696
697 TString fnewparName = "MNewImagePar";
698 MNewImagePar newpar(fnewparName);
699
700 MNewImageParCalc newimage(sourceName, fnewparName);
701 newimage.SetInput(fHilName);
702
703 MFCT1SelStandard selstandard(fHilName, fHilSrcName);
704 MContinue contstandard(&selstandard);
705
706 MFillH hfill1("MHHillas", fHilName);
707 hfill1.SetTitle("Task to plot the source independent Hillas parameters");
708
709 MFillH hfill2("MHStarMap", fHilName);
710 hfill2.SetTitle("Task to plot the star map (points along main axis of shower)");
711
712 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
713
714 MFillH hfill3("MHHillasExt", fHilSrcName);
715 hfill3.SetTitle("Task to plot the extended Hillas parameters");
716
717 MFillH hfill4("MHHillasSrc", fHilSrcName);
718 hfill4.SetTitle("Task to plot the source dependent Hillas parameters");
719
720 MFillH hfill5("MHNewImagePar", fnewparName);
721 hfill5.SetTitle("Task to plot the new image parameters");
722
723
724 //+++++++++++++++++++++++++++++++++++++++++++++++++++
725 // prepare writing of look-alike events (for gammas)
726
727 const char* mtxName = "MatrixGammas";
728 Int_t howMany = 500000;
729
730 TString outName = outPath;
731 outName += "matrix_gammas_";
732 outName += typeInput;
733 outName += ".root";
734
735
736 cout << "" << endl;
737 cout << "========================================================" << endl;
738 cout << "Matrix for (gammas)" << endl;
739 cout << "input file = " << filenamein<< endl;
740 cout << "matrix name = " << mtxName << endl;
741 cout << "max. no.of look-alike events = " << howMany << endl;
742 cout << "name of output root file = " << outName << endl;
743 cout << "" << endl;
744
745
746 // matrix limitation for look-alike events (approximate number)
747 MFEventSelector selector("", "");
748 selector.SetNumSelectEvts(howMany);
749
750 //
751 // --- setup the matrix and add it to the parlist
752 //
753 MHMatrix *pmatrix = new MHMatrix(mtxName);
754 MHMatrix& matrix = *pmatrix;
755
756 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
757 matrix.AddColumn("MSigmabar.fSigmabar");
758 matrix.AddColumn("log10(MHillas.fSize)");
759 matrix.AddColumn("MHillasSrc.fDist");
760 matrix.AddColumn("MHillas.fWidth");
761 matrix.AddColumn("MHillas.fLength");
762 matrix.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
763 matrix.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillas.fM3Long)");
764 matrix.AddColumn("MHillas.fConc");
765 matrix.AddColumn("MNewImagePar.fLeakage1");
766
767 matrix.Print("cols");
768
769 MFillH fillmatg(mtxName);
770 fillmatg.SetFilter(&flist);
771 //+++++++++++++++++++++++++++++++++++++++++++++++++++
772
773
774 if (WMC1)
775 {
776 MWriteRootFile write(outNameImage);
777
778 write.AddContainer("MRawRunHeader", "RunHeaders");
779 write.AddContainer("MTime", "Events");
780 write.AddContainer("MMcEvt", "Events");
781 write.AddContainer("MSrcPosCam", "Events");
782 write.AddContainer("MSigmabar", "Events");
783 write.AddContainer("MHillas", "Events");
784 write.AddContainer("MHillasSrc", "Events");
785 write.AddContainer("MNewImagePar", "Events");
786 }
787
788 //*****************************
789 // entries in MFilterList
790
791 flist.AddToList(&selector);
792
793
794 //*****************************
795 // entries in MParList
796
797 plist.AddToList(&tlist);
798 InitBinnings(&plist);
799
800 plist.AddToList(&source);
801 plist.AddToList(&sigbarcalc);
802
803 plist.AddToList(&hext);
804 plist.AddToList(&hilsrc);
805 plist.AddToList(&newpar);
806 plist.AddToList(&hhilext);
807
808 plist.AddToList(pmatrix);
809
810
811 //*****************************
812 // entries in MTaskList
813
814 tlist.AddToList(&read);
815 tlist.AddToList(&blind);
816 tlist.AddToList(&padthomas);
817
818 tlist.AddToList(&contbasic);
819 tlist.AddToList(&sigbarcalc);
820 tlist.AddToList(&fillsigtheta);
821 tlist.AddToList(&clean);
822
823 tlist.AddToList(&hcalc);
824 tlist.AddToList(&csrc1);
825 tlist.AddToList(&newimage);
826
827 tlist.AddToList(&flist);
828 tlist.AddToList(&fillmatg);
829
830 tlist.AddToList(&hfill1);
831 tlist.AddToList(&hfill2);
832 tlist.AddToList(&hfill3);
833 tlist.AddToList(&hfill4);
834 tlist.AddToList(&hfill5);
835
836 tlist.AddToList(&contstandard);
837 if (WMC1)
838 tlist.AddToList(&write);
839
840 //*****************************
841
842
843 //-------------------------------------------
844 // Execute event loop
845 //
846 MProgressBar bar;
847 MEvtLoop evtloop;
848 evtloop.SetParList(&plist);
849 evtloop.SetProgressBar(&bar);
850 //evtloop.Write();
851
852 //Int_t maxevents = 1000000000;
853 Int_t maxevents = 1000;
854 if ( !evtloop.Eventloop(maxevents) )
855 return;
856
857 tlist.PrintStatistics(0, kTRUE);
858 DeleteBinnings(&plist);
859
860 //-------------------------------------------
861 // Display the histograms
862 //
863 TObject *c;
864 c = ((MHSigmaTheta*)plist.FindObject("MCSigmaTheta", "MHSigmaTheta"))
865 ->DrawClone();
866 c = plist.FindObject("MHHillas")->DrawClone();
867 c = plist.FindObject("MHHillasExt")->DrawClone();
868 c = plist.FindObject("MHHillasSrc")->DrawClone();
869 c = plist.FindObject("MHNewImagePar")->DrawClone();
870 c = plist.FindObject("MHStarMap")->DrawClone();
871
872
873 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
874 //
875 // ----------------------------------------------------------
876 // Definition of the reference sample of look-alike events
877 // (this is a subsample of the original sample)
878 // ----------------------------------------------------------
879 //
880 cout << "" << endl;
881 cout << "========================================================" << endl;
882 // Select a maximum of nmaxevts events from the sample of look-alike
883 // events. They will form the reference sample.
884 Int_t nmaxevts = 2000;
885
886 // target distribution for the variable in column refcolumn (start at 0);
887 Int_t refcolumn = 0;
888 Int_t nbins = 5;
889 Float_t frombin = 0.5;
890 Float_t tobin = 1.0;
891 TH1F *thsh = new TH1F("thsh","target distribution",
892 nbins, frombin, tobin);
893 Float_t dbin = (tobin-frombin)/nbins;
894 Float_t lbin = frombin +dbin*0.5;
895 for (Int_t j=1; j<=nbins; j++)
896 {
897 thsh->Fill(lbin,1.0);
898 lbin += dbin;
899 }
900
901 cout << "" << endl;
902 cout << "========================================================" << endl;
903 cout << "Macro CT1Analysis : define reference sample for gammas" << endl;
904 cout << "Macro CT1Analysis : Parameters for reference sample : refcolumn, nmaxevts = "
905 << refcolumn << ", " << nmaxevts << endl;
906
907 matrix.EnableGraphicalOutput();
908 Bool_t matrixok = matrix.DefRefMatrix(refcolumn, *thsh, nmaxevts);
909 if ( !matrixok )
910 {
911 cout << "Macro CT1Analysis : Reference matrix for gammas cannot be defined" << endl;
912 return;
913 }
914 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
915
916 //-------------------------------------------
917 // write out look-alike events (for gammas)
918 //
919 if (WLookMC)
920 {
921 TFile *writejens = new TFile(outName, "RECREATE", "");
922 matrix.Write();
923 cout << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
924 << outName << endl;
925
926 writejens->Close();
927 delete writejens;
928 }
929
930 cout << "Macro CT1Analysis : End of Job A_MC"
931 << endl;
932 cout << "========================================================="
933 << endl;
934 }
935
936
937
938 //---------------------------------------------------------------------
939 // Job B_NN_UP
940 //============
941
942 // - read in the matrices of look-alike events for gammas and hadrons
943
944 // then read ON1.root (or MC1.root) file
945 //
946 // - calculate the hadroness for the method of nearest neighbors
947 //
948 // - update input root file, including the hadroness
949
950
951 if (JobB_NN_UP)
952 {
953 cout << "=====================================================" << endl;
954 cout << "Macro CT1Analysis : Start of Job B_NN_UP" << endl;
955
956 cout << "" << endl;
957 cout << "Macro CT1Analysis : JobB_NN_UP, WNN = "
958 << JobB_NN_UP << ", " << WNN << endl;
959
960
961 //TString typeInput = "ON";
962 TString typeInput = "MC";
963 cout << "typeInput = " << typeInput << endl;
964
965 // name of input root file
966 TString filenameData = outPath;
967 filenameData += typeInput;
968 filenameData += "1.root";
969 cout << "filenameData = " << filenameData << endl;
970
971 // name of output root file
972 TString outNameImage = outPath;
973 outNameImage += typeInput;
974 outNameImage += "2.root";
975
976 //$$$$$$$$$$$$$$$$$$$$$$$$$$
977 outNameImage = filenameData;
978 //$$$$$$$$$$$$$$$$$$$$$$$$$$
979
980 cout << "outNameImage = " << outNameImage << endl;
981
982 TString typeMatrixHadrons = "ON";
983 cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
984
985 //*************************************************************************
986 // read in matrices of look-alike events
987
988 const char* mtxName = "MatrixGammas";
989
990 TString outName = outPath;
991 outName += "matrix_gammas_";
992 outName += "MC";
993 outName += ".root";
994
995
996 cout << "" << endl;
997 cout << "========================================================" << endl;
998 cout << "Get matrix for (gammas)" << endl;
999 cout << "matrix name = " << mtxName << endl;
1000 cout << "name of root file = " << outName << endl;
1001 cout << "" << endl;
1002
1003
1004 // read in the object with the name 'mtxName' from file 'outName'
1005 //
1006 TFile *fileg = new TFile(outName);
1007
1008 MHMatrix matrixg;
1009 matrixg.Read(mtxName);
1010 matrixg.Print("cols");
1011
1012 delete fileg;
1013
1014 //*****************************************************************
1015
1016 const char* mtxName = "MatrixHadrons";
1017
1018 TString outName = outPath;
1019 outName += "matrix_hadrons_";
1020 outName += typeMatrixHadrons;
1021 outName += ".root";
1022
1023
1024 cout << "" << endl;
1025 cout << "========================================================" << endl;
1026 cout << " Get matrix for (hadrons)" << endl;
1027 cout << "matrix name = " << mtxName << endl;
1028 cout << "name of root file = " << outName << endl;
1029 cout << "" << endl;
1030
1031
1032 // read in the object with the name mtxName
1033 //
1034 TFile *fileh = new TFile(outName);
1035
1036 MHMatrix matrixh;
1037 matrixh.Read(mtxName);
1038 matrixh.Print("cols");
1039
1040 delete fileh;
1041
1042 //*************************************************************************
1043
1044
1045 MTaskList tliston;
1046 MParList pliston;
1047
1048 //-------------------------------------------
1049 // create the tasks which should be executed
1050 //
1051
1052 MReadMarsFile read("Events", filenameData);
1053 read.DisableAutoScheme();
1054
1055 TString fHilName = "MHillas";
1056 TString fHilSrcName = "MHillasSrc";
1057 TString fnewparName = "MNewImagePar";
1058
1059
1060 //.......................................................................
1061 // calculation of hadroness for method of Nearest Neighbors
1062
1063 TString hadNNName = "HadNN";
1064 MMultiDimDistCalc nncalc;
1065 nncalc.SetUseNumRows(25);
1066 nncalc.SetUseKernelMethod(kFALSE);
1067 nncalc.SetHadronnessName(hadNNName);
1068
1069 //.......................................................................
1070 // calculation of hadroness for the supercuts
1071 // (=0.25 if fullfilled, =0.75 otherwise)
1072
1073 TString hadSCName = "HadSC";
1074 MCT1SupercutsCalc sccalc(fHilName, fHilSrcName);
1075 sccalc.SetHadronnessName(hadSCName);
1076
1077 //.......................................................................
1078
1079 if (WNN)
1080 {
1081 //MWriteRootFile write(outNameImage);
1082 MWriteRootFile write(outNameImage, "UPDATE");
1083
1084 //write.AddContainer("MRawRunHeader", "RunHeaders");
1085 //write.AddContainer("MTime", "Events");
1086 //write.AddContainer("MMcEvt", "Events");
1087 //write.AddContainer("MSrcPosCam", "Events");
1088 //write.AddContainer("MSigmabar", "Events");
1089 //write.AddContainer("MHillas", "Events");
1090 //write.AddContainer("MHillasSrc", "Events");
1091 //write.AddContainer("MNewImagePar", "Events");
1092
1093 write.AddContainer(hadNNName, "Events");
1094 write.AddContainer(hadSCName, "Events");
1095 }
1096
1097 //-----------------------------------------------------------------
1098 // geometry is needed in MHHillas... classes
1099 MGeomCam *fGeom =
1100 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1101
1102
1103 Float_t hadcut = 0.40;
1104 Float_t alphacut = 100.0;
1105 MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
1106 selfinalgh.SetCuts(hadcut, alphacut);
1107 selfinalgh.SetHadronnessName(hadNNName);
1108 MContinue contfinalgh(&selfinalgh);
1109
1110 MFillH fillhadnn("hadNN[MHHadronness]", hadNNName);
1111 MFillH fillhadsc("hadSC[MHHadronness]", hadSCName);
1112
1113 Float_t hadcut = 0.40;
1114 Float_t alphacut = 20.0;
1115 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
1116 selfinal.SetCuts(hadcut, alphacut);
1117 selfinal.SetHadronnessName(hadNNName);
1118 MContinue contfinal(&selfinal);
1119
1120
1121 MFillH hfill1("MHHillas", fHilName);
1122 MFillH hfill2("MHStarMap", fHilName);
1123
1124 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
1125
1126 MFillH hfill3("MHHillasExt", fHilSrcName);
1127 MFillH hfill4("MHHillasSrc", fHilSrcName);
1128 MFillH hfill5("MHNewImagePar", fnewparName);
1129
1130
1131 //*****************************
1132 // entries in MParList
1133
1134 pliston.AddToList(&tliston);
1135 InitBinnings(&pliston);
1136
1137 pliston.AddToList(&matrixg);
1138 pliston.AddToList(&matrixh);
1139
1140 pliston.AddToList(&hhilext);
1141
1142 //*****************************
1143 // entries in MTaskList
1144
1145 tliston.AddToList(&read);
1146
1147 tliston.AddToList(&nncalc);
1148 tliston.AddToList(&sccalc);
1149 tliston.AddToList(&fillhadnn);
1150 tliston.AddToList(&fillhadsc);
1151
1152 if (WNN)
1153 tliston.AddToList(&write);
1154
1155 tliston.AddToList(&contfinalgh);
1156 tliston.AddToList(&hfill1);
1157 tliston.AddToList(&hfill2);
1158 tliston.AddToList(&hfill3);
1159 tliston.AddToList(&hfill4);
1160 tliston.AddToList(&hfill5);
1161
1162 tliston.AddToList(&contfinal);
1163
1164 //*****************************
1165
1166 //-------------------------------------------
1167 // Execute event loop
1168 //
1169 MProgressBar bar;
1170 MEvtLoop evtloop;
1171 evtloop.SetParList(&pliston);
1172 evtloop.SetProgressBar(&bar);
1173 //evtloop.Write();
1174
1175 if ( !evtloop.Eventloop() )
1176 //if ( !evtloop.Eventloop(1000) )
1177 return;
1178
1179 tliston.PrintStatistics(0, kTRUE);
1180 DeleteBinnings(&pliston);
1181
1182 //-------------------------------------------
1183 // Display the histograms
1184 //
1185 TObject *c;
1186
1187 c = pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
1188 c = pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
1189
1190 c = pliston.FindObject("MHHillas")->DrawClone();
1191 c = pliston.FindObject("MHHillasExt")->DrawClone();
1192 c = pliston.FindObject("MHHillasSrc")->DrawClone();
1193 c = pliston.FindObject("MHNewImagePar")->DrawClone();
1194 c = pliston.FindObject("MHStarMap")->DrawClone();
1195
1196
1197 cout << "Macro CT1Analysis : End of Job B_NN_UP" << endl;
1198 cout << "=======================================================" << endl;
1199 }
1200 //---------------------------------------------------------------------
1201
1202
1203
1204 //---------------------------------------------------------------------
1205 // Job B_EST_UP
1206 //============
1207
1208 // - read MC2.root file
1209 // - select g/h separation method XX
1210 // - optimize energy estimator for events passing final cuts
1211 // - write parameters of energy estimator onto file "energyest_XX.root"
1212 //
1213 // - read ON2.root and MC2.root files
1214 // - update input root file with the estimated energy
1215 // (ON_XX2.root, MC_XX2.root)
1216
1217
1218 if (JobB_EST_UP)
1219 {
1220 cout << "=====================================================" << endl;
1221 cout << "Macro CT1Analysis : Start of Job B_EST_UP" << endl;
1222
1223 cout << "" << endl;
1224 cout << "Macro CT1Analysis : JobB_EST_UP, WESTUP = "
1225 << JobB_EST_UP << ", " << WESTUP << endl;
1226
1227
1228 TString typeON = "ON";
1229 TString typeMC = "MC";
1230 TString ext = "3.root";
1231
1232 //------------------------------
1233 // name of MC file to be used for optimizing the energy estimator
1234 TString filenameOpt(outPath);
1235 filenameOpt += typeMC;
1236 filenameOpt += ext;
1237 cout << "filenameOpt = " << filenameOpt << endl;
1238
1239 //------------------------------
1240 // selection of g/h separation method
1241 // and definition of final selections
1242
1243 TString XX("NN");
1244 //TString XX("SC");
1245 //TString XX("RF");
1246 TString fhadronnessName("Had");
1247 fhadronnessName += XX;
1248 cout << "fhadronnessName = " << fhadronnessName << endl;
1249
1250 // maximum values of the hadronness and of |alpha|
1251 Float_t maxhadronness = 0.40;
1252 Float_t maxalpha = 20.0;
1253 cout << "Maximum values of hadronness and |alpha| = "
1254 << maxhadronness << ", " << maxalpha << endl;
1255
1256 // name of file containing the parameters of the energy estimator
1257 TString energyParName(outPath);
1258 energyParName += "energyest_";
1259 energyParName += XX;
1260 energyParName += ".root";
1261 cout << "energyParName = " << energyParName << endl;
1262
1263
1264 //------------------------------
1265 // name of ON file to be updated
1266 TString filenameON(outPath);
1267 filenameON += typeON;
1268 filenameON += ext;
1269 cout << "filenameON = " << filenameON << endl;
1270
1271 // name of MC file to be updated
1272 TString filenameMC(outPath);
1273 filenameMC += typeMC;
1274 filenameMC += ext;
1275 cout << "filenameMC = " << filenameMC << endl;
1276
1277 //------------------------------
1278 // name of updated ON file
1279 TString filenameONup(outPath);
1280 filenameONup += typeON;
1281 filenameONup += "_";
1282 filenameONup += XX;
1283 filenameONup += ext;
1284 cout << "filenameONup = " << filenameONup << endl;
1285
1286 // name of updated MC file
1287 TString filenameMCup(outPath);
1288 filenameMCup += typeMC;
1289 filenameMCup += "_";
1290 filenameMCup += XX;
1291 filenameMCup += ext;
1292 cout << "filenameMCup = " << filenameMCup << endl;
1293
1294 //-----------------------------------------------------------
1295
1296 TString fHilName = "MHillas";
1297 TString fHilSrcName = "MHillasSrc";
1298 TString fnewparName = "MNewImagePar";
1299
1300 //===========================================================
1301 //
1302 // Optimization of energy estimator
1303 //
1304
1305 TString inpath("");
1306 TString outpath("");
1307 Int_t howMany = 2000;
1308 CT1EEst(inpath, filenameOpt, outpath, energyParName,
1309 fHilName, fHilSrcName, fhadronnessName,
1310 howMany, maxhadronness, maxalpha);
1311
1312 //-----------------------------------------------------------
1313 //
1314 // Read in parameters of energy estimator
1315 //
1316
1317 TFile enparam(energyParName);
1318
1319 TArrayD parA(5);
1320 TArrayD parB(7);
1321 for (Int_t i=0; i<parA.GetSize(); i++)
1322 parA[i] = EnergyParams(i);
1323 for (Int_t i=0; i<parB.GetSize(); i++)
1324 parB[i] = EnergyParams( i+parA.GetSize() );
1325
1326
1327 if (WESTUP)
1328 {
1329 //========== start update ============================================
1330 //
1331 // Update ON and MC root files with the estimated energy
1332
1333 //---------------------------------------------------
1334 // Update ON data
1335 //
1336 MTaskList tliston;
1337 MParList pliston;
1338
1339 //-------------------------------------------
1340 // create the tasks which should be executed
1341 //
1342
1343 MReadMarsFile read("Events", filenameON);
1344 read.DisableAutoScheme();
1345
1346 //---------------------------
1347 // calculate estimated energy
1348
1349 MEnergyEstParam eest2(fHilName);
1350 eest2.Add(fHilSrcName);
1351
1352 eest2.SetCoeffA(parA);
1353 eest2.SetCoeffB(parB);
1354
1355 //.......................................................................
1356
1357 MWriteRootFile write(filenameONup);
1358
1359 write.AddContainer("MRawRunHeader", "RunHeaders");
1360 write.AddContainer("MTime", "Events");
1361 write.AddContainer("MMcEvt", "Events");
1362 write.AddContainer("MSrcPosCam", "Events");
1363 write.AddContainer("MSigmabar", "Events");
1364 write.AddContainer("MHillas", "Events");
1365 write.AddContainer("MHillasSrc", "Events");
1366 write.AddContainer("MNewImagePar", "Events");
1367
1368 write.AddContainer("HadNN", "Events");
1369 write.AddContainer("HadSC", "Events");
1370
1371 write.AddContainer("MEnergyEst", "Events");
1372
1373 //-----------------------------------------------------------------
1374
1375 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
1376 selfinal.SetCuts(maxhadronness, maxalpha);
1377 selfinal.SetHadronnessName(fhadronnessName);
1378 MContinue contfinal(&selfinal);
1379
1380
1381 //*****************************
1382 // entries in MParList
1383
1384 pliston.AddToList(&tliston);
1385 InitBinnings(&pliston);
1386
1387
1388 //*****************************
1389 // entries in MTaskList
1390
1391 tliston.AddToList(&read);
1392 tliston.AddToList(&eest2);
1393 tliston.AddToList(&write);
1394 tliston.AddToList(&contfinal);
1395
1396 //*****************************
1397
1398 //-------------------------------------------
1399 // Execute event loop
1400 //
1401 MProgressBar bar;
1402 MEvtLoop evtloop;
1403 evtloop.SetParList(&pliston);
1404 evtloop.SetProgressBar(&bar);
1405 //evtloop.Write();
1406
1407 Int_t maxevents = 1000000000;
1408 //Int_t maxevents = 1000;
1409 if ( !evtloop.Eventloop(maxevents) )
1410 return;
1411
1412 tliston.PrintStatistics(0, kTRUE);
1413 DeleteBinnings(&pliston);
1414
1415 //---------------------------------------------------
1416 //---------------------------------------------------
1417 // Update MC data
1418 //
1419 MTaskList tliston;
1420 MParList pliston;
1421
1422 //-------------------------------------------
1423 // create the tasks which should be executed
1424 //
1425
1426 MReadMarsFile read("Events", filenameMC);
1427 read.DisableAutoScheme();
1428
1429 //---------------------------
1430 // calculate estimated energy
1431
1432 MEnergyEstParam eest2(fHilName);
1433 eest2.Add(fHilSrcName);
1434
1435 eest2.SetCoeffA(parA);
1436 eest2.SetCoeffB(parB);
1437
1438 //.......................................................................
1439
1440 MWriteRootFile write(filenameMCup);
1441
1442 write.AddContainer("MRawRunHeader", "RunHeaders");
1443 write.AddContainer("MTime", "Events");
1444 write.AddContainer("MMcEvt", "Events");
1445 write.AddContainer("MSrcPosCam", "Events");
1446 write.AddContainer("MSigmabar", "Events");
1447 write.AddContainer("MHillas", "Events");
1448 write.AddContainer("MHillasSrc", "Events");
1449 write.AddContainer("MNewImagePar", "Events");
1450
1451 write.AddContainer("HadNN", "Events");
1452 write.AddContainer("HadSC", "Events");
1453
1454 write.AddContainer("MEnergyEst", "Events");
1455
1456 //-----------------------------------------------------------------
1457
1458 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
1459 selfinal.SetCuts(maxhadronness, maxalpha);
1460 selfinal.SetHadronnessName(fhadronnessName);
1461 MContinue contfinal(&selfinal);
1462
1463
1464 //*****************************
1465 // entries in MParList
1466
1467 pliston.AddToList(&tliston);
1468 InitBinnings(&pliston);
1469
1470
1471 //*****************************
1472 // entries in MTaskList
1473
1474 tliston.AddToList(&read);
1475 tliston.AddToList(&eest2);
1476 tliston.AddToList(&write);
1477 tliston.AddToList(&contfinal);
1478
1479 //*****************************
1480
1481 //-------------------------------------------
1482 // Execute event loop
1483 //
1484 MProgressBar bar;
1485 MEvtLoop evtloop;
1486 evtloop.SetParList(&pliston);
1487 evtloop.SetProgressBar(&bar);
1488 //evtloop.Write();
1489
1490 Int_t maxevents = 1000000000;
1491 //Int_t maxevents = 1000;
1492 if ( !evtloop.Eventloop(maxevents) )
1493 return;
1494
1495 tliston.PrintStatistics(0, kTRUE);
1496 DeleteBinnings(&pliston);
1497
1498
1499 //========== end update ============================================
1500 }
1501
1502 enparam.Close();
1503
1504 cout << "Macro CT1Analysis : End of Job B_EST_UP" << endl;
1505 cout << "=======================================================" << endl;
1506 }
1507 //---------------------------------------------------------------------
1508
1509
1510 //---------------------------------------------------------------------
1511 // Job C_NN
1512 //=========
1513
1514 // read ON1 and MC1 data files
1515 //
1516 // - calculate hadronness for method of Nearest Neighbors
1517 // - produce Neyman-Pearson plot
1518
1519 if (JobC_NN)
1520 {
1521 cout << "=====================================================" << endl;
1522 cout << "Macro CT1Analysis : Start of Job C_NN" << endl;
1523
1524 cout << "" << endl;
1525 cout << "Macro CT1Analysis : JobC_NN = " << JobC_NN << endl;
1526
1527
1528 TString typeMC = "MC";
1529 cout << "typeMC = " << typeMC << endl;
1530
1531 TString typeData = "ON";
1532 cout << "typeData = " << typeData << endl;
1533
1534 TString typeMatrixHadrons = "ON";
1535 cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
1536
1537
1538
1539 //*************************************************************************
1540 // read in matrices of look-alike events
1541
1542 const char* mtxName = "MatrixGammas";
1543
1544 TString outName = outPath;
1545 outName += "matrix_gammas_";
1546 outName += "MC";
1547 outName += ".root";
1548
1549 cout << "" << endl;
1550 cout << "========================================================" << endl;
1551 cout << "Get matrix for (gammas)" << endl;
1552 cout << "matrix name = " << mtxName << endl;
1553 cout << "name of root file = " << outName << endl;
1554 cout << "" << endl;
1555
1556
1557 // read in the object with the name 'mtxName' from file 'outName'
1558 //
1559 TFile *fileg = new TFile(outName);
1560
1561 MHMatrix matrixg;
1562 matrixg.Read(mtxName);
1563 matrixg.Print("cols");
1564
1565 delete fileg;
1566
1567
1568 //*****************************************************************
1569
1570 const char* mtxName = "MatrixHadrons";
1571
1572 TString outName = outPath;
1573 outName += "matrix_hadrons_";
1574 outName += typeMatrixHadrons;
1575 outName += ".root";
1576
1577 cout << "" << endl;
1578 cout << "========================================================" << endl;
1579 cout << " Get matrix for (hadrons)" << endl;
1580 cout << "matrix name = " << mtxName << endl;
1581 cout << "name of root file = " << outName << endl;
1582 cout << "" << endl;
1583
1584
1585 // read in the object with the name mtxName
1586 //
1587 TFile *fileh = new TFile(outName);
1588
1589 MHMatrix matrixh;
1590 matrixh.Read(mtxName);
1591 matrixh.Print("cols");
1592
1593 delete fileh;
1594
1595 //*****************************************************************
1596
1597
1598 //-----------------------------------------------------------------
1599
1600 MTaskList tliston;
1601 MParList pliston;
1602
1603 //-------------------------------------------
1604 // create the tasks which should be executed
1605 //
1606 TString filenameData = outPath;
1607 filenameData += typeData;
1608 filenameData += "1.root";
1609
1610 TString filenameMC = outPath;
1611 filenameMC += typeMC;
1612 filenameMC += "1.root";
1613
1614 cout << "filenameData = " << filenameData << endl;
1615 cout << "filenameMC = " << filenameMC << endl;
1616
1617 MReadMarsFile read("Events", filenameMC);
1618 read.AddFile(filenameData);
1619 read.DisableAutoScheme();
1620
1621 MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
1622 MFParticleId fhadron("MMcEvt", '!', kGAMMA);
1623
1624 //.......................................................................
1625 // calculate hadronnes for method of Nearest Neighbors
1626
1627 TString hadName = "HadNN";
1628 MMultiDimDistCalc nncalc;
1629 nncalc.SetUseNumRows(25);
1630 nncalc.SetUseKernelMethod(kFALSE);
1631 nncalc.SetHadronnessName(hadName);
1632
1633 //.......................................................................
1634
1635 // geometry is needed in MHHillas... classes
1636 MGeomCam *fGeom =
1637 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1638
1639 TString fHilName = "MHillas";
1640 TString fHilSrcName = "MHillasSrc";
1641 TString fnewparName = "MNewImagePar";
1642
1643 Float_t hadcut = 0.36;
1644 Float_t alphacut = 100.0;
1645 MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
1646 selfinalgh.SetCuts(hadcut, alphacut);
1647 selfinalgh.SetHadronnessName(hadName);
1648 MContinue contfinalgh(&selfinalgh);
1649
1650 MFillH fillhadnn("MHHadronness", hadName);
1651
1652 Float_t hadcut = 0.36;
1653 Float_t alphacut = 20.0;
1654 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
1655 selfinal.SetCuts(hadcut, alphacut);
1656 selfinal.SetHadronnessName(hadName);
1657 MContinue contfinal(&selfinal);
1658
1659 MFillH hfill1("MHHillas", fHilName);
1660 MFillH hfill2("MHStarMap", fHilName);
1661
1662 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
1663
1664 MFillH hfill3("MHHillasExt", fHilSrcName);
1665 MFillH hfill4("MHHillasSrc", fHilSrcName);
1666 MFillH hfill5("MHNewImagePar", fnewparName);
1667
1668
1669 //*****************************
1670 // entries in MParList
1671
1672 pliston.AddToList(&tliston);
1673 InitBinnings(&pliston);
1674
1675 pliston.AddToList(&matrixg);
1676 pliston.AddToList(&matrixh);
1677
1678 pliston.AddToList(&hhilext);
1679
1680
1681 //*****************************
1682 // entries in MTaskList
1683
1684 tliston.AddToList(&read);
1685
1686 tliston.AddToList(&fgamma);
1687 tliston.AddToList(&fhadron);
1688 tliston.AddToList(&nncalc);
1689 tliston.AddToList(&fillhadnn);
1690
1691 tliston.AddToList(&contfinalgh);
1692 tliston.AddToList(&hfill1);
1693 tliston.AddToList(&hfill2);
1694 tliston.AddToList(&hfill3);
1695 tliston.AddToList(&hfill4);
1696 tliston.AddToList(&hfill5);
1697
1698 tliston.AddToList(&contfinal);
1699
1700 //*****************************
1701
1702 //-------------------------------------------
1703 // Execute event loop
1704 //
1705 MProgressBar bar;
1706 MEvtLoop evtloop;
1707 evtloop.SetParList(&pliston);
1708 evtloop.SetProgressBar(&bar);
1709 //evtloop.Write();
1710
1711 Int_t maxevents = 1000000000;
1712 //Int_t maxevents = 35000;
1713 if ( !evtloop.Eventloop(maxevents) )
1714 return;
1715
1716 tliston.PrintStatistics(0, kTRUE);
1717 DeleteBinnings(&pliston);
1718
1719 //-------------------------------------------
1720 // Display the histograms
1721 //
1722 TObject *c;
1723 c = pliston.FindObject("MHHadronness");
1724 if (c)
1725 {
1726 c->DrawClone();
1727 c->Print();
1728 }
1729 c = pliston.FindObject("MHHillas")->DrawClone();
1730 c = pliston.FindObject("MHHillasExt")->DrawClone();
1731 c = pliston.FindObject("MHHillasSrc")->DrawClone();
1732 c = pliston.FindObject("MHNewImagePar")->DrawClone();
1733 c = pliston.FindObject("MHStarMap")->DrawClone();
1734
1735
1736
1737 cout << "Macro CT1Analysis : End of Job C_NN" << endl;
1738 cout << "===================================================" << endl;
1739 }
1740
1741
1742
1743 //---------------------------------------------------------------------
1744 // Job C_RF
1745 //=========
1746
1747 // read ON1 and MC1 data files
1748 //
1749 // - calculate hadronness for method of Random Forest
1750 // - produce Neyman-Pearson plot
1751
1752 if (JobC_RF)
1753 {
1754 cout << "=====================================================" << endl;
1755 cout << "Macro CT1Analysis : Start of Job C_RF" << endl;
1756
1757 cout << "" << endl;
1758 cout << "Macro CT1Analysis : JobC_RF = " << JobC_RF << endl;
1759
1760
1761 TString typeMC = "MC";
1762 cout << "typeMC = " << typeMC << endl;
1763
1764 TString typeData = "ON";
1765 cout << "typeData = " << typeData << endl;
1766
1767 TString typeMatrixHadrons = "ON";
1768 cout << "typeMatrixHadrons = " << typeMatrixHadrons << endl;
1769
1770 //*************************************************************************
1771 // read in matrices of look-alike events
1772
1773 const char* mtxName = "MatrixGammas";
1774
1775 TString outName = outPath;
1776 outName += "matrix_gammas_";
1777 outName += "MC";
1778 outName += ".root";
1779
1780 cout << "" << endl;
1781 cout << "========================================================" << endl;
1782 cout << "Get matrix for (gammas)" << endl;
1783 cout << "matrix name = " << mtxName << endl;
1784 cout << "name of root file = " << outName << endl;
1785 cout << "" << endl;
1786
1787
1788 // read in the object with the name 'mtxName' from file 'outName'
1789 //
1790 TFile *fileg = new TFile(outName);
1791
1792 MHMatrix matrixg;
1793 matrixg.Read(mtxName);
1794 matrixg.Print("cols");
1795
1796 delete fileg;
1797
1798
1799 //*****************************************************************
1800
1801 const char* mtxName = "MatrixHadrons";
1802
1803 TString outName = outPath;
1804 outName += "matrix_hadrons_";
1805 outName += typeMatrixHadrons;
1806 outName += ".root";
1807
1808 cout << "" << endl;
1809 cout << "========================================================" << endl;
1810 cout << " Get matrix for (hadrons)" << endl;
1811 cout << "matrix name = " << mtxName << endl;
1812 cout << "name of root file = " << outName << endl;
1813 cout << "" << endl;
1814
1815
1816 // read in the object with the name mtxName
1817 //
1818 TFile *fileh = new TFile(outName);
1819
1820 MHMatrix matrixh;
1821 matrixh.Read(mtxName);
1822 matrixh.Print("cols");
1823
1824 delete fileh;
1825
1826 //*****************************************************************
1827
1828 //-----------------------------------------------------------------
1829 // grow the trees of the random forest (event loop = tree loop)
1830
1831 cout << "Macro CT1Analysis : start growing trees" << endl;
1832
1833 MTaskList tlist2;
1834 MParList plist2;
1835 plist2.AddToList(&tlist2);
1836
1837 plist2.AddToList(&matrixg);
1838 plist2.AddToList(&matrixh);
1839
1840 MRanForestGrow rfgrow2;
1841 rfgrow2.SetNumTrees(100);
1842 rfgrow2.SetNumTry(4);
1843 rfgrow2.SetNdSize(10);
1844
1845 TString outRF = outPath;
1846 outRF += "RF.root";
1847 MWriteRootFile rfwrite2(outRF);
1848 rfwrite2.AddContainer("MRanTree", "TREE");
1849
1850 // list of tasks for the loop over the trees
1851
1852 tlist2.AddToList(&rfgrow2);
1853 tlist2.AddToList(&rfwrite2);
1854
1855 //-------------------
1856 // Execute tree loop
1857 //
1858 MEvtLoop treeloop;
1859 treeloop.SetParList(&plist2);
1860
1861 if ( !treeloop.Eventloop() )
1862 return;
1863
1864 tlist2.PrintStatistics(0, kTRUE);
1865
1866
1867 // get adresses of objects which are used in the next eventloop
1868 MRanForest *fRanForest = (MRanForest*)plist2->FindObject("MRanForest");
1869 if (!fRanForest)
1870 {
1871 *fLog << err << dbginf << "MRanForest not found... aborting." << endl;
1872 return kFALSE;
1873 }
1874
1875 MRanTree *fRanTree = (MRanTree*)plist2->FindObject("MRanTree");
1876 if (!fRanTree)
1877 {
1878 *fLog << err << dbginf << "MRanTree not found... aborting." << endl;
1879 return kFALSE;
1880 }
1881
1882 //-----------------------------------------------------------------
1883
1884 MTaskList tliston;
1885 MParList pliston;
1886
1887 //-------------------------------------------
1888 // create the tasks which should be executed
1889 //
1890 TString filenameData = outPath;
1891 filenameData += typeData;
1892 filenameData += "1.root";
1893
1894 TString filenameMC = outPath;
1895 filenameMC += typeMC;
1896 filenameMC += "1.root";
1897
1898 cout << "filenameData = " << filenameData << endl;
1899 cout << "filenameMC = " << filenameMC << endl;
1900
1901 MReadMarsFile read("Events", filenameMC);
1902 read.AddFile(filenameData);
1903 read.DisableAutoScheme();
1904
1905 MFParticleId fgamma ("MMcEvt", '=', kGAMMA);
1906 MFParticleId fhadron("MMcEvt", '!', kGAMMA);
1907
1908 //.......................................................................
1909 // calculate hadronnes for method of Random Forest
1910
1911 TString hadName = "HadRF";
1912 MRanForestCalc rfcalc;
1913 rfcalc.SetHadronnessName(hadName);
1914
1915 //.......................................................................
1916
1917 // geometry is needed in MHHillas... classes
1918 MGeomCam *fGeom =
1919 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1920
1921 TString fHilName = "MHillas";
1922 TString fHilSrcName = "MHillasSrc";
1923 TString fnewparName = "MNewImagePar";
1924
1925 Float_t hadcut = 0.36;
1926 Float_t alphacut = 100.0;
1927 MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
1928 selfinalgh.SetCuts(hadcut, alphacut);
1929 selfinalgh.SetHadronnessName(hadName);
1930 MContinue contfinalgh(&selfinalgh);
1931
1932 MFillH fillhadrf("MHHadronness", hadName);
1933 MFillH fillhrf("MHRanForest");
1934
1935 Float_t hadcut = 0.36;
1936 Float_t alphacut = 20.0;
1937 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
1938 selfinal.SetCuts(hadcut, alphacut);
1939 selfinal.SetHadronnessName(hadName);
1940 MContinue contfinal(&selfinal);
1941
1942 MFillH hfill1("MHHillas", fHilName);
1943 MFillH hfill2("MHStarMap", fHilName);
1944
1945 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
1946
1947 MFillH hfill3("MHHillasExt", fHilSrcName);
1948 MFillH hfill4("MHHillasSrc", fHilSrcName);
1949 MFillH hfill5("MHNewImagePar", fnewparName);
1950
1951
1952 //*****************************
1953 // entries in MParList
1954
1955 pliston.AddToList(&tliston);
1956 InitBinnings(&pliston);
1957
1958 pliston.AddToList(fRanForest);
1959 pliston.AddToList(fRanTree);
1960
1961 pliston.AddToList(&hhilext);
1962
1963
1964 //*****************************
1965 // entries in MTaskList
1966
1967 tliston.AddToList(&read);
1968
1969 tliston.AddToList(&fgamma);
1970 tliston.AddToList(&fhadron);
1971 tliston.AddToList(&rfcalc);
1972 tliston.AddToList(&fillhadrf);
1973 tliston.AddToList(&fillhrf);
1974
1975 tliston.AddToList(&contfinalgh);
1976 tliston.AddToList(&hfill1);
1977 tliston.AddToList(&hfill2);
1978 tliston.AddToList(&hfill3);
1979 tliston.AddToList(&hfill4);
1980 tliston.AddToList(&hfill5);
1981
1982 tliston.AddToList(&contfinal);
1983
1984 //*****************************
1985
1986 //-------------------------------------------
1987 // Execute event loop
1988 //
1989 MProgressBar bar;
1990 MEvtLoop evtloop;
1991 evtloop.SetParList(&pliston);
1992 evtloop.SetProgressBar(&bar);
1993 //evtloop.Write();
1994
1995 Int_t maxevents = 1000000000;
1996 //Int_t maxevents = 35000;
1997 if ( !evtloop.Eventloop(maxevents) )
1998 return;
1999
2000 tliston.PrintStatistics(0, kTRUE);
2001 DeleteBinnings(&pliston);
2002
2003 //-------------------------------------------
2004 // Display the histograms
2005 //
2006 TObject *c;
2007 c = pliston.FindObject("MHRanForest")->DrawClone();
2008 c = pliston.FindObject("MHHadronness");
2009 if (c)
2010 {
2011 c->DrawClone();
2012 c->Print();
2013 }
2014
2015 c = pliston.FindObject("MHHillas")->DrawClone();
2016 c = pliston.FindObject("MHHillasExt")->DrawClone();
2017 c = pliston.FindObject("MHHillasSrc")->DrawClone();
2018 c = pliston.FindObject("MHNewImagePar")->DrawClone();
2019 c = pliston.FindObject("MHStarMap")->DrawClone();
2020
2021
2022
2023 cout << "Macro CT1Analysis : End of Job C_RF" << endl;
2024 cout << "===================================================" << endl;
2025 }
2026
2027
2028 //---------------------------------------------------------------------
2029 // Job D_XX
2030 //=========
2031
2032 // - select g/h separation method XX
2033 // - read MC_XX2.root file
2034   // - calculate eff. collection area
2035 // - read ON_XX2.root file
2036 // - apply final cuts
2037 // - calculate flux
2038 // - write root file for ON data after final cuts (ON_XX3.root))
2039
2040
2041 if (JobD_XX)
2042 {
2043 cout << "=====================================================" << endl;
2044 cout << "Macro CT1Analysis : Start of Job D_XX" << endl;
2045
2046 cout << "" << endl;
2047 cout << "Macro CT1Analysis : JobD_XX, WXXD = "
2048 << JobD_XX << ", " << WXXD << endl;
2049
2050 // type of data to be analysed
2051 TString typeData = "ON";
2052 //TString typeData = "OFF";
2053 //TString typeData = "MC";
2054 cout << "typeData = " << typeData << endl;
2055
2056 TString typeMC = "MC";
2057 TString ext = "3.root";
2058 TString extout = "4.root";
2059
2060 //------------------------------
2061 // selection of g/h separation method
2062 // and definition of final selections
2063
2064 TString XX("NN");
2065 //TString XX("SC");
2066 //TString XX("RF");
2067 TString fhadronnessName("Had");
2068 fhadronnessName += XX;
2069 cout << "fhadronnessName = " << fhadronnessName << endl;
2070
2071 // maximum values of the hadronness and of |alpha|
2072 Float_t maxhadronness = 0.40;
2073 Float_t maxalpha = 20.0;
2074 cout << "Maximum values of hadronness and |alpha| = "
2075 << maxhadronness << ", " << maxalpha << endl;
2076
2077
2078 //------------------------------
2079 // name of MC file to be used for calculating the eff. collection areas
2080 TString filenameArea(outPath);
2081 filenameArea += typeMC;
2082 filenameArea += "_";
2083 filenameArea += XX;
2084 filenameArea += ext;
2085 cout << "filenameArea = " << filenameArea << endl;
2086
2087 //------------------------------
2088 // name of file containing the eff. collection areas
2089 TString collareaName(outPath);
2090 collareaName += "area_";
2091 collareaName += XX;
2092 collareaName += ".root";
2093 cout << "collareaName = " << collareaName << endl;
2094
2095 //------------------------------
2096 // name of data file to be analysed
2097 TString filenameData(outPath);
2098 filenameData += typeData;
2099 filenameData += "_";
2100 filenameData += XX;
2101 filenameData += ext;
2102 cout << "filenameData = " << filenameData << endl;
2103
2104 //------------------------------
2105 // name of output data file (after the final cuts)
2106 TString filenameDataout(outPath);
2107 filenameDataout += typeData;
2108 filenameDataout += "_";
2109 filenameDataout += XX;
2110 filenameDataout += extout;
2111 cout << "filenameDataout = " << filenameDataout << endl;
2112
2113
2114 //====================================================================
2115 cout << "-----------------------------------------------" << endl;
2116 cout << "Start calculation of effective collection areas" << endl;
2117 MParList parlist;
2118 MTaskList tasklist;
2119
2120 //---------------------------------------
2121 // Setup the tasks to be executed
2122 //
2123 MReadMarsFile reader("Events", filenameArea);
2124
2125 MFCT1SelFinal cuthadrons;
2126 cuthadrons.SetHadronnessName(fhadronnessName);
2127 cuthadrons.SetCuts(maxhadronness, maxalpha);
2128
2129 MContinue conthadrons(&cuthadrons);
2130
2131 MHMcCT1CollectionArea* collarea =
2132 new MHMcCT1CollectionArea("MHMcCT1CollectionArea","",30,2.,5.);
2133 MMcCT1CollectionAreaCalc effi;
2134
2135 //********************************
2136 // entries in MParList
2137
2138 parlist.AddToList(&tasklist);
2139 parlist.AddToList(collarea);
2140
2141 //********************************
2142 // entries in MTaskList
2143
2144 tasklist.AddToList(&reader);
2145 tasklist.AddToList(&conthadrons);
2146 tasklist.AddToList(&effi);
2147
2148 //********************************
2149
2150 //-----------------------------------------
2151 // Execute event loop
2152 //
2153 MEvtLoop magic;
2154 magic.SetParList(&parlist);
2155
2156 MProgressBar bar;
2157 magic.SetProgressBar(&bar);
2158 if (!magic.Eventloop())
2159 return;
2160
2161 tasklist.PrintStatistics(0, kTRUE);
2162
2163 // Display the histograms
2164 //
2165 parlist.FindObject("MHMcCT1CollectionArea")->DrawClone("lego");
2166
2167 //---------------------------------------------
2168 // Write histograms to a file
2169 //
2170
2171 TFile f(collareaName, "RECREATE");
2172 collarea->GetHist()->Write();
2173 collarea->GetHAll()->Write();
2174 collarea->GetHSel()->Write();
2175 f.Close();
2176
2177
2178 cout << "Calculation of effective collection areas done" << endl;
2179 cout << "-----------------------------------------------" << endl;
2180 //------------------------------------------------------------------
2181
2182
2183 //*************************************************************************
2184 //
2185 // Analyse the data
2186 //
2187
2188 MTaskList tliston;
2189 MParList pliston;
2190
2191 TString hadName("Had");
2192 hadName += XX;
2193
2194 TString fHilName = "MHillas";
2195 TString fHilSrcName = "MHillasSrc";
2196 TString fnewparName = "MNewImagePar";
2197
2198 //-------------------------------------------
2199 // create the tasks which should be executed
2200 //
2201
2202 MReadMarsFile read("Events", filenameData);
2203 read.DisableAutoScheme();
2204
2205 //.......................................................................
2206
2207 if (WXXD)
2208 {
2209 MWriteRootFile write(filenameDataout);
2210
2211 write.AddContainer("MRawRunHeader", "RunHeaders");
2212 write.AddContainer("MTime", "Events");
2213 write.AddContainer("MMcEvt", "Events");
2214 write.AddContainer("MSrcPosCam", "Events");
2215 write.AddContainer("MSigmabar", "Events");
2216 write.AddContainer("MHillas", "Events");
2217 write.AddContainer("MHillasSrc", "Events");
2218 write.AddContainer("MNewImagePar", "Events");
2219
2220 write.AddContainer("HadNN", "Events");
2221 write.AddContainer("HadSC", "Events");
2222
2223 write.AddContainer("MEnergyEst", "Events");
2224 }
2225
2226 //-----------------------------------------------------------------
2227 // geometry is needed in MHHillas... classes
2228 MGeomCam *fGeom =
2229 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
2230
2231 MFCT1SelFinal selfinalgh(fHilName, fHilSrcName);
2232 selfinalgh.SetCuts(maxhadronness, 100.0);
2233 selfinalgh.SetHadronnessName(fhadronnessName);
2234 MContinue contfinalgh(&selfinalgh);
2235
2236 MFillH fillhadnn("MHadNN[MHHadronness]", "HadNN");
2237 MFillH fillhadsc("MHadSC[MHHadronness]", "HadSC");
2238 //MFillH fillhadrf("MHadRF[MHHadronness]", "HadRF");
2239
2240 MFillH hfill1("MHHillas", fHilName);
2241 MFillH hfill2("MHStarMap", fHilName);
2242
2243 MHHillasExt hhilext("MHHillasExt", "MHHillasExt", fHilName);
2244
2245 MFillH hfill3("MHHillasExt", fHilSrcName);
2246 MFillH hfill4("MHHillasSrc", fHilSrcName);
2247 MFillH hfill5("MHNewImagePar", fnewparName);
2248
2249 MFCT1SelFinal selfinal(fHilName, fHilSrcName);
2250 selfinal.SetCuts(maxhadronness, maxalpha);
2251 selfinal.SetHadronnessName(fhadronnessName);
2252 MContinue contfinal(&selfinal);
2253
2254
2255 //*****************************
2256 // entries in MParList
2257
2258 pliston.AddToList(&tliston);
2259 InitBinnings(&pliston);
2260
2261 pliston.AddToList(&hhilext);
2262
2263 //*****************************
2264 // entries in MTaskList
2265
2266 tliston.AddToList(&read);
2267
2268 tliston.AddToList(&contfinalgh);
2269
2270 if (WXXD)
2271 tliston.AddToList(&write);
2272
2273 tliston.AddToList(&contfinal);
2274 tliston.AddToList(&fillhadnn);
2275 tliston.AddToList(&fillhadsc);
2276 //tliston.AddToList(&fillhadrf);
2277 tliston.AddToList(&hfill1);
2278 tliston.AddToList(&hfill2);
2279 tliston.AddToList(&hfill3);
2280 tliston.AddToList(&hfill4);
2281 tliston.AddToList(&hfill5);
2282
2283
2284
2285 //*****************************
2286
2287 //-------------------------------------------
2288 // Execute event loop
2289 //
2290 MProgressBar bar;
2291 MEvtLoop evtloop;
2292 evtloop.SetParList(&pliston);
2293 evtloop.SetProgressBar(&bar);
2294 //evtloop.Write();
2295
2296 Int_t maxevents = 1000000000;
2297 //Int_t maxevents = 10000;
2298 if ( !evtloop.Eventloop(maxevents) )
2299 return;
2300
2301 tliston.PrintStatistics(0, kTRUE);
2302 DeleteBinnings(&pliston);
2303
2304 //-------------------------------------------
2305 // Display the histograms
2306 //
2307 TObject *c;
2308
2309 c = pliston.FindObject("MHadNN", "MHHadronness");
2310 if (c)
2311 {
2312 c->DrawClone();
2313 c->Print();
2314 }
2315 c = pliston.FindObject("MHadRF", "MHHadronness");
2316 if (c)
2317 {
2318 c->DrawClone();
2319 c->Print();
2320 }
2321 c = pliston.FindObject("MHadSC", "MHHadronness");
2322 if (c)
2323 {
2324 c->DrawClone();
2325 c->Print();
2326 }
2327
2328 c = pliston.FindObject("MHHillas")->DrawClone();
2329 c = pliston.FindObject("MHHillasExt")->DrawClone();
2330 c = pliston.FindObject("MHHillasSrc")->DrawClone();
2331 c = pliston.FindObject("MHNewImagePar")->DrawClone();
2332 c = pliston.FindObject("MHStarMap")->DrawClone();
2333
2334
2335 cout << "Macro CT1Analysis : End of Job D_XX" << endl;
2336 cout << "=======================================================" << endl;
2337 }
2338 //---------------------------------------------------------------------
2339
2340}
2341
2342
2343
2344
2345
2346
2347
Note: See TracBrowser for help on using the repository browser.