source: trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc@ 2010

Last change on this file since 2010 was 2001, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 23.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Robert Wagner <mailto:magicsoft@rwagner.de> 10/2002
19! Author(s): Wolfgang Wittek <mailto:wittek@mppmu.mpg.de> 02/2003
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MPadSchweizer
29//
30// This task applies padding such that for a given pixel and for a given
31// Theta bin the resulting distribution of the pedestal sigma is identical
32// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta.
33//
34// The number of photons, its error and the pedestal sigmas are altered.
35// On average, the number of photons added is zero.
36//
37// The formulas used can be found in Thomas Schweizer's Thesis,
38// Section 2.2.1
39//
40// There are 2 options for the padding :
41//
42// 1) fPadFlag = 1 :
43// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
44// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
45// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
46// (fHDiffPixTheta).
47//
48// This is the preferred option as it takes into account the
49// correlations between the Sigma of a pixel and Sigmabar.
50//
51// 2) fPadFlag = 2 :
52// Generate a pedestal sigma for each pixel using the 3D-histogram
53// Theta, pixel no., Sigma (fHSigmaPixTheta).
54//
55//
56// The padding has to be done before the image cleaning because the
57// image cleaning depends on the pedestal sigmas.
58//
59// For random numbers gRandom is used.
60//
61// This implementation has been tested for CT1 data. For MAGIC some
62// modifications are necessary.
63//
64/////////////////////////////////////////////////////////////////////////////
65#include "MPadSchweizer.h"
66
67#include <stdio.h>
68
69#include <TH1.h>
70#include <TH2.h>
71#include <TH3.h>
72#include <TRandom.h>
73#include <TCanvas.h>
74
75#include "MBinning.h"
76#include "MSigmabar.h"
77#include "MMcEvt.hxx"
78#include "MLog.h"
79#include "MLogManip.h"
80#include "MParList.h"
81#include "MGeomCam.h"
82
83#include "MCerPhotPix.h"
84#include "MCerPhotEvt.h"
85
86#include "MPedestalCam.h"
87#include "MPedestalPix.h"
88
89ClassImp(MPadSchweizer);
90
91// --------------------------------------------------------------------------
92//
93// Default constructor.
94//
95MPadSchweizer::MPadSchweizer(const char *name, const char *title)
96{
97 fName = name ? name : "MPadSchweizer";
98 fTitle = title ? title : "Task for the padding (Schweizer)";
99}
100
101// --------------------------------------------------------------------------
102//
103// Destructor. STLL DOESN'T DESTRUCT EVERYTHING
104//
105MPadSchweizer::~MPadSchweizer()
106{
107 *fLog << all << "WARNING: ~MPadSchweizer called: Potential Memory Leak" << endl;
108}
109
110// --------------------------------------------------------------------------
111//
112// Set the references to the histograms to be used in the padding
113// DOCUMENTATION of the arguments MISSING
114//
115void MPadSchweizer::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff)
116{
117 fHSigmaTheta = hist2;
118 fHSigmaPixTheta = hist3;
119 fHDiffPixTheta = hist3Diff;
120
121 fHSigmaTheta->SetDirectory(NULL);
122 fHSigmaPixTheta->SetDirectory(NULL);
123 fHDiffPixTheta->SetDirectory(NULL);
124
125 Print();
126}
127
128// --------------------------------------------------------------------------
129//
130// Set the option for the padding
131//
132// There are 2 options for the padding :
133//
134// 1) fPadFlag = 1 :
135// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
136// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
137// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
138// (fHDiffPixTheta).
139//
140// This is the preferred option as it takes into account the
141// correlations between the Sigma of a pixel and Sigmabar.
142//
143// 2) fPadFlag = 2 :
144// Generate a pedestal sigma for each pixel using the 3D-histogram
145// Theta, pixel no., Sigma (fHSigmaPixTheta).
146//
147void MPadSchweizer::SetPadFlag(Int_t padflag)
148{
149 fPadFlag = padflag;
150 *fLog << "MPadSchweizer::SetPadFlag(); choose option " << fPadFlag << endl;
151}
152
153// --------------------------------------------------------------------------
154//
155// Set the pointers and prepare the histograms
156//
157Bool_t MPadSchweizer::PreProcess(MParList *pList)
158{
159 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
160 if (!fMcEvt)
161 {
162 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
163 return kFALSE;
164 }
165
166 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
167 if (!fPed)
168 {
169 *fLog << err << "MPedestalCam not found... aborting." << endl;
170 return kFALSE;
171 }
172
173 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
174 if (!fCam)
175 {
176 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
177 return kFALSE;
178 }
179
180 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
181 if (!fEvt)
182 {
183 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
184 return kFALSE;
185 }
186
187 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
188 if (!fSigmabar)
189 {
190 *fLog << err << "MSigmabar not found... aborting." << endl;
191 return kFALSE;
192 }
193
194 // Get Theta Binning
195 MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta");
196 if (!binstheta)
197 {
198 *fLog << err << "BinningTheta not found... aborting." << endl;
199 return kFALSE;
200 }
201
202 // Get Sigma Binning
203 MBinning* binssigma = (MBinning*)pList->FindObject("BinningSigmabar");
204 if (!binssigma)
205 {
206 *fLog << err << "BinningSigmabar not found... aborting." << endl;
207 return kFALSE;
208 }
209
210 // Get binning for (sigma^2-sigmabar^2)
211 MBinning* binsdiff = (MBinning*)pList->FindObject("BinningDiffsigma2");
212 if (!binsdiff)
213 {
214 *fLog << err << "BinningDiffsigma2 not found... aborting." << endl;
215 return kFALSE;
216 }
217
218 // Get binning for pixel number
219 UInt_t npix = fPed->GetSize();
220 MBinning binspix("BinningPixel");
221 MBinning* binspixel = &binspix;
222 binspixel->SetEdges(npix, -0.5, ((float)npix)-0.5 );
223
224
225 //--------------------------------------------------------------------
226 // histograms for checking the padding
227 //
228 // plot pedestal sigmas
229 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
230 100, 0.0, 5.0, 100, 0.0, 5.0);
231 fHSigmaPedestal->SetXTitle("orig. Pedestal sigma");
232 fHSigmaPedestal->SetYTitle("padded Pedestal sigma");
233
234 // plot no.of photons (before vs. after padding)
235 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
236 100, -10.0, 90.0, 100, -10, 90);
237 fHPhotons->SetXTitle("no.of photons before padding");
238 fHPhotons->SetYTitle("no.of photons after padding");
239
240 // plot of added NSB
241 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
242 fHNSB->SetXTitle("no.of added NSB photons");
243 fHNSB->SetYTitle("no.of pixels");
244
245 fHSigbarTheta = new TH2D();
246 fHSigbarTheta->SetNameTitle("fHSigmabar",
247 "2D : Sigmabar, \\Theta (after padding)");
248 MH::SetBinning(fHSigbarTheta, binstheta, binssigma);
249 fHSigbarTheta->SetXTitle("Theta");
250 fHSigbarTheta->SetYTitle("Sigmabar");
251
252 fHSigPixTh = new TH3D();
253 fHSigPixTh->SetNameTitle("fHSigPixTh","3D: \\Theta, pixel no., Sigma (after padding)");
254 MH::SetBinning(fHSigPixTh, binstheta, binspixel, binssigma);
255 fHSigPixTh->SetXTitle("\\Theta [\\circ]");
256 fHSigPixTh->SetYTitle("pixel number");
257 fHSigPixTh->SetZTitle("Sigma");
258
259 fHDiffPixTh = new TH3D();
260 fHDiffPixTh->SetNameTitle("fHDiffPixTh","3D: \\Theta, pixel no., Sigma^2-Sigmabar^2 (after padding)");
261 MH::SetBinning(fHDiffPixTh, binstheta, binspixel, binsdiff);
262 fHDiffPixTh->SetXTitle("\\Theta [\\circ]");
263 fHDiffPixTh->SetYTitle("pixel number");
264 fHDiffPixTh->SetZTitle("Sigma^2-Sigmabar^2");
265
266 //--------------------------------------------------------------------
267
268 memset(fErrors, 0, sizeof(fErrors));
269
270 return kTRUE;
271}
272
273// --------------------------------------------------------------------------
274//
275// Do the Padding
276// idealy the events to be padded should have been generated without NSB
277//
278Bool_t MPadSchweizer::Process()
279{
280 //*fLog << "Entry MPadSchweizer::Process();" << endl;
281
282 Int_t rc;
283
284 const UInt_t npix = fEvt->GetNumPixels();
285
286 Double_t SigbarOld;
287
288
289 SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
290 //*fLog << "before padding : " << endl;
291 //fSigmabar->Print("");
292
293
294 //$$$$$$$$$$$$$$$$$$$$$$$$$$
295 // to simulate the situation that before the padding the NSB and
296 // electronic noise are zero : set Sigma = 0 for all pixels
297 //for (UInt_t i=0; i<npix; i++)
298 //{
299 // MCerPhotPix &pix = fEvt->operator[](i);
300 // Int_t j = pix.GetPixId();
301
302 // MPedestalPix &ppix = fPed->operator[](j);
303 // ppix.SetMeanRms(0.0);
304 //}
305 //$$$$$$$$$$$$$$$$$$$$$$$$$$
306
307 //-------------------------------------------
308 // Calculate average sigma of the event
309 //
310 SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);
311 //fSigmabar->Print("");
312
313 if (SigbarOld > 0.0)
314 {
315 //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : "
316 // << SigbarOld << ". Stop event loop " << endl;
317 // input data should have Sigmabar = 0; stop event loop
318
319 rc = 1;
320 fErrors[rc]++;
321 return kCONTINUE;
322 }
323
324 Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
325 // *fLog << "Theta = " << Theta << endl;
326
327
328 //-------------------------------------------
329 // for the current Theta,
330 // generate a Sigmabar according to the histogram fHSigmaTheta
331 //
332 Double_t Sigmabar;
333 Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(Theta);
334
335 if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() )
336 {
337 *fLog << "MPadSchweizer::Process(); binNumber out of range : Theta, binNumber = "
338 << Theta << ", " << binNumber << "; Skip event " << endl;
339 // event cannot be padded; skip event
340
341 rc = 2;
342 fErrors[rc]++;
343 return kCONTINUE;
344 }
345 else
346 {
347 TH1D* fHSigma =
348 fHSigmaTheta->ProjectionY("", binNumber, binNumber, "");
349 if ( fHSigma->GetEntries() == 0.0 )
350 {
351 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
352 << binNumber << " has no entries; Skip event " << endl;
353 // event cannot be padded; skip event
354 delete fHSigma;
355
356 rc = 3;
357 fErrors[rc]++;
358 return kCONTINUE;
359 }
360 else
361 {
362 Sigmabar = fHSigma->GetRandom();
363
364 //*fLog << "Theta, bin number = " << Theta << ", " << binNumber
365 // << ", Sigmabar = " << Sigmabar << endl;
366 }
367 delete fHSigma;
368 }
369
370
371 //-------------------------------------------
372
373 //*fLog << "MPadSchweizer::Process(); SigbarOld, Sigmabar = "
374 // << SigbarOld << ", "<< Sigmabar << endl;
375
376 // Skip event if target Sigmabar is <= SigbarOld
377 if (Sigmabar <= SigbarOld)
378 {
379 *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : "
380 << Sigmabar << ", " << SigbarOld << ", Skip event" << endl;
381
382 rc = 4;
383 fErrors[rc]++;
384 return kCONTINUE;
385 }
386
387
388 //-------------------------------------------
389 //
390 // Calculate average number of NSB photons to be added (lambdabar)
391 // from the value of Sigmabar,
392 // - making assumptions about the average electronic noise (elNoise2) and
393 // - using a fixed value (F2excess) for the excess noise factor
394
395 Double_t elNoise2; // [photons]
396 Double_t F2excess = 1.3;
397 Double_t lambdabar; // [photons]
398
399
400 Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(Theta);
401 if (binTheta != binNumber)
402 {
403 cout << "The binnings of the 2 histograms are not identical; aborting"
404 << endl;
405 return kERROR;
406 }
407
408 // Get RMS of (Sigma^2-Sigmabar^2) in this Theta bin.
409 // The average electronic noise (to be added) has to be well above this RMS,
410 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
411 // may become negative
412 TH1D* fHNoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
413 0, 9999, "");
414 Double_t RMS = fHNoise->GetRMS(1);
415 delete fHNoise;
416 elNoise2 = TMath::Min(RMS, Sigmabar*Sigmabar - SigbarOld*SigbarOld);
417 //*fLog << "elNoise2 = " << elNoise2 << endl;
418
419 lambdabar = (Sigmabar*Sigmabar - SigbarOld*SigbarOld - elNoise2)
420 / F2excess;
421 // This value of lambdabar is the same for all pixels;
422 // note that lambdabar is normalized to the area of pixel 0
423
424
425 //---------- start loop over pixels ---------------------------------
426 // do the padding for each pixel
427 //
428 // pad only pixels - which are used (before image cleaning)
429 //
430 Double_t Sig = 0.0;
431 Double_t Sigma2 = 0.0;
432 Double_t Diff = 0.0;
433 Double_t addSig2 = 0.0;
434 Double_t elNoise2Pix = 0.0;
435
436
437 for (UInt_t i=0; i<npix; i++)
438 {
439 MCerPhotPix &pix = fEvt->operator[](i);
440 if ( !pix.IsPixelUsed() )
441 continue;
442
443 //if ( pix.GetNumPhotons() == 0.0)
444 //{
445 // *fLog << "MPadSchweizer::Process(); no.of photons is 0 for used pixel"
446 // << endl;
447 // continue;
448 //}
449
450 Int_t j = pix.GetPixId();
451 Double_t Area = fCam->GetPixRatio(j);
452
453 MPedestalPix &ppix = fPed->operator[](j);
454 Double_t oldsigma = ppix.GetMeanRms();
455
456 //---------------------------------
457 // throw the Sigma for this pixel
458 //
459 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
460
461 Int_t count;
462 Int_t OK;
463 TH1D* fHDiff;
464 TH1D* fHSig;
465
466 switch (fPadFlag)
467 {
468 case 1 :
469 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
470 fHDiff =
471 fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
472 binPixel, binPixel, "");
473 if ( fHDiff->GetEntries() == 0.0 )
474 {
475 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
476 << binTheta << " and pixel bin " << binPixel
477 << " has no entries; aborting " << endl;
478
479 rc = 5;
480 fErrors[rc]++;
481 return kCONTINUE;
482 }
483
484 count = 0;
485 OK = 0;
486 for (Int_t m=0; m<20; m++)
487 {
488 count += 1;
489 Diff = fHDiff->GetRandom();
490 // the following condition ensures that elNoise2Pix > 0.0
491 if ( (Diff+Sigmabar*Sigmabar-oldsigma*oldsigma/Area-
492 lambdabar*F2excess) > 0.0 )
493 {
494 OK = 1;
495 break;
496 }
497 }
498 if (OK == 0)
499 {
500 *fLog << "Theta, j, count, Sigmabar, Diff = " << Theta << ", "
501 << j << ", " << count << ", " << Sigmabar << ", "
502 << Diff << endl;
503 Diff = lambdabar*F2excess + oldsigma*oldsigma/Area
504 - Sigmabar*Sigmabar;
505 }
506 delete fHDiff;
507 Sigma2 = Diff + Sigmabar*Sigmabar;
508 break;
509
510 case 2 :
511 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
512 fHSig =
513 fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
514 binPixel, binPixel, "");
515 if ( fHSig->GetEntries() == 0.0 )
516 {
517 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
518 << binTheta << " and pixel bin " << binPixel
519 << " has no entries; aborting " << endl;
520
521 rc = 6;
522 fErrors[rc]++;
523 return kCONTINUE;
524 }
525
526 count = 0;
527 OK = 0;
528 for (Int_t m=0; m<20; m++)
529 {
530 count += 1;
531 Sig = fHSig->GetRandom();
532 Sigma2 = Sig*Sig/Area;
533 // the following condition ensures that elNoise2Pix > 0.0
534 if ( (Sigma2-oldsigma*oldsigma/Area-lambdabar*F2excess) > 0.0 )
535 {
536 OK = 1;
537 break;
538 }
539 }
540 if (OK == 0)
541 {
542 *fLog << "Theta, j, count, Sigmabar, Sig = " << Theta << ", "
543 << j << ", " << count << ", " << Sigmabar << ", "
544 << Sig << endl;
545 Sigma2 = lambdabar*F2excess + oldsigma*oldsigma/Area;
546 }
547 delete fHSig;
548 break;
549 }
550
551 //---------------------------------
552 // get the additional sigma^2 for this pixel (due to the padding)
553 addSig2 = Sigma2*Area - oldsigma*oldsigma;
554
555
556 //---------------------------------
557 // get the additional electronic noise for this pixel
558 elNoise2Pix = addSig2 - lambdabar*F2excess*Area;
559
560
561 //---------------------------------
562 // throw actual number of additional NSB photons (NSB)
563 // and its RMS (sigmaNSB)
564 Double_t NSB0 = gRandom->Poisson(lambdabar*Area);
565 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
566 Double_t sigmaNSB0;
567
568 if (arg >= 0.0)
569 {
570 sigmaNSB0 = sqrt( arg );
571 }
572 else
573 {
574 *fLog << "MPadSchweizer::Process(); argument of sqrt < 0.0 : "
575 << arg << endl;
576 sigmaNSB0 = 0.0000001;
577 }
578
579
580 //---------------------------------
581 // smear NSB0 according to sigmaNSB0
582 // and subtract lambdabar because of AC coupling
583 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*Area;
584
585 //---------------------------------
586
587 // add additional NSB to the number of photons
588 Double_t oldphotons = pix.GetNumPhotons();
589 Double_t newphotons = oldphotons + NSB;
590 pix.SetNumPhotons( newphotons );
591
592 fHNSB->Fill( NSB/sqrt(Area) );
593 fHPhotons->Fill( oldphotons/sqrt(Area), newphotons/sqrt(Area) );
594
595
596 // error: add sigma of padded noise quadratically
597 Double_t olderror = pix.GetErrorPhot();
598 Double_t newerror = sqrt( olderror*olderror + addSig2 );
599 pix.SetErrorPhot( newerror );
600
601
602 Double_t newsigma = sqrt( oldsigma*oldsigma + addSig2 );
603 ppix.SetMeanRms( newsigma );
604
605 fHSigmaPedestal->Fill( oldsigma, newsigma );
606 fHSigPixTh-> Fill( Theta, (Double_t) j, newsigma );
607 }
608 //---------- end of loop over pixels ---------------------------------
609
610 // Calculate Sigmabar again and crosscheck
611
612 Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt);
613 //*fLog << "after padding : " << endl;
614 //fSigmabar->Print("");
615
616 fHSigbarTheta->Fill( Theta, SigbarNew );
617
618 // this loop is only for filling the histogram fHDiffPixTh
619 for (UInt_t i=0; i<npix; i++)
620 {
621 MCerPhotPix &pix = fEvt->operator[](i);
622 if ( !pix.IsPixelUsed() )
623 continue;
624
625 //if ( pix.GetNumPhotons() == 0.0)
626 //{
627 // *fLog << "MPadSchweizer::Process(); no.of photons is 0 for used pixel"
628 // << endl;
629 // continue;
630 //}
631
632 Int_t j = pix.GetPixId();
633 Double_t Area = fCam->GetPixRatio(j);
634
635 MPedestalPix &ppix = fPed->operator[](j);
636 Double_t newsigma = ppix.GetMeanRms();
637
638 fHDiffPixTh->Fill( Theta, (Double_t) j,
639 newsigma*newsigma/Area-SigbarNew*SigbarNew);
640 }
641
642
643 //*fLog << "Exit MPadSchweizer::Process();" << endl;
644
645 rc = 0;
646 fErrors[rc]++;
647
648 return kTRUE;
649
650}
651
652// --------------------------------------------------------------------------
653//
654//
655Bool_t MPadSchweizer::PostProcess()
656{
657 if (GetNumExecutions() != 0)
658 {
659
660 *fLog << inf << endl;
661 *fLog << GetDescriptor() << " execution statistics:" << endl;
662 *fLog << dec << setfill(' ');
663 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
664 << (int)(fErrors[1]*100/GetNumExecutions())
665 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
666
667 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
668 << (int)(fErrors[2]*100/GetNumExecutions())
669 << "%) Evts skipped due to: Zenith angle out of range" << endl;
670
671 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
672 << (int)(fErrors[3]*100/GetNumExecutions())
673 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
674
675 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
676 << (int)(fErrors[4]*100/GetNumExecutions())
677 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
678
679 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
680 << (int)(fErrors[5]*100/GetNumExecutions())
681 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
682
683 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
684 << (int)(fErrors[6]*100/GetNumExecutions())
685 << "%) Evts skipped due to: No data for generating Sigma" << endl;
686
687 *fLog << " " << fErrors[0] << " ("
688 << (int)(fErrors[0]*100/GetNumExecutions())
689 << "%) Evts survived the padding!" << endl;
690 *fLog << endl;
691
692 }
693
694 //---------------------------------------------------------------
695 TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200));
696 c.Divide(3, 4);
697 gROOT->SetSelectedPad(NULL);
698
699
700 c.cd(1);
701 fHSigmaTheta->SetDirectory(NULL);
702 fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)");
703 fHSigmaTheta->DrawCopy();
704 fHSigmaTheta->SetBit(kCanDelete);
705
706 //c.cd(1);
707 //fHSigmaPixTheta->DrawCopy();
708 //fHSigmaPixTheta->SetBit(kCanDelete);
709
710 c.cd(4);
711 fHSigbarTheta->SetDirectory(NULL);
712 fHSigbarTheta->DrawCopy();
713 fHSigbarTheta->SetBit(kCanDelete);
714
715
716 c.cd(7);
717 fHNSB->SetDirectory(NULL);
718 fHNSB->DrawCopy();
719 fHNSB->SetBit(kCanDelete);
720
721
722 //--------------------------------------------------------------------
723 // draw the 3D histogram : Theta, pixel, Sigma^2-Sigmabar^2
724
725 TH2D *l;
726
727 c.cd(2);
728 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx");
729 l->SetDirectory(NULL);
730 l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
731 l->SetXTitle("\\Theta [\\circ]");
732 l->SetYTitle("Sigma^2-Sigmabar^2");
733
734 l->DrawCopy("box");
735 l->SetBit(kCanDelete);;
736
737 c.cd(5);
738 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy");
739 l->SetDirectory(NULL);
740 l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
741 l->SetXTitle("pixel");
742 l->SetYTitle("Sigma^2-Sigmabar^2");
743
744 l->DrawCopy("box");
745 l->SetBit(kCanDelete);;
746
747 c.cd(8);
748 ((TH2*)fHDiffPixTh)->DrawCopy();
749
750
751 //--------------------------------------------------------------------
752 // draw the 3D histogram : Theta, pixel, Sigma
753
754 TH2D *k;
755
756 c.cd(3);
757 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx");
758 k->SetDirectory(NULL);
759 k->SetTitle("Sigma vs. \\Theta (all pixels)");
760 k->SetXTitle("\\Theta [\\circ]");
761 k->SetYTitle("Sigma");
762
763 k->DrawCopy("box");
764 k->SetBit(kCanDelete);
765
766 c.cd(6);
767 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy");
768 k->SetDirectory(NULL);
769 k->SetTitle("Sigma vs. pixel number (all \\Theta)");
770 k->SetXTitle("pixel");
771 k->SetYTitle("Sigma");
772
773 k->DrawCopy("box");
774 k->SetBit(kCanDelete);;
775
776 c.cd(9);
777 ((TH2*)fHSigPixTh)->DrawCopy();
778
779 //--------------------------------------------------------------------
780
781 c.cd(10);
782 fHSigmaPedestal->SetDirectory(NULL);
783 fHSigmaPedestal->DrawCopy();
784 fHSigmaPedestal->SetBit(kCanDelete);
785
786 c.cd(11);
787 fHPhotons->SetDirectory(NULL);
788 fHPhotons->DrawCopy();
789 fHPhotons->SetBit(kCanDelete);
790
791 //--------------------------------------------------------------------
792
793 return kTRUE;
794}
795
796// --------------------------------------------------------------------------
797
Note: See TracBrowser for help on using the repository browser.