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

Last change on this file since 2015 was 2012, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 23.4 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=0;
283
284 const UInt_t npix = fEvt->GetNumPixels();
285
286 fSigmabar->Calc(*fCam, *fPed, *fEvt);
287 //*fLog << "before padding : " << endl;
288 //fSigmabar->Print("");
289
290
291 //$$$$$$$$$$$$$$$$$$$$$$$$$$
292 // to simulate the situation that before the padding the NSB and
293 // electronic noise are zero : set Sigma = 0 for all pixels
294 //for (UInt_t i=0; i<npix; i++)
295 //{
296 // MCerPhotPix &pix = fEvt->operator[](i);
297 // Int_t j = pix.GetPixId();
298
299 // MPedestalPix &ppix = fPed->operator[](j);
300 // ppix.SetMeanRms(0.0);
301 //}
302 //$$$$$$$$$$$$$$$$$$$$$$$$$$
303
304 //-------------------------------------------
305 // Calculate average sigma of the event
306 //
307 Double_t sigbarold = fSigmabar->Calc(*fCam, *fPed, *fEvt);
308 //fSigmabar->Print("");
309
310 if (sigbarold > 0)
311 {
312 //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : "
313 // << SigbarOld << ". Stop event loop " << endl;
314 // input data should have Sigmabar = 0; stop event loop
315
316 rc = 1;
317 fErrors[rc]++;
318 return kCONTINUE;
319 }
320
321 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
322 // *fLog << "Theta = " << Theta << endl;
323
324
325 //-------------------------------------------
326 // for the current Theta,
327 // generate a Sigmabar according to the histogram fHSigmaTheta
328 //
329 Double_t sigmabar=0;
330 Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta);
331
332 if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() )
333 {
334 *fLog << "MPadSchweizer::Process(); binNumber out of range : Theta, binNumber = "
335 << theta << ", " << binNumber << "; Skip event " << endl;
336 // event cannot be padded; skip event
337
338 rc = 2;
339 fErrors[rc]++;
340 return kCONTINUE;
341 }
342 else
343 {
344 TH1D* hsigma =
345 fHSigmaTheta->ProjectionY("", binNumber, binNumber, "");
346 if ( hsigma->GetEntries() == 0 )
347 {
348 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
349 << binNumber << " has no entries; Skip event " << endl;
350 // event cannot be padded; skip event
351 delete hsigma;
352
353 rc = 3;
354 fErrors[rc]++;
355 return kCONTINUE;
356 }
357 else
358 {
359 sigmabar = hsigma->GetRandom();
360
361 //*fLog << "Theta, bin number = " << Theta << ", " << binNumber
362 // << ", Sigmabar = " << Sigmabar << endl;
363 }
364 delete hsigma;
365 }
366
367
368 //-------------------------------------------
369
370 //*fLog << "MPadSchweizer::Process(); SigbarOld, Sigmabar = "
371 // << SigbarOld << ", "<< Sigmabar << endl;
372
373 // Skip event if target Sigmabar is <= SigbarOld
374 if (sigmabar <= sigbarold)
375 {
376 *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : "
377 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
378
379 rc = 4;
380 fErrors[rc]++;
381 return kCONTINUE;
382 }
383
384
385 //-------------------------------------------
386 //
387 // Calculate average number of NSB photons to be added (lambdabar)
388 // from the value of Sigmabar,
389 // - making assumptions about the average electronic noise (elNoise2) and
390 // - using a fixed value (F2excess) for the excess noise factor
391
392 Double_t F2excess = 1.3;
393
394 const Double_t sigmabar2 = sigmabar*sigmabar;
395
396 Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta);
397 if (binTheta != binNumber)
398 {
399 cout << "The binnings of the 2 histograms are not identical; aborting"
400 << endl;
401 return kERROR;
402 }
403
404 // Get RMS of (Sigma^2-Sigmabar^2) in this Theta bin.
405 // The average electronic noise (to be added) has to be well above this RMS,
406 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
407 // may become negative
408 TH1D* hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
409 0, 9999, "");
410
411 const Double_t elNoise2 = TMath::Min(hnoise->GetRMS(1),
412 sigmabar2 - sigbarold*sigbarold); // [photons]
413
414 delete hnoise;
415
416 //*fLog << "elNoise2 = " << elNoise2 << endl;
417
418 const Double_t lambdabar = (sigmabar2 - sigbarold*sigbarold - elNoise2)
419 / F2excess; // [photons]
420
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 //---------- start loop over pixels ---------------------------------
425 // do the padding for each pixel
426 //
427 // pad only pixels - which are used (before image cleaning)
428 //
429 Double_t sig = 0;
430 Double_t sigma2 = 0;
431 Double_t diff = 0;
432 Double_t addSig2 = 0;
433 Double_t elNoise2Pix = 0;
434
435
436 for (UInt_t i=0; i<npix; i++)
437 {
438 MCerPhotPix &pix = (*fEvt)[i];
439 if ( !pix.IsPixelUsed() )
440 continue;
441
442 //if ( pix.GetNumPhotons() == 0.0)
443 //{
444 // *fLog << "MPadSchweizer::Process(); no.of photons is 0 for used pixel"
445 // << endl;
446 // continue;
447 //}
448
449 Int_t j = pix.GetPixId();
450 Double_t area = fCam->GetPixRatio(j);
451
452 MPedestalPix &ppix = (*fPed)[j];
453 Double_t oldsigma = ppix.GetMeanRms();
454
455 //---------------------------------
456 // throw the Sigma for this pixel
457 //
458 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
459
460 Int_t count;
461 Bool_t ok=kFALSE;
462 TH1D* hDiff;
463 TH1D* hSig;
464
465 switch (fPadFlag)
466 {
467 case 1 :
468 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
469 hDiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
470 binPixel, binPixel, "");
471 if ( hDiff->GetEntries() == 0 )
472 {
473 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
474 << binTheta << " and pixel bin " << binPixel
475 << " has no entries; aborting " << endl;
476
477 rc = 5;
478 fErrors[rc]++;
479 return kCONTINUE;
480 }
481
482 count = 0;
483 for (Int_t m=0; m<20; m++)
484 {
485 count += 1;
486 diff = hDiff->GetRandom();
487 // the following condition ensures that elNoise2Pix > 0.0
488 if ( (diff+sigmabar2-oldsigma*oldsigma/area-
489 lambdabar*F2excess) > 0 )
490 {
491 ok = kTRUE;
492 break;
493 }
494 }
495 if (!ok)
496 {
497 *fLog << "Theta, j, count, Sigmabar, Diff = " << theta << ", "
498 << j << ", " << count << ", " << sigmabar << ", "
499 << diff << endl;
500 diff = lambdabar*F2excess + oldsigma*oldsigma/area
501 - sigmabar2;
502 }
503 delete hDiff;
504 sigma2 = diff + sigmabar2;
505 break;
506
507 case 2 :
508 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
509 hSig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
510 binPixel, binPixel, "");
511 if ( hSig->GetEntries() == 0 )
512 {
513 *fLog << "MPadSchweizer::Process(); projection for Theta bin "
514 << binTheta << " and pixel bin " << binPixel
515 << " has no entries; aborting " << endl;
516
517 rc = 6;
518 fErrors[rc]++;
519 return kCONTINUE;
520 }
521
522 count = 0;
523 for (Int_t m=0; m<20; m++)
524 {
525 count += 1;
526 sig = hSig->GetRandom();
527 sigma2 = sig*sig/area;
528 // the following condition ensures that elNoise2Pix > 0.0
529 if ( (sigma2-oldsigma*oldsigma/area-lambdabar*F2excess) > 0.0 )
530 {
531 ok = kTRUE;
532 break;
533 }
534 }
535 if (!ok)
536 {
537 *fLog << "Theta, j, count, Sigmabar, Sig = " << theta << ", "
538 << j << ", " << count << ", " << sigmabar << ", "
539 << sig << endl;
540 sigma2 = lambdabar*F2excess + oldsigma*oldsigma/area;
541 }
542 delete hSig;
543 break;
544 }
545
546 //---------------------------------
547 // get the additional sigma^2 for this pixel (due to the padding)
548 addSig2 = sigma2*area - oldsigma*oldsigma;
549
550
551 //---------------------------------
552 // get the additional electronic noise for this pixel
553 elNoise2Pix = addSig2 - lambdabar*F2excess*area;
554
555
556 //---------------------------------
557 // throw actual number of additional NSB photons (NSB)
558 // and its RMS (sigmaNSB)
559 Double_t NSB0 = gRandom->Poisson(lambdabar*area);
560 Double_t arg = NSB0*(F2excess-1) + elNoise2Pix;
561 Double_t sigmaNSB0;
562
563 if (arg >= 0)
564 {
565 sigmaNSB0 = sqrt( arg );
566 }
567 else
568 {
569 *fLog << "MPadSchweizer::Process(); argument of sqrt < 0 : "
570 << arg << endl;
571 sigmaNSB0 = 0.0000001;
572 }
573
574
575 //---------------------------------
576 // smear NSB0 according to sigmaNSB0
577 // and subtract lambdabar because of AC coupling
578 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*area;
579
580 //---------------------------------
581
582 // add additional NSB to the number of photons
583 Double_t oldphotons = pix.GetNumPhotons();
584 Double_t newphotons = oldphotons + NSB;
585 pix.SetNumPhotons( newphotons );
586
587 const Double_t areasqrt = sqrt(area);
588
589 fHNSB->Fill( NSB/areasqrt );
590 fHPhotons->Fill( oldphotons/areasqrt, newphotons/areasqrt );
591
592
593 // error: add sigma of padded noise quadratically
594 Double_t olderror = pix.GetErrorPhot();
595 Double_t newerror = sqrt( olderror*olderror + addSig2 );
596 pix.SetErrorPhot( newerror );
597
598
599 Double_t newsigma = sqrt( oldsigma*oldsigma + addSig2 );
600 ppix.SetMeanRms( newsigma );
601
602 fHSigmaPedestal->Fill( oldsigma, newsigma );
603 fHSigPixTh-> Fill( theta, (Double_t) j, newsigma );
604 }
605 //---------- end of loop over pixels ---------------------------------
606
607 // Calculate Sigmabar again and crosscheck
608
609 Double_t sigbarnew = fSigmabar->Calc(*fCam, *fPed, *fEvt);
610 //*fLog << "after padding : " << endl;
611 //fSigmabar->Print("");
612
613 fHSigbarTheta->Fill( theta, sigbarnew );
614
615 // this loop is only for filling the histogram fHDiffPixTh
616 for (UInt_t i=0; i<npix; i++)
617 {
618 MCerPhotPix &pix = (*fEvt)[i];
619 if ( !pix.IsPixelUsed() )
620 continue;
621
622 //if ( pix.GetNumPhotons() == 0.0)
623 //{
624 // *fLog << "MPadSchweizer::Process(); no.of photons is 0 for used pixel"
625 // << endl;
626 // continue;
627 //}
628
629 Int_t j = pix.GetPixId();
630 Double_t area = fCam->GetPixRatio(j);
631
632 MPedestalPix &ppix = (*fPed)[j];
633 Double_t newsigma = ppix.GetMeanRms();
634
635 fHDiffPixTh->Fill( theta, (Double_t) j,
636 newsigma*newsigma/area-sigbarnew*sigbarnew);
637 }
638
639
640 //*fLog << "Exit MPadSchweizer::Process();" << endl;
641
642 rc = 0;
643 fErrors[rc]++;
644
645 return kTRUE;
646
647}
648
649// --------------------------------------------------------------------------
650//
651//
652Bool_t MPadSchweizer::PostProcess()
653{
654 if (GetNumExecutions() != 0)
655 {
656
657 *fLog << inf << endl;
658 *fLog << GetDescriptor() << " execution statistics:" << endl;
659 *fLog << dec << setfill(' ');
660 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
661 << (int)(fErrors[1]*100/GetNumExecutions())
662 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
663
664 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
665 << (int)(fErrors[2]*100/GetNumExecutions())
666 << "%) Evts skipped due to: Zenith angle out of range" << endl;
667
668 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
669 << (int)(fErrors[3]*100/GetNumExecutions())
670 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
671
672 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
673 << (int)(fErrors[4]*100/GetNumExecutions())
674 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
675
676 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
677 << (int)(fErrors[5]*100/GetNumExecutions())
678 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
679
680 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
681 << (int)(fErrors[6]*100/GetNumExecutions())
682 << "%) Evts skipped due to: No data for generating Sigma" << endl;
683
684 *fLog << " " << fErrors[0] << " ("
685 << (int)(fErrors[0]*100/GetNumExecutions())
686 << "%) Evts survived the padding!" << endl;
687 *fLog << endl;
688
689 }
690
691 //---------------------------------------------------------------
692 TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200));
693 c.Divide(3, 4);
694 gROOT->SetSelectedPad(NULL);
695
696
697 c.cd(1);
698 fHSigmaTheta->SetDirectory(NULL);
699 fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)");
700 fHSigmaTheta->DrawCopy();
701 fHSigmaTheta->SetBit(kCanDelete);
702
703 //c.cd(1);
704 //fHSigmaPixTheta->DrawCopy();
705 //fHSigmaPixTheta->SetBit(kCanDelete);
706
707 c.cd(4);
708 fHSigbarTheta->SetDirectory(NULL);
709 fHSigbarTheta->DrawCopy();
710 fHSigbarTheta->SetBit(kCanDelete);
711
712
713 c.cd(7);
714 fHNSB->SetDirectory(NULL);
715 fHNSB->DrawCopy();
716 fHNSB->SetBit(kCanDelete);
717
718
719 //--------------------------------------------------------------------
720 // draw the 3D histogram : Theta, pixel, Sigma^2-Sigmabar^2
721
722 TH2D *l;
723
724 c.cd(2);
725 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx");
726 l->SetDirectory(NULL);
727 l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
728 l->SetXTitle("\\Theta [\\circ]");
729 l->SetYTitle("Sigma^2-Sigmabar^2");
730
731 l->DrawCopy("box");
732 l->SetBit(kCanDelete);;
733
734 c.cd(5);
735 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy");
736 l->SetDirectory(NULL);
737 l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
738 l->SetXTitle("pixel");
739 l->SetYTitle("Sigma^2-Sigmabar^2");
740
741 l->DrawCopy("box");
742 l->SetBit(kCanDelete);;
743
744 c.cd(8);
745 ((TH2*)fHDiffPixTh)->DrawCopy();
746
747
748 //--------------------------------------------------------------------
749 // draw the 3D histogram : Theta, pixel, Sigma
750
751 TH2D *k;
752
753 c.cd(3);
754 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx");
755 k->SetDirectory(NULL);
756 k->SetTitle("Sigma vs. \\Theta (all pixels)");
757 k->SetXTitle("\\Theta [\\circ]");
758 k->SetYTitle("Sigma");
759
760 k->DrawCopy("box");
761 k->SetBit(kCanDelete);
762
763 c.cd(6);
764 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy");
765 k->SetDirectory(NULL);
766 k->SetTitle("Sigma vs. pixel number (all \\Theta)");
767 k->SetXTitle("pixel");
768 k->SetYTitle("Sigma");
769
770 k->DrawCopy("box");
771 k->SetBit(kCanDelete);;
772
773 c.cd(9);
774 ((TH2*)fHSigPixTh)->DrawCopy();
775
776 //--------------------------------------------------------------------
777
778 c.cd(10);
779 fHSigmaPedestal->SetDirectory(NULL);
780 fHSigmaPedestal->DrawCopy();
781 fHSigmaPedestal->SetBit(kCanDelete);
782
783 c.cd(11);
784 fHPhotons->SetDirectory(NULL);
785 fHPhotons->DrawCopy();
786 fHPhotons->SetBit(kCanDelete);
787
788 //--------------------------------------------------------------------
789
790 return kTRUE;
791}
792
793// --------------------------------------------------------------------------
794
Note: See TracBrowser for help on using the repository browser.