source: trunk/MagicSoft/Mars/manalysis/MCT1PadONOFF.cc@ 2194

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