source: trunk/MagicSoft/Mars/manalysisct1/MCT1PadONOFF.cc@ 8139

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