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

Last change on this file since 1767 was 1767, checked in by wittek, 22 years ago
*** empty log message ***
File size: 22.2 KB
Line 
1void AnalyseCT1()
2{
3 const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc";
4 const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc";
5 //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6";
6 //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc";
7 const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc";
8 const char *filename;
9
10 Bool_t Data = kFALSE; // loop over data ?
11 Bool_t WLook = kFALSE; // write out look-alike events for hadrons ?
12 Bool_t WHist = kFALSE; // write out histograms for padding ?
13
14 Bool_t MC = kTRUE; // loop over MC gamma data ?
15 Bool_t RHist = kTRUE; // read in histograms for padding ?
16 Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ?
17
18
19 cout << "" << endl;
20 cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ", "
21 << WHist << endl;
22
23 cout << " RHist, MC = " << RHist << ", "
24 << MC << endl;
25 cout << "" << endl;
26
27 //---------------------------------------------------------------------
28 // start of section for ON data
29 // (in the CT1 analysis the ON data are also used as a sample representing
30 // the hadrons for the g/h separation)
31
32 // read ON data file
33
34 // - to produce the 2D-histogram "sigmabar versus Theta"
35 // (to be used for the padding of the MC gamma data)
36
37 // - to write a file of look-alike events (for hadrons)
38 // (to be used later together with the corresponding MC gamma file
39 // for the g/h separation in the analysis of the ON data)
40
41 // - to write a file of ON events (after the standard cuts,
42 // before the g/h separation)
43 // (to be used together with the corresponding MC gamma file
44 // for the optimization of the g/h separation)
45
46 if (Data)
47 {
48 cout << "=====================================================" << endl;
49 cout << "Macro AnalyseCT1 : Start of section for ON data" << endl;
50
51 MTaskList tliston;
52 MParList pliston;
53 pliston.AddToList(&tliston);
54
55
56 //::::::::::::::::::::::::::::::::::::::::::
57
58 MBinning binssize("BinningSize");
59 binssize.SetEdgesLog(50, 10, 1.0e5);
60 pliston.AddToList(&binssize);
61
62 MBinning binsdistc("BinningDist");
63 binsdistc.SetEdges(50, 0, 1.4);
64 pliston.AddToList(&binsdistc);
65
66 MBinning binswidth("BinningWidth");
67 binswidth.SetEdges(50, 0, 1.0);
68 pliston.AddToList(&binswidth);
69
70 MBinning binslength("BinningLength");
71 binslength.SetEdges(50, 0, 1.0);
72 pliston.AddToList(&binslength);
73
74 MBinning binsalpha("BinningAlpha");
75 binsalpha.SetEdges(100, -100, 100);
76 pliston.AddToList(&binsalpha);
77
78 MBinning binsasym("BinningAsym");
79 binsasym.SetEdges(50, -1.5, 1.5);
80 pliston.AddToList(&binsasym);
81
82 MBinning binsht("BinningHeadTail");
83 binsht.SetEdges(50, -1.5, 1.5);
84 pliston.AddToList(&binsht);
85
86 MBinning binsm3l("BinningM3Long");
87 binsm3l.SetEdges(50, -1.5, 1.5);
88 pliston.AddToList(&binsm3l);
89
90 MBinning binsm3t("BinningM3Trans");
91 binsm3t.SetEdges(50, -1.5, 1.5);
92 pliston.AddToList(&binsm3t);
93
94
95 //.....
96 MBinning binsb("BinningSigmabar");
97 binsb.SetEdges( 100, 0.0, 5.0);
98 pliston.AddToList(&binsb);
99
100 MBinning binth("BinningTheta");
101 binth.SetEdges( 70, -0.5, 69.5);
102 pliston.AddToList(&binth);
103
104 MBinning binsdiff("BinningDiffsigma2");
105 binsdiff.SetEdges(100, -5.0, 20.0);
106 pliston.AddToList(&binsdiff);
107 //::::::::::::::::::::::::::::::::::::::::::
108
109
110 //-------------------------------------------
111 // create the tasks which should be executed
112 //
113
114 filename = onfile;
115 readname = "ReadCT1data";
116 MCT1ReadPreProc read(filename, readname);
117
118 MSelBasic selbasic;
119
120 MBlindPixelCalc blind;
121 blind.SetUseInterpolation();
122 blind.SetUseBlindPixels();
123 // blind.SetUseCetralPixel();
124
125 // create an object of class MSigmabar,
126 // because class MHSigmaTheta will look for it
127 MSigmabar sigbar;
128 pliston.AddToList(&sigbar);
129 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt");
130 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta");
131
132 MImgCleanStd clean;
133
134 // calculate also extended image parameters
135 TString fHilName = "Hillas";
136 MHillasExt hext(fHilName);
137 pliston.AddToList(&hext);
138 MHillasExt *fHillas = &hext;
139
140 // name of output container for MHillasCalc
141 // = name of MHillas object
142 MHillasCalc hcalc(fHilName);
143
144 // name of output container for MHillasSrcCalc
145 // = name of MHillasSrc object
146 TString fHilSrcName = "HillasSrc";
147 MHillasSrc hilsrc(fHilSrcName);
148 MHillasSrc *fHillasSrc = &hilsrc;
149 pliston.AddToList(&hilsrc);
150 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
151
152 MSelStandard selstand(fHillas, fHillasSrc);
153
154 MFillH hfill1("MHHillas", fHilName);
155 MFillH hfill2("MHStarMap", fHilName);
156
157 TString nxt = "HillasExt";
158 MHHillasExt hhilext(nxt, "", fHilName);
159 pliston.AddToList(&hhilext);
160 TString namext = nxt;
161 namext += "[MHHillasExt]";
162 MFillH hfill3(namext, fHilSrcName);
163
164 MFillH hfill4("MHHillasSrc", fHilSrcName);
165
166 //+++++++++++++++++++++++++++++++++++++++++++++++++++
167 // prepare writing of look-alike events (for hadrons)
168
169 const char* mtxName = "MatrixHadrons";
170 Int_t howMany = 30000;
171 char* outName = "matrix_hadrons.root";
172
173 cout << "" << endl;
174 cout << "========================================================" << endl;
175 cout << "Matrix for (hadrons)" << endl;
176 cout << "input file = " << filename<< endl;
177 cout << "matrix name = " << mtxName << endl;
178 cout << "max. no.of look-alike events = " << howMany << endl;
179 cout << "name of output root file = " << outName << endl;
180 cout << "" << endl;
181
182
183 // matrix limitation for look-alike events (approximate number)
184 MFEventSelector selector("", "", readname);
185 selector.SetNumSelectEvts(howMany);
186 MFilterList flist;
187 flist.AddToList(&selector);
188
189 //
190 // --- setup the matrix and add it to the parlist
191 //
192 MHMatrix *pmatrix = new MHMatrix(mtxName);
193 MHMatrix& matrix = *pmatrix;
194 matrix.AddColumn("Hillas.fWidth");
195 matrix.AddColumn("Hillas.fLength");
196 matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
197 matrix.AddColumn("abs(Hillas.fAsym)");
198 matrix.AddColumn("abs(Hillas.fM3Long)");
199 matrix.AddColumn("abs(Hillas.fM3Trans)");
200 matrix.AddColumn("abs(HillasSrc.fHeadTail)");
201 matrix.AddColumn("Hillas.fConc");
202 matrix.AddColumn("Hillas.fConc1");
203 matrix.AddColumn("HillasSrc.fDist");
204 matrix.AddColumn("log10(Hillas.fSize)");
205 matrix.AddColumn("HillasSrc.fAlpha");
206 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
207 //matrix.AddColumn("MMcEvt.fTheta");
208 pliston.AddToList(pmatrix);
209
210 MFillH fillmatg(mtxName);
211 fillmatg.SetFilter(&flist);
212
213 //+++++++++++++++++++++++++++++++++++++++++++++++++++
214
215 MSelFinal selfinal(fHillas, fHillasSrc);
216
217 //*****************************
218 // tasks to be executed
219
220 tliston.AddToList(&read);
221
222 tliston.AddToList(&selbasic);
223 tliston.AddToList(&blind);
224 tliston.AddToList(&fillsigtheta);
225 tliston.AddToList(&clean);
226
227 tliston.AddToList(&hcalc);
228 tliston.AddToList(&csrc1);
229
230 tliston.AddToList(&selstand);
231 tliston.AddToList(&hfill1);
232 tliston.AddToList(&hfill2);
233 tliston.AddToList(&hfill3);
234 tliston.AddToList(&hfill4);
235
236 tliston.AddToList(&flist);
237 tliston.AddToList(&fillmatg);
238
239 //tliston.AddToList(&selfinal);
240 //*****************************
241
242 //-------------------------------------------
243 // Execute event loop
244 //
245 MEvtLoop evtloop;
246 evtloop.SetParList(&pliston);
247
248 Int_t maxevents = 1000000000;
249 if ( !evtloop.Eventloop(maxevents) )
250 return;
251
252 tliston.PrintStatistics(0, kTRUE);
253
254 //-------------------------------------------
255 // Display the histograms
256 //
257 pliston.FindObject("MHHillas")->DrawClone();
258 pliston.FindObject("MHHillasSrc")->DrawClone();
259
260 //pliston.FindObject("MHHillasExt")->DrawClone();
261 pliston.FindObject(nxt)->DrawClone();
262
263 pliston.FindObject("MHStarMap")->DrawClone();
264
265
266 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
267 //
268 // ----------------------------------------------------------
269 // Definition of the reference sample of look-alike events
270 // (this is a subsample of the original sample)
271 // ----------------------------------------------------------
272 //
273 cout << "" << endl;
274 cout << "========================================================" << endl;
275 // Select a maximum of nmaxevts events from the sample of look-alike
276 // events. They will form the reference sample.
277 Int_t nmaxevts = 2000;
278
279 // target distribution for the variable in column refcolumn;
280 // set refcolumn negative if distribution is not to be changed
281 Int_t refcolumn = 12;
282 Int_t nbins = 10;
283 Float_t frombin = 0.5;
284 Float_t tobin = 1.0;
285 TH1F *thsh = new TH1F("thsh","target distribution",
286 nbins, frombin, tobin);
287 thsh->SetXTitle("cos( \\Theta)");
288 thsh->SetYTitle("Counts");
289 Float_t dbin = (tobin-frombin)/nbins;
290 Float_t lbin = frombin +dbin*0.5;
291 for (Int_t j=1; j<=nbins; j++)
292 {
293 thsh->Fill(lbin,1.0);
294 lbin += dbin;
295 }
296
297 cout << "" << endl;
298 cout << "========================================================" << endl;
299 cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl;
300 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
301 << refcolumn << ", " << nmaxevts << endl;
302
303 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
304 {
305 cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl;
306 return;
307 }
308 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
309
310 //-------------------------------------------
311 // write out look-alike events (for hadrons)
312 //
313 if (WLook)
314 {
315 TFile *writejens = new TFile(outName, "RECREATE", "");
316 matrix.Write();
317 cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file "
318 << outName << endl;
319
320 writejens->Close();
321 delete writejens;
322 }
323
324 //-------------------------------------------
325 // Write histograms onto a file
326 if (WHist)
327 {
328 MHSigmaTheta *sigtheta =
329 (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
330 if (!sigtheta)
331 {
332 *fLog << "Object 'SigmaTheta' not found" << endl;
333 return;
334 }
335 TH2D *fHSigTh = sigtheta->GetSigmaTheta();
336 TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta();
337 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta();
338
339 TFile outfile("SigmaTheta.root", "RECREATE");
340 fHSigTh->Write();
341 fHSigPixTh->Write();
342 fHDifPixTh->Write();
343 outfile.Close();
344
345 cout << "File 'SigmaTheta.root' was written out" << endl;
346 }
347
348
349 cout << "Macro AnalyseCT1 : End of section for ON data" << endl;
350 cout << "===================================================" << endl;
351 }
352 //---------------------------------------------------------------------
353
354 //---------------------------------------------------------------------
355 // start of section for MC gamma data
356
357 // read MC gamma data
358
359 // - to pad them
360 // (using the 2D-histogram "sigmabar versus Theta" of the OFF data)
361
362 // - to write a file of look-alike events (for gammas)
363 // (to be used later together with the corresponding hadron file
364 // for the g/h separation in the analysis of the ON data)
365
366 // - to write a file of padded MC gamma events
367 // (after the standard cuts, before the g/h separation)
368 // (to be used together with the corresponding hadron file
369 // for the optimization of the g/h separation)
370
371 // - to write a file of padded MC gamma events (after the final cuts)
372 // (to be used for the optimization of the energy estimator
373 // and for calculating the effective collection areas)
374
375 if (MC)
376 {
377 cout << "============================================================="
378 << endl;
379 cout << "Macro : AnalyseCT1 : Start of section for MC gamma data"
380 << endl;
381
382 MTaskList tlist;
383 MParList plist;
384
385 plist.AddToList(&tlist);
386
387 //------------------------------------
388 // Get the histograms "2D-ThetaSigmabar"
389 // and "3D-ThetaPixSigma"
390 // and "3D-ThetaPixDiff"
391 if (RHist)
392 {
393 cout << "Reading in file 'SigmaTheta.root' " << endl;
394
395 TFile *infile = new TFile("SigmaTheta.root");
396 infile->ls();
397
398 TH2D *fHSigmaTheta =
399 (TH2D*) gROOT.FindObject("2D-ThetaSigmabar");
400 if (!fHSigmaTheta)
401 {
402 *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl;
403 return;
404 }
405 cout << "Object '2D-ThetaSigmabar' was read in" << endl;
406
407 TH3D *fHSigmaPixTheta =
408 (TH3D*) gROOT.FindObject("3D-ThetaPixSigma");
409 if (!fHSigmaPixTheta)
410 {
411 *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl;
412 return;
413 }
414 cout << "Object '3D-ThetaPixSigma' was read in" << endl;
415
416 TH3D *fHDiffPixTheta =
417 (TH3D*) gROOT.FindObject("3D-ThetaPixDiff");
418 if (!fHSigmaPixTheta)
419 {
420 *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl;
421 return;
422 }
423 cout << "Object '3D-ThetaPixDiff' was read in" << endl;
424 }
425 else
426 {
427 MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta");
428 if (!sigtheta)
429 {
430 cout << "Object with name 'SigmaTheta' not found" << endl;
431 return;
432 }
433
434 TH2D *fHSigmaTheta = sigtheta->GetSigmaTheta();
435 TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta();
436 TH3D *fHDiffPixTheta = sigtheta->GetDiffPixTheta();
437 }
438 //------------------------------------
439
440
441
442 //::::::::::::::::::::::::::::::::::::::::::
443
444 MBinning binssize("BinningSize");
445 binssize.SetEdgesLog(50, 10, 1.0e5);
446 plist.AddToList(&binssize);
447
448 MBinning binsdistc("BinningDist");
449 binsdistc.SetEdges(50, 0, 1.4);
450 plist.AddToList(&binsdistc);
451
452 MBinning binswidth("BinningWidth");
453 binswidth.SetEdges(50, 0, 1.0);
454 plist.AddToList(&binswidth);
455
456 MBinning binslength("BinningLength");
457 binslength.SetEdges(50, 0, 1.0);
458 plist.AddToList(&binslength);
459
460 MBinning binsalpha("BinningAlpha");
461 binsalpha.SetEdges(100, -100, 100);
462 plist.AddToList(&binsalpha);
463
464 MBinning binsasym("BinningAsym");
465 binsasym.SetEdges(50, -1.5, 1.5);
466 plist.AddToList(&binsasym);
467
468 MBinning binsht("BinningHeadTail");
469 binsht.SetEdges(50, -1.5, 1.5);
470 plist.AddToList(&binsht);
471
472 MBinning binsm3l("BinningM3Long");
473 binsm3l.SetEdges(50, -1.5, 1.5);
474 plist.AddToList(&binsm3l);
475
476 MBinning binsm3t("BinningM3Trans");
477 binsm3t.SetEdges(50, -1.5, 1.5);
478 plist.AddToList(&binsm3t);
479
480
481 //.....
482 MBinning binsb("BinningSigmabar");
483 binsb.SetEdges( 100, 0.0, 5.0);
484 plist.AddToList(&binsb);
485
486 MBinning binth("BinningTheta");
487 binth.SetEdges( 70, -0.5, 69.5);
488 plist.AddToList(&binth);
489
490 MBinning binsdiff("BinningDiffsigma2");
491 binsdiff.SetEdges(100, -5.0, 20.0);
492 plist.AddToList(&binsdiff);
493
494 //::::::::::::::::::::::::::::::::::::::::::
495
496 //-------------------------------------------
497 // create the tasks which should be executed
498 //
499 filename = mcfile;
500 readname = "ReadCT1MC";
501 MCT1ReadPreProc read(filename, readname);
502
503 MSelBasic selbasic;
504
505 MBlindPixelCalc blind;
506 blind.SetUseInterpolation();
507 blind.SetUseBlindPixels();
508 // blind.SetUseCetralPixel();
509
510 // There are 2 options for Thomas Schweizer's padding
511 // fPadFlag = 1 get Sigmabar from fHSigmaTheta
512 // and Sigma from fHDiffPixTheta
513 // fPadFlag = 2 get Sigma from fHSigmaPixTheta
514
515 MPadSchweizer padthomas("MPadSchweizer","Task for the padding (Schweizer)",
516 fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta);
517 padthomas.SetPadFlag(1);
518
519 //...........................................
520
521 MImgCleanStd clean; //(0, 0);
522
523 // calculate also extended image parameters
524 TString fHilName = "Hillas";
525 MHillasExt hext(fHilName);
526 plist.AddToList(&hext);
527 MHillasExt *fHillas = &hext;
528
529 // name of output container for MHillasCalc
530 // = name of MHillas object
531 MHillasCalc hcalc(fHilName);
532
533 // name of output container for MHillasSrcCalc
534 // = name of MHillasSrc object
535 TString fHilSrcName = "HillasSrc";
536 MHillasSrc hilsrc(fHilSrcName);
537 MHillasSrc *fHillasSrc = &hilsrc;
538 plist.AddToList(&hilsrc);
539 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName);
540
541 MSelStandard selstand(fHillas, fHillasSrc);
542
543 MFillH hfill1("MHHillas", fHilName);
544 MFillH hfill2("MHStarMap", fHilName);
545
546 TString nxt = "HillasExt";
547 MHHillasExt hhilext(nxt, "", fHilName);
548 plist.AddToList(&hhilext);
549 TString namext = nxt;
550 namext += "[MHHillasExt]";
551 MFillH hfill3(namext, fHilSrcName);
552
553 MFillH hfill4("MHHillasSrc", fHilSrcName);
554
555 //+++++++++++++++++++++++++++++++++++++++++++++++++++
556 // prepare writing of look-alike events (for gammas)
557
558 const char* mtxName = "MatrixGammas";
559 Int_t howMany = 30000;
560 char* outName = "matrix_gammas.root";
561
562 cout << "" << endl;
563 cout << "========================================================" << endl;
564 cout << "Matrix for (gammas)" << endl;
565 cout << "input file = " << filename<< endl;
566 cout << "matrix name = " << mtxName << endl;
567 cout << "max. no.of look-alike events = " << howMany << endl;
568 cout << "name of output root file = " << outName << endl;
569 cout << "" << endl;
570
571
572 // matrix limitation for look-alike events (approximate number)
573 MFEventSelector selector("", "", readname);
574 selector.SetNumSelectEvts(howMany);
575 MFilterList flist;
576 flist.AddToList(&selector);
577
578 //
579 // --- setup the matrix and add it to the parlist
580 //
581 MHMatrix *pmatrix = new MHMatrix(mtxName);
582 MHMatrix& matrix = *pmatrix;
583 matrix.AddColumn("Hillas.fWidth");
584 matrix.AddColumn("Hillas.fLength");
585 matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize");
586 matrix.AddColumn("abs(Hillas.fAsym)");
587 matrix.AddColumn("abs(Hillas.fM3Long)");
588 matrix.AddColumn("abs(Hillas.fM3Trans)");
589 matrix.AddColumn("abs(HillasSrc.fHeadTail)");
590 matrix.AddColumn("Hillas.fConc");
591 matrix.AddColumn("Hillas.fConc1");
592 matrix.AddColumn("HillasSrc.fDist");
593 matrix.AddColumn("log10(Hillas.fSize)");
594 matrix.AddColumn("HillasSrc.fAlpha");
595 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)");
596 //matrix.AddColumn("MMcEvt.fTheta");
597 plist.AddToList(pmatrix);
598
599 MFillH fillmatg(mtxName);
600 fillmatg.SetFilter(&flist);
601
602 //+++++++++++++++++++++++++++++++++++++++++++++++++++
603
604
605 MSelFinal selfinal(fHillas, fHillasSrc);
606
607 //*****************************
608 // tasks to be executed
609
610 tlist.AddToList(&read);
611
612 tlist.AddToList(&selbasic);
613 tlist.AddToList(&blind);
614 tlist.AddToList(&padthomas);
615 tlist.AddToList(&clean);
616
617 tlist.AddToList(&hcalc);
618 tlist.AddToList(&csrc1);
619
620 tlist.AddToList(&selstand);
621 tlist.AddToList(&hfill1);
622 tlist.AddToList(&hfill2);
623 tlist.AddToList(&hfill3);
624 tlist.AddToList(&hfill4);
625
626 tlist.AddToList(&flist);
627 tlist.AddToList(&fillmatg);
628
629 //tlist.AddToList(&selfinal);
630 //*****************************
631
632
633 //-------------------------------------------
634 // Execute event loop
635 //
636 MEvtLoop evtloop;
637 evtloop.SetParList(&plist);
638
639 Int_t maxevents = 1000000000;
640 if ( !evtloop.Eventloop(maxevents) )
641 return;
642
643 tlist.PrintStatistics(0, kTRUE);
644
645 //-------------------------------------------
646 // Display the histograms
647 //
648 plist.FindObject("MHHillas")->DrawClone();
649 plist.FindObject("MHHillasSrc")->DrawClone();
650
651 //plist.FindObject("MHHillasExt")->DrawClone();
652 plist.FindObject(nxt)->DrawClone();
653
654 plist.FindObject("MHStarMap")->DrawClone();
655
656
657 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
658 //
659 // ----------------------------------------------------------
660 // Definition of the reference sample of look-alike events
661 // (this is a subsample of the original sample)
662 // ----------------------------------------------------------
663 //
664 cout << "" << endl;
665 cout << "========================================================" << endl;
666 // Select a maximum of nmaxevts events from the sample of look-alike
667 // events. They will form the reference sample.
668 Int_t nmaxevts = 2000;
669
670 // target distribution for the variable in column refcolumn;
671 // set refcolumn negative if distribution is not to be changed
672 Int_t refcolumn = 12;
673 Int_t nbins = 10;
674 Float_t frombin = 0.5;
675 Float_t tobin = 1.0;
676 TH1F *thsh = new TH1F("thsh","target distribution",
677 nbins, frombin, tobin);
678 Float_t dbin = (tobin-frombin)/nbins;
679 Float_t lbin = frombin +dbin*0.5;
680 for (Int_t j=1; j<=nbins; j++)
681 {
682 thsh->Fill(lbin,1.0);
683 lbin += dbin;
684 }
685
686 cout << "" << endl;
687 cout << "========================================================" << endl;
688 cout << "Macro AnalyseCT1 : define reference sample for gammas" << endl;
689 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = "
690 << refcolumn << ", " << nmaxevts << endl;
691
692 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) )
693 {
694 cout << "Macro AnalyseCT1 : Reference matrix for gammas cannot be defined" << endl;
695 return;
696 }
697 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
698
699 //-------------------------------------------
700 // write out look-alike events (for gammas)
701 //
702 if (WLookMC)
703 {
704 TFile *writejens = new TFile(outName, "RECREATE", "");
705 matrix.Write();
706 cout << "Macro AnalyseCT1 : matrix of look-alike events for gammas written onto file "
707 << outName << endl;
708
709 writejens->Close();
710 delete writejens;
711 }
712
713 cout << "Macro AnalyseCT1 : End of section for MC gamma data"
714 << endl;
715 cout << "========================================================="
716 << endl;
717 }
718}
719
720
721
722
Note: See TracBrowser for help on using the repository browser.