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

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