source: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc@ 2155

Last change on this file since 2155 was 2155, checked in by wittek, 21 years ago
*** empty log message ***
File size: 37.7 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 <mailto:wittek@mppmu.mpg.de> 06/2003
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MPadONOFF
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// 2) fPadFlag = 2 :
51// Generate a pedestal sigma for each pixel using the 3D-histogram
52// Theta, pixel no., Sigma (fHSigmaPixTheta).
53//
54//
55// The padding has to be done before the image cleaning because the
56// image cleaning depends on the pedestal sigmas.
57//
58// For random numbers gRandom is used.
59//
60// This implementation has been tested for CT1 data. For MAGIC some
61// modifications are necessary.
62//
63/////////////////////////////////////////////////////////////////////////////
64#include "MPadONOFF.h"
65
66#include <stdio.h>
67
68#include <TH1.h>
69#include <TH2.h>
70#include <TH3.h>
71#include <TRandom.h>
72#include <TCanvas.h>
73#include <TFile.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#include "MBlindPixels.h"
89
90#include "MRead.h"
91#include "MFilterList.h"
92#include "MTaskList.h"
93#include "MBlindPixelCalc.h"
94#include "MHBlindPixels.h"
95#include "MFillH.h"
96#include "MHSigmaTheta.h"
97#include "MEvtLoop.h"
98
99
100ClassImp(MPadONOFF);
101
102
103// --------------------------------------------------------------------------
104//
105// Default constructor.
106//
107MPadONOFF::MPadONOFF(const char *name, const char *title)
108{
109 fName = name ? name : "MPadONOFF";
110 fTitle = title ? title : "Task for the ON-OFF padding";
111
112 fPadFlag = 1;
113 *fLog << "MPadONOFF: fPadFlag = " << fPadFlag << endl;
114
115 fType = "";
116
117 fHSigmaTheta = NULL;
118 fHSigmaPixTheta = NULL;
119 fHDiffPixTheta = NULL;
120
121 fHgON = NULL;
122 fHgOFF = NULL;
123
124 fHBlindPixIdTheta = NULL;
125 fHBlindPixNTheta = NULL;
126
127 fHSigmaPedestal = NULL;
128 fHPhotons = NULL;
129 fHNSB = NULL;
130}
131
132// --------------------------------------------------------------------------
133//
134// Destructor.
135//
136MPadONOFF::~MPadONOFF()
137{
138 if (fHSigmaTheta != NULL) delete fHSigmaTheta;
139
140 if (fHSigmaPedestal != NULL) delete fHSigmaPedestal;
141 if (fHPhotons != NULL) delete fHPhotons;
142 if (fHNSB != NULL) delete fHNSB;
143
144 if (fInfile != NULL) delete fInfile;
145}
146
147// --------------------------------------------------------------------------
148//
149// Merge the distributions of ON and OFF data
150//
151// fHSigmaTheta 2D-histogram (Theta, sigmabar)
152// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma)
153// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2)
154// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id)
155// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels )
156//
157// and define the target distributions for the padding
158//
159Bool_t MPadONOFF::MergeHistograms(TH2D *sigthon, TH2D *sigthoff,
160 TH3D *sigpixthon, TH3D *sigpixthoff,
161 TH3D *diffpixthon, TH3D *diffpixthoff,
162 TH2D *blindidthon, TH2D *blindidthoff,
163 TH2D *blindnthon, TH2D *blindnthoff)
164{
165 //-------------------------------------------------------------
166 // merge the distributions of ON and OFF events
167 // to obtain the target distribution fHSigmaTheta
168 //
169
170 fHSigmaTheta = new TH2D( (TH2D&)*sigthon );
171 fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");
172
173 Int_t nbinstheta = sigthon->GetNbinsX();
174 Int_t n = sigthon->GetNbinsY();
175 Double_t eps = 1.e-10;
176
177 fHgON = new TH3D;
178 fHgON->SetBins(nbinstheta, 0.5, 0.5+(Double_t)nbinstheta,
179 n, 0.5, 0.5+(Double_t)n,
180 n, 0.5, 0.5+(Double_t)n);
181 fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");
182
183 fHgOFF = new TH3D;
184 fHgOFF->SetBins(nbinstheta, 0.5, 0.5+(Double_t)nbinstheta,
185 n, 0.5, 0.5+(Double_t)n,
186 n, 0.5, 0.5+(Double_t)n);
187 fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF");
188
189 //............ loop over Theta bins ...........................
190
191 // hista is the normalized 1D histogram of sigmabar for ON data
192 // histb is the normalized 1D histogram of sigmabar for OFF data
193 // histc is the difference ON-OFF
194
195 // at the beginning, histap is a copy of hista
196 // at the end, it will be the 1D histogram for ON data after padding
197
198 // at the beginning, histbp is a copy of histb
199 // at the end, it will be the 1D histogram for OFF data after padding
200
201 // at the beginning, histcp is a copy of histc
202 // at the end, it should be zero
203
204 for (Int_t l=1; l<=nbinstheta; l++)
205 {
206 TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, "");
207 Stat_t suma = hista->Integral();
208 hista->Scale(1./suma);
209
210 TH1D *histap = new TH1D( (const TH1D&)*hista );
211
212 TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, "");
213 Stat_t sumb = histb->Integral();
214 histb->Scale(1./sumb);
215
216 TH1D *histbp = new TH1D( (const TH1D&)*histb );
217
218 TH1D *histc = new TH1D( (const TH1D&)*hista );
219 histc->Add(histb, -1.0);
220
221 TH1D *histcp = new TH1D( (const TH1D&)*histc );
222
223
224 // calculate matrix g
225
226 // fHg[k][j] (if <0.0) tells how many ON events in bin k should be padded
227 // from sigma_k to sigma_j
228
229
230 // fHg[k][j] (if >0.0) tells how many OFF events in bin k should be padded
231 // from sigma_k to sigma_j
232
233 //-------- start j loop ------------------------------------------------
234 // loop over bins in histc, starting from the end
235 Double_t v, s, w, t, x, u, a, b, c, arg;
236
237 for (Int_t j=n; j >= 1; j--)
238 {
239 v = histcp->GetBinContent(j);
240 if ( fabs(v) < eps ) continue;
241 if (v >= 0.0)
242 s = 1.0;
243 else
244 s = -1.0;
245
246 //................ start k loop ......................................
247 // look for a bin k which may compensate the content of bin j
248 for (Int_t k=j-1; k >= 1; k--)
249 {
250 w = histcp->GetBinContent(k);
251 if ( fabs(w) < eps ) continue;
252 if (w >= 0.0)
253 t = 1.0;
254 else
255 t = -1.0;
256
257 if (s==t) continue;
258
259 x = v + w;
260 if (x >= 0.0)
261 u = 1.0;
262 else
263 u = -1.0;
264
265 if (u == s)
266 {
267 arg = -w;
268
269 if (arg <=0.0)
270 {
271 fHgON->SetBinContent (l, k, j, -arg);
272 fHgOFF->SetBinContent(l, k, j, 0.0);
273 }
274 else
275 {
276 fHgON->SetBinContent (l, k, j, 0.0);
277 fHgOFF->SetBinContent(l, k, j, arg);
278 }
279
280 *fLog << "k, j, arg = " << k << ", " << j << ", " << arg << endl;
281
282 // g(k-1, j-1) = arg;
283 //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", "
284 // << arg << endl;
285
286 //......................................
287 // this is for checking the procedure
288 if (arg < 0.0)
289 {
290 a = histap->GetBinContent(k);
291 histap->SetBinContent(k, a+arg);
292 a = histap->GetBinContent(j);
293 histap->SetBinContent(j, a-arg);
294 }
295 else
296 {
297 b = histbp->GetBinContent(k);
298 histbp->SetBinContent(k, b-arg);
299 b = histbp->GetBinContent(j);
300 histbp->SetBinContent(j, b+arg);
301 }
302 //......................................
303
304 histcp->SetBinContent(k, 0.0);
305 histcp->SetBinContent(j, x);
306
307 //......................................
308 // redefine v
309 v = histcp->GetBinContent(j);
310 if ( fabs(v) < eps ) break;
311 if (v >= 0.0)
312 s = 1.0;
313 else
314 s = -1.0;
315 //......................................
316
317 continue;
318 }
319
320 arg = v;
321
322 if (arg <=0.0)
323 {
324 fHgON->SetBinContent (l, k, j, -arg);
325 fHgOFF->SetBinContent(l, k, j, 0.0);
326 }
327 else
328 {
329 fHgON->SetBinContent (l, k, j, 0.0);
330 fHgOFF->SetBinContent(l, k, j, arg);
331 }
332
333 *fLog << "k, j, arg = " << k << ", " << j << ", " << arg << endl;
334 //g(k-1, j-1) = arg;
335 //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", "
336 // << arg << endl;
337
338 //......................................
339 // this is for checking the procedure
340 if (arg < 0.0)
341 {
342 a = histap->GetBinContent(k);
343 histap->SetBinContent(k, a+arg);
344 a = histap->GetBinContent(j);
345 histap->SetBinContent(j, a-arg);
346 }
347 else
348 {
349 b = histbp->GetBinContent(k);
350
351 histbp->SetBinContent(k, b-arg);
352 b = histbp->GetBinContent(j);
353 histbp->SetBinContent(j, b+arg);
354 }
355 //......................................
356
357 histcp->SetBinContent(k, x);
358 histcp->SetBinContent(j, 0.0);
359
360 break;
361 }
362 //................ end k loop ......................................
363 }
364 //-------- end j loop ------------------------------------------------
365
366 // check results for this Theta bin
367 for (Int_t j=1; j<=n; j++)
368 {
369 a = histap->GetBinContent(j);
370 b = histbp->GetBinContent(j);
371 c = histcp->GetBinContent(j);
372
373 if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps )
374 *fLog << "MPadONOFF::Read; inconsistency in results; a, b, c = "
375 << a << ", " << b << ", " << c << endl;
376 }
377
378
379 // fill target distribution SigmaTheta
380 // for this Theta bin
381 //
382 for (Int_t j=1; j<=n; j++)
383 {
384 a = histap->GetBinContent(j);
385 fHSigmaTheta->SetBinContent(l, j, a);
386 }
387
388 delete hista;
389 delete histb;
390 delete histc;
391
392 delete histap;
393 delete histbp;
394 delete histcp;
395 }
396 //............ end of loop over Theta bins ....................
397
398
399 // target distributions
400 // SigmaPixTheta and DiffPixTheta
401 // BlindPixIdTheta and BlindPixNTheta
402 // are taken from the ON data
403
404 fHSigmaPixTheta = sigpixthon;
405 fHDiffPixTheta = diffpixthon;
406
407 fHBlindPixIdTheta = blindidthon;
408 fHBlindPixNTheta = blindnthon;
409
410 //--------------------------------------------
411
412 fHSigmaTheta->SetDirectory(NULL);
413 fHSigmaPixTheta->SetDirectory(NULL);
414 fHDiffPixTheta->SetDirectory(NULL);
415 fHBlindPixIdTheta->SetDirectory(NULL);
416 fHBlindPixNTheta->SetDirectory(NULL);
417 fHBlindPixNTheta->SetDirectory(NULL);
418
419 fHgON->SetDirectory(NULL);
420 fHgOFF->SetDirectory(NULL);
421
422
423 return kTRUE;
424}
425
426
427// --------------------------------------------------------------------------
428//
429// Read target distributions from a file
430//
431//
432Bool_t MPadONOFF::ReadTargetDist(const char* namefilein)
433{
434 *fLog << "Read padding histograms from file " << namefilein << endl;
435
436 fInfile = new TFile(namefilein);
437 fInfile->ls();
438
439 //------------------------------------
440
441 fHSigmaTheta =
442 (TH2D*) gROOT->FindObject("2D-ThetaSigmabar");
443 if (!fHSigmaTheta)
444 {
445 *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl;
446 return kFALSE;
447 }
448 *fLog << "Object '2D-ThetaSigmabar' was read in" << endl;
449
450 *fLog << "ReadTargetDist : fHSigmaTheta = " << fHSigmaTheta << endl;
451
452 fHSigmaPixTheta =
453 (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
454 if (!fHSigmaPixTheta)
455 {
456 *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl;
457 return kFALSE;
458 }
459 *fLog << "Object '3D-ThetaPixSigma' was read in" << endl;
460
461 fHDiffPixTheta =
462 (TH3D*) gROOT->FindObject("3D-ThetaPixDiff");
463 if (!fHDiffPixTheta)
464 {
465 *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl;
466 return kFALSE;
467 }
468 *fLog << "Object '3D-ThetaPixDiff' was read in" << endl;
469
470 fHgON =
471 (TH3D*) gROOT->FindObject("3D-PaddingMatrixON");
472 if (!fHgON)
473 {
474 *fLog << "Object '3D-PaddingMatrixON' not found on root file" << endl;
475 return kFALSE;
476 }
477 *fLog << "Object '3D-PaddingMatrixON' was read in" << endl;
478
479 fHgOFF =
480 (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF");
481 if (!fHgOFF)
482 {
483 *fLog << "Object '3D-PaddingMatrixOFF' not found on root file" << endl;
484 return kFALSE;
485 }
486 *fLog << "Object '3D-PaddingMatrixOFF' was read in" << endl;
487
488
489 fHBlindPixIdTheta =
490 (TH2D*) gROOT->FindObject("2D-IdBlindPixels");
491 if (!fHBlindPixIdTheta)
492 {
493 *fLog << "Object '2D-IdBlindPixels' not found on root file" << endl;
494 return kFALSE;
495 }
496 *fLog << "Object '2D-IdBlindPixels' was read in" << endl;
497
498 fHBlindPixNTheta =
499 (TH2D*) gROOT->FindObject("2D-NBlindPixels");
500 if (!fHBlindPixNTheta)
501 {
502 *fLog << "Object '2D-NBlindPixels' not found on root file" << endl;
503 return kFALSE;
504 }
505 *fLog << "Object '2D-NBlindPixels' was read in" << endl;
506
507 *fLog << "ReadTargetDist : fHBlindPixNTheta = " << fHBlindPixNTheta << endl;
508
509 //------------------------------------
510
511 return kTRUE;
512}
513
514
515// --------------------------------------------------------------------------
516//
517// Write target distributions onto a file
518//
519//
520Bool_t MPadONOFF::WriteTargetDist(const char* namefileout)
521{
522 *fLog << "Write padding histograms onto file " << namefileout << endl;
523
524 TFile outfile(namefileout, "RECREATE");
525
526 fHBlindPixNTheta->Write();
527 fHBlindPixIdTheta->Write();
528
529 fHSigmaTheta->Write();
530 fHSigmaPixTheta->Write();
531 fHDiffPixTheta->Write();
532
533 fHgON->Write();
534 fHgOFF->Write();
535
536 *fLog << "MPadONOFF::WriteTargetDist; target histograms written onto file "
537 << namefileout << endl;
538
539 return kTRUE;
540}
541
542// --------------------------------------------------------------------------
543//
544// Set type of data to be padded
545//
546// this is not necessary if the type of the data can be recognized
547// directly from the input
548//
549//
550void MPadONOFF::SetDataType(const char* type)
551{
552 fType = type;
553 *fLog << "MPadONOFF::SetDataType(); type of data to be padded : "
554 << fType << endl;
555
556 return;
557}
558
559
560// --------------------------------------------------------------------------
561//
562// Set the option for the padding
563//
564// There are 2 options for the padding :
565//
566// 1) fPadFlag = 1 :
567// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
568// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
569// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
570// (fHDiffPixTheta).
571//
572// This is the preferred option as it takes into account the
573// correlations between the Sigma of a pixel and Sigmabar.
574//
575// 2) fPadFlag = 2 :
576// Generate a pedestal sigma for each pixel using the 3D-histogram
577// Theta, pixel no., Sigma (fHSigmaPixTheta).
578//
579void MPadONOFF::SetPadFlag(Int_t padflag)
580{
581 fPadFlag = padflag;
582 *fLog << "MPadONOFF::SetPadFlag(); choose option " << fPadFlag << endl;
583}
584
585// --------------------------------------------------------------------------
586//
587// Set the pointers and prepare the histograms
588//
589Bool_t MPadONOFF::PreProcess(MParList *pList)
590{
591 if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta ||
592 !fHBlindPixIdTheta || !fHBlindPixNTheta ||
593 !fHgON || !fHgOFF)
594 {
595 *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl;
596 return kFALSE;
597 }
598
599 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
600 if (!fMcEvt)
601 {
602 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
603 return kFALSE;
604 }
605
606 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
607 if (!fPed)
608 {
609 *fLog << err << "MPedestalCam not found... aborting." << endl;
610 return kFALSE;
611 }
612
613 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
614 if (!fCam)
615 {
616 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
617 return kFALSE;
618 }
619
620 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
621 if (!fEvt)
622 {
623 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
624 return kFALSE;
625 }
626
627 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
628 if (!fSigmabar)
629 {
630 *fLog << err << "MSigmabar not found... aborting." << endl;
631 return kFALSE;
632 }
633
634 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
635 if (!fBlinds)
636 {
637 *fLog << err << "MBlindPixels not found... aborting." << endl;
638 return kFALSE;
639 }
640
641 if (fType !="ON" && fType !="OFF" && fType !="MC")
642 {
643 *fLog << err << "Type of data to be padded not defined... aborting." << endl;
644 return kFALSE;
645 }
646
647
648 //--------------------------------------------------------------------
649 // histograms for checking the padding
650 //
651 // plot pedestal sigmas
652 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
653 100, 0.0, 5.0, 100, 0.0, 5.0);
654 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
655 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
656
657 // plot no.of photons (before vs. after padding)
658 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
659 100, -10.0, 90.0, 100, -10, 90);
660 fHPhotons->SetXTitle("no.of photons before padding");
661 fHPhotons->SetYTitle("no.of photons after padding");
662
663 // plot of added NSB
664 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
665 fHNSB->SetXTitle("no.of added NSB photons");
666 fHNSB->SetYTitle("no.of pixels");
667
668
669 //--------------------------------------------------------------------
670
671 memset(fErrors, 0, sizeof(fErrors));
672
673 return kTRUE;
674}
675
676// --------------------------------------------------------------------------
677//
678// Do the Padding
679// idealy the events to be padded should have been generated without NSB
680//
681Bool_t MPadONOFF::Process()
682{
683 //*fLog << "Entry MPadONOFF::Process();" << endl;
684
685 //------------------------------------------------
686 // add pixels to MCerPhotEvt which are not yet in;
687 // set their number of photons equal to zero
688
689 UInt_t imaxnumpix = fCam->GetNumPixels();
690
691 for (UInt_t i=0; i<imaxnumpix; i++)
692 {
693 Bool_t alreadythere = kFALSE;
694 UInt_t npix = fEvt->GetNumPixels();
695 for (UInt_t j=0; j<npix; j++)
696 {
697 MCerPhotPix &pix = (*fEvt)[j];
698 UInt_t id = pix.GetPixId();
699 if (i==id)
700 {
701 alreadythere = kTRUE;
702 break;
703 }
704 }
705 if (alreadythere)
706 continue;
707
708 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetMeanRms());
709 }
710
711
712
713 //-----------------------------------------
714 Int_t rc=0;
715
716 const UInt_t npix = fEvt->GetNumPixels();
717
718 //fSigmabar->Calc(*fCam, *fPed, *fEvt);
719 //*fLog << "before padding : " << endl;
720 //fSigmabar->Print("");
721
722
723 //$$$$$$$$$$$$$$$$$$$$$$$$$$
724 // to simulate the situation that before the padding the NSB and
725 // electronic noise are zero : set Sigma = 0 for all pixels
726 //for (UInt_t i=0; i<npix; i++)
727 //{
728 // MCerPhotPix &pix = fEvt->operator[](i);
729 // Int_t j = pix.GetPixId();
730
731 // MPedestalPix &ppix = fPed->operator[](j);
732 // ppix.SetMeanRms(0.0);
733 //}
734 //$$$$$$$$$$$$$$$$$$$$$$$$$$
735
736 //-------------------------------------------
737 // Calculate average sigma of the event
738 //
739 Double_t sigbarold = fSigmabar->Calc(*fCam, *fPed, *fEvt);
740 Double_t sigbarold2 = sigbarold*sigbarold;
741 //fSigmabar->Print("");
742
743 // for MC data : expect sigmabar to be zero before the padding
744 if (fType == "MC")
745 {
746 if (sigbarold > 0)
747 {
748 //*fLog << "MPadONOFF::Process(); sigmabar of event to be padded is > 0 : "
749 // << sigbarold << ". Stop event loop " << endl;
750 // input data should have sigmabar = 0; stop event loop
751
752 rc = 1;
753 fErrors[rc]++;
754 return kCONTINUE;
755 }
756 }
757
758 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
759 // *fLog << "theta = " << theta << endl;
760
761 Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
762 if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() )
763 {
764 //*fLog << "MPadONOFF::Process(); binNumber out of range : theta, binTheta = "
765 // << theta << ", " << binTheta << "; Skip event " << endl;
766 // event cannot be padded; skip event
767
768 rc = 2;
769 fErrors[rc]++;
770 return kCONTINUE;
771 }
772
773 //-------------------------------------------
774 // for the current theta,
775 // generate blind pixels according to the histograms
776 // fHBlindPixNTheta and fHBlindPixIDTheta
777 //
778
779
780
781 // numBlind is the number of blind pixels in this event
782 TH1D *nblind;
783 UInt_t numBlind;
784
785 nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, "");
786 if ( nblind->GetEntries() == 0.0 )
787 {
788 //*fLog << "MPadONOFF::Process(); projection for Theta bin "
789 // << binTheta << " has no entries; Skip event " << endl;
790 // event cannot be padded; skip event
791 delete nblind;
792
793 rc = 7;
794 fErrors[rc]++;
795 return kCONTINUE;
796 }
797 else
798 {
799 numBlind = (Int_t) (nblind->GetRandom()+0.5);
800 }
801 delete nblind;
802
803
804 // throw the Id of numBlind different pixels in this event
805 TH1D *hblind;
806 UInt_t idBlind;
807 UInt_t listId[npix];
808 UInt_t nlist = 0;
809 Bool_t equal;
810
811 hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, "");
812 if ( hblind->GetEntries() > 0.0 )
813 for (UInt_t i=0; i<numBlind; i++)
814 {
815 while(1)
816 {
817 idBlind = (Int_t) (hblind->GetRandom()+0.5);
818 equal = kFALSE;
819 for (UInt_t j=0; j<nlist; j++)
820 if (idBlind == listId[j])
821 {
822 equal = kTRUE;
823 break;
824 }
825 if (!equal) break;
826 }
827 listId[nlist] = idBlind;
828 nlist++;
829
830 fBlinds->SetPixelBlind(idBlind);
831 //*fLog << "idBlind = " << idBlind << endl;
832 }
833 fBlinds->SetReadyToSave();
834
835 delete hblind;
836
837
838 //******************************************************************
839 // has the event to be padded ?
840 // if yes : to which sigmabar should it be padded ?
841 //
842
843 Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold);
844
845 Double_t prob;
846 TH1D *hpad = NULL;
847 if (fType == "ON")
848 {
849 hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");
850 prob = hpad->Integral();
851 }
852 else if (fType == "OFF")
853 {
854 hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, "");
855 prob = hpad->Integral();
856 }
857 else
858 prob = 1.0;
859
860 if ( fType != "MC" &&
861 (prob <= 0.0 || gRandom->Uniform() > prob) )
862 {
863 delete hpad;
864 // event should not be padded
865 rc = 8;
866 fErrors[rc]++;
867 return kTRUE;
868 }
869 // event should be padded
870
871 //-------------------------------------------
872 // for the current theta, generate a sigmabar
873 //
874 // for ON/OFF data according to the matrix fHg
875 // for MC data according to the histogram fHSigmaTheta
876 //
877 Double_t sigmabar=0;
878 if (fType == "ON" || fType == "OFF")
879 {
880 sigmabar = hpad->GetRandom();
881 delete hpad;
882 }
883
884 else if (fType == "MC")
885 {
886 Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta);
887
888 if (binTheta != bincheck)
889 {
890 cout << "The binnings of the 2 histograms are not identical; aborting"
891 << endl;
892 return kERROR;
893 }
894
895 TH1D *hsigma;
896
897 hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, "");
898 if ( hsigma->GetEntries() == 0.0 )
899 {
900 *fLog << "MPadONOFF::Process(); projection for Theta bin "
901 << binTheta << " has no entries; Skip event " << endl;
902 // event cannot be padded; skip event
903 delete hsigma;
904
905 rc = 3;
906 fErrors[rc]++;
907 return kCONTINUE;
908 }
909 else
910 {
911 sigmabar = hsigma->GetRandom();
912 //*fLog << "Theta, bin number = " << theta << ", " << binTheta
913 // << ", sigmabar = " << sigmabar << endl
914 }
915 delete hsigma;
916 }
917
918 const Double_t sigmabar2 = sigmabar*sigmabar;
919
920 //-------------------------------------------
921
922 //*fLog << "MPadONOFF::Process(); sigbarold, sigmabar = "
923 // << sigbarold << ", "<< sigmabar << endl;
924
925 // Skip event if target sigmabar is <= sigbarold
926 if (sigmabar <= sigbarold)
927 {
928 *fLog << "MPadONOFF::Process(); target sigmabar is less than sigbarold : "
929 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
930
931 rc = 4;
932 fErrors[rc]++;
933 return kCONTINUE;
934 }
935
936
937 //-------------------------------------------
938 //
939 // Calculate average number of NSB photons to be added (lambdabar)
940 // from the value of sigmabar,
941 // - making assumptions about the average electronic noise (elNoise2) and
942 // - using a fixed value (F2excess) for the excess noise factor
943
944 Double_t elNoise2; // [photons]
945 Double_t F2excess = 1.3;
946 Double_t lambdabar; // [photons]
947
948
949
950 Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta);
951 if (binTheta != bincheck)
952 {
953 cout << "The binnings of the 2 histograms are not identical; aborting"
954 << endl;
955 return kERROR;
956 }
957
958 // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
959 // The average electronic noise (to be added) has to be well above this RMS,
960 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
961 // may become negative
962
963 TH1D *hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
964 0, 9999, "");
965 Double_t RMS = hnoise->GetRMS(1);
966 delete hnoise;
967
968 elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2);
969 //*fLog << "elNoise2 = " << elNoise2 << endl;
970
971 lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photons]
972
973 // This value of lambdabar is the same for all pixels;
974 // note that lambdabar is normalized to the area of pixel 0
975
976 //---------- start loop over pixels ---------------------------------
977 // do the padding for each pixel
978 //
979 // pad only pixels - which are used (before image cleaning)
980 //
981 Double_t sig = 0;
982 Double_t sigma2 = 0;
983 Double_t diff = 0;
984 Double_t addSig2 = 0;
985 Double_t elNoise2Pix = 0;
986
987
988 for (UInt_t i=0; i<npix; i++)
989 {
990 MCerPhotPix &pix = (*fEvt)[i];
991 if ( !pix.IsPixelUsed() )
992 continue;
993
994 //if ( pix.GetNumPhotons() == 0.0)
995 //{
996 // *fLog << "MPadONOFF::Process(); no.of photons is 0 for used pixel"
997 // << endl;
998 // continue;
999 //}
1000
1001 Int_t j = pix.GetPixId();
1002
1003 Double_t ratioArea = fCam->GetPixRatio(j);
1004
1005 MPedestalPix &ppix = (*fPed)[j];
1006 Double_t oldsigma = ppix.GetMeanRms();
1007 Double_t oldsigma2 = oldsigma*oldsigma;
1008
1009 //---------------------------------
1010 // throw the Sigma for this pixel
1011 //
1012 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
1013
1014 Int_t count;
1015 Bool_t ok;
1016
1017 TH1D *hdiff;
1018 TH1D *hsig;
1019
1020 switch (fPadFlag)
1021 {
1022 case 1 :
1023 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
1024
1025 hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
1026 binPixel, binPixel, "");
1027 if ( hdiff->GetEntries() == 0 )
1028 {
1029 *fLog << "MPadONOFF::Process(); projection for Theta bin "
1030 << binTheta << " and pixel bin " << binPixel
1031 << " has no entries; aborting " << endl;
1032 delete hdiff;
1033
1034 rc = 5;
1035 fErrors[rc]++;
1036 return kCONTINUE;
1037 }
1038
1039 count = 0;
1040 ok = kFALSE;
1041 for (Int_t m=0; m<20; m++)
1042 {
1043 count += 1;
1044 diff = hdiff->GetRandom();
1045 // the following condition ensures that elNoise2Pix > 0.0
1046
1047 if ( (diff + sigmabar2 - oldsigma2/ratioArea
1048 - lambdabar*F2excess) > 0.0 )
1049 {
1050 ok = kTRUE;
1051 break;
1052 }
1053 }
1054 if (!ok)
1055 {
1056
1057 *fLog << "theta, j, count, sigmabar, diff = " << theta << ", "
1058 << j << ", " << count << ", " << sigmabar << ", "
1059 << diff << endl;
1060 diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
1061 }
1062 delete hdiff;
1063 sigma2 = diff + sigmabar2;
1064 break;
1065
1066 case 2 :
1067 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
1068
1069 hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
1070 binPixel, binPixel, "");
1071 if ( hsig->GetEntries() == 0 )
1072 {
1073 *fLog << "MPadONOFF::Process(); projection for Theta bin "
1074 << binTheta << " and pixel bin " << binPixel
1075 << " has no entries; aborting " << endl;
1076 delete hsig;
1077
1078 rc = 6;
1079 fErrors[rc]++;
1080 return kCONTINUE;
1081 }
1082
1083 count = 0;
1084 ok = kFALSE;
1085 for (Int_t m=0; m<20; m++)
1086 {
1087 count += 1;
1088
1089 sig = hsig->GetRandom();
1090 sigma2 = sig*sig/ratioArea;
1091 // the following condition ensures that elNoise2Pix > 0.0
1092
1093 if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 )
1094 {
1095 ok = kTRUE;
1096 break;
1097 }
1098 }
1099 if (!ok)
1100 {
1101
1102 *fLog << "theta, j, count, sigmabar, sig = " << theta << ", "
1103 << j << ", " << count << ", " << sigmabar << ", "
1104 << sig << endl;
1105 sigma2 = lambdabar*F2excess + oldsigma2/ratioArea;
1106 }
1107 delete hsig;
1108 break;
1109 }
1110
1111 //---------------------------------
1112 // get the additional sigma^2 for this pixel (due to the padding)
1113
1114 addSig2 = sigma2*ratioArea - oldsigma2;
1115
1116
1117 //---------------------------------
1118 // get the additional electronic noise for this pixel
1119
1120 elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;
1121
1122
1123 //---------------------------------
1124 // throw actual number of additional NSB photons (NSB)
1125 // and its RMS (sigmaNSB)
1126
1127 Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
1128 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
1129 Double_t sigmaNSB0;
1130
1131 if (arg >= 0)
1132 {
1133 sigmaNSB0 = sqrt( arg );
1134 }
1135 else
1136 {
1137 *fLog << "MPadONOFF::Process(); argument of sqrt < 0 : "
1138 << arg << endl;
1139 sigmaNSB0 = 0.0000001;
1140 }
1141
1142
1143 //---------------------------------
1144 // smear NSB0 according to sigmaNSB0
1145 // and subtract lambdabar because of AC coupling
1146
1147 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
1148
1149 //---------------------------------
1150
1151 // add additional NSB to the number of photons
1152 Double_t oldphotons = pix.GetNumPhotons();
1153 Double_t newphotons = oldphotons + NSB;
1154 pix.SetNumPhotons( newphotons );
1155
1156
1157 fHNSB->Fill( NSB/sqrt(ratioArea) );
1158 fHPhotons->Fill( oldphotons/sqrt(ratioArea), newphotons/sqrt(ratioArea) );
1159
1160
1161 // error: add sigma of padded noise quadratically
1162 Double_t olderror = pix.GetErrorPhot();
1163 Double_t newerror = sqrt( olderror*olderror + addSig2 );
1164 pix.SetErrorPhot( newerror );
1165
1166
1167 Double_t newsigma = sqrt( oldsigma2 + addSig2 );
1168 ppix.SetMeanRms( newsigma );
1169
1170 fHSigmaPedestal->Fill( oldsigma, newsigma );
1171 }
1172 //---------- end of loop over pixels ---------------------------------
1173
1174 // Calculate sigmabar again and crosscheck
1175
1176 //fSigmabar->Calc(*fCam, *fPed, *fEvt);
1177 //*fLog << "after padding : " << endl;
1178 //fSigmabar->Print("");
1179
1180 rc = 0;
1181 fErrors[rc]++;
1182 return kTRUE;
1183 //******************************************************************
1184}
1185
1186// --------------------------------------------------------------------------
1187//
1188//
1189Bool_t MPadONOFF::PostProcess()
1190{
1191 if (GetNumExecutions() != 0)
1192 {
1193
1194 *fLog << inf << endl;
1195 *fLog << GetDescriptor() << " execution statistics:" << endl;
1196 *fLog << dec << setfill(' ');
1197 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
1198 << (int)(fErrors[1]*100/GetNumExecutions())
1199 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
1200
1201 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
1202 << (int)(fErrors[2]*100/GetNumExecutions())
1203 << "%) Evts skipped due to: Zenith angle out of range" << endl;
1204
1205 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
1206 << (int)(fErrors[3]*100/GetNumExecutions())
1207 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
1208
1209 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
1210 << (int)(fErrors[4]*100/GetNumExecutions())
1211 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
1212
1213 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
1214 << (int)(fErrors[5]*100/GetNumExecutions())
1215 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
1216
1217 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
1218 << (int)(fErrors[6]*100/GetNumExecutions())
1219 << "%) Evts skipped due to: No data for generating Sigma" << endl;
1220
1221 *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3)
1222 << (int)(fErrors[7]*100/GetNumExecutions())
1223 << "%) Evts skipped due to: No data for generating Blind pixels" << endl;
1224
1225 *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3)
1226 << (int)(fErrors[8]*100/GetNumExecutions())
1227 << "%) Evts didn't have to be padded" << endl;
1228
1229 *fLog << " " << fErrors[0] << " ("
1230 << (int)(fErrors[0]*100/GetNumExecutions())
1231 << "%) Evts were successfully padded!" << endl;
1232 *fLog << endl;
1233
1234 }
1235
1236 //---------------------------------------------------------------
1237 TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200));
1238 c.Divide(3, 4);
1239 gROOT->SetSelectedPad(NULL);
1240
1241 c.cd(1);
1242 fHNSB->SetDirectory(NULL);
1243 fHNSB->DrawCopy();
1244 fHNSB->SetBit(kCanDelete);
1245
1246 c.cd(2);
1247 fHSigmaPedestal->SetDirectory(NULL);
1248 fHSigmaPedestal->DrawCopy();
1249 fHSigmaPedestal->SetBit(kCanDelete);
1250
1251 c.cd(3);
1252 fHPhotons->SetDirectory(NULL);
1253 fHPhotons->DrawCopy();
1254 fHPhotons->SetBit(kCanDelete);
1255
1256 //--------------------------------------------------------------------
1257
1258
1259 c.cd(4);
1260 fHSigmaTheta->SetDirectory(NULL);
1261 fHSigmaTheta->SetTitle("(Input) 2D : Sigmabar, \\Theta");
1262 fHSigmaTheta->DrawCopy();
1263 fHSigmaTheta->SetBit(kCanDelete);
1264
1265
1266 //--------------------------------------------------------------------
1267
1268
1269 c.cd(7);
1270 fHBlindPixNTheta->SetDirectory(NULL);
1271 fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta");
1272 fHBlindPixNTheta->DrawCopy();
1273 fHBlindPixNTheta->SetBit(kCanDelete);
1274
1275 //--------------------------------------------------------------------
1276
1277
1278 c.cd(10);
1279 fHBlindPixIdTheta->SetDirectory(NULL);
1280 fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta");
1281 fHBlindPixIdTheta->DrawCopy();
1282 fHBlindPixIdTheta->SetBit(kCanDelete);
1283
1284
1285 //--------------------------------------------------------------------
1286 // draw the 3D histogram (input): Theta, pixel, Sigma^2-Sigmabar^2
1287
1288 c.cd(5);
1289 TH2D *l1;
1290 l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
1291 l1->SetDirectory(NULL);
1292 l1->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
1293 l1->SetXTitle("\\Theta [\\circ]");
1294 l1->SetYTitle("Sigma^2-Sigmabar^2");
1295
1296 l1->DrawCopy("box");
1297 l1->SetBit(kCanDelete);;
1298
1299 c.cd(8);
1300 TH2D *l2;
1301 l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
1302 l2->SetDirectory(NULL);
1303 l2->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
1304 l2->SetXTitle("pixel");
1305 l2->SetYTitle("Sigma^2-Sigmabar^2");
1306
1307 l2->DrawCopy("box");
1308 l2->SetBit(kCanDelete);;
1309
1310 //--------------------------------------------------------------------
1311 // draw the 3D histogram (input): Theta, pixel, Sigma
1312
1313 c.cd(6);
1314 TH2D *k1;
1315 k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
1316 k1->SetDirectory(NULL);
1317 k1->SetTitle("(Input) Sigma vs. \\Theta (all pixels)");
1318 k1->SetXTitle("\\Theta [\\circ]");
1319 k1->SetYTitle("Sigma");
1320
1321 k1->DrawCopy("box");
1322 k1->SetBit(kCanDelete);
1323
1324 c.cd(9);
1325 TH2D *k2;
1326 k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
1327 k2->SetDirectory(NULL);
1328 k2->SetTitle("(Input) Sigma vs. pixel number (all \\Theta)");
1329 k2->SetXTitle("pixel");
1330 k2->SetYTitle("Sigma");
1331
1332 k2->DrawCopy("box");
1333 k2->SetBit(kCanDelete);;
1334
1335
1336 //--------------------------------------------------------------------
1337
1338 c.cd(11);
1339 TH2D *m1;
1340 m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
1341 m1->SetDirectory(NULL);
1342 m1->SetTitle("(Input) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
1343 m1->SetXTitle("Sigmabar-old (bin)");
1344 m1->SetYTitle("Sigmabar-new (bin)");
1345
1346 m1->DrawCopy("box");
1347 m1->SetBit(kCanDelete);;
1348
1349 c.cd(12);
1350 TH2D *m2;
1351 m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
1352 m2->SetDirectory(NULL);
1353 m2->SetTitle("(Input) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
1354 m2->SetXTitle("Sigmabar-old (bin)");
1355 m2->SetYTitle("Sigmabar-new (bin)");
1356
1357 m2->DrawCopy("box");
1358 m2->SetBit(kCanDelete);;
1359
1360
1361 return kTRUE;
1362}
1363
1364// --------------------------------------------------------------------------
1365
1366
1367
Note: See TracBrowser for help on using the repository browser.