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

Last change on this file since 2717 was 2663, checked in by wittek, 22 years ago
*** empty log message ***
File size: 52.2 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 << "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 << "----------------------------------------------------------------------------------" << endl;
178 *fLog << all << "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 << "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 << "Read; 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 << "The target distributions for the padding have been created"
645 << endl;
646 *fLog << all << "----------------------------------------------------------"
647 << endl;
648 //--------------------------------------------
649
650 fHSigmaTheta->SetDirectory(NULL);
651 fHSigmaThetaON->SetDirectory(NULL);
652 fHSigmaThetaOFF->SetDirectory(NULL);
653
654 fHSigmaPixTheta->SetDirectory(NULL);
655 fHSigmaPixThetaON->SetDirectory(NULL);
656 fHSigmaPixThetaOFF->SetDirectory(NULL);
657
658 fHDiffPixTheta->SetDirectory(NULL);
659 fHDiffPixThetaON->SetDirectory(NULL);
660 fHDiffPixThetaOFF->SetDirectory(NULL);
661
662 fHBlindPixIdTheta->SetDirectory(NULL);
663 fHBlindPixNTheta->SetDirectory(NULL);
664
665
666 fHgON->SetDirectory(NULL);
667 fHgOFF->SetDirectory(NULL);
668
669
670 return kTRUE;
671}
672
673
674// --------------------------------------------------------------------------
675//
676// Read target distributions from a file
677//
678//
679Bool_t MPadONOFF::ReadTargetDist(const char* namefilein)
680{
681 *fLog << all << "Read padding histograms from file " << namefilein << endl;
682
683 fInfile = new TFile(namefilein);
684 fInfile->ls();
685
686 //------------------------------------
687
688 fHSigmaTheta =
689 (TH2D*) gROOT->FindObject("2D-ThetaSigmabar");
690 if (!fHSigmaTheta)
691 {
692 *fLog << all << "Object '2D-ThetaSigmabar' not found on root file"
693 << endl;
694 return kFALSE;
695 }
696 *fLog << all << "Object '2D-ThetaSigmabar' was read in" << endl;
697
698 fHSigmaThetaON =
699 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON");
700 if (!fHSigmaThetaON)
701 {
702 *fLog << all << "Object '2D-ThetaSigmabarON' not found on root file"
703 << endl;
704 return kFALSE;
705 }
706 *fLog << all << "Object '2D-ThetaSigmabarON' was read in" << endl;
707
708 fHSigmaThetaOFF =
709 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF");
710 if (!fHSigmaThetaOFF)
711 {
712 *fLog << "Object '2D-ThetaSigmabarOFF' not found on root file"
713 << endl;
714 return kFALSE;
715 }
716 *fLog << all << "Object '2D-ThetaSigmabarOFF' was read in" << endl;
717
718 fHSigmaPixTheta =
719 (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
720 if (!fHSigmaPixTheta)
721 {
722 *fLog << all << "Object '3D-ThetaPixSigma' not found on root file"
723 << endl;
724 return kFALSE;
725 }
726 *fLog << all << "Object '3D-ThetaPixSigma' was read in" << endl;
727
728 fHSigmaPixThetaON =
729 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON");
730 if (!fHSigmaPixThetaON)
731 {
732 *fLog << all << "Object '3D-ThetaPixSigmaON' not found on root file"
733 << endl;
734 return kFALSE;
735 }
736 *fLog << all << "Object '3D-ThetaPixSigmaON' was read in" << endl;
737
738 fHSigmaPixThetaOFF =
739 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF");
740 if (!fHSigmaPixThetaOFF)
741 {
742 *fLog << all << "Object '3D-ThetaPixSigmaOFF' not found on root file"
743 << endl;
744 return kFALSE;
745 }
746 *fLog << "Object '3D-ThetaPixSigmaOFF' was read in" << endl;
747
748 fHDiffPixTheta =
749 (TH3D*) gROOT->FindObject("3D-ThetaPixDiff");
750 if (!fHDiffPixTheta)
751 {
752 *fLog << all << "Object '3D-ThetaPixDiff' not found on root file"
753 << endl;
754 return kFALSE;
755 }
756 *fLog << all << "Object '3D-ThetaPixDiff' was read in" << endl;
757
758 fHDiffPixThetaON =
759 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON");
760 if (!fHDiffPixThetaON)
761 {
762 *fLog << all << "Object '3D-ThetaPixDiffON' not found on root file"
763 << endl;
764 return kFALSE;
765 }
766 *fLog << "Object '3D-ThetaPixDiffON' was read in" << endl;
767
768 fHDiffPixThetaOFF =
769 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF");
770 if (!fHDiffPixThetaOFF)
771 {
772 *fLog << all << "Object '3D-ThetaPixDiffOFF' not found on root file"
773 << endl;
774 return kFALSE;
775 }
776 *fLog << all << "Object '3D-ThetaPixDiffOFF' was read in" << endl;
777
778 fHgON =
779 (TH3D*) gROOT->FindObject("3D-PaddingMatrixON");
780 if (!fHgON)
781 {
782 *fLog << all << "Object '3D-PaddingMatrixON' not found on root file"
783 << endl;
784 return kFALSE;
785 }
786 *fLog << all << "Object '3D-PaddingMatrixON' was read in" << endl;
787
788 fHgOFF =
789 (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF");
790 if (!fHgOFF)
791 {
792 *fLog << all << "Object '3D-PaddingMatrixOFF' not found on root file"
793 << endl;
794 return kFALSE;
795 }
796 *fLog << all << "Object '3D-PaddingMatrixOFF' was read in" << endl;
797
798
799 fHBlindPixIdTheta =
800 (TH2D*) gROOT->FindObject("2D-ThetaBlindId");
801 if (!fHBlindPixIdTheta)
802 {
803 *fLog << all << "Object '2D-ThetaBlindId' not found on root file"
804 << endl;
805 return kFALSE;
806 }
807 *fLog << all << "Object '2D-ThetaBlindId' was read in" << endl;
808
809 fHBlindPixNTheta =
810 (TH2D*) gROOT->FindObject("2D-ThetaBlindN");
811 if (!fHBlindPixNTheta)
812 {
813 *fLog << all << "Object '2D-ThetaBlindN' not found on root file"
814 << endl;
815 return kFALSE;
816 }
817 *fLog << all << "Object '2D-ThetaBlindN' was read in" << endl;
818
819
820 //------------------------------------
821
822 return kTRUE;
823}
824
825
826// --------------------------------------------------------------------------
827//
828// Write target distributions onto a file
829//
830//
831Bool_t MPadONOFF::WriteTargetDist(const char* namefileout)
832{
833 *fLog << all << "Write padding histograms onto file " << namefileout << endl;
834
835 TFile outfile(namefileout, "RECREATE");
836
837 fHBlindPixNTheta->Write();
838 fHBlindPixIdTheta->Write();
839
840 fHSigmaTheta->Write();
841 fHSigmaThetaON->Write();
842 fHSigmaThetaOFF->Write();
843
844 fHSigmaPixTheta->Write();
845 fHSigmaPixThetaON->Write();
846 fHSigmaPixThetaOFF->Write();
847
848 fHDiffPixTheta->Write();
849 fHDiffPixThetaON->Write();
850 fHDiffPixThetaOFF->Write();
851
852 fHgON->Write();
853 fHgOFF->Write();
854
855 *fLog << all << "WriteTargetDist; target histograms written onto file "
856 << namefileout << endl;
857
858 return kTRUE;
859}
860
861// --------------------------------------------------------------------------
862//
863// Set type of data to be padded
864//
865// this is not necessary if the type of the data can be recognized
866// directly from the input
867//
868//
869void MPadONOFF::SetDataType(const char* type)
870{
871 fType = type;
872 *fLog << all << "SetDataType(); type of data to be padded : "
873 << fType << endl;
874
875 return;
876}
877
878
879// --------------------------------------------------------------------------
880//
881// Set the option for the padding
882//
883// There are 2 options for the padding :
884//
885// 1) fPadFlag = 1 :
886// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta
887// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using
888// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2
889// (fHDiffPixTheta).
890//
891// This is the preferred option as it takes into account the
892// correlations between the Sigma of a pixel and Sigmabar.
893//
894// NOT YET TESTED :
895// 2) fPadFlag = 2 :
896// Generate a pedestal sigma for each pixel using the 3D-histogram
897// Theta, pixel no., Sigma (fHSigmaPixTheta).
898//
899void MPadONOFF::SetPadFlag(Int_t padflag)
900{
901 fPadFlag = padflag;
902 *fLog << all << "SetPadFlag(); choose option " << fPadFlag << endl;
903}
904
905// --------------------------------------------------------------------------
906//
907// Set the pointers and prepare the histograms
908//
909Int_t MPadONOFF::PreProcess(MParList *pList)
910{
911 if ( !fHSigmaTheta || !fHSigmaThetaON || !fHSigmaThetaOFF ||
912 !fHSigmaPixTheta || !fHSigmaPixThetaON || !fHSigmaPixThetaOFF ||
913 !fHDiffPixTheta || !fHDiffPixThetaON || !fHDiffPixThetaOFF ||
914 !fHBlindPixIdTheta || !fHBlindPixNTheta ||
915 !fHgON || !fHgOFF)
916 {
917 *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl;
918 return kFALSE;
919 }
920
921 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
922 if (!fMcEvt)
923 {
924 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
925 return kFALSE;
926 }
927
928 fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
929 if (!fPed)
930 {
931 *fLog << err << "MPedestalCam not found... aborting." << endl;
932 return kFALSE;
933 }
934
935 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
936 if (!fCam)
937 {
938 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
939 return kFALSE;
940 }
941
942 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
943 if (!fEvt)
944 {
945 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
946 return kFALSE;
947 }
948
949 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
950 if (!fSigmabar)
951 {
952 *fLog << err << "MSigmabar not found... aborting." << endl;
953 return kFALSE;
954 }
955
956 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
957 if (!fBlinds)
958 {
959 *fLog << err << "MBlindPixels not found... aborting." << endl;
960 return kFALSE;
961 }
962
963 if (fType !="ON" && fType !="OFF" && fType !="MC")
964 {
965 *fLog << err << "Type of data to be padded not defined... aborting." << endl;
966 return kFALSE;
967 }
968
969
970 //--------------------------------------------------------------------
971 // histograms for checking the padding
972 //
973 // plot pedestal sigmas
974 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
975 100, 0.0, 5.0, 100, 0.0, 5.0);
976 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
977 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
978
979 // plot no.of photons (before vs. after padding)
980 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
981 100, -10.0, 90.0, 100, -10, 90);
982 fHPhotons->SetXTitle("no.of photons before padding");
983 fHPhotons->SetYTitle("no.of photons after padding");
984
985 // plot of added NSB
986 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
987 fHNSB->SetXTitle("no.of added NSB photons");
988 fHNSB->SetYTitle("no.of pixels");
989
990
991 //--------------------------------------------------------------------
992
993 fIter = 20;
994
995 memset(fErrors, 0, sizeof(fErrors));
996
997 return kTRUE;
998}
999
1000// --------------------------------------------------------------------------
1001//
1002// Do the Padding
1003// idealy the events to be padded should have been generated without NSB
1004//
1005Int_t MPadONOFF::Process()
1006{
1007 *fLog << all << "Entry Process();" << endl;
1008
1009 // this is the conversion factor from the number of photons
1010 // to the number of photo electrons
1011 // later to be taken from MCalibration
1012 // fPEpePhoton = PW * LG * QE * 1D
1013 Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9;
1014
1015 //------------------------------------------------
1016 // add pixels to MCerPhotEvt which are not yet in;
1017 // set their number of photons equal to zero
1018
1019 UInt_t imaxnumpix = fCam->GetNumPixels();
1020
1021 for (UInt_t i=0; i<imaxnumpix; i++)
1022 {
1023 Bool_t alreadythere = kFALSE;
1024 UInt_t npix = fEvt->GetNumPixels();
1025 for (UInt_t j=0; j<npix; j++)
1026 {
1027 MCerPhotPix &pix = (*fEvt)[j];
1028 UInt_t id = pix.GetPixId();
1029 if (i==id)
1030 {
1031 alreadythere = kTRUE;
1032 break;
1033 }
1034 }
1035 if (alreadythere)
1036 continue;
1037
1038 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetPedestalRms());
1039 }
1040
1041
1042
1043 //-----------------------------------------
1044 Int_t rc=0;
1045 Int_t rw=0;
1046
1047 const UInt_t npix = fEvt->GetNumPixels();
1048
1049 //fSigmabar->Calc(*fCam, *fPed, *fEvt);
1050 // *fLog << all << "before padding : " << endl;
1051 //fSigmabar->Print("");
1052
1053
1054 //$$$$$$$$$$$$$$$$$$$$$$$$$$
1055 // to simulate the situation that before the padding the NSB and
1056 // electronic noise are zero : set Sigma = 0 for all pixels
1057 //for (UInt_t i=0; i<npix; i++)
1058 //{
1059 // MCerPhotPix &pix = fEvt->operator[](i);
1060 // Int_t j = pix.GetPixId();
1061
1062 // MPedestalPix &ppix = fPed->operator[](j);
1063 // ppix.SetMeanRms(0.0);
1064 //}
1065 //$$$$$$$$$$$$$$$$$$$$$$$$$$
1066
1067 //-------------------------------------------
1068 // Calculate average sigma of the event
1069 //
1070 Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt);
1071 *fLog << all << "before padding : " << endl;
1072 fSigmabar->Print("");
1073 Double_t sigbarold = sigbarold_phot * fPEperPhoton;
1074
1075
1076 Double_t sigbarold2 = sigbarold*sigbarold;
1077
1078
1079 // for MC data : expect sigmabar to be zero before the padding
1080 if (fType == "MC")
1081 {
1082 if (sigbarold_phot > 0)
1083 {
1084 *fLog << warn << "Process(); sigmabar of event to be padded is > 0 : "
1085 << sigbarold_phot << ". Stop event loop " << endl;
1086 // input data should have sigmabar = 0; stop event loop
1087
1088 rc = 1;
1089 fErrors[rc]++;
1090 return kCONTINUE;
1091 }
1092 }
1093
1094 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
1095 *fLog << all << "theta = " << theta << endl;
1096
1097 Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
1098 if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() )
1099 {
1100 *fLog << warn << "Process(); binNumber out of range : theta, binTheta = "
1101 << theta << ", " << binTheta << "; Skip event " << endl;
1102 // event cannot be padded; skip event
1103
1104 rc = 2;
1105 fErrors[rc]++;
1106 return kCONTINUE;
1107 }
1108
1109 //-------------------------------------------
1110 // get number of events in this theta bin
1111 // and number of events in this sigbarold_phot bin
1112
1113 Double_t nTheta;
1114 Double_t nSigma;
1115 if (fType == "ON")
1116 {
1117 TH1D *hn;
1118
1119 hn = fHSigmaThetaON->ProjectionY("", binTheta, binTheta, "");
1120 nTheta = hn->Integral();
1121 Int_t binSigma = hn->FindBin(sigbarold_phot);
1122 nSigma = hn->GetBinContent(binSigma);
1123
1124 *fLog << all
1125 << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1126 << theta << ", " << sigbarold_phot << ", " << binTheta << ", "
1127 << binSigma << ", " << nTheta << ", " << nSigma << endl;
1128
1129 delete hn;
1130 }
1131
1132 else if (fType == "OFF")
1133 {
1134 TH1D *hn;
1135
1136 hn = fHSigmaThetaOFF->ProjectionY("", binTheta, binTheta, "");
1137 nTheta = hn->Integral();
1138 Int_t binSigma = hn->FindBin(sigbarold_phot);
1139 nSigma = hn->GetBinContent(binSigma);
1140
1141 *fLog << all
1142 << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1143 << theta << ", " << sigbarold_phot << ", " << binTheta << ", "
1144 << binSigma << ", " << nTheta << ", " << nSigma << endl;
1145
1146 delete hn;
1147 }
1148
1149 else
1150 {
1151 nTheta = 0.0;
1152 nSigma = 0.0;
1153 }
1154
1155 //-------------------------------------------
1156 // for the current theta,
1157 // generate blind pixels according to the histograms
1158 // fHBlindPixNTheta and fHBlindPixIDTheta
1159 // do this only for MC data
1160
1161
1162 if (fType == "MC")
1163 {
1164
1165 // numBlind is the number of blind pixels in this event
1166 TH1D *nblind;
1167 UInt_t numBlind;
1168
1169 nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, "");
1170 if ( nblind->Integral() == 0.0 )
1171 {
1172 *fLog << warn << "Process(); projection for Theta bin "
1173 << binTheta << " has no entries; Skip event " << endl;
1174 // event cannot be padded; skip event
1175 delete nblind;
1176
1177 rc = 7;
1178 fErrors[rc]++;
1179 return kCONTINUE;
1180 }
1181 else
1182 {
1183 numBlind = (Int_t) (nblind->GetRandom()+0.5);
1184 }
1185 delete nblind;
1186
1187 //warn Code commented out due no compilation errors on Alpha OSF (tgb)
1188
1189 // throw the Id of numBlind different pixels in this event
1190 TH1D *hblind;
1191 UInt_t idBlind;
1192 UInt_t listId[npix];
1193 UInt_t nlist = 0;
1194 Bool_t equal;
1195
1196 hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, "");
1197 if ( hblind->Integral() > 0.0 )
1198 for (UInt_t i=0; i<numBlind; i++)
1199 {
1200 while(1)
1201 {
1202 idBlind = (Int_t) (hblind->GetRandom()+0.5);
1203 equal = kFALSE;
1204 for (UInt_t j=0; j<nlist; j++)
1205 if (idBlind == listId[j])
1206 {
1207 equal = kTRUE;
1208 break;
1209 }
1210 if (!equal) break;
1211 }
1212 listId[nlist] = idBlind;
1213 nlist++;
1214
1215 fBlinds->SetPixelBlind(idBlind);
1216 *fLog << all << "idBlind = " << idBlind << endl;
1217 }
1218 fBlinds->SetReadyToSave();
1219
1220 delete hblind;
1221
1222 }
1223
1224 //******************************************************************
1225 // has the event to be padded ?
1226 // if yes : to which sigmabar should it be padded ?
1227 //
1228
1229 Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot);
1230 *fLog << all << "binSig, sigbarold_phot = " << binSig << ", "
1231 << sigbarold_phot << endl;
1232
1233 Double_t prob;
1234 TH1D *hpad = NULL;
1235 if (fType == "ON")
1236 {
1237 hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");
1238
1239 //Int_t nb = hpad->GetNbinsX();
1240 //for (Int_t i=1; i<=nb; i++)
1241 //{
1242 // if (hpad->GetBinContent(i) != 0.0)
1243 // *fLog << all << "i, fHgON = " << i << ", "
1244 // << hpad->GetBinContent(i) << endl;
1245 //}
1246
1247 if (nSigma != 0.0)
1248 prob = hpad->Integral() * nTheta/nSigma;
1249 else
1250 prob = 0.0;
1251
1252 *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma
1253 << ", " << prob << endl;
1254 }
1255 else if (fType == "OFF")
1256 {
1257 hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, "");
1258 if (nSigma != 0.0)
1259 prob = hpad->Integral() * nTheta/nSigma;
1260 else
1261 prob = 0.0;
1262 }
1263 else
1264 prob = 1.0;
1265
1266 if ( fType != "MC" &&
1267 (prob <= 0.0 || gRandom->Uniform() > prob) )
1268 {
1269 delete hpad;
1270 // event should not be padded
1271 *fLog << all << "event is not padded" << endl;
1272
1273 rc = 8;
1274 fErrors[rc]++;
1275 return kTRUE;
1276 }
1277 // event should be padded
1278 *fLog << all << "event is padded" << endl;
1279
1280
1281 //-------------------------------------------
1282 // for the current theta, generate a sigmabar
1283 //
1284 // for ON/OFF data according to the matrix fHgON/OFF
1285 // for MC data according to the histogram fHSigmaTheta
1286 //
1287 Double_t sigmabar_phot = 0;
1288 Double_t sigmabar = 0;
1289 if (fType == "ON" || fType == "OFF")
1290 {
1291 //Int_t nbinsx = hpad->GetNbinsX();
1292 //for (Int_t i=1; i<=nbinsx; i++)
1293 //{
1294 // if (hpad->GetBinContent(i) != 0.0)
1295 // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i)
1296 // << endl;
1297 //}
1298
1299 sigmabar_phot = hpad->GetRandom();
1300 sigmabar = sigmabar_phot * fPEperPhoton;
1301
1302 *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl;
1303
1304 delete hpad;
1305 }
1306
1307 else if (fType == "MC")
1308 {
1309 Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta);
1310
1311 if (binTheta != bincheck)
1312 {
1313 cout << "The binnings of the 2 histograms are not identical; aborting"
1314 << endl;
1315 return kERROR;
1316 }
1317
1318 TH1D *hsigma;
1319
1320 hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, "");
1321 if ( hsigma->Integral() == 0.0 )
1322 {
1323 *fLog << warn << "Process(); projection for Theta bin "
1324 << binTheta << " has no entries; Skip event " << endl;
1325 // event cannot be padded; skip event
1326 delete hsigma;
1327
1328 rc = 3;
1329 fErrors[rc]++;
1330 return kCONTINUE;
1331 }
1332 else
1333 {
1334 sigmabar_phot = hsigma->GetRandom();
1335 sigmabar = sigmabar_phot * fPEperPhoton;
1336
1337 *fLog << all << "Theta, bin number = " << theta << ", " << binTheta
1338 << ", sigmabar_phot = " << sigmabar_phot << endl;
1339 }
1340 delete hsigma;
1341 }
1342
1343 const Double_t sigmabar2 = sigmabar*sigmabar;
1344
1345 //-------------------------------------------
1346
1347 *fLog << all << "Process(); sigbarold, sigmabar = "
1348 << sigbarold << ", "<< sigmabar << endl;
1349
1350 // Skip event if target sigmabar is <= sigbarold
1351 if (sigmabar <= sigbarold)
1352 {
1353 *fLog << all << "Process(); target sigmabar is less than sigbarold : "
1354 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
1355
1356 rc = 4;
1357 fErrors[rc]++;
1358 return kCONTINUE;
1359 }
1360
1361
1362 //-------------------------------------------
1363 //
1364 // Calculate average number of NSB photo electrons to be added (lambdabar)
1365 // from the value of sigmabar,
1366 // - making assumptions about the average electronic noise (elNoise2) and
1367 // - using a fixed value (F2excess) for the excess noise factor
1368
1369 Double_t elNoise2; // [photons]
1370 Double_t F2excess = 1.3;
1371 Double_t lambdabar; // [photons]
1372
1373
1374
1375 Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta);
1376 if (binTheta != bincheck)
1377 {
1378 cout << "The binnings of the 2 histograms are not identical; aborting"
1379 << endl;
1380 return kERROR;
1381 }
1382
1383 // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
1384 // The average electronic noise (to be added) has to be well above this RMS,
1385 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
1386 // may become negative
1387
1388 TH1D *hnoise;
1389 if (fType == "MC")
1390 hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1391 else if (fType == "ON")
1392 hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1393 else if (fType == "OFF")
1394 hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1395 else
1396 {
1397 *fLog << err << "Process; illegal data type... aborting" << endl;
1398 return kERROR;
1399 }
1400
1401 Double_t RMS_phot = hnoise->GetRMS(1);
1402 Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
1403 delete hnoise;
1404
1405 elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2);
1406 *fLog << all << "elNoise2 = " << elNoise2 << endl;
1407
1408 lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photo electrons]
1409
1410 // This value of lambdabar is the same for all pixels;
1411 // note that lambdabar is normalized to the area of pixel 0
1412
1413 //---------- start loop over pixels ---------------------------------
1414 // do the padding for each pixel
1415 //
1416 // pad only pixels - which are used (before image cleaning)
1417 //
1418 Double_t sig_phot = 0;
1419 Double_t sig = 0;
1420
1421 Double_t sigma2 = 0;
1422
1423 Double_t diff_phot = 0;
1424 Double_t diff = 0;
1425
1426 Double_t addSig2_phot= 0;
1427 Double_t addSig2 = 0;
1428
1429 Double_t elNoise2Pix = 0;
1430
1431
1432 for (UInt_t i=0; i<npix; i++)
1433 {
1434 MCerPhotPix &pix = (*fEvt)[i];
1435 if ( !pix.IsPixelUsed() )
1436 continue;
1437
1438 //if ( pix.GetNumPhotons() == 0.0)
1439 //{
1440 // *fLog << warn << "Process(); no.of photons is 0 for used pixel"
1441 // << endl;
1442 // continue;
1443 //}
1444
1445 Int_t j = pix.GetPixId();
1446
1447 // GetPixRatio returns (area of pixel 0 / area of current pixel)
1448 Double_t ratioArea = 1.0 / fCam->GetPixRatio(j);
1449
1450 MPedestalPix &ppix = (*fPed)[j];
1451 Double_t oldsigma_phot = ppix.GetPedestalRms();
1452 Double_t oldsigma = oldsigma_phot * fPEperPhoton;
1453 Double_t oldsigma2 = oldsigma*oldsigma;
1454
1455 //---------------------------------
1456 // throw the Sigma for this pixel
1457 //
1458 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
1459
1460 Int_t count;
1461 Bool_t ok;
1462
1463 TH1D *hdiff;
1464 TH1D *hsig;
1465
1466 switch (fPadFlag)
1467 {
1468 case 1 :
1469 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
1470
1471 if (fType == "MC")
1472 hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta,
1473 binPixel, binPixel, "");
1474 else if (fType == "ON")
1475 hdiff = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta,
1476 binPixel, binPixel, "");
1477 else
1478 hdiff = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta,
1479 binPixel, binPixel, "");
1480
1481 if ( hdiff->Integral() == 0 )
1482 {
1483 *fLog << warn << "Process(); projection for Theta bin "
1484 << binTheta << " and pixel bin " << binPixel
1485 << " has no entries; aborting " << endl;
1486 delete hdiff;
1487
1488 rc = 5;
1489 fErrors[rc]++;
1490 return kCONTINUE;
1491 }
1492
1493 count = 0;
1494 ok = kFALSE;
1495 for (Int_t m=0; m<fIter; m++)
1496 {
1497 count += 1;
1498 diff_phot = hdiff->GetRandom();
1499 diff = diff_phot * fPEperPhoton * fPEperPhoton;
1500
1501 // the following condition ensures that elNoise2Pix > 0.0
1502 if ( (diff + sigmabar2 - oldsigma2/ratioArea
1503 - lambdabar*F2excess) > 0.0 )
1504 {
1505 ok = kTRUE;
1506 break;
1507 }
1508 }
1509 if (!ok)
1510 {
1511 *fLog << all << "theta, j, count, sigmabar, diff = " << theta
1512 << ", "
1513 << j << ", " << count << ", " << sigmabar << ", "
1514 << diff << endl;
1515 diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
1516
1517 rw = 1;
1518 fWarnings[rw]++;
1519 }
1520 else
1521 {
1522 rw = 0;
1523 fWarnings[rw]++;
1524 }
1525
1526 delete hdiff;
1527 sigma2 = diff + sigmabar2;
1528 break;
1529
1530 case 2 :
1531 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta
1532
1533 if (fType == "MC")
1534 hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta,
1535 binPixel, binPixel, "");
1536 else if (fType == "ON")
1537 hsig = fHSigmaPixThetaOFF->ProjectionZ("", binTheta, binTheta,
1538 binPixel, binPixel, "");
1539 else
1540 hsig = fHSigmaPixThetaON->ProjectionZ("", binTheta, binTheta,
1541 binPixel, binPixel, "");
1542
1543 if ( hsig->Integral() == 0 )
1544 {
1545 *fLog << warn << "Process(); projection for Theta bin "
1546 << binTheta << " and pixel bin " << binPixel
1547 << " has no entries; aborting " << endl;
1548 delete hsig;
1549
1550 rc = 6;
1551 fErrors[rc]++;
1552 return kCONTINUE;
1553 }
1554
1555 count = 0;
1556 ok = kFALSE;
1557 for (Int_t m=0; m<fIter; m++)
1558 {
1559 count += 1;
1560
1561 sig_phot = hsig->GetRandom();
1562 sig = sig_phot * fPEperPhoton;
1563
1564 sigma2 = sig*sig/ratioArea;
1565 // the following condition ensures that elNoise2Pix > 0.0
1566
1567 if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 )
1568 {
1569 ok = kTRUE;
1570 break;
1571 }
1572 }
1573 if (!ok)
1574 {
1575 *fLog << all << "theta, j, count, sigmabar, sig = " << theta
1576 << ", "
1577 << j << ", " << count << ", " << sigmabar << ", "
1578 << sig << endl;
1579 sigma2 = lambdabar*F2excess + oldsigma2/ratioArea;
1580
1581 rw = 1;
1582 fWarnings[rw]++;
1583 }
1584 else
1585 {
1586 rw = 0;
1587 fWarnings[rw]++;
1588 }
1589
1590 delete hsig;
1591 break;
1592 }
1593
1594 //---------------------------------
1595 // get the additional sigma^2 for this pixel (due to the padding)
1596
1597 addSig2 = sigma2*ratioArea - oldsigma2;
1598 addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
1599
1600 //---------------------------------
1601 // get the additional electronic noise for this pixel
1602
1603 elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;
1604
1605
1606 //---------------------------------
1607 // throw actual number of additional NSB photo electrons (NSB)
1608 // and its RMS (sigmaNSB)
1609
1610 Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
1611 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
1612 Double_t sigmaNSB0;
1613
1614 if (arg >= 0.0)
1615 {
1616 sigmaNSB0 = sqrt( arg );
1617 }
1618 else
1619 {
1620 sigmaNSB0 = 0.0000001;
1621 if (arg < -1.e-10)
1622 {
1623 *fLog << err << "Process(); argument of sqrt < 0 : "
1624 << arg << endl;
1625 }
1626 }
1627
1628
1629 //---------------------------------
1630 // smear NSB0 according to sigmaNSB0
1631 // and subtract lambdabar because of AC coupling
1632
1633 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
1634 Double_t NSB_phot = NSB / fPEperPhoton;
1635
1636 //---------------------------------
1637
1638 // add additional NSB to the number of photons
1639 Double_t oldphotons_phot = pix.GetNumPhotons();
1640 Double_t oldphotons = oldphotons_phot * fPEperPhoton;
1641 Double_t newphotons = oldphotons + NSB;
1642 Double_t newphotons_phot = newphotons / fPEperPhoton;;
1643 pix.SetNumPhotons( newphotons_phot );
1644
1645
1646 fHNSB->Fill( NSB_phot/ratioArea );
1647 fHPhotons->Fill( oldphotons_phot/ratioArea,
1648 newphotons_phot/ratioArea );
1649
1650
1651 // error: add sigma of padded noise quadratically
1652 Double_t olderror_phot = pix.GetErrorPhot();
1653 Double_t newerror_phot =
1654 sqrt( olderror_phot*olderror_phot + addSig2_phot );
1655 pix.SetErrorPhot( newerror_phot );
1656
1657
1658 Double_t newsigma = sqrt( oldsigma2 + addSig2 );
1659 Double_t newsigma_phot = newsigma / fPEperPhoton;;
1660 ppix.SetPedestalRms( newsigma_phot );
1661
1662 fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot );
1663 }
1664 //---------- end of loop over pixels ---------------------------------
1665
1666 // Calculate sigmabar again and crosscheck
1667
1668 fSigmabar->Calc(*fCam, *fPed, *fEvt);
1669 *fLog << all << "after padding : " << endl;
1670 fSigmabar->Print("");
1671
1672 rc = 0;
1673 fErrors[rc]++;
1674 return kTRUE;
1675 //******************************************************************
1676}
1677
1678// --------------------------------------------------------------------------
1679//
1680//
1681Int_t MPadONOFF::PostProcess()
1682{
1683 if (GetNumExecutions() != 0)
1684 {
1685
1686 *fLog << inf << endl;
1687 *fLog << GetDescriptor() << " execution statistics:" << endl;
1688 *fLog << dec << setfill(' ');
1689
1690 int temp;
1691 if (fWarnings[0] != 0.0)
1692 temp = (int)(fWarnings[1]*100/fWarnings[0]);
1693 else
1694 temp = 100;
1695
1696 *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3)
1697 << temp
1698 << "%) Warning : iteration to find acceptable sigma failed"
1699 << ", fIter = " << fIter << endl;
1700
1701 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
1702 << (int)(fErrors[1]*100/GetNumExecutions())
1703 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
1704
1705 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
1706 << (int)(fErrors[2]*100/GetNumExecutions())
1707 << "%) Evts skipped due to: Zenith angle out of range" << endl;
1708
1709 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
1710 << (int)(fErrors[3]*100/GetNumExecutions())
1711 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
1712
1713 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
1714 << (int)(fErrors[4]*100/GetNumExecutions())
1715 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
1716
1717 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
1718 << (int)(fErrors[5]*100/GetNumExecutions())
1719 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
1720
1721 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
1722 << (int)(fErrors[6]*100/GetNumExecutions())
1723 << "%) Evts skipped due to: No data for generating Sigma" << endl;
1724
1725 *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3)
1726 << (int)(fErrors[7]*100/GetNumExecutions())
1727 << "%) Evts skipped due to: No data for generating Blind pixels" << endl;
1728
1729 *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3)
1730 << (int)(fErrors[8]*100/GetNumExecutions())
1731 << "%) Evts didn't have to be padded" << endl;
1732
1733 *fLog << " " << fErrors[0] << " ("
1734 << (int)(fErrors[0]*100/GetNumExecutions())
1735 << "%) Evts were successfully padded!" << endl;
1736 *fLog << endl;
1737
1738 }
1739
1740 //---------------------------------------------------------------
1741 TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200));
1742 c.Divide(3, 4);
1743 gROOT->SetSelectedPad(NULL);
1744
1745 c.cd(1);
1746 fHNSB->SetDirectory(NULL);
1747 fHNSB->DrawCopy();
1748 fHNSB->SetBit(kCanDelete);
1749
1750 c.cd(2);
1751 fHSigmaPedestal->SetDirectory(NULL);
1752 fHSigmaPedestal->DrawCopy();
1753 fHSigmaPedestal->SetBit(kCanDelete);
1754
1755 c.cd(3);
1756 fHPhotons->SetDirectory(NULL);
1757 fHPhotons->DrawCopy();
1758 fHPhotons->SetBit(kCanDelete);
1759
1760 //--------------------------------------------------------------------
1761
1762
1763 c.cd(4);
1764 fHSigmaTheta->SetDirectory(NULL);
1765 fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta");
1766 fHSigmaTheta->DrawCopy();
1767 fHSigmaTheta->SetBit(kCanDelete);
1768
1769
1770 //--------------------------------------------------------------------
1771
1772
1773 c.cd(7);
1774 fHBlindPixNTheta->SetDirectory(NULL);
1775 fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta");
1776 fHBlindPixNTheta->DrawCopy();
1777 fHBlindPixNTheta->SetBit(kCanDelete);
1778
1779 //--------------------------------------------------------------------
1780
1781
1782 c.cd(10);
1783 fHBlindPixIdTheta->SetDirectory(NULL);
1784 fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta");
1785 fHBlindPixIdTheta->DrawCopy();
1786 fHBlindPixIdTheta->SetBit(kCanDelete);
1787
1788
1789 //--------------------------------------------------------------------
1790 // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2
1791
1792 c.cd(5);
1793 TH2D *l1;
1794 l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
1795 l1->SetDirectory(NULL);
1796 l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
1797 l1->SetXTitle("\\Theta [\\circ]");
1798 l1->SetYTitle("Sigma^2-Sigmabar^2");
1799
1800 l1->DrawCopy("box");
1801 l1->SetBit(kCanDelete);;
1802
1803 c.cd(8);
1804 TH2D *l2;
1805 l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
1806 l2->SetDirectory(NULL);
1807 l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
1808 l2->SetXTitle("pixel");
1809 l2->SetYTitle("Sigma^2-Sigmabar^2");
1810
1811 l2->DrawCopy("box");
1812 l2->SetBit(kCanDelete);;
1813
1814 //--------------------------------------------------------------------
1815 // draw the 3D histogram (target): Theta, pixel, Sigma
1816
1817 c.cd(6);
1818 TH2D *k1;
1819 k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
1820 k1->SetDirectory(NULL);
1821 k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)");
1822 k1->SetXTitle("\\Theta [\\circ]");
1823 k1->SetYTitle("Sigma");
1824
1825 k1->DrawCopy("box");
1826 k1->SetBit(kCanDelete);
1827
1828 c.cd(9);
1829 TH2D *k2;
1830 k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
1831 k2->SetDirectory(NULL);
1832 k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)");
1833 k2->SetXTitle("pixel");
1834 k2->SetYTitle("Sigma");
1835
1836 k2->DrawCopy("box");
1837 k2->SetBit(kCanDelete);;
1838
1839
1840 //--------------------------------------------------------------------
1841
1842 c.cd(11);
1843 TH2D *m1;
1844 m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
1845 m1->SetDirectory(NULL);
1846 m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
1847 m1->SetXTitle("Sigmabar-old");
1848 m1->SetYTitle("Sigmabar-new");
1849
1850 m1->DrawCopy("box");
1851 m1->SetBit(kCanDelete);;
1852
1853 c.cd(12);
1854 TH2D *m2;
1855 m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
1856 m2->SetDirectory(NULL);
1857 m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
1858 m2->SetXTitle("Sigmabar-old");
1859 m2->SetYTitle("Sigmabar-new");
1860
1861 m2->DrawCopy("box");
1862 m2->SetBit(kCanDelete);;
1863
1864
1865 return kTRUE;
1866}
1867
1868// --------------------------------------------------------------------------
1869
1870
1871
1872
1873
1874
Note: See TracBrowser for help on using the repository browser.