source: trunk/MagicSoft/Mars/manalysis/MPadON.cc@ 2773

Last change on this file since 2773 was 2746, checked in by wittek, 22 years ago
*** empty log message ***
File size: 26.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): Wolfgang Wittek, 12/2003 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MPadON
28//
29// This task applies padding such that for a given pixel and for a given
30// Theta bin the resulting distribution of the pedestal sigma is identical
31// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta.
32//
33// The number of photons, its error and the pedestal sigmas are altered.
34// On average, the number of photons added is zero.
35//
36// The formulas used can be found in Thomas Schweizer's Thesis,
37// Section 2.2.1
38//
39// There are 2 options for the padding :
40//
41// 1) fPadFlag = 1 :
42// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
43// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
44// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
45// (fHDiffPixTheta).
46//
47// This is the preferred option as it takes into account the
48// correlations between the Sigma of a pixel and Sigmabar.
49//
50// THE FOLLOWING OPTION HAS NOT BEEN TESTED :
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/////////////////////////////////////////////////////////////////////////////
62#include "MPadON.h"
63
64#include <stdio.h>
65
66#include <TH1.h>
67#include <TH2.h>
68#include <TH3.h>
69#include <TRandom.h>
70#include <TCanvas.h>
71
72#include "MBinning.h"
73#include "MSigmabar.h"
74#include "MMcEvt.hxx"
75#include "MLog.h"
76#include "MLogManip.h"
77#include "MParList.h"
78#include "MGeomCam.h"
79
80#include "MCerPhotPix.h"
81#include "MCerPhotEvt.h"
82
83#include "MPedestalCam.h"
84#include "MPedestalPix.h"
85#include "MBlindPixels.h"
86
87ClassImp(MPadON);
88
89using namespace std;
90
91// --------------------------------------------------------------------------
92//
93// Default constructor.
94//
95MPadON::MPadON(const char *name, const char *title)
96{
97 fName = name ? name : "MPadON";
98 fTitle = title ? title : "Task for the padding (Schweizer)";
99
100 fHSigmaTheta = NULL;
101 fHSigmaPixTheta = NULL;
102 fHDiffPixTheta = NULL;
103 fHBlindPixIdTheta = NULL;
104 fHBlindPixNTheta = NULL;
105
106 fHSigmaPedestal = NULL;
107 fHPhotons = NULL;
108 fHNSB = NULL;
109}
110
111// --------------------------------------------------------------------------
112//
113// Destructor.
114//
115MPadON::~MPadON()
116{
117 if (fHSigmaPedestal != NULL) delete fHSigmaPedestal;
118 if (fHPhotons != NULL) delete fHPhotons;
119 if (fHNSB != NULL) delete fHNSB;
120}
121
122// --------------------------------------------------------------------------
123//
124// Set the references to the histograms to be used in the padding
125//
126// fHSigmaTheta 2D-histogram (Theta, sigmabar)
127// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma)
128// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2)
129// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id)
130// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels )
131//
132//
133void MPadON::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff,
134 TH2D *hist2Pix, TH2D *hist2PixN)
135{
136 fHSigmaTheta = hist2;
137 fHSigmaPixTheta = hist3;
138 fHDiffPixTheta = hist3Diff;
139 fHBlindPixIdTheta = hist2Pix;
140 fHBlindPixNTheta = hist2PixN;
141
142 fHSigmaTheta->SetDirectory(NULL);
143 fHSigmaPixTheta->SetDirectory(NULL);
144 fHDiffPixTheta->SetDirectory(NULL);
145 fHBlindPixIdTheta->SetDirectory(NULL);
146 fHBlindPixNTheta->SetDirectory(NULL);
147
148 Print();
149}
150
151// --------------------------------------------------------------------------
152//
153// Set the option for the padding
154//
155// There are 2 options for the padding :
156//
157// 1) fPadFlag = 1 :
158// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
159// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
160// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
161// (fHDiffPixTheta).
162//
163// This is the preferred option as it takes into account the
164// correlations between the Sigma of a pixel and Sigmabar.
165//
166// NOT YET TESTED :
167// 2) fPadFlag = 2 :
168// Generate a pedestal sigma for each pixel using the 3D-histogram
169// Theta, pixel no., Sigma (fHSigmaPixTheta).
170//
171void MPadON::SetPadFlag(Int_t padflag)
172{
173 fPadFlag = padflag;
174 *fLog << inf << "MPadON::SetPadFlag(); choose option " << fPadFlag << endl;
175}
176
177// --------------------------------------------------------------------------
178//
179// Set the pointers and prepare the histograms
180//
181Int_t MPadON::PreProcess(MParList *pList)
182{
183 if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta ||
184 !fHBlindPixIdTheta || !fHBlindPixNTheta)
185 {
186 *fLog << err
187 << "MPadON : At least one of the histograms needed for the padding is not defined ... aborting."
188 << endl;
189 return kFALSE;
190 }
191
192 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
193 if (!fMcEvt)
194 {
195 *fLog << err << dbginf << "MPadON : MMcEvt not found... aborting."
196 << endl;
197 return kFALSE;
198 }
199
200 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
201 if (!fPed)
202 {
203 *fLog << err << "MPadON : MPedestalCam not found... aborting." << endl;
204 return kFALSE;
205 }
206
207 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
208 if (!fCam)
209 {
210 *fLog << err
211 << "MPadON : MGeomCam not found (no geometry information available)... aborting."
212 << endl;
213 return kFALSE;
214 }
215
216 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
217 if (!fEvt)
218 {
219 *fLog << err << "MPadON : MCerPhotEvt not found... aborting." << endl;
220 return kFALSE;
221 }
222
223 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
224 if (!fSigmabar)
225 {
226 *fLog << err << "MPadON : MSigmabar not found... aborting." << endl;
227 return kFALSE;
228 }
229
230 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
231 if (!fBlinds)
232 {
233 *fLog << err << "MPadON : MBlindPixels not found... aborting." << endl;
234 return kFALSE;
235 }
236
237
238 //--------------------------------------------------------------------
239 // histograms for checking the padding
240 //
241 // plot pedestal sigmas
242 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
243 100, 0.0, 5.0, 100, 0.0, 5.0);
244 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
245 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
246
247 // plot no.of photons (before vs. after padding)
248 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
249 100, -10.0, 90.0, 100, -10, 90);
250 fHPhotons->SetXTitle("no.of photons before padding");
251 fHPhotons->SetYTitle("no.of photons after padding");
252
253 // plot of added NSB
254 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
255 fHNSB->SetXTitle("no.of added NSB photons");
256 fHNSB->SetYTitle("no.of pixels");
257
258
259 //--------------------------------------------------------------------
260
261 fIter = 20;
262
263 memset(fErrors, 0, sizeof(fErrors));
264
265 return kTRUE;
266}
267
268// --------------------------------------------------------------------------
269//
270// Do the Padding
271// idealy the events to be padded should have been generated without NSB
272//
273Int_t MPadON::Process()
274{
275 *fLog << all << "Entry MPadON::Process();" << endl;
276
277 // this is the conversion factor from the number of photons
278 // to the number of photo electrons
279 // later to be taken from MCalibration
280 // fPEperPhoton = PW * LG * QE * 1D
281 Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9;
282
283 //------------------------------------------------
284 // add pixels to MCerPhotEvt which are not yet in;
285 // set their numnber of photons equal to zero
286
287 UInt_t imaxnumpix = fCam->GetNumPixels();
288
289 for (UInt_t i=0; i<imaxnumpix; i++)
290 {
291 Bool_t alreadythere = kFALSE;
292 UInt_t npix = fEvt->GetNumPixels();
293 for (UInt_t j=0; j<npix; j++)
294 {
295 MCerPhotPix &pix = (*fEvt)[j];
296 UInt_t id = pix.GetPixId();
297 if (i==id)
298 {
299 alreadythere = kTRUE;
300 break;
301 }
302 }
303 if (alreadythere)
304 continue;
305
306 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetPedestalRms());
307 }
308
309
310
311 //-----------------------------------------
312 Int_t rc=0;
313
314 const UInt_t npix = fEvt->GetNumPixels();
315
316 //fSigmabar->Calc(*fCam, *fPed, *fEvt);
317 //*fLog << all << "before padding : " << endl;
318 //fSigmabar->Print("");
319
320
321 //$$$$$$$$$$$$$$$$$$$$$$$$$$
322 // to simulate the situation that before the padding the NSB and
323 // electronic noise are zero : set Sigma = 0 for all pixels
324 //for (UInt_t i=0; i<npix; i++)
325 //{
326 // MCerPhotPix &pix = fEvt->operator[](i);
327 // Int_t j = pix.GetPixId();
328
329 // MPedestalPix &ppix = fPed->operator[](j);
330 // ppix.SetMeanRms(0.0);
331 //}
332 //$$$$$$$$$$$$$$$$$$$$$$$$$$
333
334 //-------------------------------------------
335 // Calculate average sigma of the event
336 //
337 Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt);
338 *fLog << all << "before padding : " << endl;
339 fSigmabar->Print("");
340 Double_t sigbarold = sigbarold_phot * fPEperPhoton;
341
342 Double_t sigbarold2 = sigbarold*sigbarold;
343
344
345 // for MC data : expect sigmabar to be zero before the padding
346 if (sigbarold_phot > 0)
347 {
348 *fLog << err
349 << "MPadON::Process(); sigmabar of event to be padded is > 0 : "
350 << sigbarold_phot << ". Stop event loop " << endl;
351 // input data should have sigmabar = 0; stop event loop
352
353 rc = 1;
354 fErrors[rc]++;
355 return kERROR;
356 }
357
358 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
359 *fLog << all << "theta = " << theta << endl;
360
361
362 //-------------------------------------------
363 // for the current theta,
364 // generate blind pixels according to the histograms
365 // fHBlindPixNTheta and fHBlindPixIDTheta
366 //
367
368
369 Int_t binPix = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
370
371 if ( binPix < 1 || binPix > fHBlindPixNTheta->GetNbinsX() )
372 {
373 *fLog << warn
374 << "MPadON::Process(); binNumber out of range : theta, binPix = "
375 << theta << ", " << binPix << "; Skip event " << endl;
376 // event cannot be padded; skip event
377
378 rc = 2;
379 fErrors[rc]++;
380 return kCONTINUE;
381 }
382
383 // numBlind is the number of blind pixels in this event
384 TH1D *nblind;
385 UInt_t numBlind;
386
387 nblind = fHBlindPixNTheta->ProjectionY("", binPix, binPix, "");
388 if ( nblind->GetEntries() == 0.0 )
389 {
390 *fLog << warn << "MPadON::Process(); projection for Theta bin "
391 << binPix << " has no entries; Skip event " << endl;
392 // event cannot be padded; skip event
393 delete nblind;
394
395 rc = 7;
396 fErrors[rc]++;
397 return kCONTINUE;
398 }
399 else
400 {
401 numBlind = (Int_t) (nblind->GetRandom()+0.5);
402 }
403 delete nblind;
404
405
406 // throw the Id of numBlind different pixels in this event
407 TH1D *hblind;
408 UInt_t idBlind;
409 UInt_t listId[npix];
410 UInt_t nlist = 0;
411 Bool_t equal;
412
413 hblind = fHBlindPixIdTheta->ProjectionY("", binPix, binPix, "");
414 if ( hblind->GetEntries() > 0.0 )
415 for (UInt_t i=0; i<numBlind; i++)
416 {
417 while(1)
418 {
419 idBlind = (Int_t) (hblind->GetRandom()+0.5);
420 equal = kFALSE;
421 for (UInt_t j=0; j<nlist; j++)
422 if (idBlind == listId[j])
423 {
424 equal = kTRUE;
425 break;
426 }
427 if (!equal) break;
428 }
429 listId[nlist] = idBlind;
430 nlist++;
431
432 fBlinds->SetPixelBlind(idBlind);
433 *fLog << "idBlind = " << idBlind << endl;
434 }
435 fBlinds->SetReadyToSave();
436
437 delete hblind;
438
439
440 //-------------------------------------------
441 // for the current theta,
442 // generate a sigmabar according to the histogram fHSigmaTheta
443 //
444 Double_t sigmabar_phot = 0;
445 Double_t sigmabar = 0;
446 Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta);
447
448 if (binPix != binNumber)
449 {
450 *fLog << err
451 << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting"
452 << endl;
453 return kERROR;
454 }
455
456 TH1D *hsigma;
457
458 hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, "");
459 if ( hsigma->GetEntries() == 0.0 )
460 {
461 *fLog << warn << "MPadON::Process(); projection for Theta bin "
462 << binNumber << " has no entries; Skip event " << endl;
463 // event cannot be padded; skip event
464 delete hsigma;
465
466 rc = 3;
467 fErrors[rc]++;
468 return kCONTINUE;
469 }
470 else
471 {
472 sigmabar_phot = hsigma->GetRandom();
473 sigmabar = sigmabar_phot * fPEperPhoton;
474 *fLog << all << "Theta, bin number = " << theta << ", " << binNumber
475 << ", sigmabar_phot = " << sigmabar_phot << endl;
476 }
477 delete hsigma;
478
479 const Double_t sigmabar2 = sigmabar*sigmabar;
480
481 //-------------------------------------------
482
483 *fLog << all << "MPadON::Process(); sigbarold, sigmabar = "
484 << sigbarold << ", "<< sigmabar << endl;
485
486 // Skip event if target sigmabar is <= sigbarold
487 if (sigmabar <= sigbarold)
488 {
489 *fLog << all
490 << "MPadON::Process(); target sigmabar is less than sigbarold : "
491 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
492
493 rc = 4;
494 fErrors[rc]++;
495 return kCONTINUE;
496 }
497
498
499 //-------------------------------------------
500 //
501 // Calculate average number of NSB photons to be added (lambdabar)
502 // from the value of sigmabar,
503 // - making assumptions about the average electronic noise (elNoise2) and
504 // - using a fixed value (F2excess) for the excess noise factor
505
506 Double_t elNoise2; // [photo electrons]
507 Double_t F2excess = 1.3;
508 Double_t lambdabar; // [photo electrons]
509
510
511
512 Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta);
513 if (binTheta != binNumber)
514 {
515 *fLog << err
516 << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting"
517 << endl;
518 return kERROR;
519 }
520
521 // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
522 // The average electronic noise (to be added) has to be well above this RMS,
523 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
524 // may become negative
525
526 TH1D *hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
527 0, 9999, "");
528 Double_t RMS_phot = hnoise->GetRMS(1);
529 Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
530 delete hnoise;
531
532 elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2);
533 *fLog << all << "elNoise2 = " << elNoise2 << endl;
534
535 lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess;
536
537 // This value of lambdabar is the same for all pixels;
538 // note that lambdabar is normalized to the area of pixel 0
539
540 //---------- start loop over pixels ---------------------------------
541 // do the padding for each pixel
542 //
543 // pad only pixels - which are used (before image cleaning)
544 //
545 Double_t sig_phot = 0;
546 Double_t sig = 0;
547
548 Double_t sigma2 = 0;
549
550 Double_t diff_phot = 0;
551 Double_t diff = 0;
552
553 Double_t addSig2_phot= 0;
554 Double_t addSig2 = 0;
555
556 Double_t elNoise2Pix = 0;
557
558
559 for (UInt_t i=0; i<npix; i++)
560 {
561 MCerPhotPix &pix = (*fEvt)[i];
562 if ( !pix.IsPixelUsed() )
563 continue;
564
565 //if ( pix.GetNumPhotons() == 0.0)
566 //{
567 // *fLog << warn
568 // << "MPadON::Process(); no.of photons is 0 for used pixel"
569 // << endl;
570 // continue;
571 //}
572
573 Int_t j = pix.GetPixId();
574
575 // GetPixRatio returns (area of pixel 0 / area of current pixel)
576 Double_t ratioArea = 1.0 / fCam->GetPixRatio(j);
577
578 MPedestalPix &ppix = (*fPed)[j];
579 Double_t oldsigma_phot = ppix.GetPedestalRms();
580 Double_t oldsigma = oldsigma_phot * fPEperPhoton;
581 Double_t oldsigma2 = oldsigma*oldsigma;
582
583 //---------------------------------
584 // throw the Sigma for this pixel
585 //
586 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
587
588 Int_t count;
589 Bool_t ok;
590
591 TH1D *hdiff;
592 TH1D *hsig;
593
594 switch (fPadFlag)
595 {
596 case 1 :
597 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
598
599 hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
600 binPixel, binPixel, "");
601 if ( hdiff->GetEntries() == 0 )
602 {
603 *fLog << warn << "MPadON::Process(); projection for Theta bin "
604 << binTheta << " and pixel bin " << binPixel
605 << " has no entries; aborting " << endl;
606 delete hdiff;
607
608 rc = 5;
609 fErrors[rc]++;
610 return kCONTINUE;
611 }
612
613 count = 0;
614 ok = kFALSE;
615 for (Int_t m=0; m<fIter; m++)
616 {
617 count += 1;
618 diff_phot = hdiff->GetRandom();
619 diff = diff_phot * fPEperPhoton * fPEperPhoton;
620
621 // the following condition ensures that elNoise2Pix > 0.0
622 if ( (diff + sigmabar2 - oldsigma2/ratioArea
623 - lambdabar*F2excess) > 0.0 )
624 {
625 ok = kTRUE;
626 break;
627 }
628 }
629 if (!ok)
630 {
631 *fLog << all << "theta, j, count, sigmabar, diff = " << theta << ", "
632 << j << ", " << count << ", " << sigmabar << ", "
633 << diff << endl;
634 diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
635 }
636 delete hdiff;
637 sigma2 = diff + sigmabar2;
638 break;
639
640 case 2 :
641 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
642
643 hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
644 binPixel, binPixel, "");
645 if ( hsig->GetEntries() == 0 )
646 {
647 *fLog << err << "MPadON::Process(); projection for Theta bin "
648 << binTheta << " and pixel bin " << binPixel
649 << " has no entries; aborting " << endl;
650 delete hsig;
651
652 rc = 6;
653 fErrors[rc]++;
654 return kCONTINUE;
655 }
656
657 count = 0;
658 ok = kFALSE;
659 for (Int_t m=0; m<fIter; m++)
660 {
661 count += 1;
662
663 sig_phot = hsig->GetRandom();
664 sig = sig_phot * fPEperPhoton;
665
666 sigma2 = sig*sig/ratioArea;
667 // the following condition ensures that elNoise2Pix > 0.0
668
669 if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 )
670 {
671 ok = kTRUE;
672 break;
673 }
674 }
675 if (!ok)
676 {
677
678 *fLog << "theta, j, count, sigmabar, sig = " << theta << ", "
679 << j << ", " << count << ", " << sigmabar << ", "
680 << sig << endl;
681 sigma2 = lambdabar*F2excess + oldsigma2/ratioArea;
682 }
683 delete hsig;
684 break;
685 }
686
687 //---------------------------------
688 // get the additional sigma^2 for this pixel (due to the padding)
689
690 addSig2 = sigma2*ratioArea - oldsigma2;
691 addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
692
693 //---------------------------------
694 // get the additional electronic noise for this pixel
695
696 elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;
697
698
699 //---------------------------------
700 // throw actual number of additional NSB photons (NSB)
701 // and its RMS (sigmaNSB)
702
703 Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
704 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
705 Double_t sigmaNSB0;
706
707 if (arg >= 0)
708 {
709 sigmaNSB0 = sqrt( arg );
710 }
711 else
712 {
713 sigmaNSB0 = 0.0000001;
714 if (arg < -1.e-10)
715 {
716 *fLog << warn << "MPadON::Process(); argument of sqrt < 0 : "
717 << arg << endl;
718 }
719 }
720
721
722 //---------------------------------
723 // smear NSB0 according to sigmaNSB0
724 // and subtract lambdabar because of AC coupling
725
726 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
727 Double_t NSB_phot = NSB / fPEperPhoton;
728
729 //---------------------------------
730
731 // add additional NSB to the number of photons
732 Double_t oldphotons_phot = pix.GetNumPhotons();
733 Double_t oldphotons = oldphotons_phot * fPEperPhoton;
734 Double_t newphotons = oldphotons + NSB;
735 Double_t newphotons_phot = newphotons / fPEperPhoton;
736 pix.SetNumPhotons(newphotons_phot);
737
738
739 fHNSB->Fill( NSB_phot/sqrt(ratioArea) );
740 fHPhotons->Fill( oldphotons_phot/ratioArea,
741 newphotons_phot/ratioArea );
742
743
744 // error: add sigma of padded noise quadratically
745 Double_t olderror_phot = pix.GetErrorPhot();
746 Double_t newerror_phot =
747 sqrt(olderror_phot*olderror_phot + addSig2_phot);
748 pix.SetErrorPhot(newerror_phot);
749
750
751 Double_t newsigma = sqrt(oldsigma2 + addSig2);
752 Double_t newsigma_phot = newsigma / fPEperPhoton;
753 ppix.SetPedestalRms(newsigma_phot);
754
755 fHSigmaPedestal->Fill(oldsigma_phot, newsigma_phot);
756 }
757 //---------- end of loop over pixels ---------------------------------
758
759 // Calculate sigmabar again and crosscheck
760
761 fSigmabar->Calc(*fCam, *fPed, *fEvt);
762 *fLog << all << "MPadON::Process(); after padding : " << endl;
763 fSigmabar->Print("");
764
765
766 *fLog << all << "Exit MPadON::Process();" << endl;
767
768 rc = 0;
769 fErrors[rc]++;
770
771 return kTRUE;
772
773}
774
775// --------------------------------------------------------------------------
776//
777//
778Int_t MPadON::PostProcess()
779{
780 if (GetNumExecutions() != 0)
781 {
782
783 *fLog << inf << endl;
784 *fLog << GetDescriptor() << " execution statistics:" << endl;
785 *fLog << dec << setfill(' ');
786 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
787 << (int)(fErrors[1]*100/GetNumExecutions())
788 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
789
790 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
791 << (int)(fErrors[2]*100/GetNumExecutions())
792 << "%) Evts skipped due to: Zenith angle out of range" << endl;
793
794 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
795 << (int)(fErrors[3]*100/GetNumExecutions())
796 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
797
798 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
799 << (int)(fErrors[4]*100/GetNumExecutions())
800 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
801
802 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
803 << (int)(fErrors[5]*100/GetNumExecutions())
804 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
805
806 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
807 << (int)(fErrors[6]*100/GetNumExecutions())
808 << "%) Evts skipped due to: No data for generating Sigma" << endl;
809
810 *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3)
811 << (int)(fErrors[7]*100/GetNumExecutions())
812 << "%) Evts skipped due to: No data for generating Blind pixels" << endl;
813
814 *fLog << " " << fErrors[0] << " ("
815 << (int)(fErrors[0]*100/GetNumExecutions())
816 << "%) Evts survived the padding!" << endl;
817 *fLog << endl;
818
819 }
820
821 //---------------------------------------------------------------
822 TCanvas &c = *(MH::MakeDefCanvas("PadON", "", 900, 1200));
823 c.Divide(3, 4);
824 gROOT->SetSelectedPad(NULL);
825
826 c.cd(1);
827 fHNSB->SetDirectory(NULL);
828 fHNSB->DrawCopy();
829 fHNSB->SetBit(kCanDelete);
830
831 c.cd(2);
832 fHSigmaPedestal->SetDirectory(NULL);
833 fHSigmaPedestal->DrawCopy();
834 fHSigmaPedestal->SetBit(kCanDelete);
835
836 c.cd(3);
837 fHPhotons->SetDirectory(NULL);
838 fHPhotons->DrawCopy();
839 fHPhotons->SetBit(kCanDelete);
840
841 //--------------------------------------------------------------------
842
843
844 c.cd(4);
845 fHSigmaTheta->SetDirectory(NULL);
846 fHSigmaTheta->SetTitle("(Input) 2D : Sigmabar, \\Theta");
847 fHSigmaTheta->DrawCopy();
848 fHSigmaTheta->SetBit(kCanDelete);
849
850 //--------------------------------------------------------------------
851
852
853 c.cd(7);
854 fHBlindPixNTheta->SetDirectory(NULL);
855 fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta");
856 fHBlindPixNTheta->DrawCopy();
857 fHBlindPixNTheta->SetBit(kCanDelete);
858
859 //--------------------------------------------------------------------
860
861
862 c.cd(10);
863 fHBlindPixIdTheta->SetDirectory(NULL);
864 fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta");
865 fHBlindPixIdTheta->DrawCopy();
866 fHBlindPixIdTheta->SetBit(kCanDelete);
867
868
869 //--------------------------------------------------------------------
870 // draw the 3D histogram (input): Theta, pixel, Sigma^2-Sigmabar^2
871
872 c.cd(5);
873 TH2D *l1;
874 l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
875 l1->SetDirectory(NULL);
876 l1->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
877 l1->SetXTitle("\\Theta [\\circ]");
878 l1->SetYTitle("Sigma^2-Sigmabar^2");
879
880 l1->DrawCopy("box");
881 l1->SetBit(kCanDelete);;
882
883 c.cd(8);
884 TH2D *l2;
885 l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
886 l2->SetDirectory(NULL);
887 l2->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
888 l2->SetXTitle("pixel");
889 l2->SetYTitle("Sigma^2-Sigmabar^2");
890
891 l2->DrawCopy("box");
892 l2->SetBit(kCanDelete);;
893
894 //--------------------------------------------------------------------
895 // draw the 3D histogram (input): Theta, pixel, Sigma
896
897 c.cd(6);
898 TH2D *k1;
899 k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
900 k1->SetDirectory(NULL);
901 k1->SetTitle("(Input) Sigma vs. \\Theta (all pixels)");
902 k1->SetXTitle("\\Theta [\\circ]");
903 k1->SetYTitle("Sigma");
904
905 k1->DrawCopy("box");
906 k1->SetBit(kCanDelete);
907
908 c.cd(9);
909 TH2D *k2;
910 k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
911 k2->SetDirectory(NULL);
912 k2->SetTitle("(Input) Sigma vs. pixel number (all \\Theta)");
913 k2->SetXTitle("pixel");
914 k2->SetYTitle("Sigma");
915
916 k2->DrawCopy("box");
917 k2->SetBit(kCanDelete);;
918
919
920 //--------------------------------------------------------------------
921
922
923 return kTRUE;
924}
925
926// --------------------------------------------------------------------------
927
928
929
930
931
932
933
934
935
936
937
938
939
940
Note: See TracBrowser for help on using the repository browser.