source: trunk/MagicSoft/Mars/macros/AnalyseCT1.C@ 1835

Last change on this file since 1835 was 1809, checked in by wittek, 22 years ago
*** empty log message ***
File size: 57.9 KB
Line 
1void AnalyseCT1()
2{
3 //-----------------------------------------------
4 //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc";
5 //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc";
6 //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6";
7 //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc";
8 //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
9
10 //-----------------------------------------------
11 const char *offfile = "~/Mars200203/CT1data/offdata.preproc";
12
13 const char *onfile = "~/Mars200203/CT1data/mkn421_on.preproc";
14 //const char *onfile = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6";
15
16 //const char *mcfile = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5
17 //const char *mcfile = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6";
18 const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc";
19
20 //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6";
21 //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303";
22
23 //-----------------------------------------------
24 //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc";
25 //const char *onfile = "~magican/home/ct1test/PreprocData/mkn421_on.preproc";
26 //const char *mcfile = "~magican/home/ct1test/PreprocData/mc_gammas.preproc";
27
28
29 const char *filename;
30
31 //************************************************************************
32 // Job 1 : read OFF data
33 // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1);
34 // generate hadron matrix for g/h separation)
35
36 Bool_t Job1 = kFALSE;
37 Bool_t WHistOFF = kTRUE; // write out histogram sigmabar vs. Theta ?
38 Bool_t WLookOFF = kTRUE; // write out look-alike events for hadrons ?
39 Bool_t WOFF1 = kTRUE; // write out root file OFF1 ?
40
41
42 // Job 2 : read ON data
43 // (generate sigmabar vs. Theta plot; write root file for ON data (ON1))
44
45 Bool_t Job2 = kFALSE;
46 Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ?
47 Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ?
48 Bool_t WON1 = kFALSE; // write out root file ON1 ?
49
50 // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data
51 // (do padding; write root file for MC gammas (MC1);
52 // generate gamma matrix for g/h separation)
53
54 Bool_t Job3 = kFALSE;
55 Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ?
56 Bool_t WMC1 = kTRUE; // write out root file MC1 ?
57
58
59
60 // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix
61 // (produce the Neyman-Pearson plot)
62 Bool_t Job4 = kTRUE;
63
64
65 // Job 5 : read MC1 file, read hadron and gamma matrix
66 // (do g/h separation; write root file for MC gammas after final cuts (MC2))
67 Bool_t Job5 = kFALSE;
68 Bool_t WMC2 = kFALSE; // write out root file MC2 ?
69
70
71 // Job 6 : read ON data, read hadron and gamma matrix
72 // (do g/h separation; write root file for ON data after final cuts (ON2))
73 Bool_t Job6 = kFALSE;
74 Bool_t WON2 = kTRUE; // write out root file ON2 ?
75 //************************************************************************
76
77
78
79
80 //---------------------------------------------------------------------
81 // Job 1
82 //======
83
84 // read OFF data file
85
86 // - produce the 2D-histogram "sigmabar versus Theta" for OFF data
87 // (to be used for the padding of the MC gamma data)
88
89 // - write a file of look-alike events (hadron matrix)
90 // (to be used later together with the corresponding MC gamma file
91 // for the g/h separation in the analysis of the ON data)
92
93 // - write a file of OFF events (OFF1)
94 // (after the standard cuts, before the g/h separation)
95 // (to be used together with the corresponding MC gamma file
96 // for the optimization of the g/h separation)
97
98 if (Job1)
99 {
100 cout << "=====================================================" << endl;
101 cout << "Macro AnalyseCT1 : Start of Job 1" << endl;
102
103 cout << "" << endl;
104 cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = "
105 << Job1 << ", " << WHistOFF << ", " << WLookOFF << ", "
106 << WOFF1 << endl;
107
108
109 TString typeData = "OFF";
110 cout << "typeData = " << typeData << endl;
111
112 MTaskList tliston;
113 MParList pliston;
114 pliston.AddToList(&tliston);
115
116
117 //::::::::::::::::::::::::::::::::::::::::::
118
119 MBinning binssize("BinningSize");
120 binssize.SetEdgesLog(50, 10, 1.0e5);
121 pliston.AddToList(&binssize);
122
123 MBinning binsdistc("BinningDist");
124 binsdistc.SetEdges(50, 0, 1.4);
125 pliston.AddToList(&binsdistc);
126
127 MBinning binswidth("BinningWidth");
128 binswidth.SetEdges(50, 0, 1.0);
129 pliston.AddToList(&binswidth);
130
131 MBinning binslength("BinningLength");
132 binslength.SetEdges(50, 0, 1.0);
133 pliston.AddToList(&binslength);
134
135 MBinning binsalpha("BinningAlpha");
136 binsalpha.SetEdges(100, -100, 100);
137 pliston.AddToList(&binsalpha);
138
139 MBinning binsasym("BinningAsym");
140 binsasym.SetEdges(50, -1.5, 1.5);
141 pliston.AddToList(&binsasym);
142
143 MBinning binsht("BinningHeadTail");
144 binsht.SetEdges(50, -1.5, 1.5);
145 pliston.AddToList(&binsht);
146
147 MBinning binsm3l("BinningM3Long");
148 binsm3l.SetEdges(50, -1.5, 1.5);
149 pliston.AddToList(&binsm3l);
150
151 MBinning binsm3t("BinningM3Trans");
152 binsm3t.SetEdges(50, -1.5, 1.5);
153 pliston.AddToList(&binsm3t);
154
155
156 //.....
157 MBinning binsb("BinningSigmabar");
158 binsb.SetEdges( 100, 0.0, 5.0);
159 pliston.AddToList(&binsb);
160
161 MBinning binth("BinningTheta");
162 binth.SetEdges( 70, -0.5, 69.5);
163 pliston.AddToList(&binth);
164
165 MBinning binsdiff("BinningDiffsigma2");
166 binsdiff.SetEdges(100, -5.0, 20.0);
167 pliston.AddToList(&binsdiff);
168 //::::::::::::::::::::::::::::::::::::::::::
169
170
171 //-------------------------------------------
172 // create the tasks which should be executed
173 //
174
175 filename = offfile;
176 RName = "ReadCT1data";
177 MCT1ReadPreProc read(filename, RName);
178
179 MSelBasic selbasic;
180
181 MBlindPixelCalc blind;
182 blind.SetUseInterpolation();
183 blind.SetUseBlindPixels();
184 // blind.SetUseCetralPixel();
185
186 // create an object of class MSigmabar,
187 // because class MHSigmaTheta will look for it
188 MSigmabar sigbar;
189 pliston.AddToList(&sigbar);
190 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
191 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
192
193 MImgCleanStd clean;
194
195 // calculate also extended image parameters
196 TString fHilName = "Hillas";
197 MHillasExt hext(fHilName);
198 pliston.AddToList(&hext);
199 MHillasExt *fHillas = &hext;
200
201 // name of output container for MHillasCalc
202 // = name of MHillas object
203 MHillasCalc hcalc(fHilName);
204
205 // name of output container for MHillasSrcCalc
206 // = name of MHillasSrc object
207 TString fHilSrcName = "HillasSrc";
208 MHillasSrc hilsrc(fHilSrcName);
209 MHillasSrc *fHillasSrc = &hilsrc;
210 pliston.AddToList(&hilsrc);
211 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
212
213 MSelStandard selstand(fHilName, fHilSrcName);
214
215 MFillH hfill1("MHHillas", fHilName);
216 MFillH hfill2("MHStarMap", fHilName);
217
218 TString nxt = "HillasExt";
219 MHHillasExt hhilext(nxt, "", fHilName);
220 pliston.AddToList(&hhilext);
221 TString namext = nxt;
222 namext += "[MHHillasExt]";
223 MFillH hfill3(namext, fHilSrcName);
224
225 MFillH hfill4("MHHillasSrc", fHilSrcName);
226
227 //+++++++++++++++++++++++++++++++++++++++++++++++++++
228 // prepare writing of look-alike events (for hadrons)
229
230 const char* mtxName = "MatrixHadrons";
231 Int_t howMany = 50000;
232 TString outName = "matrix_hadrons_";
233 outName += typeData;
234 outName += ".root";
235
236 cout << "" << endl;
237 cout << "========================================================" << endl;
238 cout << "Matrix for (hadrons)" << endl;
239 cout << "input file = " << filename<< endl;
240 cout << "matrix name = " << mtxName << endl;
241 cout << "max. no.of look-alike events = " << howMany << endl;
242 cout << "name of output root file = " << outName << endl;
243 cout << "" << endl;
244
245
246 // matrix limitation for look-alike events (approximate number)
247 MFEventSelector selector("", "", RName);
248 selector.SetNumSelectEvts(howMany);
249 MFilterList flist;
250 flist.AddToList(&selector);
251
252 //
253 // --- setup the matrix and add it to the parlist
254 //
255 MHMatrix *pmatrix = new MHMatrix(mtxName);
256 MHMatrix& matrix = *pmatrix;
257
258 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
259 matrix.AddColumn("MSigmabar.fSigmabar");
260 matrix.AddColumn("log10(Hillas.fSize)");
261 matrix.AddColumn("HillasSrc.fDist");
262 matrix.AddColumn("Hillas.fWidth");
263 matrix.AddColumn("Hillas.fLength");
264 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
265 matrix.AddColumn("abs(Hillas.fM3Long)");
266 matrix.AddColumn("Hillas.fConc");
267 matrix.AddColumn("Hillas.fConc1");
268
269 pliston.AddToList(pmatrix);
270
271 MFillH fillmatg(mtxName);
272 fillmatg.SetFilter(&flist);
273
274
275 if (WOFF1)
276 {
277 TString outNameImage = typeData;
278 outNameImage += "1.root";
279 MWriteRootFile write(outNameImage);
280 write.AddContainer("MSigmabar", "Events");
281 write.AddContainer("Hillas", "Events");
282 write.AddContainer("MMcEvt", "Events");
283 write.AddContainer("HillasSrc", "Events");
284 write.AddContainer("MRawRunHeader", "RunHeaders");
285 //write.AddContainer("MMcRunHeader","RunHeaders");
286 write.AddContainer("MSrcPosCam", "RunHeaders");
287 }
288
289 //+++++++++++++++++++++++++++++++++++++++++++++++++++
290
291
292 //*****************************
293 // tasks to be executed
294
295 tliston.AddToList(&read);
296
297 tliston.AddToList(&selbasic);
298 tliston.AddToList(&blind);
299 tliston.AddToList(&fillsigtheta);
300 tliston.AddToList(&clean);
301
302 tliston.AddToList(&hcalc);
303 tliston.AddToList(&csrc1);
304
305 tliston.AddToList(&selstand);
306 if (WOFF1)
307 tliston.AddToList(&write);
308
309 tliston.AddToList(&flist);
310 tliston.AddToList(&fillmatg);
311
312 tliston.AddToList(&hfill1);
313 tliston.AddToList(&hfill2);
314 tliston.AddToList(&hfill3);
315 tliston.AddToList(&hfill4);
316
317 //*****************************
318
319 //-------------------------------------------
320 // Execute event loop
321 //
322 MProgressBar bar;
323 MEvtLoop evtloop;
324 evtloop.SetParList(&pliston);
325 evtloop.SetProgressBar(&bar);
326
327 Int_t maxevents = 1000000000;
328 //Int_t maxevents = 1000;
329 if ( !evtloop.Eventloop(maxevents) )
330 return;
331
332 tliston.PrintStatistics(0, kTRUE);
333
334 //-------------------------------------------
335 // Display the histograms
336 //
337 TObject *c;
338 c = pliston.FindObject("MHHillas")->DrawClone();
339
340 c = pliston.FindObject("MHHillasSrc")->DrawClone();
341
342 //pliston.FindObject("MHHillasExt")->DrawClone();
343 c = pliston.FindObject(nxt)->DrawClone();
344
345 c = pliston.FindObject("MHStarMap")->DrawClone();
346
347
348 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
349 //
350 // ----------------------------------------------------------
351 // Definition of the reference sample of look-alike events
352 // (this is a subsample of the original sample)
353 // ----------------------------------------------------------
354 //
355 cout << "" << endl;
356 cout << "========================================================" << endl;
357 // Select a maximum of nmaxevts events from the sample of look-alike
358 // events. They will form the reference sample.
359 Int_t nmaxevts = 2000;
360
361 // target distribution for the variable in column refcolumn (start at 1);
362 // set refcolumn negative if distribution is not to be changed
363 Int_t refcolumn = 1;
364 Int_t nbins = 10;
365 Float_t frombin = 0.5;
366 Float_t tobin = 1.0;
367 TH1F *thsh = new TH1F("thsh","target distribution",
368 nbins, frombin, tobin);
369 thsh->SetXTitle("cos( \\Theta)");
370 thsh->SetYTitle("Counts");
371 Float_t dbin = (tobin-frombin)/nbins;
372 Float_t lbin = frombin +dbin*0.5;
373 for (Int_t j=1; j<=nbins; j++)
374 {
375 thsh->Fill(lbin,1.0);
376 lbin += dbin;
377 }
378
379 cout << "" << endl;
380 cout << "========================================================" << endl;
381 cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl;
382 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
383 << refcolumn << ", " << nmaxevts << endl;
384
385 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
386 {
387 cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
388 return;
389 }
390 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
391
392 //-------------------------------------------
393 // write out look-alike events (for hadrons)
394 //
395 if (WLookOFF)
396 {
397 TFile *writejens = new TFile(outName, "RECREATE", "");
398 matrix.Write();
399 cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file "
400 << outName << endl;
401
402 writejens->Close();
403 delete writejens;
404 }
405
406 //-------------------------------------------
407 // Write histograms onto a file
408 if (WHistOFF)
409 {
410 MHSigmaTheta *sigtheta =
411 (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
412 if (!sigtheta)
413 {
414 cout << "Object 'SigmaTheta' not found" << endl;
415 return;
416 }
417 TH2D *fHSigTh = sigtheta->GetSigmaTheta();
418 TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
419 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
420
421 TString outNameSigTh = "SigmaTheta_";
422 outNameSigTh += typeData;
423 outNameSigTh += ".root";
424
425 TFile outfile(outNameSigTh, "RECREATE");
426 fHSigTh->Write();
427 fHSigPixTh->Write();
428 fHDifPixTh->Write();
429 outfile.Close();
430
431 cout << "File " << outNameSigTh << " was written out" << endl;
432 }
433
434
435 cout << "Macro AnalyseCT1 : End of Job 1" << endl;
436 cout << "===================================================" << endl;
437 }
438 //---------------------------------------------------------------------
439
440 //---------------------------------------------------------------------
441 // Job 2
442 //======
443 // read ON data file
444
445 // - produce the 2D-histogram "sigmabar versus Theta" for ON data
446 // (to be used for the padding of the MC gamma data)
447
448 // - write a file of look-alike events (hadron matrix)
449 // (to be used later together with the corresponding MC gamma file
450 // for the g/h separation in the analysis of the ON data)
451
452 // - write a file of ON events (ON1)
453 // (after the standard cuts, before the g/h separation)
454 // (to be used together with the corresponding MC gamma file
455 // for the optimization of the g/h separation)
456
457 if (Job2)
458 {
459 cout << "=====================================================" << endl;
460 cout << "Macro AnalyseCT1 : Start of Job 2" << endl;
461
462 cout << "" << endl;
463 cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = "
464 << Job2 << ", " << WHistON << ", " << WLookON << ", "
465 << WON1 << endl;
466
467 TString typeData = "ON";
468 cout << "typeData = " << typeData << endl;
469
470 MTaskList tliston;
471 MParList pliston;
472 pliston.AddToList(&tliston);
473
474
475 //::::::::::::::::::::::::::::::::::::::::::
476
477 MBinning binssize("BinningSize");
478 binssize.SetEdgesLog(50, 10, 1.0e5);
479 pliston.AddToList(&binssize);
480
481 MBinning binsdistc("BinningDist");
482 binsdistc.SetEdges(50, 0, 1.4);
483 pliston.AddToList(&binsdistc);
484
485 MBinning binswidth("BinningWidth");
486 binswidth.SetEdges(50, 0, 1.0);
487 pliston.AddToList(&binswidth);
488
489 MBinning binslength("BinningLength");
490 binslength.SetEdges(50, 0, 1.0);
491 pliston.AddToList(&binslength);
492
493 MBinning binsalpha("BinningAlpha");
494 binsalpha.SetEdges(100, -100, 100);
495 pliston.AddToList(&binsalpha);
496
497 MBinning binsasym("BinningAsym");
498 binsasym.SetEdges(50, -1.5, 1.5);
499 pliston.AddToList(&binsasym);
500
501 MBinning binsht("BinningHeadTail");
502 binsht.SetEdges(50, -1.5, 1.5);
503 pliston.AddToList(&binsht);
504
505 MBinning binsm3l("BinningM3Long");
506 binsm3l.SetEdges(50, -1.5, 1.5);
507 pliston.AddToList(&binsm3l);
508
509 MBinning binsm3t("BinningM3Trans");
510 binsm3t.SetEdges(50, -1.5, 1.5);
511 pliston.AddToList(&binsm3t);
512
513
514 //.....
515 MBinning binsb("BinningSigmabar");
516 binsb.SetEdges( 100, 0.0, 5.0);
517 pliston.AddToList(&binsb);
518
519 MBinning binth("BinningTheta");
520 binth.SetEdges( 70, -0.5, 69.5);
521 pliston.AddToList(&binth);
522
523 MBinning binsdiff("BinningDiffsigma2");
524 binsdiff.SetEdges(100, -5.0, 20.0);
525 pliston.AddToList(&binsdiff);
526 //::::::::::::::::::::::::::::::::::::::::::
527
528
529 //-------------------------------------------
530 // create the tasks which should be executed
531 //
532
533 filename = onfile;
534 RName = "ReadCT1data";
535 MCT1ReadPreProc read(filename, RName);
536
537 MSelBasic selbasic;
538
539 MBlindPixelCalc blind;
540 blind.SetUseInterpolation();
541 blind.SetUseBlindPixels();
542 // blind.SetUseCetralPixel();
543
544 // create an object of class MSigmabar,
545 // because class MHSigmaTheta will look for it
546 MSigmabar sigbar;
547 pliston.AddToList(&sigbar);
548 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
549 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
550
551 MImgCleanStd clean;
552
553 // calculate also extended image parameters
554 TString fHilName = "Hillas";
555 MHillasExt hext(fHilName);
556 pliston.AddToList(&hext);
557 MHillasExt *fHillas = &hext;
558
559 // name of output container for MHillasCalc
560 // = name of MHillas object
561 MHillasCalc hcalc(fHilName);
562
563 // name of output container for MHillasSrcCalc
564 // = name of MHillasSrc object
565 TString fHilSrcName = "HillasSrc";
566 MHillasSrc hilsrc(fHilSrcName);
567 MHillasSrc *fHillasSrc = &hilsrc;
568 pliston.AddToList(&hilsrc);
569 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
570
571 MSelStandard selstand(fHilName, fHilSrcName);
572
573 MFillH hfill1("MHHillas", fHilName);
574 MFillH hfill2("MHStarMap", fHilName);
575
576 TString nxt = "HillasExt";
577 MHHillasExt hhilext(nxt, "", fHilName);
578 pliston.AddToList(&hhilext);
579 TString namext = nxt;
580 namext += "[MHHillasExt]";
581 MFillH hfill3(namext, fHilSrcName);
582
583 MFillH hfill4("MHHillasSrc", fHilSrcName);
584
585 //+++++++++++++++++++++++++++++++++++++++++++++++++++
586 // prepare writing of look-alike events (for hadrons)
587
588 const char* mtxName = "MatrixHadrons";
589 Int_t howMany = 50000;
590 TString outName = "matrix_hadrons_";
591 outName += typeData;
592 outName += ".root";
593
594
595 cout << "" << endl;
596 cout << "========================================================" << endl;
597 cout << "Matrix for (hadrons)" << endl;
598 cout << "input file = " << filename<< endl;
599 cout << "matrix name = " << mtxName << endl;
600 cout << "max. no.of look-alike events = " << howMany << endl;
601 cout << "name of output root file = " << outName << endl;
602 cout << "" << endl;
603
604
605 // matrix limitation for look-alike events (approximate number)
606 MFEventSelector selector("", "", RName);
607 selector.SetNumSelectEvts(howMany);
608 MFilterList flist;
609 flist.AddToList(&selector);
610
611 //
612 // --- setup the matrix and add it to the parlist
613 //
614 MHMatrix *pmatrix = new MHMatrix(mtxName);
615 MHMatrix& matrix = *pmatrix;
616
617 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
618 matrix.AddColumn("MSigmabar.fSigmabar");
619 matrix.AddColumn("log10(Hillas.fSize)");
620 matrix.AddColumn("HillasSrc.fDist");
621 matrix.AddColumn("Hillas.fWidth");
622 matrix.AddColumn("Hillas.fLength");
623 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
624 matrix.AddColumn("abs(Hillas.fM3Long)");
625 matrix.AddColumn("Hillas.fConc");
626 matrix.AddColumn("Hillas.fConc1");
627
628 pliston.AddToList(pmatrix);
629
630 MFillH fillmatg(mtxName);
631 fillmatg.SetFilter(&flist);
632
633
634 if (WON1)
635 {
636 TString outNameImage = typeData;
637 outNameImage += "1.root";
638 MWriteRootFile write(outNameImage);
639 write.AddContainer("MSigmabar", "Events");
640 write.AddContainer("Hillas", "Events");
641 write.AddContainer("MMcEvt", "Events");
642 write.AddContainer("HillasSrc", "Events");
643 write.AddContainer("MRawRunHeader", "RunHeaders");
644 //write.AddContainer("MMcRunHeader","RunHeaders");
645 write.AddContainer("MSrcPosCam", "RunHeaders");
646 }
647
648 //+++++++++++++++++++++++++++++++++++++++++++++++++++
649
650
651 //*****************************
652 // tasks to be executed
653
654 tliston.AddToList(&read);
655
656 tliston.AddToList(&selbasic);
657 tliston.AddToList(&blind);
658 tliston.AddToList(&fillsigtheta);
659 tliston.AddToList(&clean);
660
661 tliston.AddToList(&hcalc);
662 tliston.AddToList(&csrc1);
663
664 tliston.AddToList(&selstand);
665 if (WON1)
666 tliston.AddToList(&write);
667
668 tliston.AddToList(&flist);
669 tliston.AddToList(&fillmatg);
670
671 tliston.AddToList(&hfill1);
672 tliston.AddToList(&hfill2);
673 tliston.AddToList(&hfill3);
674 tliston.AddToList(&hfill4);
675
676
677 //*****************************
678
679 //-------------------------------------------
680 // Execute event loop
681 //
682 MProgressBar bar;
683 MEvtLoop evtloop;
684 evtloop.SetParList(&pliston);
685 evtloop.SetProgressBar(&bar);
686
687 Int_t maxevents = 1000000000;
688 //Int_t maxevents = 1000;
689 if ( !evtloop.Eventloop(maxevents) )
690 return;
691
692 tliston.PrintStatistics(0, kTRUE);
693
694 //-------------------------------------------
695 // Display the histograms
696 //
697 TObject *c;
698 c = pliston.FindObject("MHHillas")->DrawClone();
699
700 c = pliston.FindObject("MHHillasSrc")->DrawClone();
701
702 //pliston.FindObject("MHHillasExt")->DrawClone();
703 c = pliston.FindObject(nxt)->DrawClone();
704
705 c = pliston.FindObject("MHStarMap")->DrawClone();
706
707
708 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
709 //
710 // ----------------------------------------------------------
711 // Definition of the reference sample of look-alike events
712 // (this is a subsample of the original sample)
713 // ----------------------------------------------------------
714 //
715 cout << "" << endl;
716 cout << "========================================================" << endl;
717 // Select a maximum of nmaxevts events from the sample of look-alike
718 // events. They will form the reference sample.
719 Int_t nmaxevts = 2000;
720
721 // target distribution for the variable in column refcolumn (start at 1);
722 // set refcolumn negative if distribution is not to be changed
723 Int_t refcolumn = 1;
724 Int_t nbins = 10;
725 Float_t frombin = 0.5;
726 Float_t tobin = 1.0;
727 TH1F *thsh = new TH1F("thsh","target distribution",
728 nbins, frombin, tobin);
729 thsh->SetXTitle("cos( \\Theta)");
730 thsh->SetYTitle("Counts");
731 Float_t dbin = (tobin-frombin)/nbins;
732 Float_t lbin = frombin +dbin*0.5;
733 for (Int_t j=1; j<=nbins; j++)
734 {
735 thsh->Fill(lbin,1.0);
736 lbin += dbin;
737 }
738
739 cout << "" << endl;
740 cout << "========================================================" << endl;
741 cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl;
742 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
743 << refcolumn << ", " << nmaxevts << endl;
744
745 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
746 {
747 cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
748 return;
749 }
750 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
751
752 //-------------------------------------------
753 // write out look-alike events (for hadrons)
754 //
755 if (WLookON)
756 {
757 TFile *writejens = new TFile(outName, "RECREATE", "");
758 matrix.Write();
759 cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file "
760 << outName << endl;
761
762 writejens->Close();
763 delete writejens;
764 }
765
766 //-------------------------------------------
767 // Write histograms onto a file
768 if (WHistON)
769 {
770 MHSigmaTheta *sigtheta =
771 (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
772 if (!sigtheta)
773 {
774 cout << "Object 'SigmaTheta' not found" << endl;
775 return;
776 }
777 TH2D *fHSigTh = sigtheta->GetSigmaTheta();
778 TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
779 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
780
781 TString outNameSigTh = "SigmaTheta_";
782 outNameSigTh += typeData;
783 outNameSigTh += ".root";
784
785 TFile outfile(outNameSigTh, "RECREATE");
786 fHSigTh->Write();
787 fHSigPixTh->Write();
788 fHDifPixTh->Write();
789 outfile.Close();
790
791 cout << "File " << outNameSigTh << " was written out" << endl;
792 }
793
794
795 cout << "Macro AnalyseCT1 : End of Job 2" << endl;
796 cout << "===================================================" << endl;
797 }
798
799
800 //---------------------------------------------------------------------
801 // Job 3
802 //======
803
804 // read MC gamma data
805
806 // - to pad them
807 // (using the 2D-histogram "sigmabar versus Theta" of the OFF data)
808
809 // - to write a file of look-alike events (gammas matrix)
810 // (to be used later together with the corresponding hadron file
811 // for the g/h separation in the analysis of the ON data)
812
813 // - to write a file of padded MC gamma events (MC1)
814 // (after the standard cuts, before the g/h separation)
815 // (to be used together with the corresponding hadron file
816 // for the optimization of the g/h separation)
817
818
819 if (Job3)
820 {
821 cout << "=====================================================" << endl;
822 cout << "Macro AnalyseCT1 : Start of Job 3" << endl;
823
824 cout << "" << endl;
825 cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = "
826 << Job3 << ", " << WLookMC << ", " << WMC1 << endl;
827
828
829 TString typeMC = "MC";
830 cout << "typeMC = " << typeMC << endl;
831
832 // use for padding sigmabar vs. Theta from ON or OFF data
833 TString typeData = "ON";
834 //TString typeData = "OFF";
835 cout << "typeData = " << typeData << endl;
836
837 //------------------------------------
838 // Get the histograms "2D-ThetaSigmabar"
839 // and "3D-ThetaPixSigma"
840 // and "3D-ThetaPixDiff"
841
842
843 TString outNameSigTh = "SigmaTheta_";
844 outNameSigTh += typeData;
845 outNameSigTh += ".root";
846
847
848 cout << "Reading in file " << outNameSigTh << endl;
849
850 TFile *infile = new TFile(outNameSigTh);
851 infile->ls();
852
853 TH2D *fHSigmaTheta =
854 (TH2D*) gROOT.FindObject("2D-ThetaSigmabar");
855 if (!fHSigmaTheta)
856 {
857 cout << "Object '2D-ThetaSigmabar' not found on root file" << endl;
858 return;
859 }
860 cout << "Object '2D-ThetaSigmabar' was read in" << endl;
861
862 TH3D *fHSigmaPixTheta =
863 (TH3D*) gROOT.FindObject("3D-ThetaPixSigma");
864 if (!fHSigmaPixTheta)
865 {
866 cout << "Object '3D-ThetaPixSigma' not found on root file" << endl;
867 return;
868 }
869 cout << "Object '3D-ThetaPixSigma' was read in" << endl;
870
871 TH3D *fHDiffPixTheta =
872 (TH3D*) gROOT.FindObject("3D-ThetaPixDiff");
873 if (!fHSigmaPixTheta)
874 {
875 cout << "Object '3D-ThetaPixDiff' not found on root file" << endl;
876 return;
877 }
878 cout << "Object '3D-ThetaPixDiff' was read in" << endl;
879
880 //------------------------------------
881
882 MTaskList tlist;
883 MParList plist;
884
885 plist.AddToList(&tlist);
886
887
888 //::::::::::::::::::::::::::::::::::::::::::
889
890 MBinning binssize("BinningSize");
891 binssize.SetEdgesLog(50, 10, 1.0e5);
892 plist.AddToList(&binssize);
893
894 MBinning binsdistc("BinningDist");
895 binsdistc.SetEdges(50, 0, 1.4);
896 plist.AddToList(&binsdistc);
897
898 MBinning binswidth("BinningWidth");
899 binswidth.SetEdges(50, 0, 1.0);
900 plist.AddToList(&binswidth);
901
902 MBinning binslength("BinningLength");
903 binslength.SetEdges(50, 0, 1.0);
904 plist.AddToList(&binslength);
905
906 MBinning binsalpha("BinningAlpha");
907 binsalpha.SetEdges(100, -100, 100);
908 plist.AddToList(&binsalpha);
909
910 MBinning binsasym("BinningAsym");
911 binsasym.SetEdges(50, -1.5, 1.5);
912 plist.AddToList(&binsasym);
913
914 MBinning binsht("BinningHeadTail");
915 binsht.SetEdges(50, -1.5, 1.5);
916 plist.AddToList(&binsht);
917
918 MBinning binsm3l("BinningM3Long");
919 binsm3l.SetEdges(50, -1.5, 1.5);
920 plist.AddToList(&binsm3l);
921
922 MBinning binsm3t("BinningM3Trans");
923 binsm3t.SetEdges(50, -1.5, 1.5);
924 plist.AddToList(&binsm3t);
925
926
927 //.....
928 MBinning binsb("BinningSigmabar");
929 binsb.SetEdges( 100, 0.0, 5.0);
930 plist.AddToList(&binsb);
931
932 MBinning binth("BinningTheta");
933 binth.SetEdges( 70, -0.5, 69.5);
934 plist.AddToList(&binth);
935
936 MBinning binsdiff("BinningDiffsigma2");
937 binsdiff.SetEdges(100, -5.0, 20.0);
938 plist.AddToList(&binsdiff);
939
940 //::::::::::::::::::::::::::::::::::::::::::
941
942 //-------------------------------------------
943 // create the tasks which should be executed
944 //
945 filename = mcfile;
946 readname = "ReadCT1MC";
947 MCT1ReadPreProc read(filename, readname);
948
949 MSelBasic selbasic;
950
951 MBlindPixelCalc blind;
952 blind.SetUseInterpolation();
953 blind.SetUseBlindPixels();
954 // blind.SetUseCetralPixel();
955
956 // There are 2 options for Thomas Schweizer's padding
957 // fPadFlag = 1 get Sigmabar from fHSigmaTheta
958 // and Sigma from fHDiffPixTheta
959 // fPadFlag = 2 get Sigma from fHSigmaPixTheta
960
961 MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)",
962 fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta);
963 padthomas.SetPadFlag(1);
964
965 //...........................................
966
967 MImgCleanStd clean; //(0, 0);
968
969 // calculate also extended image parameters
970 TString fHilName = "Hillas";
971 MHillasExt hext(fHilName);
972 plist.AddToList(&hext);
973 MHillasExt *fHillas = &hext;
974
975 // name of output container for MHillasCalc
976 // = name of MHillas object
977 MHillasCalc hcalc(fHilName);
978
979 // name of output container for MHillasSrcCalc
980 // = name of MHillasSrc object
981 TString fHilSrcName = "HillasSrc";
982 MHillasSrc hilsrc(fHilSrcName);
983 MHillasSrc *fHillasSrc = &hilsrc;
984 plist.AddToList(&hilsrc);
985 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
986
987 MSelStandard selstand(fHilName, fHilSrcName);
988
989 MFillH hfill1("MHHillas", fHilName);
990 MFillH hfill2("MHStarMap", fHilName);
991
992 TString nxt = "HillasExt";
993 MHHillasExt hhilext(nxt, "", fHilName);
994 plist.AddToList(&hhilext);
995 TString namext = nxt;
996 namext += "[MHHillasExt]";
997 MFillH hfill3(namext, fHilSrcName);
998
999 MFillH hfill4("MHHillasSrc", fHilSrcName);
1000
1001 //+++++++++++++++++++++++++++++++++++++++++++++++++++
1002 // prepare writing of look-alike events (for gammas)
1003
1004 const char* mtxName = "MatrixGammas";
1005 Int_t howMany = 50000;
1006
1007 TString outName = "matrix_gammas_";
1008 outName += typeMC;
1009 outName += "_";
1010 outName += typeData;
1011 outName += ".root";
1012
1013
1014 cout << "" << endl;
1015 cout << "========================================================" << endl;
1016 cout << "Matrix for (gammas)" << endl;
1017 cout << "input file = " << filename<< endl;
1018 cout << "matrix name = " << mtxName << endl;
1019 cout << "max. no.of look-alike events = " << howMany << endl;
1020 cout << "name of output root file = " << outName << endl;
1021 cout << "" << endl;
1022
1023
1024 // matrix limitation for look-alike events (approximate number)
1025 MFEventSelector selector("", "", readname);
1026 selector.SetNumSelectEvts(howMany);
1027 MFilterList flist;
1028 flist.AddToList(&selector);
1029
1030 //
1031 // --- setup the matrix and add it to the parlist
1032 //
1033 MHMatrix *pmatrix = new MHMatrix(mtxName);
1034 MHMatrix& matrix = *pmatrix;
1035
1036 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
1037 matrix.AddColumn("MSigmabar.fSigmabar");
1038 matrix.AddColumn("log10(Hillas.fSize)");
1039 matrix.AddColumn("HillasSrc.fDist");
1040 matrix.AddColumn("Hillas.fWidth");
1041 matrix.AddColumn("Hillas.fLength");
1042 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)");
1043 matrix.AddColumn("abs(Hillas.fM3Long)");
1044 matrix.AddColumn("Hillas.fConc");
1045 matrix.AddColumn("Hillas.fConc1");
1046
1047 plist.AddToList(pmatrix);
1048
1049 MFillH fillmatg(mtxName);
1050 fillmatg.SetFilter(&flist);
1051
1052
1053 if (WMC1)
1054 {
1055 TString outNameImage = typeMC;
1056 outNameImage += "_";
1057 outNameImage += typeData;
1058 outNameImage += "1.root";
1059 MWriteRootFile write(outNameImage);
1060 write.AddContainer("MSigmabar", "Events");
1061 write.AddContainer("Hillas", "Events");
1062 write.AddContainer("MMcEvt", "Events");
1063 write.AddContainer("HillasSrc", "Events");
1064 write.AddContainer("MRawRunHeader", "RunHeaders");
1065 //write.AddContainer("MMcRunHeader","RunHeaders");
1066 write.AddContainer("MSrcPosCam", "RunHeaders");
1067 }
1068
1069
1070 //+++++++++++++++++++++++++++++++++++++++++++++++++++
1071
1072
1073 //*****************************
1074 // tasks to be executed
1075
1076 tlist.AddToList(&read);
1077
1078 tlist.AddToList(&selbasic);
1079 tlist.AddToList(&blind);
1080 tlist.AddToList(&padthomas);
1081 tlist.AddToList(&clean);
1082
1083 tlist.AddToList(&hcalc);
1084 tlist.AddToList(&csrc1);
1085
1086 tlist.AddToList(&selstand);
1087 if (WMC1)
1088 tlist.AddToList(&write);
1089
1090 tlist.AddToList(&flist);
1091 tlist.AddToList(&fillmatg);
1092
1093 tlist.AddToList(&hfill1);
1094 tlist.AddToList(&hfill2);
1095 tlist.AddToList(&hfill3);
1096 tlist.AddToList(&hfill4);
1097
1098
1099 //*****************************
1100
1101
1102 //-------------------------------------------
1103 // Execute event loop
1104 //
1105 MProgressBar bar;
1106 MEvtLoop evtloop;
1107 evtloop.SetParList(&plist);
1108 evtloop.SetProgressBar(&bar);
1109
1110 Int_t maxevents = 1000000000;
1111 //Int_t maxevents = 100;
1112 if ( !evtloop.Eventloop(maxevents) )
1113 return;
1114
1115 tlist.PrintStatistics(0, kTRUE);
1116
1117 //-------------------------------------------
1118 // Display the histograms
1119 //
1120 plist.FindObject("MHHillas")->DrawClone();
1121 plist.FindObject("MHHillasSrc")->DrawClone();
1122
1123 //plist.FindObject("MHHillasExt")->DrawClone();
1124 plist.FindObject(nxt)->DrawClone();
1125
1126 plist.FindObject("MHStarMap")->DrawClone();
1127
1128
1129 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1130 //
1131 // ----------------------------------------------------------
1132 // Definition of the reference sample of look-alike events
1133 // (this is a subsample of the original sample)
1134 // ----------------------------------------------------------
1135 //
1136 cout << "" << endl;
1137 cout << "========================================================" << endl;
1138 // Select a maximum of nmaxevts events from the sample of look-alike
1139 // events. They will form the reference sample.
1140 Int_t nmaxevts = 2000;
1141
1142 // target distribution for the variable in column refcolumn;
1143 // set refcolumn negative if distribution is not to be changed
1144 Int_t refcolumn = 1;
1145 Int_t nbins = 10;
1146 Float_t frombin = 0.5;
1147 Float_t tobin = 1.0;
1148 TH1F *thsh = new TH1F("thsh","target distribution",
1149 nbins, frombin, tobin);
1150 Float_t dbin = (tobin-frombin)/nbins;
1151 Float_t lbin = frombin +dbin*0.5;
1152 for (Int_t j=1; j<=nbins; j++)
1153 {
1154 thsh->Fill(lbin,1.0);
1155 lbin += dbin;
1156 }
1157
1158 cout << "" << endl;
1159 cout << "========================================================" << endl;
1160 cout << "Macro AnalyseCT1 : define reference sample for gammas" << endl;
1161 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
1162 << refcolumn << ", " << nmaxevts << endl;
1163
1164 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
1165 {
1166 cout << "Macro AnalyseCT1 : Reference matrix for gammas cannot be defined" << endl;
1167 return;
1168 }
1169 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1170
1171 //-------------------------------------------
1172 // write out look-alike events (for gammas)
1173 //
1174 if (WLookMC)
1175 {
1176 TFile *writejens = new TFile(outName, "RECREATE", "");
1177 matrix.Write();
1178 cout << "Macro AnalyseCT1 : matrix of look-alike events for gammas written onto file "
1179 << outName << endl;
1180
1181 writejens->Close();
1182 delete writejens;
1183 }
1184
1185 cout << "Macro AnalyseCT1 : End of Job 3"
1186 << endl;
1187 cout << "========================================================="
1188 << endl;
1189 }
1190
1191
1192
1193 //---------------------------------------------------------------------
1194 // Job 4
1195 //======
1196
1197 // read OFF1 and MC1 data files
1198 //
1199 // - produce Neyman-Pearson plot
1200
1201 if (Job4)
1202 {
1203 cout << "=====================================================" << endl;
1204 cout << "Macro AnalyseCT1 : Start of Job 4" << endl;
1205
1206 cout << "" << endl;
1207 cout << "Macro AnalyseCT1 : Job4 = " << Job4 << endl;
1208
1209
1210 TString typeMC = "MC";
1211 cout << "typeMC = " << typeMC << endl;
1212
1213 // use hadron matrix from ON or OFF data
1214 TString typeData = "ON";
1215 //TString typeData = "OFF";
1216 cout << "typeData = " << typeData << endl;
1217
1218 //*************************************************************************
1219 // read in matrices of look-alike events
1220
1221 const char* mtxName = "MatrixGammas";
1222
1223
1224 TString outName = "matrix_gammas_";
1225 outName += typeMC;
1226 outName += "_";
1227 outName += typeData;
1228 outName += ".root";
1229
1230 cout << "" << endl;
1231 cout << "========================================================" << endl;
1232 cout << "Get matrix for (gammas)" << endl;
1233 cout << "matrix name = " << mtxName << endl;
1234 cout << "name of root file = " << outName << endl;
1235 cout << "" << endl;
1236
1237
1238 // read in the object with the name 'mtxName' from file 'outName'
1239 //
1240 TFile *fileg = new TFile(outName);
1241
1242 MHMatrix matrixg;
1243 matrixg.Read(mtxName);
1244 pmatrixg = &matrixg;
1245
1246 delete fileg;
1247
1248
1249 //*****************************************************************
1250
1251 const char* mtxName = "MatrixHadrons";
1252
1253 TString outName = "matrix_hadrons_";
1254 outName += typeData;
1255 outName += ".root";
1256
1257 cout << "" << endl;
1258 cout << "========================================================" << endl;
1259 cout << " Get matrix for (hadrons)" << endl;
1260 cout << "matrix name = " << mtxName << endl;
1261 cout << "name of root file = " << outName << endl;
1262 cout << "" << endl;
1263
1264
1265 // read in the object with the name mtxName
1266 //
1267 TFile *fileh = new TFile(outName);
1268
1269 MHMatrix matrixh;
1270 matrixh.Read(mtxName);
1271 pmatrixh = &matrixh;
1272
1273 delete fileh;
1274
1275 //*****************************************************************
1276
1277 MTaskList tliston;
1278 MParList pliston;
1279 pliston.AddToList(&tliston);
1280
1281 pliston.AddToList(pmatrixg);
1282 pliston.AddToList(pmatrixh);
1283
1284 //::::::::::::::::::::::::::::::::::::::::::
1285
1286 MBinning binssize("BinningSize");
1287 binssize.SetEdgesLog(50, 10, 1.0e5);
1288 pliston.AddToList(&binssize);
1289
1290 MBinning binsdistc("BinningDist");
1291 binsdistc.SetEdges(50, 0, 1.4);
1292 pliston.AddToList(&binsdistc);
1293
1294 MBinning binswidth("BinningWidth");
1295 binswidth.SetEdges(50, 0, 1.0);
1296 pliston.AddToList(&binswidth);
1297
1298 MBinning binslength("BinningLength");
1299 binslength.SetEdges(50, 0, 1.0);
1300 pliston.AddToList(&binslength);
1301
1302 MBinning binsalpha("BinningAlpha");
1303 binsalpha.SetEdges(100, -100, 100);
1304 pliston.AddToList(&binsalpha);
1305
1306 MBinning binsasym("BinningAsym");
1307 binsasym.SetEdges(50, -1.5, 1.5);
1308 pliston.AddToList(&binsasym);
1309
1310 MBinning binsht("BinningHeadTail");
1311 binsht.SetEdges(50, -1.5, 1.5);
1312 pliston.AddToList(&binsht);
1313
1314 MBinning binsm3l("BinningM3Long");
1315 binsm3l.SetEdges(50, -1.5, 1.5);
1316 pliston.AddToList(&binsm3l);
1317
1318 MBinning binsm3t("BinningM3Trans");
1319 binsm3t.SetEdges(50, -1.5, 1.5);
1320 pliston.AddToList(&binsm3t);
1321
1322
1323 //.....
1324 MBinning binsb("BinningSigmabar");
1325 binsb.SetEdges( 100, 0.0, 5.0);
1326 pliston.AddToList(&binsb);
1327
1328 MBinning binth("BinningTheta");
1329 binth.SetEdges( 70, -0.5, 69.5);
1330 pliston.AddToList(&binth);
1331
1332 MBinning binsdiff("BinningDiffsigma2");
1333 binsdiff.SetEdges(100, -5.0, 20.0);
1334 pliston.AddToList(&binsdiff);
1335 //::::::::::::::::::::::::::::::::::::::::::
1336
1337
1338 //-------------------------------------------
1339 // create the tasks which should be executed
1340 //
1341
1342 TString filenameData = typeData;
1343 filenameData += "1.root";
1344
1345 TString filenameMC = typeMC;
1346 filenameMC += "_";
1347 filenameMC += typeData;
1348 filenameMC += "1.root";
1349
1350
1351 cout << "filenameData = " << filenameData << endl;
1352 cout << "filenameMC = " << filenameMC << endl;
1353
1354 MReadMarsFile read("Events", filenameMC);
1355 read.AddFile(filenameData);
1356 read.DisableAutoScheme();
1357
1358 //.......................................................................
1359 // g/h separation
1360
1361 MMultiDimDistCalc ghsep;
1362 ghsep.SetUseNumRows(25);
1363 ghsep.SetUseKernelMethod(kFALSE);
1364
1365 //.......................................................................
1366
1367 // geometry is needed in MHHillas... classes
1368 MGeomCam *fGeom =
1369 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1370
1371 TString fHilName = "Hillas";
1372 TString fHilSrcName = "HillasSrc";
1373
1374 Float_t hadcut = 0.2;
1375 MSelFinal selfinalgh(fHilName, fHilSrcName);
1376 selfinalgh.SetHadronnessCut(hadcut);
1377 selfinalgh.SetAlphaCut(100.0);
1378
1379 MFillH fillh("MHHadronness");
1380
1381 MSelFinal selfinal(fHilName, fHilSrcName);
1382 selfinal.SetHadronnessCut(hadcut);
1383 selfinal.SetAlphaCut(20.0);
1384
1385 MFillH hfill1("MHHillas", fHilName);
1386 MFillH hfill2("MHStarMap", fHilName);
1387
1388 TString nxt = "HillasExt";
1389 MHHillasExt hhilext(nxt, "", fHilName);
1390 pliston.AddToList(&hhilext);
1391 TString namext = nxt;
1392 namext += "[MHHillasExt]";
1393 MFillH hfill3(namext, fHilSrcName);
1394
1395 MFillH hfill4("MHHillasSrc", fHilSrcName);
1396
1397
1398 //*****************************
1399 // tasks to be executed
1400
1401 tliston.AddToList(&read);
1402
1403 tliston.AddToList(&ghsep);
1404 tliston.AddToList(&fillh);
1405
1406 tliston.AddToList(&selfinalgh);
1407 tliston.AddToList(&hfill1);
1408 tliston.AddToList(&hfill2);
1409 tliston.AddToList(&hfill3);
1410 tliston.AddToList(&hfill4);
1411
1412 tliston.AddToList(&selfinal);
1413
1414 //*****************************
1415
1416 //-------------------------------------------
1417 // Execute event loop
1418 //
1419 MProgressBar bar;
1420 MEvtLoop evtloop;
1421 evtloop.SetParList(&pliston);
1422 evtloop.SetProgressBar(&bar);
1423
1424 Int_t maxevents = 1000000000;
1425 //Int_t maxevents = 1000;
1426 if ( !evtloop.Eventloop(maxevents) )
1427 return;
1428
1429 tliston.PrintStatistics(0, kTRUE);
1430
1431
1432 //-------------------------------------------
1433 // Display the histograms
1434 //
1435 TObject *c;
1436
1437 c = pliston.FindObject("MHHadronness")->DrawClone();
1438 c->Print();
1439
1440 //c = pliston.FindObject("MHHillas")->DrawClone();
1441
1442 c = pliston.FindObject("MHHillasSrc")->DrawClone();
1443
1444 //pliston.FindObject("MHHillasExt")->DrawClone();
1445 //c = pliston.FindObject(nxt)->DrawClone();
1446
1447 c = pliston.FindObject("MHStarMap")->DrawClone();
1448
1449
1450
1451 cout << "Macro AnalyseCT1 : End of Job 4" << endl;
1452 cout << "===================================================" << endl;
1453 }
1454
1455
1456 //---------------------------------------------------------------------
1457 // Job 5
1458 //======
1459
1460 // - read in the matrices of look-alike events for gammas and hadrons
1461
1462 // then read MC1 data file
1463 //
1464 // - perform the g/h separation
1465 // - apply the final cuts
1466 //
1467 // - write a file of MC gamma events (MC2)
1468 // (after the g/h separation and after the final cuts)
1469
1470 if (Job5)
1471 {
1472 cout << "=====================================================" << endl;
1473 cout << "Macro AnalyseCT1 : Start of Job 5" << endl;
1474
1475 cout << "" << endl;
1476 cout << "Macro AnalyseCT1 : Job5, WMC2 = "
1477 << Job5 << ", " << WMC2 << endl;
1478
1479 TString typeInput = "MC";
1480 cout << "typeInput = " << typeInput << endl;
1481
1482 TString typeMC = "MC";
1483 cout << "typeMC = " << typeMC << endl;
1484
1485 // use hadron matrix from ON or OFF data
1486 TString typeData = "ON";
1487 //TString typeData = "OFF";
1488 cout << "typeData = " << typeData << endl;
1489
1490 //*************************************************************************
1491 // read in matrices of look-alike events
1492
1493 const char* mtxName = "MatrixGammas";
1494
1495 TString outName = "matrix_gammas_";
1496 outName += typeMC;
1497 outName += "_";
1498 outName += typeData;
1499 outName += ".root";
1500
1501
1502 cout << "" << endl;
1503 cout << "========================================================" << endl;
1504 cout << "Get matrix for (gammas)" << endl;
1505 cout << "matrix name = " << mtxName << endl;
1506 cout << "name of root file = " << outName << endl;
1507 cout << "" << endl;
1508
1509
1510 // read in the object with the name 'mtxName' from file 'outName'
1511 //
1512 TFile *fileg = new TFile(outName);
1513
1514 MHMatrix matrixg;
1515 matrixg.Read(mtxName);
1516 pmatrixg = &matrixg;
1517
1518 delete fileg;
1519
1520 //*****************************************************************
1521
1522 const char* mtxName = "MatrixHadrons";
1523
1524 TString outName = "matrix_hadrons_";
1525 outName += typeData;
1526 outName += ".root";
1527
1528
1529 cout << "" << endl;
1530 cout << "========================================================" << endl;
1531 cout << " Get matrix for (hadrons)" << endl;
1532 cout << "matrix name = " << mtxName << endl;
1533 cout << "name of root file = " << outName << endl;
1534 cout << "" << endl;
1535
1536
1537 // read in the object with the name mtxName
1538 //
1539 TFile *fileh = new TFile(outName);
1540
1541 MHMatrix matrixh;
1542 matrixh.Read(mtxName);
1543 pmatrixh = &matrixh;
1544
1545 delete fileh;
1546
1547 //*************************************************************************
1548
1549
1550 MTaskList tliston;
1551 MParList pliston;
1552 pliston.AddToList(&tliston);
1553
1554 pliston.AddToList(pmatrixg);
1555 pliston.AddToList(pmatrixh);
1556
1557 //::::::::::::::::::::::::::::::::::::::::::
1558
1559 MBinning binssize("BinningSize");
1560 binssize.SetEdgesLog(50, 10, 1.0e5);
1561 pliston.AddToList(&binssize);
1562
1563 MBinning binsdistc("BinningDist");
1564 binsdistc.SetEdges(50, 0, 1.4);
1565 pliston.AddToList(&binsdistc);
1566
1567 MBinning binswidth("BinningWidth");
1568 binswidth.SetEdges(50, 0, 1.0);
1569 pliston.AddToList(&binswidth);
1570
1571 MBinning binslength("BinningLength");
1572 binslength.SetEdges(50, 0, 1.0);
1573 pliston.AddToList(&binslength);
1574
1575 MBinning binsalpha("BinningAlpha");
1576 binsalpha.SetEdges(100, -100, 100);
1577 pliston.AddToList(&binsalpha);
1578
1579 MBinning binsasym("BinningAsym");
1580 binsasym.SetEdges(50, -1.5, 1.5);
1581 pliston.AddToList(&binsasym);
1582
1583 MBinning binsht("BinningHeadTail");
1584 binsht.SetEdges(50, -1.5, 1.5);
1585 pliston.AddToList(&binsht);
1586
1587 MBinning binsm3l("BinningM3Long");
1588 binsm3l.SetEdges(50, -1.5, 1.5);
1589 pliston.AddToList(&binsm3l);
1590
1591 MBinning binsm3t("BinningM3Trans");
1592 binsm3t.SetEdges(50, -1.5, 1.5);
1593 pliston.AddToList(&binsm3t);
1594
1595
1596 //.....
1597 MBinning binsb("BinningSigmabar");
1598 binsb.SetEdges( 100, 0.0, 5.0);
1599 pliston.AddToList(&binsb);
1600
1601 MBinning binth("BinningTheta");
1602 binth.SetEdges( 70, -0.5, 69.5);
1603 pliston.AddToList(&binth);
1604
1605 MBinning binsdiff("BinningDiffsigma2");
1606 binsdiff.SetEdges(100, -5.0, 20.0);
1607 pliston.AddToList(&binsdiff);
1608 //::::::::::::::::::::::::::::::::::::::::::
1609
1610
1611 //-------------------------------------------
1612 // create the tasks which should be executed
1613 //
1614
1615 TString filenameMC = typeInput;
1616 filenameMC += "_";
1617 filenameMC += typeData;
1618 filenameMC += "1.root";
1619
1620 readname = "ReadCT1MCdata";
1621 MReadMarsFile read("Events", filenameMC);
1622 read.DisableAutoScheme();
1623
1624
1625 //.......................................................................
1626 // g/h separation
1627
1628 MMultiDimDistCalc ghsep;
1629 ghsep.SetUseNumRows(25);
1630 ghsep.SetUseKernelMethod(kFALSE);
1631
1632
1633
1634 //.......................................................................
1635
1636 // geometry is needed in MHHillas... classes
1637 MGeomCam *fGeom =
1638 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1639
1640 TString fHilName = "Hillas";
1641 TString fHilSrcName = "HillasSrc";
1642
1643 Float_t hadcut = 0.2;
1644 MSelFinal selfinalgh(fHilName, fHilSrcName);
1645 selfinalgh.SetHadronnessCut(hadcut);
1646 selfinalgh.SetAlphaCut(100.0);
1647
1648 MFillH fillh("MHHadronness");
1649
1650 MSelFinal selfinal(fHilName, fHilSrcName);
1651 selfinal.SetHadronnessCut(hadcut);
1652 selfinal.SetAlphaCut(20.0);
1653
1654 if (WMC2)
1655 {
1656 TString outNameImage = typeInput;
1657 outNameImage += "_";
1658 outNameImage += typeData;
1659 outNameImage += "2.root";
1660 MWriteRootFile write(outNameImage);
1661 write.AddContainer("MHadronness", "Events");
1662 write.AddContainer("MSigmabar", "Events");
1663 write.AddContainer("Hillas", "Events");
1664 write.AddContainer("MMcEvt", "Events");
1665 write.AddContainer("HillasSrc", "Events");
1666 write.AddContainer("MRawRunHeader", "RunHeaders");
1667 //write.AddContainer("MMcRunHeader","RunHeaders");
1668 write.AddContainer("MSrcPosCam", "RunHeaders");
1669 }
1670
1671
1672 MFillH hfill1("MHHillas", fHilName);
1673 MFillH hfill2("MHStarMap", fHilName);
1674
1675 TString nxt = "HillasExt";
1676 MHHillasExt hhilext(nxt, "", fHilName);
1677 pliston.AddToList(&hhilext);
1678 TString namext = nxt;
1679 namext += "[MHHillasExt]";
1680 MFillH hfill3(namext, fHilSrcName);
1681
1682 MFillH hfill4("MHHillasSrc", fHilSrcName);
1683
1684
1685
1686 //*****************************
1687 // tasks to be executed
1688
1689 tliston.AddToList(&read);
1690
1691 tliston.AddToList(&ghsep);
1692 tliston.AddToList(&fillh);
1693
1694 tliston.AddToList(&selfinalgh);
1695 tliston.AddToList(&hfill1);
1696 tliston.AddToList(&hfill2);
1697 tliston.AddToList(&hfill3);
1698 tliston.AddToList(&hfill4);
1699
1700 tliston.AddToList(&selfinal);
1701 if (WMC2)
1702 tliston.AddToList(&write);
1703
1704 //*****************************
1705
1706 //-------------------------------------------
1707 // Execute event loop
1708 //
1709 MProgressBar bar;
1710 MEvtLoop evtloop;
1711 evtloop.SetParList(&pliston);
1712 evtloop.SetProgressBar(&bar);
1713
1714 Int_t maxevents = 1000000000;
1715 //Int_t maxevents = 1000;
1716 if ( !evtloop.Eventloop(maxevents) )
1717 return;
1718
1719 tliston.PrintStatistics(0, kTRUE);
1720
1721
1722 //-------------------------------------------
1723 // Display the histograms
1724 //
1725 TObject *c;
1726
1727 c = pliston.FindObject("MHHadronness")->DrawClone();
1728 c->Print();
1729
1730 c = pliston.FindObject("MHHillas")->DrawClone();
1731
1732 c = pliston.FindObject("MHHillasSrc")->DrawClone();
1733
1734 //pliston.FindObject("MHHillasExt")->DrawClone();
1735 c = pliston.FindObject(nxt)->DrawClone();
1736
1737 c = pliston.FindObject("MHStarMap")->DrawClone();
1738
1739
1740
1741 cout << "Macro AnalyseCT1 : End of Job 5" << endl;
1742 cout << "===================================================" << endl;
1743 }
1744
1745
1746 //---------------------------------------------------------------------
1747 // Job 6
1748 //======
1749
1750 // - read in the matrices of look-alike events for gammas and hadrons
1751
1752 // then read ON1 data file
1753 //
1754 // - perform the g/h separation
1755 // - apply the final cuts
1756 //
1757 // - write a file of ON events (ON2)
1758 // (after the g/h separation and after the final cuts)
1759 // (to be used for the flux calculation)
1760 //
1761 // - do the flux calculation
1762
1763 if (Job6)
1764 {
1765 cout << "=====================================================" << endl;
1766 cout << "Macro AnalyseCT1 : Start of Job 6" << endl;
1767
1768 cout << "" << endl;
1769 cout << "Macro AnalyseCT1 : Job6, WON2 = "
1770 << Job6 << ", " << WON2 << endl;
1771
1772 TString typeInput = "ON";
1773 cout << "typeInput = " << typeInput << endl;
1774
1775 TString typeMC = "MC";
1776 cout << "typeMC = " << typeMC << endl;
1777
1778 // use hadron matrix from ON or OFF data
1779 TString typeData = "ON";
1780 //TString typeData = "OFF";
1781 cout << "typeData = " << typeData << endl;
1782
1783
1784 //*************************************************************************
1785 // read in matrices of look-alike events
1786
1787 const char* mtxName = "MatrixGammas";
1788
1789 TString outName = "matrix_gammas_";
1790 outName += typeMC;
1791 outName += "_";
1792 outName += typeData;
1793 outName += ".root";
1794
1795
1796 cout << "" << endl;
1797 cout << "========================================================" << endl;
1798 cout << "Get matrix for (gammas)" << endl;
1799 cout << "matrix name = " << mtxName << endl;
1800 cout << "name of root file = " << outName << endl;
1801 cout << "" << endl;
1802
1803
1804 // read in the object with the name 'mtxName' from file 'outName'
1805 //
1806 TFile *fileg = new TFile(outName);
1807
1808 MHMatrix matrixg;
1809 matrixg.Read(mtxName);
1810 pmatrixg = &matrixg;
1811
1812 delete fileg;
1813
1814 //*****************************************************************
1815
1816 const char* mtxName = "MatrixHadrons";
1817
1818 TString outName = "matrix_hadrons_";
1819 outName += typeData;
1820 outName += ".root";
1821
1822
1823 cout << "" << endl;
1824 cout << "========================================================" << endl;
1825 cout << " Get matrix for (hadrons)" << endl;
1826 cout << "matrix name = " << mtxName << endl;
1827 cout << "name of root file = " << outName << endl;
1828 cout << "" << endl;
1829
1830
1831 // read in the object with the name mtxName
1832 //
1833 TFile *fileh = new TFile(outName);
1834
1835 MHMatrix matrixh;
1836 matrixh.Read(mtxName);
1837 pmatrixh = &matrixh;
1838
1839 delete fileh;
1840
1841 //*************************************************************************
1842
1843
1844 MTaskList tliston;
1845 MParList pliston;
1846 pliston.AddToList(&tliston);
1847
1848 pliston.AddToList(pmatrixg);
1849 pliston.AddToList(pmatrixh);
1850
1851 //::::::::::::::::::::::::::::::::::::::::::
1852
1853 MBinning binssize("BinningSize");
1854 binssize.SetEdgesLog(50, 10, 1.0e5);
1855 pliston.AddToList(&binssize);
1856
1857 MBinning binsdistc("BinningDist");
1858 binsdistc.SetEdges(50, 0, 1.4);
1859 pliston.AddToList(&binsdistc);
1860
1861 MBinning binswidth("BinningWidth");
1862 binswidth.SetEdges(50, 0, 1.0);
1863 pliston.AddToList(&binswidth);
1864
1865 MBinning binslength("BinningLength");
1866 binslength.SetEdges(50, 0, 1.0);
1867 pliston.AddToList(&binslength);
1868
1869 MBinning binsalpha("BinningAlpha");
1870 binsalpha.SetEdges(100, -100, 100);
1871 pliston.AddToList(&binsalpha);
1872
1873 MBinning binsasym("BinningAsym");
1874 binsasym.SetEdges(50, -1.5, 1.5);
1875 pliston.AddToList(&binsasym);
1876
1877 MBinning binsht("BinningHeadTail");
1878 binsht.SetEdges(50, -1.5, 1.5);
1879 pliston.AddToList(&binsht);
1880
1881 MBinning binsm3l("BinningM3Long");
1882 binsm3l.SetEdges(50, -1.5, 1.5);
1883 pliston.AddToList(&binsm3l);
1884
1885 MBinning binsm3t("BinningM3Trans");
1886 binsm3t.SetEdges(50, -1.5, 1.5);
1887 pliston.AddToList(&binsm3t);
1888
1889
1890 //.....
1891 MBinning binsb("BinningSigmabar");
1892 binsb.SetEdges( 100, 0.0, 5.0);
1893 pliston.AddToList(&binsb);
1894
1895 MBinning binth("BinningTheta");
1896 binth.SetEdges( 70, -0.5, 69.5);
1897 pliston.AddToList(&binth);
1898
1899 MBinning binsdiff("BinningDiffsigma2");
1900 binsdiff.SetEdges(100, -5.0, 20.0);
1901 pliston.AddToList(&binsdiff);
1902 //::::::::::::::::::::::::::::::::::::::::::
1903
1904
1905 //-------------------------------------------
1906 // create the tasks which should be executed
1907 //
1908
1909 TString filenameData = typeInput;
1910 filenameData += "1.root";
1911
1912 readname = "ReadCT1ONdata";
1913 MReadMarsFile read("Events", filenameData);
1914 read.DisableAutoScheme();
1915
1916
1917 //.......................................................................
1918 // g/h separation
1919
1920 MMultiDimDistCalc ghsep;
1921 ghsep.SetUseNumRows(25);
1922 ghsep.SetUseKernelMethod(kFALSE);
1923
1924
1925
1926 //.......................................................................
1927
1928 // geometry is needed in MHHillas... classes
1929 MGeomCam *fGeom =
1930 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
1931
1932 TString fHilName = "Hillas";
1933 TString fHilSrcName = "HillasSrc";
1934
1935 Float_t hadcut = 0.2;
1936 MSelFinal selfinalgh(fHilName, fHilSrcName);
1937 selfinalgh.SetHadronnessCut(hadcut);
1938 selfinalgh.SetAlphaCut(100.0);
1939
1940 MFillH fillh("MHHadronness");
1941
1942 MSelFinal selfinal(fHilName, fHilSrcName);
1943 selfinal.SetHadronnessCut(hadcut);
1944 selfinal.SetAlphaCut(20.0);
1945
1946 if (WON2)
1947 {
1948 TString outNameImage = typeInput;
1949 outNameImage += "_";
1950 outNameImage += typeData;
1951 outNameImage += "2.root";
1952 MWriteRootFile write(outNameImage);
1953 write.AddContainer("MHadronness", "Events");
1954 write.AddContainer("MSigmabar", "Events");
1955 write.AddContainer("Hillas", "Events");
1956 write.AddContainer("MMcEvt", "Events");
1957 write.AddContainer("HillasSrc", "Events");
1958 write.AddContainer("MRawRunHeader", "RunHeaders");
1959 //write.AddContainer("MMcRunHeader","RunHeaders");
1960 write.AddContainer("MSrcPosCam", "RunHeaders");
1961 }
1962
1963
1964 MFillH hfill1("MHHillas", fHilName);
1965 MFillH hfill2("MHStarMap", fHilName);
1966
1967 TString nxt = "HillasExt";
1968 MHHillasExt hhilext(nxt, "", fHilName);
1969 pliston.AddToList(&hhilext);
1970 TString namext = nxt;
1971 namext += "[MHHillasExt]";
1972 MFillH hfill3(namext, fHilSrcName);
1973
1974 MFillH hfill4("MHHillasSrc", fHilSrcName);
1975
1976
1977
1978 //*****************************
1979 // tasks to be executed
1980
1981 tliston.AddToList(&read);
1982
1983 tliston.AddToList(&ghsep);
1984 tliston.AddToList(&fillh);
1985
1986 tliston.AddToList(&selfinalgh);
1987
1988 tliston.AddToList(&hfill1);
1989 tliston.AddToList(&hfill2);
1990 tliston.AddToList(&hfill3);
1991 tliston.AddToList(&hfill4);
1992
1993 tliston.AddToList(&selfinal);
1994 if (WON2)
1995 tliston.AddToList(&write);
1996
1997
1998
1999 //*****************************
2000
2001 //-------------------------------------------
2002 // Execute event loop
2003 //
2004 MProgressBar bar;
2005 MEvtLoop evtloop;
2006 evtloop.SetParList(&pliston);
2007 evtloop.SetProgressBar(&bar);
2008
2009 Int_t maxevents = 1000000000;
2010 //Int_t maxevents = 1000;
2011 if ( !evtloop.Eventloop(maxevents) )
2012 return;
2013
2014 tliston.PrintStatistics(0, kTRUE);
2015
2016
2017 //-------------------------------------------
2018 // Display the histograms
2019 //
2020 TObject *c;
2021
2022 c = pliston.FindObject("MHHadronness")->DrawClone();
2023 c->Print();
2024
2025 //c = pliston.FindObject("MHHillas")->DrawClone();
2026
2027 c = pliston.FindObject("MHHillasSrc")->DrawClone();
2028
2029 ////pliston.FindObject("MHHillasExt")->DrawClone();
2030 //c = pliston.FindObject(nxt)->DrawClone();
2031
2032 c = pliston.FindObject("MHStarMap")->DrawClone();
2033
2034
2035
2036 cout << "Macro AnalyseCT1 : End of Job 6" << endl;
2037 cout << "=======================================================" << endl;
2038 }
2039 //---------------------------------------------------------------------
2040}
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
Note: See TracBrowser for help on using the repository browser.