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

Last change on this file since 2756 was 2746, checked in by wittek, 21 years ago
*** empty log message ***
File size: 53.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, 06/2003 <mailto:wittek@mppmu.mpg.de>
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// 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 "MPadONOFF.h"
63
64#include <math.h>
65#include <stdio.h>
66
67#include <TH1.h>
68#include <TH2.h>
69#include <TH3.h>
70#include <TRandom.h>
71#include <TCanvas.h>
72#include <TFile.h>
73
74#include "MBinning.h"
75#include "MSigmabar.h"
76#include "MMcEvt.hxx"
77#include "MLog.h"
78#include "MLogManip.h"
79#include "MParList.h"
80#include "MGeomCam.h"
81
82#include "MCerPhotPix.h"
83#include "MCerPhotEvt.h"
84
85#include "MPedestalCam.h"
86#include "MPedestalPix.h"
87#include "MBlindPixels.h"
88
89#include "MRead.h"
90#include "MFilterList.h"
91#include "MTaskList.h"
92#include "MBlindPixelCalc.h"
93#include "MHBlindPixels.h"
94#include "MFillH.h"
95#include "MHSigmaTheta.h"
96#include "MEvtLoop.h"
97
98ClassImp(MPadONOFF);
99
100using namespace std;
101
102// --------------------------------------------------------------------------
103//
104// Default constructor.
105//
106MPadONOFF::MPadONOFF(const char *name, const char *title)
107{
108 fName = name ? name : "MPadONOFF";
109 fTitle = title ? title : "Task for the ON-OFF padding";
110
111 fPadFlag = 1;
112 *fLog << all << "MadONOFF : fPadFlag = " << fPadFlag << endl;
113
114 fType = "";
115
116 fHSigmaTheta = NULL;
117 fHSigmaThetaON = NULL;
118 fHSigmaThetaOFF = NULL;
119
120 fHSigmaPixTheta = NULL;
121 fHSigmaPixThetaON = NULL;
122 fHSigmaPixThetaOFF = NULL;
123
124 fHDiffPixTheta = NULL;
125 fHDiffPixThetaON = NULL;
126 fHDiffPixThetaOFF = NULL;
127
128 fHgON = NULL;
129 fHgOFF = NULL;
130
131 fHBlindPixIdTheta = NULL;
132 fHBlindPixNTheta = NULL;
133
134 fHSigmaPedestal = NULL;
135 fHPhotons = NULL;
136 fHNSB = NULL;
137}
138
139// --------------------------------------------------------------------------
140//
141// Destructor.
142//
143MPadONOFF::~MPadONOFF()
144{
145 if (fHBlindPixIdTheta != NULL) delete fHBlindPixIdTheta;
146 if (fHBlindPixNTheta != NULL) delete fHBlindPixNTheta;
147
148 if (fHSigmaTheta != NULL) delete fHSigmaTheta;
149 if (fHSigmaPixTheta != NULL) delete fHSigmaPixTheta;
150 if (fHDiffPixTheta != NULL) delete fHDiffPixTheta;
151
152 if (fHSigmaPedestal != NULL) delete fHSigmaPedestal;
153 if (fHPhotons != NULL) delete fHPhotons;
154 if (fHNSB != NULL) delete fHNSB;
155
156 if (fInfile != NULL) delete fInfile;
157}
158
159// --------------------------------------------------------------------------
160//
161// Merge the distributions of ON and OFF data
162//
163// fHSigmaTheta 2D-histogram (Theta, sigmabar)
164// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma)
165// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2)
166// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id)
167// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels )
168//
169// and define the target distributions for the padding
170//
171Bool_t MPadONOFF::MergeHistograms(TH2D *sigthon, TH2D *sigthoff,
172 TH3D *sigpixthon, TH3D *sigpixthoff,
173 TH3D *diffpixthon, TH3D *diffpixthoff,
174 TH2D *blindidthon, TH2D *blindidthoff,
175 TH2D *blindnthon, TH2D *blindnthoff)
176{
177 *fLog << all << "----------------------------------------------------------------------------------" << endl;
178 *fLog << all << "MPadONOFF::MergeHistograms(); Merge the ON and OFF histograms to obtain the target distributions for the padding"
179 << endl;
180
181 //-------------------------------------------------------------
182 // merge the distributions of ON and OFF events
183 // to obtain the target distribution fHSigmaTheta
184 //
185
186 Double_t eps = 1.e-10;
187
188 fHSigmaTheta = new TH2D( (TH2D&)*sigthon );
189 fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");
190
191 // get binning for fHgON and fHgOFF from sigthon
192 // binning in Theta
193 TAxis *ax = sigthon->GetXaxis();
194 Int_t nbinstheta = ax->GetNbins();
195 TArrayD edgesx;
196 edgesx.Set(nbinstheta+1);
197 for (Int_t i=0; i<=nbinstheta; i++)
198 {
199 edgesx[i] = ax->GetBinLowEdge(i+1);
200 *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i]
201 << endl;
202 }
203 MBinning binth;
204 binth.SetEdges(edgesx);
205
206 // binning in sigmabar
207 TAxis *ay = sigthon->GetYaxis();
208 Int_t nbinssig = ay->GetNbins();
209 TArrayD edgesy;
210 edgesy.Set(nbinssig+1);
211 for (Int_t i=0; i<=nbinssig; i++)
212 {
213 edgesy[i] = ay->GetBinLowEdge(i+1);
214 *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i]
215 << endl;
216 }
217 MBinning binsg;
218 binsg.SetEdges(edgesy);
219
220
221 fHgON = new TH3D;
222 MH::SetBinning(fHgON, &binth, &binsg, &binsg);
223 fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");
224
225 fHgOFF = new TH3D;
226 MH::SetBinning(fHgOFF, &binth, &binsg, &binsg);
227 fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF");
228
229 //............ loop over Theta bins ...........................
230
231 // hista is the normalized 1D histogram of sigmabar for ON data
232 // histb is the normalized 1D histogram of sigmabar for OFF data
233 // histc is the difference ON-OFF
234
235 // at the beginning, histap is a copy of hista
236 // at the end, it will be the 1D histogram for ON data after padding
237
238 // at the beginning, histbp is a copy of histb
239 // at the end, it will be the 1D histogram for OFF data after padding
240
241 // at the beginning, histcp is a copy of histc
242 // at the end, it should be zero
243
244 *fLog << all << "MPadONOFF::MergeHistograms(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
245 for (Int_t l=1; l<=nbinstheta; l++)
246 {
247 TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, "");
248 Stat_t suma = hista->Integral();
249 hista->Scale(1./suma);
250
251 TH1D *histap = new TH1D( (const TH1D&)*hista );
252
253 TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, "");
254 Stat_t sumb = histb->Integral();
255 histb->Scale(1./sumb);
256
257 TH1D *histbp = new TH1D( (const TH1D&)*histb );
258
259 TH1D *histc = new TH1D( (const TH1D&)*hista );
260 histc->Add(histb, -1.0);
261
262 TH1D *histcp = new TH1D( (const TH1D&)*histc );
263
264
265 // calculate matrix g
266
267 // fHgON[k][j] tells which fraction of ON events should be padded
268 // from sigma_k to sigma_j
269
270
271 // fHgOFF[k][j] tells which fraction of OFF events should be padded
272 // from sigma_k to sigma_j
273
274 //-------- start j loop ------------------------------------------------
275 // loop over bins in histc, starting from the end
276 Double_t v, s, w, t, x, u, a, b, c, arg;
277
278 for (Int_t j=nbinssig; j >= 1; j--)
279 {
280 v = histcp->GetBinContent(j);
281 if ( fabs(v) < eps ) continue;
282 if (v >= 0.0)
283 s = 1.0;
284 else
285 s = -1.0;
286
287 //................ start k loop ......................................
288 // look for a bin k which may compensate the content of bin j
289 for (Int_t k=j-1; k >= 1; k--)
290 {
291 w = histcp->GetBinContent(k);
292 if ( fabs(w) < eps ) continue;
293 if (w >= 0.0)
294 t = 1.0;
295 else
296 t = -1.0;
297
298 if (s==t) continue;
299
300 x = v + w;
301 if (x >= 0.0)
302 u = 1.0;
303 else
304 u = -1.0;
305
306 if (u == s)
307 {
308 arg = -w;
309
310 if (arg <=0.0)
311 {
312 fHgON->SetBinContent (l, k, j, -arg);
313 fHgOFF->SetBinContent(l, k, j, 0.0);
314 }
315 else
316 {
317 fHgON->SetBinContent (l, k, j, 0.0);
318 fHgOFF->SetBinContent(l, k, j, arg);
319 }
320
321
322 *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j
323 << ", " << arg << endl;
324
325 // g(k-1, j-1) = arg;
326 //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", "
327 // << arg << endl;
328
329 //......................................
330 // this is for checking the procedure
331 if (arg < 0.0)
332 {
333 a = histap->GetBinContent(k);
334 histap->SetBinContent(k, a+arg);
335 a = histap->GetBinContent(j);
336 histap->SetBinContent(j, a-arg);
337 }
338 else
339 {
340 b = histbp->GetBinContent(k);
341 histbp->SetBinContent(k, b-arg);
342 b = histbp->GetBinContent(j);
343 histbp->SetBinContent(j, b+arg);
344 }
345 //......................................
346
347 histcp->SetBinContent(k, 0.0);
348 histcp->SetBinContent(j, x);
349
350 //......................................
351 // redefine v
352 v = histcp->GetBinContent(j);
353 if ( fabs(v) < eps ) break;
354 if (v >= 0.0)
355 s = 1.0;
356 else
357 s = -1.0;
358 //......................................
359
360 continue;
361 }
362
363 arg = v;
364
365 if (arg <=0.0)
366 {
367 fHgON->SetBinContent (l, k, j, -arg);
368 fHgOFF->SetBinContent(l, k, j, 0.0);
369 }
370 else
371 {
372 fHgON->SetBinContent (l, k, j, 0.0);
373 fHgOFF->SetBinContent(l, k, j, arg);
374 }
375
376 *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j
377 << ", " << arg << endl;
378
379 //g(k-1, j-1) = arg;
380 //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", "
381 // << arg << endl;
382
383 //......................................
384 // this is for checking the procedure
385 if (arg < 0.0)
386 {
387 a = histap->GetBinContent(k);
388 histap->SetBinContent(k, a+arg);
389 a = histap->GetBinContent(j);
390 histap->SetBinContent(j, a-arg);
391 }
392 else
393 {
394 b = histbp->GetBinContent(k);
395
396 histbp->SetBinContent(k, b-arg);
397 b = histbp->GetBinContent(j);
398 histbp->SetBinContent(j, b+arg);
399 }
400 //......................................
401
402 histcp->SetBinContent(k, x);
403 histcp->SetBinContent(j, 0.0);
404
405 break;
406 }
407 //................ end k loop ......................................
408 }
409 //-------- end j loop ------------------------------------------------
410
411 // check results for this Theta bin
412 for (Int_t j=1; j<=nbinssig; j++)
413 {
414 a = histap->GetBinContent(j);
415 b = histbp->GetBinContent(j);
416 c = histcp->GetBinContent(j);
417
418 if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps )
419 *fLog << err << "MPadONOFF::MergeHistograms(); inconsistency in results; a, b, c = "
420 << a << ", " << b << ", " << c << endl;
421 }
422
423
424 // fill target distribution SigmaTheta
425 // for this Theta bin
426 //
427 for (Int_t j=1; j<=nbinssig; j++)
428 {
429 a = histap->GetBinContent(j);
430 fHSigmaTheta->SetBinContent(l, j, a);
431 }
432
433 delete hista;
434 delete histb;
435 delete histc;
436
437 delete histap;
438 delete histbp;
439 delete histcp;
440 }
441 //............ end of loop over Theta bins ....................
442
443
444 // target distributions for MC
445 // SigmaPixTheta and DiffPixTheta
446 // BlindPixIdTheta and BlindPixNTheta
447 // are calculated as averages of the ON and OFF distributions
448
449 // Note : in all histograms the pedestal RMS is assumed to be given in units of "number of photons"
450
451 fHSigmaThetaON = sigthon;
452 fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (target)");
453
454 fHSigmaThetaOFF = sigthoff;
455 fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (target)");
456
457
458 fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon );
459 fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
460
461 fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon );
462 fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
463
464 fHSigmaPixTheta = new TH3D( (TH3D&)*sigpixthon );
465 fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");
466
467 fHSigmaPixThetaON = sigpixthon;
468 fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (target)");
469
470 fHSigmaPixThetaOFF = sigpixthoff;
471 fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (target)");
472
473 fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon );
474 fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");
475
476 fHDiffPixThetaON = diffpixthon;
477 fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (target)");
478
479 fHDiffPixThetaOFF = diffpixthoff;
480 fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (target)");
481
482
483 for (Int_t j=1; j<=nbinstheta; j++)
484 {
485 // fraction of ON/OFF events to be padded to a sigmabar
486 // of the OFF/ON events : fracON, fracOFF
487
488 Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, "");
489 Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, "");
490
491 TAxis *ax;
492 TAxis *ay;
493 TAxis *az;
494
495 Double_t entON;
496 Double_t entOFF;
497
498 Double_t normON;
499 Double_t normOFF;
500
501 // set ranges for 2D-projections of 3D-histograms
502 ax = sigpixthon->GetXaxis();
503 ax->SetRange(j, j);
504 ax = sigpixthoff->GetXaxis();
505 ax->SetRange(j, j);
506
507 ax = diffpixthon->GetXaxis();
508 ax->SetRange(j, j);
509 ax = diffpixthoff->GetXaxis();
510 ax->SetRange(j, j);
511
512 TH2D *hist;
513 TH2D *histOFF;
514
515 TH1D *hist1;
516 TH1D *hist1OFF;
517
518 // weights for ON and OFF distrubtions when
519 // calculating the weighted average
520 Double_t facON = 0.5 * (1. - fracON + fracOFF);
521 Double_t facOFF = 0.5 * (1. + fracON - fracOFF);
522
523 *fLog << all << fracON << ", " << fracOFF << ", "
524 << facON << ", " << facOFF << endl;
525 //-------------------------------------------
526 ay = blindnthon->GetYaxis();
527 Int_t nbinsn = ay->GetNbins();
528
529 hist1 = (TH1D*)blindnthon->ProjectionY ("", j, j, "");
530 hist1->SetName("dummy");
531 hist1OFF = (TH1D*)blindnthoff->ProjectionY("", j, j, "");
532
533 // number of events in this theta bin
534 entON = hist1->Integral();
535 entOFF = hist1OFF->Integral();
536
537 *fLog << all << "BlindPixNTheta : j, entON, entOFF = " << j << ", "
538 << entON << ", " << entOFF << endl;
539
540 normON = entON<=0.0 ? 0.0 : 1.0/entON;
541 normOFF = entOFF<=0.0 ? 0.0 : 1.0/entOFF;
542
543 hist1->Scale(normON);
544 hist1OFF->Scale(normOFF);
545
546 // weighted average of ON and OFF distributions
547 hist1->Scale(facON);
548 hist1->Add(hist1OFF, facOFF);
549
550 for (Int_t k=1; k<=nbinsn; k++)
551 {
552 Double_t cont = hist1->GetBinContent(k);
553 fHBlindPixNTheta->SetBinContent(j, k, cont);
554 }
555
556 delete hist1;
557 delete hist1OFF;
558
559 //-------------------------------------------
560 ay = blindidthon->GetYaxis();
561 Int_t nbinsid = ay->GetNbins();
562
563 hist1 = (TH1D*)blindidthon->ProjectionY ("", j, j, "");
564 hist1->SetName("dummy");
565 hist1OFF = (TH1D*)blindidthoff->ProjectionY("", j, j, "");
566
567 hist1->Scale(normON);
568 hist1OFF->Scale(normOFF);
569
570 // weighted average of ON and OFF distributions
571 hist1->Scale(facON);
572 hist1->Add(hist1OFF, facOFF);
573
574 for (Int_t k=1; k<=nbinsid; k++)
575 {
576 Double_t cont = hist1->GetBinContent(k);
577 fHBlindPixIdTheta->SetBinContent(j, k, cont);
578 }
579
580 delete hist1;
581 delete hist1OFF;
582
583 //-------------------------------------------
584 ay = sigpixthon->GetYaxis();
585 Int_t nbinspix = ay->GetNbins();
586
587 az = sigpixthon->GetZaxis();
588 Int_t nbinssigma = az->GetNbins();
589
590 hist = (TH2D*)sigpixthon->Project3D("zy");
591 hist->SetName("dummy");
592 histOFF = (TH2D*)sigpixthoff->Project3D("zy");
593
594 hist->Scale(normON);
595 histOFF->Scale(normOFF);
596
597 // weighted average of ON and OFF distributions
598
599 hist->Scale(facON);
600 hist->Add(histOFF, facOFF);
601
602 for (Int_t k=1; k<=nbinspix; k++)
603 for (Int_t l=1; l<=nbinssigma; l++)
604 {
605 Double_t cont = hist->GetBinContent(k,l);
606 fHSigmaPixTheta->SetBinContent(j, k, l, cont);
607 }
608
609 delete hist;
610 delete histOFF;
611
612 //-------------------------------------------
613 ay = diffpixthon->GetYaxis();
614 Int_t nbinspixel = ay->GetNbins();
615
616 az = diffpixthon->GetZaxis();
617 Int_t nbinsdiff = az->GetNbins();
618
619 hist = (TH2D*)diffpixthon->Project3D("zy");
620 hist->SetName("dummy");
621 histOFF = (TH2D*)diffpixthoff->Project3D("zy");
622
623 hist->Scale(normON);
624 histOFF->Scale(normOFF);
625
626 // weighted average of ON and OFF distributions
627 hist->Scale(facON);
628 hist->Add(histOFF, facOFF);
629
630 for (Int_t k=1; k<=nbinspixel; k++)
631 for (Int_t l=1; l<=nbinsdiff; l++)
632 {
633 Double_t cont = hist->GetBinContent(k,l);
634 fHDiffPixTheta->SetBinContent(j, k, l, cont);
635 }
636
637 delete hist;
638 delete histOFF;
639
640 //-------------------------------------------
641 }
642
643
644 *fLog << all
645 << "MPadONOFF::MergeHistograms(); The target distributions for the padding have been created"
646 << endl;
647 *fLog << all << "----------------------------------------------------------"
648 << endl;
649 //--------------------------------------------
650
651 fHSigmaTheta->SetDirectory(NULL);
652 fHSigmaThetaON->SetDirectory(NULL);
653 fHSigmaThetaOFF->SetDirectory(NULL);
654
655 fHSigmaPixTheta->SetDirectory(NULL);
656 fHSigmaPixThetaON->SetDirectory(NULL);
657 fHSigmaPixThetaOFF->SetDirectory(NULL);
658
659 fHDiffPixTheta->SetDirectory(NULL);
660 fHDiffPixThetaON->SetDirectory(NULL);
661 fHDiffPixThetaOFF->SetDirectory(NULL);
662
663 fHBlindPixIdTheta->SetDirectory(NULL);
664 fHBlindPixNTheta->SetDirectory(NULL);
665
666
667 fHgON->SetDirectory(NULL);
668 fHgOFF->SetDirectory(NULL);
669
670
671 return kTRUE;
672}
673
674
675// --------------------------------------------------------------------------
676//
677// Read target distributions from a file
678//
679//
680Bool_t MPadONOFF::ReadTargetDist(const char* namefilein)
681{
682 *fLog << all << "MPadONOFF : Read padding histograms from file "
683 << namefilein << endl;
684
685 fInfile = new TFile(namefilein);
686 fInfile->ls();
687
688 //------------------------------------
689
690 fHSigmaTheta =
691 (TH2D*) gROOT->FindObject("2D-ThetaSigmabar");
692 if (!fHSigmaTheta)
693 {
694 *fLog << all
695 << "MPadONOFF : Object '2D-ThetaSigmabar' not found on root file"
696 << endl;
697 return kFALSE;
698 }
699 *fLog << all
700 << "MPadONOFF : Object '2D-ThetaSigmabar' was read in" << endl;
701
702 fHSigmaThetaON =
703 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON");
704 if (!fHSigmaThetaON)
705 {
706 *fLog << all
707 << "MPadONOFF : Object '2D-ThetaSigmabarON' not found on root file"
708 << endl;
709 return kFALSE;
710 }
711 *fLog << all
712 << "MPadONOFF : Object '2D-ThetaSigmabarON' was read in" << endl;
713
714 fHSigmaThetaOFF =
715 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF");
716 if (!fHSigmaThetaOFF)
717 {
718 *fLog << all
719 << "MPadONOFF : Object '2D-ThetaSigmabarOFF' not found on root file"
720 << endl;
721 return kFALSE;
722 }
723 *fLog << all
724 << "MPadONOFF :Object '2D-ThetaSigmabarOFF' was read in" << endl;
725
726 fHSigmaPixTheta =
727 (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
728 if (!fHSigmaPixTheta)
729 {
730 *fLog << all
731 << "MPadONOFF : Object '3D-ThetaPixSigma' not found on root file"
732 << endl;
733 return kFALSE;
734 }
735 *fLog << all
736 << "MPadONOFF : Object '3D-ThetaPixSigma' was read in" << endl;
737
738 fHSigmaPixThetaON =
739 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON");
740 if (!fHSigmaPixThetaON)
741 {
742 *fLog << all
743 << "MPadONOFF : Object '3D-ThetaPixSigmaON' not found on root file"
744 << endl;
745 return kFALSE;
746 }
747 *fLog << all
748 << "MPadONOFF : Object '3D-ThetaPixSigmaON' was read in" << endl;
749
750 fHSigmaPixThetaOFF =
751 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF");
752 if (!fHSigmaPixThetaOFF)
753 {
754 *fLog << all
755 << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' not found on root file"
756 << endl;
757 return kFALSE;
758 }
759 *fLog << all
760 << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' was read in" << endl;
761
762 fHDiffPixTheta =
763 (TH3D*) gROOT->FindObject("3D-ThetaPixDiff");
764 if (!fHDiffPixTheta)
765 {
766 *fLog << all
767 << "MPadONOFF : Object '3D-ThetaPixDiff' not found on root file"
768 << endl;
769 return kFALSE;
770 }
771 *fLog << all
772 << "MPadONOFF : Object '3D-ThetaPixDiff' was read in" << endl;
773
774 fHDiffPixThetaON =
775 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON");
776 if (!fHDiffPixThetaON)
777 {
778 *fLog << all
779 << "MPadONOFF : Object '3D-ThetaPixDiffON' not found on root file"
780 << endl;
781 return kFALSE;
782 }
783 *fLog << all
784 << "MPadONOFF : Object '3D-ThetaPixDiffON' was read in" << endl;
785
786 fHDiffPixThetaOFF =
787 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF");
788 if (!fHDiffPixThetaOFF)
789 {
790 *fLog << all
791 << "MPadONOFF : Object '3D-ThetaPixDiffOFF' not found on root file"
792 << endl;
793 return kFALSE;
794 }
795 *fLog << all
796 << "MPadONOFF : Object '3D-ThetaPixDiffOFF' was read in" << endl;
797
798 fHgON =
799 (TH3D*) gROOT->FindObject("3D-PaddingMatrixON");
800 if (!fHgON)
801 {
802 *fLog << all
803 << "MPadONOFF : Object '3D-PaddingMatrixON' not found on root file"
804 << endl;
805 return kFALSE;
806 }
807 *fLog << all
808 << "MPadONOFF : Object '3D-PaddingMatrixON' was read in" << endl;
809
810 fHgOFF =
811 (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF");
812 if (!fHgOFF)
813 {
814 *fLog << all
815 << "MPadONOFF : Object '3D-PaddingMatrixOFF' not found on root file"
816 << endl;
817 return kFALSE;
818 }
819 *fLog << all
820 << "MPadONOFF : Object '3D-PaddingMatrixOFF' was read in" << endl;
821
822
823 fHBlindPixIdTheta =
824 (TH2D*) gROOT->FindObject("2D-ThetaBlindId");
825 if (!fHBlindPixIdTheta)
826 {
827 *fLog << all
828 << "MPadONOFF : Object '2D-ThetaBlindId' not found on root file"
829 << endl;
830 return kFALSE;
831 }
832 *fLog << all
833 << "MPadONOFF : Object '2D-ThetaBlindId' was read in" << endl;
834
835 fHBlindPixNTheta =
836 (TH2D*) gROOT->FindObject("2D-ThetaBlindN");
837 if (!fHBlindPixNTheta)
838 {
839 *fLog << all
840 << "MPadONOFF : Object '2D-ThetaBlindN' not found on root file"
841 << endl;
842 return kFALSE;
843 }
844 *fLog << all
845 << "MPadONOFF : Object '2D-ThetaBlindN' was read in" << endl;
846
847
848 //------------------------------------
849
850 return kTRUE;
851}
852
853
854// --------------------------------------------------------------------------
855//
856// Write target distributions onto a file
857//
858//
859Bool_t MPadONOFF::WriteTargetDist(const char* namefileout)
860{
861 *fLog << all << "MPadONOFF : Write padding histograms onto file "
862 << namefileout << endl;
863
864 TFile outfile(namefileout, "RECREATE");
865
866 fHBlindPixNTheta->Write();
867 fHBlindPixIdTheta->Write();
868
869 fHSigmaTheta->Write();
870 fHSigmaThetaON->Write();
871 fHSigmaThetaOFF->Write();
872
873 fHSigmaPixTheta->Write();
874 fHSigmaPixThetaON->Write();
875 fHSigmaPixThetaOFF->Write();
876
877 fHDiffPixTheta->Write();
878 fHDiffPixThetaON->Write();
879 fHDiffPixThetaOFF->Write();
880
881 fHgON->Write();
882 fHgOFF->Write();
883
884 *fLog << all
885 << "MPadONOFF::WriteTargetDist(); target histograms written onto file "
886 << namefileout << endl;
887
888 return kTRUE;
889}
890
891// --------------------------------------------------------------------------
892//
893// Set type of data to be padded
894//
895// this is not necessary if the type of the data can be recognized
896// directly from the input
897//
898//
899void MPadONOFF::SetDataType(const char* type)
900{
901 fType = type;
902 *fLog << all << "MPadONOFF::SetDataType(); type of data to be padded : "
903 << fType << endl;
904
905 return;
906}
907
908
909// --------------------------------------------------------------------------
910//
911// Set the option for the padding
912//
913// There are 2 options for the padding :
914//
915// 1) fPadFlag = 1 :
916// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
917// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
918// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
919// (fHDiffPixTheta).
920//
921// This is the preferred option as it takes into account the
922// correlations between the Sigma of a pixel and Sigmabar.
923//
924// NOT YET TESTED :
925// 2) fPadFlag = 2 :
926// Generate a pedestal sigma for each pixel using the 3D-histogram
927// Theta, pixel no., Sigma (fHSigmaPixTheta).
928//
929void MPadONOFF::SetPadFlag(Int_t padflag)
930{
931 fPadFlag = padflag;
932 *fLog << all << "MPadONOFF::SetPadFlag(); choose option " << fPadFlag
933 << endl;
934}
935
936// --------------------------------------------------------------------------
937//
938// Set the pointers and prepare the histograms
939//
940Int_t MPadONOFF::PreProcess(MParList *pList)
941{
942 if ( !fHSigmaTheta || !fHSigmaThetaON || !fHSigmaThetaOFF ||
943 !fHSigmaPixTheta || !fHSigmaPixThetaON || !fHSigmaPixThetaOFF ||
944 !fHDiffPixTheta || !fHDiffPixThetaON || !fHDiffPixThetaOFF ||
945 !fHBlindPixIdTheta || !fHBlindPixNTheta ||
946 !fHgON || !fHgOFF)
947 {
948 *fLog << err
949 << "MPadONOFF : At least one of the histograms needed for the padding is not defined ... aborting."
950 << endl;
951 return kFALSE;
952 }
953
954 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
955 if (!fMcEvt)
956 {
957 *fLog << err << dbginf << "MPadONOFF : MMcEvt not found... aborting."
958 << endl;
959 return kFALSE;
960 }
961
962 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
963 if (!fPed)
964 {
965 *fLog << err << "MPadONOFF : MPedestalCam not found... aborting."
966 << endl;
967 return kFALSE;
968 }
969
970 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
971 if (!fCam)
972 {
973 *fLog << err
974 << "MPadONOFF : MGeomCam not found (no geometry information available)... aborting."
975 << endl;
976 return kFALSE;
977 }
978
979 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
980 if (!fEvt)
981 {
982 *fLog << err << "MPadONOFF : MCerPhotEvt not found... aborting."
983 << endl;
984 return kFALSE;
985 }
986
987 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
988 if (!fSigmabar)
989 {
990 *fLog << err << "MPadONOFF : MSigmabar not found... aborting." << endl;
991 return kFALSE;
992 }
993
994 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
995 if (!fBlinds)
996 {
997 *fLog << err << "MPadONOFF : MBlindPixels not found... aborting."
998 << endl;
999 return kFALSE;
1000 }
1001
1002 if (fType !="ON" && fType !="OFF" && fType !="MC")
1003 {
1004 *fLog << err
1005 << "MPadONOFF : Type of data to be padded not defined... aborting."
1006 << endl;
1007 return kFALSE;
1008 }
1009
1010
1011 //--------------------------------------------------------------------
1012 // histograms for checking the padding
1013 //
1014 // plot pedestal sigmas
1015 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
1016 100, 0.0, 5.0, 100, 0.0, 5.0);
1017 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
1018 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
1019
1020 // plot no.of photons (before vs. after padding)
1021 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
1022 100, -10.0, 90.0, 100, -10, 90);
1023 fHPhotons->SetXTitle("no.of photons before padding");
1024 fHPhotons->SetYTitle("no.of photons after padding");
1025
1026 // plot of added NSB
1027 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
1028 fHNSB->SetXTitle("no.of added NSB photons");
1029 fHNSB->SetYTitle("no.of pixels");
1030
1031
1032 //--------------------------------------------------------------------
1033
1034 fIter = 20;
1035
1036 memset(fErrors, 0, sizeof(fErrors));
1037
1038 return kTRUE;
1039}
1040
1041// --------------------------------------------------------------------------
1042//
1043// Do the Padding
1044// idealy the events to be padded should have been generated without NSB
1045//
1046Int_t MPadONOFF::Process()
1047{
1048 *fLog << all << "Entry MPadONOFF::Process();" << endl;
1049
1050 // this is the conversion factor from the number of photons
1051 // to the number of photo electrons
1052 // later to be taken from MCalibration
1053 // fPEperPhoton = PW * LG * QE * 1D
1054 Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9;
1055
1056 //------------------------------------------------
1057 // add pixels to MCerPhotEvt which are not yet in;
1058 // set their number of photons equal to zero
1059
1060 UInt_t imaxnumpix = fCam->GetNumPixels();
1061
1062 for (UInt_t i=0; i<imaxnumpix; i++)
1063 {
1064 Bool_t alreadythere = kFALSE;
1065 UInt_t npix = fEvt->GetNumPixels();
1066 for (UInt_t j=0; j<npix; j++)
1067 {
1068 MCerPhotPix &pix = (*fEvt)[j];
1069 UInt_t id = pix.GetPixId();
1070 if (i==id)
1071 {
1072 alreadythere = kTRUE;
1073 break;
1074 }
1075 }
1076 if (alreadythere)
1077 continue;
1078
1079 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetPedestalRms());
1080 }
1081
1082
1083
1084 //-----------------------------------------
1085 Int_t rc=0;
1086 Int_t rw=0;
1087
1088 const UInt_t npix = fEvt->GetNumPixels();
1089
1090 //fSigmabar->Calc(*fCam, *fPed, *fEvt);
1091 // *fLog << all << "before padding : " << endl;
1092 //fSigmabar->Print("");
1093
1094
1095 //$$$$$$$$$$$$$$$$$$$$$$$$$$
1096 // to simulate the situation that before the padding the NSB and
1097 // electronic noise are zero : set Sigma = 0 for all pixels
1098 //for (UInt_t i=0; i<npix; i++)
1099 //{
1100 // MCerPhotPix &pix = fEvt->operator[](i);
1101 // Int_t j = pix.GetPixId();
1102
1103 // MPedestalPix &ppix = fPed->operator[](j);
1104 // ppix.SetMeanRms(0.0);
1105 //}
1106 //$$$$$$$$$$$$$$$$$$$$$$$$$$
1107
1108 //-------------------------------------------
1109 // Calculate average sigma of the event
1110 //
1111 Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt);
1112 *fLog << all << "MPadONOFF::Process(); before padding : " << endl;
1113 fSigmabar->Print("");
1114 Double_t sigbarold = sigbarold_phot * fPEperPhoton;
1115
1116
1117 Double_t sigbarold2 = sigbarold*sigbarold;
1118
1119
1120 // for MC data : expect sigmabar to be zero before the padding
1121 if (fType == "MC")
1122 {
1123 if (sigbarold_phot > 0)
1124 {
1125 *fLog << err
1126 << "MPadONOFF::Process(); sigmabar of event to be padded is > 0 : "
1127 << sigbarold_phot << ". Stop event loop " << endl;
1128 // input data should have sigmabar = 0; stop event loop
1129
1130 rc = 1;
1131 fErrors[rc]++;
1132 return kERROR;
1133 }
1134 }
1135
1136 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
1137 *fLog << all << "theta = " << theta << endl;
1138
1139 Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
1140 if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() )
1141 {
1142 *fLog << warn
1143 << "MPadONOFF::Process(); binNumber out of range : theta, binTheta = "
1144 << theta << ", " << binTheta << "; Skip event " << endl;
1145 // event cannot be padded; skip event
1146
1147 rc = 2;
1148 fErrors[rc]++;
1149 return kCONTINUE;
1150 }
1151
1152 //-------------------------------------------
1153 // get number of events in this theta bin
1154 // and number of events in this sigbarold_phot bin
1155
1156 Double_t nTheta;
1157 Double_t nSigma;
1158 if (fType == "ON")
1159 {
1160 TH1D *hn;
1161
1162 hn = fHSigmaThetaON->ProjectionY("", binTheta, binTheta, "");
1163 nTheta = hn->Integral();
1164 Int_t binSigma = hn->FindBin(sigbarold_phot);
1165 nSigma = hn->GetBinContent(binSigma);
1166
1167 *fLog << all
1168 << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1169 << theta << ", " << sigbarold_phot << ", " << binTheta << ", "
1170 << binSigma << ", " << nTheta << ", " << nSigma << endl;
1171
1172 delete hn;
1173 }
1174
1175 else if (fType == "OFF")
1176 {
1177 TH1D *hn;
1178
1179 hn = fHSigmaThetaOFF->ProjectionY("", binTheta, binTheta, "");
1180 nTheta = hn->Integral();
1181 Int_t binSigma = hn->FindBin(sigbarold_phot);
1182 nSigma = hn->GetBinContent(binSigma);
1183
1184 *fLog << all
1185 << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1186 << theta << ", " << sigbarold_phot << ", " << binTheta << ", "
1187 << binSigma << ", " << nTheta << ", " << nSigma << endl;
1188
1189 delete hn;
1190 }
1191
1192 else
1193 {
1194 nTheta = 0.0;
1195 nSigma = 0.0;
1196 }
1197
1198 //-------------------------------------------
1199 // for the current theta,
1200 // generate blind pixels according to the histograms
1201 // fHBlindPixNTheta and fHBlindPixIDTheta
1202 // do this only for MC data
1203
1204
1205 if (fType == "MC")
1206 {
1207
1208 // numBlind is the number of blind pixels in this event
1209 TH1D *nblind;
1210 UInt_t numBlind;
1211
1212 nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, "");
1213 if ( nblind->Integral() == 0.0 )
1214 {
1215 *fLog << warn << "MPadONOFF::Process(); projection for Theta bin "
1216 << binTheta << " has no entries; Skip event " << endl;
1217 // event cannot be padded; skip event
1218 delete nblind;
1219
1220 rc = 7;
1221 fErrors[rc]++;
1222 return kCONTINUE;
1223 }
1224 else
1225 {
1226 numBlind = (Int_t) (nblind->GetRandom()+0.5);
1227 }
1228 delete nblind;
1229
1230 //warn Code commented out due no compilation errors on Alpha OSF (tgb)
1231
1232 // throw the Id of numBlind different pixels in this event
1233 TH1D *hblind;
1234 UInt_t idBlind;
1235 UInt_t listId[npix];
1236 UInt_t nlist = 0;
1237 Bool_t equal;
1238
1239 hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, "");
1240 if ( hblind->Integral() > 0.0 )
1241 for (UInt_t i=0; i<numBlind; i++)
1242 {
1243 while(1)
1244 {
1245 idBlind = (Int_t) (hblind->GetRandom()+0.5);
1246 equal = kFALSE;
1247 for (UInt_t j=0; j<nlist; j++)
1248 if (idBlind == listId[j])
1249 {
1250 equal = kTRUE;
1251 break;
1252 }
1253 if (!equal) break;
1254 }
1255 listId[nlist] = idBlind;
1256 nlist++;
1257
1258 fBlinds->SetPixelBlind(idBlind);
1259 *fLog << all << "idBlind = " << idBlind << endl;
1260 }
1261 fBlinds->SetReadyToSave();
1262
1263 delete hblind;
1264
1265 }
1266
1267 //******************************************************************
1268 // has the event to be padded ?
1269 // if yes : to which sigmabar should it be padded ?
1270 //
1271
1272 Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot);
1273 *fLog << all << "binSig, sigbarold_phot = " << binSig << ", "
1274 << sigbarold_phot << endl;
1275
1276 Double_t prob;
1277 TH1D *hpad = NULL;
1278 if (fType == "ON")
1279 {
1280 hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");
1281
1282 //Int_t nb = hpad->GetNbinsX();
1283 //for (Int_t i=1; i<=nb; i++)
1284 //{
1285 // if (hpad->GetBinContent(i) != 0.0)
1286 // *fLog << all << "i, fHgON = " << i << ", "
1287 // << hpad->GetBinContent(i) << endl;
1288 //}
1289
1290 if (nSigma != 0.0)
1291 prob = hpad->Integral() * nTheta/nSigma;
1292 else
1293 prob = 0.0;
1294
1295 *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma
1296 << ", " << prob << endl;
1297 }
1298 else if (fType == "OFF")
1299 {
1300 hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, "");
1301 if (nSigma != 0.0)
1302 prob = hpad->Integral() * nTheta/nSigma;
1303 else
1304 prob = 0.0;
1305 }
1306 else
1307 prob = 1.0;
1308
1309 if ( fType != "MC" &&
1310 (prob <= 0.0 || gRandom->Uniform() > prob) )
1311 {
1312 delete hpad;
1313 // event should not be padded
1314 *fLog << all << "event is not padded" << endl;
1315
1316 rc = 8;
1317 fErrors[rc]++;
1318 return kTRUE;
1319 }
1320 // event should be padded
1321 *fLog << all << "event is padded" << endl;
1322
1323
1324 //-------------------------------------------
1325 // for the current theta, generate a sigmabar
1326 //
1327 // for ON/OFF data according to the matrix fHgON/OFF
1328 // for MC data according to the histogram fHSigmaTheta
1329 //
1330 Double_t sigmabar_phot = 0;
1331 Double_t sigmabar = 0;
1332 if (fType == "ON" || fType == "OFF")
1333 {
1334 //Int_t nbinsx = hpad->GetNbinsX();
1335 //for (Int_t i=1; i<=nbinsx; i++)
1336 //{
1337 // if (hpad->GetBinContent(i) != 0.0)
1338 // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i)
1339 // << endl;
1340 //}
1341
1342 sigmabar_phot = hpad->GetRandom();
1343 sigmabar = sigmabar_phot * fPEperPhoton;
1344
1345 *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl;
1346
1347 delete hpad;
1348 }
1349
1350 else if (fType == "MC")
1351 {
1352 Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta);
1353
1354 if (binTheta != bincheck)
1355 {
1356 cout << "The binnings of the 2 histograms are not identical; aborting"
1357 << endl;
1358 return kERROR;
1359 }
1360
1361 TH1D *hsigma;
1362
1363 hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, "");
1364 if ( hsigma->GetEntries() == 0.0 )
1365 {
1366 *fLog << warn << "MPadONOFF::Process(); projection for Theta bin "
1367 << binTheta << " has no entries; Skip event " << endl;
1368 // event cannot be padded; skip event
1369 delete hsigma;
1370
1371 rc = 3;
1372 fErrors[rc]++;
1373 return kCONTINUE;
1374 }
1375 else
1376 {
1377 sigmabar_phot = hsigma->GetRandom();
1378 sigmabar = sigmabar_phot * fPEperPhoton;
1379
1380 *fLog << all << "Theta, bin number = " << theta << ", " << binTheta
1381 << ", sigmabar_phot = " << sigmabar_phot << endl;
1382 }
1383 delete hsigma;
1384 }
1385
1386 const Double_t sigmabar2 = sigmabar*sigmabar;
1387
1388 //-------------------------------------------
1389
1390 *fLog << all << "MPadONOFF::Process(); sigbarold, sigmabar = "
1391 << sigbarold << ", "<< sigmabar << endl;
1392
1393 // Skip event if target sigmabar is <= sigbarold
1394 if (sigmabar <= sigbarold)
1395 {
1396 *fLog << all << "MPadONOFF::Process(); target sigmabar is less than sigbarold : "
1397 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
1398
1399 rc = 4;
1400 fErrors[rc]++;
1401 return kCONTINUE;
1402 }
1403
1404
1405 //-------------------------------------------
1406 //
1407 // Calculate average number of NSB photo electrons to be added (lambdabar)
1408 // from the value of sigmabar,
1409 // - making assumptions about the average electronic noise (elNoise2) and
1410 // - using a fixed value (F2excess) for the excess noise factor
1411
1412 Double_t elNoise2; // [photo electrons]
1413 Double_t F2excess = 1.3;
1414 Double_t lambdabar; // [photo electrons]
1415
1416
1417
1418 Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta);
1419 if (binTheta != bincheck)
1420 {
1421 *fLog << err
1422 << "MPadONOFF::Process(); The binnings of the 2 histograms are not identical; aborting"
1423 << endl;
1424 return kERROR;
1425 }
1426
1427 // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
1428 // The average electronic noise (to be added) has to be well above this RMS,
1429 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
1430 // may become negative
1431
1432 TH1D *hnoise;
1433 if (fType == "MC")
1434 hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1435 else if (fType == "ON")
1436 hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1437 else if (fType == "OFF")
1438 hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1439 else
1440 {
1441 *fLog << err << "MPadONOFF::Process(); illegal data type... aborting"
1442 << endl;
1443 return kERROR;
1444 }
1445
1446 Double_t RMS_phot = hnoise->GetRMS(1);
1447 Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
1448 delete hnoise;
1449
1450 elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2);
1451 *fLog << all << "elNoise2 = " << elNoise2 << endl;
1452
1453 lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess;
1454
1455 // This value of lambdabar is the same for all pixels;
1456 // note that lambdabar is normalized to the area of pixel 0
1457
1458 //---------- start loop over pixels ---------------------------------
1459 // do the padding for each pixel
1460 //
1461 // pad only pixels - which are used (before image cleaning)
1462 //
1463 Double_t sig_phot = 0;
1464 Double_t sig = 0;
1465
1466 Double_t sigma2 = 0;
1467
1468 Double_t diff_phot = 0;
1469 Double_t diff = 0;
1470
1471 Double_t addSig2_phot= 0;
1472 Double_t addSig2 = 0;
1473
1474 Double_t elNoise2Pix = 0;
1475
1476
1477 for (UInt_t i=0; i<npix; i++)
1478 {
1479 MCerPhotPix &pix = (*fEvt)[i];
1480 if ( !pix.IsPixelUsed() )
1481 continue;
1482
1483 //if ( pix.GetNumPhotons() == 0.0)
1484 //{
1485 // *fLog << warn
1486 // << "MPadONOFF::Process(); no.of photons is 0 for used pixel"
1487 // << endl;
1488 // continue;
1489 //}
1490
1491 Int_t j = pix.GetPixId();
1492
1493 // GetPixRatio returns (area of pixel 0 / area of current pixel)
1494 Double_t ratioArea = 1.0 / fCam->GetPixRatio(j);
1495
1496 MPedestalPix &ppix = (*fPed)[j];
1497 Double_t oldsigma_phot = ppix.GetPedestalRms();
1498 Double_t oldsigma = oldsigma_phot * fPEperPhoton;
1499 Double_t oldsigma2 = oldsigma*oldsigma;
1500
1501 //---------------------------------
1502 // throw the Sigma for this pixel
1503 //
1504 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
1505
1506 Int_t count;
1507 Bool_t ok;
1508
1509 TH1D *hdiff;
1510 TH1D *hsig;
1511
1512 switch (fPadFlag)
1513 {
1514 case 1 :
1515 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
1516
1517 if (fType == "MC")
1518 hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
1519 binPixel, binPixel, "");
1520 else if (fType == "ON")
1521 hdiff = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta,
1522 binPixel, binPixel, "");
1523 else
1524 hdiff = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta,
1525 binPixel, binPixel, "");
1526
1527 if ( hdiff->GetEntries() == 0 )
1528 {
1529 *fLog << warn << "MPadONOFF::Process(); projection for Theta bin "
1530 << binTheta << " and pixel bin " << binPixel
1531 << " has no entries; aborting " << endl;
1532 delete hdiff;
1533
1534 rc = 5;
1535 fErrors[rc]++;
1536 return kCONTINUE;
1537 }
1538
1539 count = 0;
1540 ok = kFALSE;
1541 for (Int_t m=0; m<fIter; m++)
1542 {
1543 count += 1;
1544 diff_phot = hdiff->GetRandom();
1545 diff = diff_phot * fPEperPhoton * fPEperPhoton;
1546
1547 // the following condition ensures that elNoise2Pix > 0.0
1548 if ( (diff + sigmabar2 - oldsigma2/ratioArea
1549 - lambdabar*F2excess) > 0.0 )
1550 {
1551 ok = kTRUE;
1552 break;
1553 }
1554 }
1555 if (!ok)
1556 {
1557 *fLog << all << "theta, j, count, sigmabar, diff = " << theta
1558 << ", "
1559 << j << ", " << count << ", " << sigmabar << ", "
1560 << diff << endl;
1561 diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
1562
1563 rw = 1;
1564 fWarnings[rw]++;
1565 }
1566 else
1567 {
1568 rw = 0;
1569 fWarnings[rw]++;
1570 }
1571
1572 delete hdiff;
1573 sigma2 = diff + sigmabar2;
1574 break;
1575
1576 case 2 :
1577 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
1578
1579 if (fType == "MC")
1580 hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
1581 binPixel, binPixel, "");
1582 else if (fType == "ON")
1583 hsig = fHSigmaPixThetaOFF->ProjectionZ("", binTheta, binTheta,
1584 binPixel, binPixel, "");
1585 else
1586 hsig = fHSigmaPixThetaON->ProjectionZ("", binTheta, binTheta,
1587 binPixel, binPixel, "");
1588
1589 if ( hsig->GetEntries() == 0 )
1590 {
1591 *fLog << warn << "MPadONOFF::Process(); projection for Theta bin "
1592 << binTheta << " and pixel bin " << binPixel
1593 << " has no entries; aborting " << endl;
1594 delete hsig;
1595
1596 rc = 6;
1597 fErrors[rc]++;
1598 return kCONTINUE;
1599 }
1600
1601 count = 0;
1602 ok = kFALSE;
1603 for (Int_t m=0; m<fIter; m++)
1604 {
1605 count += 1;
1606
1607 sig_phot = hsig->GetRandom();
1608 sig = sig_phot * fPEperPhoton;
1609
1610 sigma2 = sig*sig/ratioArea;
1611 // the following condition ensures that elNoise2Pix > 0.0
1612
1613 if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 )
1614 {
1615 ok = kTRUE;
1616 break;
1617 }
1618 }
1619 if (!ok)
1620 {
1621 *fLog << all << "theta, j, count, sigmabar, sig = " << theta
1622 << ", "
1623 << j << ", " << count << ", " << sigmabar << ", "
1624 << sig << endl;
1625 sigma2 = lambdabar*F2excess + oldsigma2/ratioArea;
1626
1627 rw = 1;
1628 fWarnings[rw]++;
1629 }
1630 else
1631 {
1632 rw = 0;
1633 fWarnings[rw]++;
1634 }
1635
1636 delete hsig;
1637 break;
1638 }
1639
1640 //---------------------------------
1641 // get the additional sigma^2 for this pixel (due to the padding)
1642
1643 addSig2 = sigma2*ratioArea - oldsigma2;
1644 addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
1645
1646 //---------------------------------
1647 // get the additional electronic noise for this pixel
1648
1649 elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;
1650
1651
1652 //---------------------------------
1653 // throw actual number of additional NSB photo electrons (NSB)
1654 // and its RMS (sigmaNSB)
1655
1656 Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
1657 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
1658 Double_t sigmaNSB0;
1659
1660 if (arg >= 0.0)
1661 {
1662 sigmaNSB0 = sqrt( arg );
1663 }
1664 else
1665 {
1666 sigmaNSB0 = 0.0000001;
1667 if (arg < -1.e-10)
1668 {
1669 *fLog << warn << "MPadONOFF::Process(); argument of sqrt < 0 : "
1670 << arg << endl;
1671 }
1672 }
1673
1674
1675 //---------------------------------
1676 // smear NSB0 according to sigmaNSB0
1677 // and subtract lambdabar because of AC coupling
1678
1679 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
1680 Double_t NSB_phot = NSB / fPEperPhoton;
1681
1682 //---------------------------------
1683
1684 // add additional NSB to the number of photons
1685 Double_t oldphotons_phot = pix.GetNumPhotons();
1686 Double_t oldphotons = oldphotons_phot * fPEperPhoton;
1687 Double_t newphotons = oldphotons + NSB;
1688 Double_t newphotons_phot = newphotons / fPEperPhoton;
1689 pix.SetNumPhotons( newphotons_phot );
1690
1691
1692 fHNSB->Fill( NSB_phot/ratioArea );
1693 fHPhotons->Fill( oldphotons_phot/ratioArea,
1694 newphotons_phot/ratioArea );
1695
1696
1697 // error: add sigma of padded noise quadratically
1698 Double_t olderror_phot = pix.GetErrorPhot();
1699 Double_t newerror_phot =
1700 sqrt( olderror_phot*olderror_phot + addSig2_phot );
1701 pix.SetErrorPhot( newerror_phot );
1702
1703
1704 Double_t newsigma = sqrt( oldsigma2 + addSig2 );
1705 Double_t newsigma_phot = newsigma / fPEperPhoton;
1706 ppix.SetPedestalRms( newsigma_phot );
1707
1708 fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot );
1709 }
1710 //---------- end of loop over pixels ---------------------------------
1711
1712 // Calculate sigmabar again and crosscheck
1713
1714 fSigmabar->Calc(*fCam, *fPed, *fEvt);
1715 *fLog << all << "MPadONOFF::Process(); after padding : " << endl;
1716 fSigmabar->Print("");
1717
1718 *fLog << all << "Exit MPadONOFF::Process();" << endl;
1719
1720 rc = 0;
1721 fErrors[rc]++;
1722 return kTRUE;
1723 //******************************************************************
1724}
1725
1726// --------------------------------------------------------------------------
1727//
1728//
1729Int_t MPadONOFF::PostProcess()
1730{
1731 if (GetNumExecutions() != 0)
1732 {
1733
1734 *fLog << inf << endl;
1735 *fLog << GetDescriptor() << " execution statistics:" << endl;
1736 *fLog << dec << setfill(' ');
1737
1738 int temp;
1739 if (fWarnings[0] != 0.0)
1740 temp = (int)(fWarnings[1]*100/fWarnings[0]);
1741 else
1742 temp = 100;
1743
1744 *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3)
1745 << temp
1746 << "%) Warning : iteration to find acceptable sigma failed"
1747 << ", fIter = " << fIter << endl;
1748
1749 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
1750 << (int)(fErrors[1]*100/GetNumExecutions())
1751 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
1752
1753 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
1754 << (int)(fErrors[2]*100/GetNumExecutions())
1755 << "%) Evts skipped due to: Zenith angle out of range" << endl;
1756
1757 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
1758 << (int)(fErrors[3]*100/GetNumExecutions())
1759 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
1760
1761 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
1762 << (int)(fErrors[4]*100/GetNumExecutions())
1763 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
1764
1765 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
1766 << (int)(fErrors[5]*100/GetNumExecutions())
1767 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
1768
1769 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
1770 << (int)(fErrors[6]*100/GetNumExecutions())
1771 << "%) Evts skipped due to: No data for generating Sigma" << endl;
1772
1773 *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3)
1774 << (int)(fErrors[7]*100/GetNumExecutions())
1775 << "%) Evts skipped due to: No data for generating Blind pixels" << endl;
1776
1777 *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3)
1778 << (int)(fErrors[8]*100/GetNumExecutions())
1779 << "%) Evts didn't have to be padded" << endl;
1780
1781 *fLog << " " << fErrors[0] << " ("
1782 << (int)(fErrors[0]*100/GetNumExecutions())
1783 << "%) Evts were successfully padded!" << endl;
1784 *fLog << endl;
1785
1786 }
1787
1788 //---------------------------------------------------------------
1789 TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200));
1790 c.Divide(3, 4);
1791 gROOT->SetSelectedPad(NULL);
1792
1793 c.cd(1);
1794 fHNSB->SetDirectory(NULL);
1795 fHNSB->DrawCopy();
1796 fHNSB->SetBit(kCanDelete);
1797
1798 c.cd(2);
1799 fHSigmaPedestal->SetDirectory(NULL);
1800 fHSigmaPedestal->DrawCopy();
1801 fHSigmaPedestal->SetBit(kCanDelete);
1802
1803 c.cd(3);
1804 fHPhotons->SetDirectory(NULL);
1805 fHPhotons->DrawCopy();
1806 fHPhotons->SetBit(kCanDelete);
1807
1808 //--------------------------------------------------------------------
1809
1810
1811 c.cd(4);
1812 fHSigmaTheta->SetDirectory(NULL);
1813 fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta");
1814 fHSigmaTheta->DrawCopy();
1815 fHSigmaTheta->SetBit(kCanDelete);
1816
1817
1818 //--------------------------------------------------------------------
1819
1820
1821 c.cd(7);
1822 fHBlindPixNTheta->SetDirectory(NULL);
1823 fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta");
1824 fHBlindPixNTheta->DrawCopy();
1825 fHBlindPixNTheta->SetBit(kCanDelete);
1826
1827 //--------------------------------------------------------------------
1828
1829
1830 c.cd(10);
1831 fHBlindPixIdTheta->SetDirectory(NULL);
1832 fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta");
1833 fHBlindPixIdTheta->DrawCopy();
1834 fHBlindPixIdTheta->SetBit(kCanDelete);
1835
1836
1837 //--------------------------------------------------------------------
1838 // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2
1839
1840 c.cd(5);
1841 TH2D *l1;
1842 l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
1843 l1->SetDirectory(NULL);
1844 l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
1845 l1->SetXTitle("\\Theta [\\circ]");
1846 l1->SetYTitle("Sigma^2-Sigmabar^2");
1847
1848 l1->DrawCopy("box");
1849 l1->SetBit(kCanDelete);;
1850
1851 c.cd(8);
1852 TH2D *l2;
1853 l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
1854 l2->SetDirectory(NULL);
1855 l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
1856 l2->SetXTitle("pixel");
1857 l2->SetYTitle("Sigma^2-Sigmabar^2");
1858
1859 l2->DrawCopy("box");
1860 l2->SetBit(kCanDelete);;
1861
1862 //--------------------------------------------------------------------
1863 // draw the 3D histogram (target): Theta, pixel, Sigma
1864
1865 c.cd(6);
1866 TH2D *k1;
1867 k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
1868 k1->SetDirectory(NULL);
1869 k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)");
1870 k1->SetXTitle("\\Theta [\\circ]");
1871 k1->SetYTitle("Sigma");
1872
1873 k1->DrawCopy("box");
1874 k1->SetBit(kCanDelete);
1875
1876 c.cd(9);
1877 TH2D *k2;
1878 k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
1879 k2->SetDirectory(NULL);
1880 k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)");
1881 k2->SetXTitle("pixel");
1882 k2->SetYTitle("Sigma");
1883
1884 k2->DrawCopy("box");
1885 k2->SetBit(kCanDelete);;
1886
1887
1888 //--------------------------------------------------------------------
1889
1890 c.cd(11);
1891 TH2D *m1;
1892 m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
1893 m1->SetDirectory(NULL);
1894 m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
1895 m1->SetXTitle("Sigmabar-old");
1896 m1->SetYTitle("Sigmabar-new");
1897
1898 m1->DrawCopy("box");
1899 m1->SetBit(kCanDelete);;
1900
1901 c.cd(12);
1902 TH2D *m2;
1903 m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
1904 m2->SetDirectory(NULL);
1905 m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
1906 m2->SetXTitle("Sigmabar-old");
1907 m2->SetYTitle("Sigmabar-new");
1908
1909 m2->DrawCopy("box");
1910 m2->SetBit(kCanDelete);;
1911
1912
1913 return kTRUE;
1914}
1915
1916// --------------------------------------------------------------------------
1917
1918
1919
1920
1921
1922
Note: See TracBrowser for help on using the repository browser.