source: trunk/MagicSoft/Mars/manalysis/MPad.cc@ 3907

Last change on this file since 3907 was 2798, checked in by wittek, 22 years ago
*** empty log message ***
File size: 69.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// MPad
28//
29// The aim of the padding is to make the noise level (NSB + electronic noise)
30// equal for the MC, ON and OFF data samples
31//
32// The target noise level is determined by 'merging' the (sigmabar vs. Theta)
33// distributions of MC, ON and OFF data.
34//
35// The prescription on which fraction of events has to padded from
36// bin k of sigmabar to bin j is stored in the matrices
37// fHgMC, fHgON and fHgOFF.
38//
39// By the padding the number of photons, its error and the pedestal sigmas
40// are altered. On average, the number of photons added is zero.
41//
42// The formulas used can be found in Thomas Schweizer's Thesis,
43// Section 2.2.1
44//
45// The padding is done as follows :
46//
47// Increase the sigmabar according to the matrices fHg.... Then generate
48// a pedestal sigma for each pixel using the 3D-histogram Theta, pixel no.,
49// Sigma^2-Sigmabar^2 (fHDiffPixTheta).
50//
51// The padding has to be done before the image cleaning because the
52// image cleaning depends on the pedestal sigmas.
53//
54/////////////////////////////////////////////////////////////////////////////
55#include "MPad.h"
56
57#include <math.h>
58#include <stdio.h>
59
60#include <TH1.h>
61#include <TH2.h>
62#include <TH3.h>
63#include <TRandom.h>
64#include <TCanvas.h>
65#include <TFile.h>
66
67#include "MBinning.h"
68#include "MSigmabar.h"
69#include "MMcEvt.hxx"
70#include "MLog.h"
71#include "MLogManip.h"
72#include "MParList.h"
73#include "MGeomCam.h"
74
75#include "MCerPhotPix.h"
76#include "MCerPhotEvt.h"
77
78#include "MPedPhotCam.h"
79#include "MPedPhotPix.h"
80
81#include "MBlindPixels.h"
82
83#include "MRead.h"
84#include "MFilterList.h"
85#include "MTaskList.h"
86#include "MBlindPixelCalc.h"
87#include "MHBlindPixels.h"
88#include "MFillH.h"
89#include "MHSigmaTheta.h"
90#include "MEvtLoop.h"
91
92ClassImp(MPad);
93
94using namespace std;
95
96// --------------------------------------------------------------------------
97//
98// Default constructor.
99//
100MPad::MPad(const char *name, const char *title)
101{
102 fName = name ? name : "MPad";
103 fTitle = title ? title : "Task for the padding";
104
105 fType = "";
106
107 fHSigmaThetaMC = NULL;
108 fHSigmaThetaON = NULL;
109 fHSigmaThetaOFF = NULL;
110
111 fHDiffPixThetaMC = NULL;
112 fHDiffPixThetaON = NULL;
113 fHDiffPixThetaOFF = NULL;
114
115 fHSigmaPixThetaMC = NULL;
116 fHSigmaPixThetaON = NULL;
117 fHSigmaPixThetaOFF = NULL;
118
119 fHgMC = NULL;
120 fHgON = NULL;
121 fHgOFF = NULL;
122
123 fHBlindPixIdThetaMC = NULL;
124 fHBlindPixIdThetaON = NULL;
125 fHBlindPixIdThetaOFF = NULL;
126
127 fHBlindPixNThetaMC = NULL;
128 fHBlindPixNThetaON = NULL;
129 fHBlindPixNThetaOFF = NULL;
130
131 fHSigmaPedestal = NULL;
132 fHPhotons = NULL;
133 fHNSB = NULL;
134}
135
136// --------------------------------------------------------------------------
137//
138// Destructor.
139//
140MPad::~MPad()
141{
142 delete fHgMC;
143 delete fHgON;
144 delete fHgOFF;
145
146 delete fHSigmaTheta;
147 delete fHSigmaPixTheta;
148 delete fHDiffPixTheta;
149 delete fHBlindPixNTheta;
150 delete fHBlindPixIdTheta;
151
152 delete fHSigmaPedestal;
153 delete fHPhotons;
154 delete fHNSB;
155
156 delete fInfile;
157}
158// --------------------------------------------------------------------------
159//
160// Merge the sigmabar distributions of 2 samples (samples A and B)
161//
162// input : the original distributions for samples A and B (hista, histb)
163//
164// output : the prescription which fraction of events should be padded
165// from the sigmabar bin k to bin j (fHgMC, fHgON, fHgOFF)
166//
167
168Bool_t MPad::Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap,
169 TH2D *fHga, TH2D *fHgb, Int_t nbinssig )
170{
171 // hista is the normalized 1D histogram of sigmabar for sample A
172 // histb is the normalized 1D histogram of sigmabar for sample B
173 // histc is the difference A-B
174
175 // at the beginning, histap is a copy of hista
176 // at the end, it will be the 1D histogram for sample A after padding
177
178 // at the beginning, histbp is a copy of histb
179 // at the end, it will be the 1D histogram for sample B after padding
180
181 // at the beginning, histcp is a copy of histc
182 // at the end, it should be the difference histap-histbp,
183 // which should be zero
184
185 // fHga[k][j] tells which fraction of events from sample A should be padded
186 // from sigma_k to sigma_j
187
188
189 // fHgb[k][j] tells which fraction of events from sample B should be padded
190 // from sigma_k to sigma_j
191
192 Double_t eps = 1.e-10;
193
194 TH1D *histbp = new TH1D( (const TH1D&)*histb );
195
196 TH1D *histc = new TH1D( (const TH1D&)*hista );
197 histc->Add(histb, -1.0);
198
199 TH1D *histcp = new TH1D( (const TH1D&)*histc );
200
201
202 //-------- start j loop ------------------------------------------------
203 // loop over bins in histc,
204 // starting from the end (i.e. at the largest sigmabar)
205 Double_t v, s, w, t, x, u, a, b, arg;
206
207 for (Int_t j=nbinssig; j >= 1; j--)
208 {
209 v = histcp->GetBinContent(j);
210 if ( fabs(v) < eps ) continue;
211 if (v >= 0.0)
212 s = 1.0;
213 else
214 s = -1.0;
215
216 //................ start k loop ......................................
217 // look for a bin k which may compensate the content of bin j
218 for (Int_t k=j-1; k >= 1; k--)
219 {
220 w = histcp->GetBinContent(k);
221 if ( fabs(w) < eps ) continue;
222 if (w >= 0.0)
223 t = 1.0;
224 else
225 t = -1.0;
226
227 // if s==t bin k cannot compensate, go to next k bin
228 if (s == t) continue;
229
230 x = v + w;
231 if (x >= 0.0)
232 u = 1.0;
233 else
234 u = -1.0;
235
236 // if u==s bin k will partly compensate : pad the whole fraction
237 // w from bin k to bin j
238 if (u == s)
239 arg = -w;
240
241 // if u!=s bin k would overcompensate : pad only the fraction
242 // v from bin k to bin j
243 else
244 arg = v;
245
246 if (arg <=0.0)
247 {
248 fHga->SetBinContent(k, j, -arg);
249 fHgb->SetBinContent(k, j, 0.0);
250 }
251 else
252 {
253 fHga->SetBinContent(k, j, 0.0);
254 fHgb->SetBinContent(k, j, arg);
255 }
256
257 *fLog << all << "k, j, arg = " << k << ", " << j
258 << ", " << arg << endl;
259
260 //......................................
261 // this is for checking the procedure
262 if (arg < 0.0)
263 {
264 a = histap->GetBinContent(k);
265 histap->SetBinContent(k, a+arg);
266 a = histap->GetBinContent(j);
267 histap->SetBinContent(j, a-arg);
268 }
269 else
270 {
271 b = histbp->GetBinContent(k);
272 histbp->SetBinContent(k, b-arg);
273 b = histbp->GetBinContent(j);
274 histbp->SetBinContent(j, b+arg);
275 }
276 //......................................
277
278 if (u == s)
279 {
280 histcp->SetBinContent(k, 0.0);
281 histcp->SetBinContent(j, x);
282
283 v = histcp->GetBinContent(j);
284 if ( fabs(v) < eps ) break;
285
286 // bin j was only partly compensated : go to next k bin
287 continue;
288 }
289 else
290 {
291 histcp->SetBinContent(k, x);
292 histcp->SetBinContent(j, 0.0);
293
294 // bin j was completely compensated : go to next j bin
295 break;
296 }
297
298 }
299 //................ end k loop ......................................
300 }
301 //-------- end j loop ------------------------------------------------
302
303 // check results
304 for (Int_t j=1; j<=nbinssig; j++)
305 {
306 Double_t a = histap->GetBinContent(j);
307 Double_t b = histbp->GetBinContent(j);
308 Double_t c = histcp->GetBinContent(j);
309
310 if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps )
311 *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = "
312 << j << ", " << a << ", " << b << ", " << c << endl;
313 }
314
315 //---------------------------------------------------------------
316 TCanvas &c = *(MH::MakeDefCanvas("Merging", "", 600, 600));
317 c.Divide(2, 2);
318 gROOT->SetSelectedPad(NULL);
319
320 c.cd(1);
321 hista->SetDirectory(NULL);
322 hista->DrawCopy();
323 hista->SetBit(kCanDelete);
324
325 c.cd(2);
326 histb->SetDirectory(NULL);
327 histb->DrawCopy();
328 histb->SetBit(kCanDelete);
329
330 c.cd(3);
331 histap->SetDirectory(NULL);
332 histap->DrawCopy();
333 histap->SetBit(kCanDelete);
334
335 //--------------------------------------------------------------------
336
337 delete histc;
338 delete histbp;
339 delete histcp;
340
341 return kTRUE;
342}
343
344
345// --------------------------------------------------------------------------
346//
347// Merge the distributions
348// fHSigmaTheta 2D-histogram (Theta, sigmabar)
349//
350// ===> of ON, OFF and MC data <==============
351//
352// and define the matrices fHgMC, fHgON and fHgOFF
353//
354// These matrices tell which fraction of events should be padded
355// from bin k of sigmabar to bin j
356//
357
358Bool_t MPad::MergeONOFFMC(
359 TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc,
360 TH2D *blindidthmc, TH2D *blindnthmc,
361 TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon,
362 TH2D *blindidthon, TH2D *blindnthon,
363 TH2D *sigthoff, TH3D *diffpixthoff, TH3D *sigmapixthoff,
364 TH2D *blindidthoff, TH2D *blindnthoff)
365{
366 *fLog << all << "----------------------------------------------------------------------------------" << endl;
367 *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding"
368 << endl;
369
370 fHSigmaTheta = new TH2D( (TH2D&)*sigthon );
371 fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");
372
373 fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon );
374 fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");
375
376 fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon );
377 fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");
378
379 fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon );
380 fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
381
382 fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon );
383 fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
384
385 //--------------------------
386 fHSigmaThetaMC = sigthmc;
387 fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
388 fHSigmaThetaON = sigthon;
389 fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
390 fHSigmaThetaOFF = sigthoff;
391 fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)");
392
393 //--------------------------
394 fHBlindPixNThetaMC = blindnthmc;
395 fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
396 fHBlindPixNThetaON = blindnthon;
397 fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
398 fHBlindPixNThetaOFF = blindnthoff;
399 fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)");
400
401 //--------------------------
402 fHBlindPixIdThetaMC = blindidthmc;
403 fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
404 fHBlindPixIdThetaON = blindidthon;
405 fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)");
406 fHBlindPixIdThetaOFF = blindidthoff;
407 fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)");
408
409 //--------------------------
410 fHDiffPixThetaMC = diffpixthmc;
411 fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)");
412 fHDiffPixThetaON = diffpixthon;
413 fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)");
414 fHDiffPixThetaOFF = diffpixthoff;
415 fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)");
416
417 //--------------------------
418 fHSigmaPixThetaMC = sigmapixthmc;
419 fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)");
420 fHSigmaPixThetaON = sigmapixthon;
421 fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)");
422 fHSigmaPixThetaOFF = sigmapixthoff;
423 fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)");
424 //--------------------------
425
426
427 // get binning for fHgON, fHgOFF and fHgMC from sigthon
428 // binning in Theta
429 TAxis *ax = sigthon->GetXaxis();
430 Int_t nbinstheta = ax->GetNbins();
431 TArrayD edgesx;
432 edgesx.Set(nbinstheta+1);
433 for (Int_t i=0; i<=nbinstheta; i++)
434 {
435 edgesx[i] = ax->GetBinLowEdge(i+1);
436 *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i]
437 << endl;
438 }
439 MBinning binth;
440 binth.SetEdges(edgesx);
441
442 // binning in sigmabar
443 TAxis *ay = sigthon->GetYaxis();
444 Int_t nbinssig = ay->GetNbins();
445 TArrayD edgesy;
446 edgesy.Set(nbinssig+1);
447 for (Int_t i=0; i<=nbinssig; i++)
448 {
449 edgesy[i] = ay->GetBinLowEdge(i+1);
450 *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i]
451 << endl;
452 }
453 MBinning binsg;
454 binsg.SetEdges(edgesy);
455
456
457 fHgMC = new TH3D;
458 MH::SetBinning(fHgMC, &binth, &binsg, &binsg);
459 fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC");
460
461 fHgON = new TH3D;
462 MH::SetBinning(fHgON, &binth, &binsg, &binsg);
463 fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");
464
465 fHgOFF = new TH3D;
466 MH::SetBinning(fHgOFF, &binth, &binsg, &binsg);
467 fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF");
468
469 //............ loop over Theta bins ...........................
470
471
472
473 *fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
474 for (Int_t l=1; l<=nbinstheta; l++)
475 {
476 TH2D * fHga = new TH2D;
477 MH::SetBinning(fHga, &binsg, &binsg);
478 TH2D * fHgb = new TH2D;
479 MH::SetBinning(fHgb, &binsg, &binsg);
480
481 //-------------------------------------------
482 // merge ON and OFF distributions
483 // input : hista is the normalized 1D distr. of sigmabar for ON data
484 // histb is the normalized 1D distr. of sigmabar for OFF data
485 // output : histap will be the merged distribution (ON-OFF)
486 // fHga(k,j) will tell which fraction of ON events should be
487 // padded from bin k to bin j
488 // fHgb(k,j) will tell which fraction of OFF events should be
489 // padded from bin k to bin j
490
491 TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, "");
492 Stat_t suma = hista->Integral();
493 hista->Scale(1./suma);
494
495 TH1D *histap = new TH1D( (const TH1D&)*hista );
496
497 TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, "");
498 Stat_t sumb = histb->Integral();
499 histb->Scale(1./sumb);
500
501 Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
502 delete hista;
503 delete histb;
504
505 // fill fHgON and fHgOFF
506 for (Int_t k=1; k<=nbinssig; k++)
507 for (Int_t j=1; j<=nbinssig; j++)
508 {
509 Double_t a = fHga->GetBinContent(k,j);
510 fHga->SetBinContent (k, j, 0.0);
511 fHgON->SetBinContent(l, k, j, a);
512
513 Double_t b = fHgb->GetBinContent(k,j);
514 fHgb->SetBinContent (k, j, 0.0);
515 fHgOFF->SetBinContent(l, k, j, b);
516 }
517
518 //-------------------------------------------
519 // now merge ON-OFF and MC distributions
520 // input : hista is the result of the merging of ON and OFF distributions
521 // histb is the normalized 1D distr. of sigmabar for MC data
522 // output : histap will be the merged distribution (target distribution)
523 // fHga(k,j) will tell which fraction of ON, OFF events should be
524 // padded from bin k to bin j
525 // fHgb(k,j) will tell which fraction of MC events should be
526 // padded from bin k to bin j
527
528 hista = new TH1D( (const TH1D&)*histap);
529
530 histb = sigthmc->ProjectionY("sigmc_y", l, l, "");
531 sumb = histb->Integral();
532 histb->Scale(1./sumb);
533
534 Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
535 delete hista;
536 delete histb;
537
538 // update fHgON and fHgOFF and fill fHgMC
539 for (Int_t k=1; k<=nbinssig; k++)
540 for (Int_t j=1; j<=nbinssig; j++)
541 {
542 Double_t a = fHga->GetBinContent (k,j);
543
544 Double_t aON = fHgON->GetBinContent(l,k,j);
545 fHgON->SetBinContent(l, k, j, aON+a);
546
547 Double_t aOFF = fHgOFF->GetBinContent(l,k,j);
548 fHgOFF->SetBinContent(l, k, j, aOFF+a);
549
550 Double_t b = fHgb->GetBinContent(k,j);
551 fHgMC->SetBinContent(l, k, j, b);
552 }
553
554 // fill target distribution fHSigmaTheta
555 // (only for plotting, it will not be used in the padding)
556 for (Int_t j=1; j<=nbinssig; j++)
557 {
558 Double_t a = histap->GetBinContent(j);
559 fHSigmaTheta->SetBinContent(l, j, a);
560 }
561
562 delete fHga;
563 delete fHgb;
564
565 delete histap;
566 }
567 //............ end of loop over Theta bins ....................
568
569
570 // Define the target distributions
571 // fHDiffPixTheta, fHSigmaPixTheta
572 // (they are calculated as
573 // averages of the ON and OFF distributions)
574 // fHBlindPixNTheta, fHBlindPixIdTheta
575 // (they are calculated as
576 // the sum of the ON and OFF distributions)
577 // (fHDiffPixTheta will be used in the padding of MC events)
578
579 //............ start of new loop over Theta bins ....................
580 for (Int_t j=1; j<=nbinstheta; j++)
581 {
582 // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC
583
584 Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, "");
585 Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, "");
586 Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, "");
587
588 Double_t normON;
589 Double_t normOFF;
590
591 // set ranges for 2D-projections of 3D-histograms
592 TAxis *ax = diffpixthon->GetXaxis();
593 ax->SetRange(j, j);
594 ax = diffpixthoff->GetXaxis();
595 ax->SetRange(j, j);
596
597 TH2D *hist;
598 TH2D *histOFF;
599
600 TH1D *hist1;
601 TH1D *hist1OFF;
602
603 // weights for ON and OFF distrubtions when
604 // calculating the weighted average
605 Double_t facON = 0.5 * (1. - fracON + fracOFF);
606 Double_t facOFF = 0.5 * (1. + fracON - fracOFF);
607
608 *fLog << all << "Theta bin j : fracON, fracOFF, facON, facOFF = "
609 << j << ", " << fracON << ", " << fracOFF << ", "
610 << fracMC << ", " << facON << ", " << facOFF << endl;
611
612 //------------------------------------------------------------------
613 // define target distribution 'sigma^2-sigmavar^2 vs. pixel number'
614 ay = diffpixthon->GetYaxis();
615 Int_t nbinspixel = ay->GetNbins();
616
617 TAxis *az = diffpixthon->GetZaxis();
618 Int_t nbinsdiff = az->GetNbins();
619
620 hist = (TH2D*)diffpixthon->Project3D("zy");
621 hist->SetName("dummy");
622 histOFF = (TH2D*)diffpixthoff->Project3D("zy");
623
624 normON = hist->Integral();
625 normOFF = histOFF->Integral();
626 if (normON != 0.0) hist->Scale(1.0/normON);
627 if (normOFF != 0.0) histOFF->Scale(1.0/normOFF);
628
629 // weighted average of ON and OFF distributions
630 hist->Scale(facON);
631 hist->Add(histOFF, facOFF);
632
633 for (Int_t k=1; k<=nbinspixel; k++)
634 for (Int_t l=1; l<=nbinsdiff; l++)
635 {
636 Double_t cont = hist->GetBinContent(k,l);
637 fHDiffPixTheta->SetBinContent(j, k, l, cont);
638 }
639
640 delete hist;
641 delete histOFF;
642
643 //------------------------------------------------------------------
644 // define target distribution 'sigma vs. pixel number'
645 ay = sigmapixthon->GetYaxis();
646 nbinspixel = ay->GetNbins();
647
648 az = sigmapixthon->GetZaxis();
649 Int_t nbinssigma = az->GetNbins();
650
651 hist = (TH2D*)sigmapixthon->Project3D("zy");
652 hist->SetName("dummy");
653 histOFF = (TH2D*)sigmapixthoff->Project3D("zy");
654
655 normON = hist->Integral();
656 normOFF = histOFF->Integral();
657 if (normON != 0.0) hist->Scale(1.0/normON);
658 if (normOFF != 0.0) histOFF->Scale(1.0/normOFF);
659
660 // weighted average of ON and OFF distributions
661 hist->Scale(facON);
662 hist->Add(histOFF, facOFF);
663
664 for (Int_t k=1; k<=nbinspixel; k++)
665 for (Int_t l=1; l<=nbinssigma; l++)
666 {
667 Double_t cont = hist->GetBinContent(k,l);
668 fHSigmaPixTheta->SetBinContent(j, k, l, cont);
669 }
670
671 delete hist;
672 delete histOFF;
673
674
675 //------------------------------------------------------------------
676 // define target distribution 'number of blind pixels per event'
677 ay = blindnthon->GetYaxis();
678 Int_t nbinsn = ay->GetNbins();
679
680 hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, "");
681 hist1->SetName("dummy");
682 hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, "");
683
684 normON = hist1->Integral();
685 normOFF = hist1OFF->Integral();
686 if (normON != 0.0) hist1->Scale(1.0/normON);
687 if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
688
689 // sum of ON and OFF distributions
690 hist1->Add(hist1OFF, 1.0);
691 hist1->Scale( 1.0/hist1->Integral() );
692
693 for (Int_t k=1; k<=nbinsn; k++)
694 {
695 Double_t cont = hist1->GetBinContent(k);
696 fHBlindPixNTheta->SetBinContent(j, k, cont);
697 }
698
699 delete hist1;
700 delete hist1OFF;
701
702 //------------------------------------------------------------------
703 // define target distribution 'id of blind pixel'
704 ay = blindidthon->GetYaxis();
705 Int_t nbinsid = ay->GetNbins();
706
707 hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
708 hist1->SetName("dummy");
709 hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, "");
710
711 // divide by the number of events (from fHBlindPixNTheta)
712 if (normON != 0.0) hist1->Scale(1.0/normON);
713 if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);
714
715 // sum of ON and OFF distributions
716 hist1->Add(hist1OFF, 1.0);
717
718 for (Int_t k=1; k<=nbinsid; k++)
719 {
720 Double_t cont = hist1->GetBinContent(k);
721 fHBlindPixIdTheta->SetBinContent(j, k, cont);
722 }
723
724 delete hist1;
725 delete hist1OFF;
726 }
727 //............ end of new loop over Theta bins ....................
728
729 *fLog << all
730 << "MPad::MergeONOFFMC(); The target distributions for the padding have been created"
731 << endl;
732 *fLog << all << "----------------------------------------------------------"
733 << endl;
734 //--------------------------------------------
735
736 fHSigmaThetaMC->SetDirectory(NULL);
737 fHSigmaThetaON->SetDirectory(NULL);
738 fHSigmaThetaOFF->SetDirectory(NULL);
739
740 fHDiffPixThetaMC->SetDirectory(NULL);
741 fHDiffPixThetaON->SetDirectory(NULL);
742 fHDiffPixThetaOFF->SetDirectory(NULL);
743
744 fHSigmaPixThetaMC->SetDirectory(NULL);
745 fHSigmaPixThetaON->SetDirectory(NULL);
746 fHSigmaPixThetaOFF->SetDirectory(NULL);
747
748 fHBlindPixIdThetaMC->SetDirectory(NULL);
749 fHBlindPixIdThetaON->SetDirectory(NULL);
750 fHBlindPixIdThetaOFF->SetDirectory(NULL);
751
752 fHBlindPixNThetaMC->SetDirectory(NULL);
753 fHBlindPixNThetaON->SetDirectory(NULL);
754 fHBlindPixNThetaOFF->SetDirectory(NULL);
755
756 fHgMC->SetDirectory(NULL);
757 fHgON->SetDirectory(NULL);
758 fHgOFF->SetDirectory(NULL);
759
760 return kTRUE;
761}
762
763// --------------------------------------------------------------------------
764//
765// Merge the distributions
766// fHSigmaTheta 2D-histogram (Theta, sigmabar)
767//
768// ===> of ON and MC data <==============
769//
770// and define the matrices fHgMC and fHgON
771//
772// These matrices tell which fraction of events should be padded
773// from bin k of sigmabar to bin j
774//
775
776Bool_t MPad::MergeONMC(
777 TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc,
778 TH2D *blindidthmc, TH2D *blindnthmc,
779 TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon,
780 TH2D *blindidthon, TH2D *blindnthon)
781{
782 *fLog << all << "----------------------------------------------------------------------------------" << endl;
783 *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding"
784 << endl;
785
786 fHSigmaTheta = new TH2D( (TH2D&)*sigthon );
787 fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");
788
789 fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon );
790 fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");
791
792 fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon );
793 fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");
794
795 fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon );
796 fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");
797
798 fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon );
799 fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");
800
801 //--------------------------
802 fHSigmaThetaMC = sigthmc;
803 fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
804 fHSigmaThetaON = sigthon;
805 fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
806
807 //--------------------------
808 fHBlindPixNThetaMC = blindnthmc;
809 fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
810 fHBlindPixNThetaON = blindnthon;
811 fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
812
813 //--------------------------
814 fHBlindPixIdThetaMC = blindidthmc;
815 fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
816 fHBlindPixIdThetaON = blindidthon;
817 fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)");
818
819 //--------------------------
820 fHDiffPixThetaMC = diffpixthmc;
821 fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)");
822 fHDiffPixThetaON = diffpixthon;
823 fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)");
824
825 //--------------------------
826 fHSigmaPixThetaMC = sigmapixthmc;
827 fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)");
828 fHSigmaPixThetaON = sigmapixthon;
829 fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)");
830
831 //--------------------------
832
833
834 // get binning for fHgON and fHgMC from sigthon
835 // binning in Theta
836 TAxis *ax = sigthon->GetXaxis();
837 Int_t nbinstheta = ax->GetNbins();
838 TArrayD edgesx;
839 edgesx.Set(nbinstheta+1);
840 for (Int_t i=0; i<=nbinstheta; i++)
841 {
842 edgesx[i] = ax->GetBinLowEdge(i+1);
843 *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i]
844 << endl;
845 }
846 MBinning binth;
847 binth.SetEdges(edgesx);
848
849 // binning in sigmabar
850 TAxis *ay = sigthon->GetYaxis();
851 Int_t nbinssig = ay->GetNbins();
852 TArrayD edgesy;
853 edgesy.Set(nbinssig+1);
854 for (Int_t i=0; i<=nbinssig; i++)
855 {
856 edgesy[i] = ay->GetBinLowEdge(i+1);
857 *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i]
858 << endl;
859 }
860 MBinning binsg;
861 binsg.SetEdges(edgesy);
862
863
864 fHgMC = new TH3D;
865 MH::SetBinning(fHgMC, &binth, &binsg, &binsg);
866 fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC");
867
868 fHgON = new TH3D;
869 MH::SetBinning(fHgON, &binth, &binsg, &binsg);
870 fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");
871
872
873 //............ loop over Theta bins ...........................
874
875
876
877 *fLog << all << "MPad::MergeONMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
878 for (Int_t l=1; l<=nbinstheta; l++)
879 {
880 TH2D * fHga = new TH2D;
881 MH::SetBinning(fHga, &binsg, &binsg);
882 TH2D * fHgb = new TH2D;
883 MH::SetBinning(fHgb, &binsg, &binsg);
884
885 //-------------------------------------------
886 // merge ON and MC distributions
887 // input : hista is the normalized 1D distr. of sigmabar for ON data
888 // histb is the normalized 1D distr. of sigmabar for MC data
889 // output : histap will be the merged distribution (ON-MC)
890 // fHga(k,j) will tell which fraction of ON events should be
891 // padded from bin k to bin j
892 // fHgb(k,j) will tell which fraction of MC events should be
893 // padded from bin k to bin j
894
895 TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, "");
896 Stat_t suma = hista->Integral();
897 hista->Scale(1./suma);
898
899 TH1D *histap = new TH1D( (const TH1D&)*hista );
900
901 TH1D *histb = sigthmc->ProjectionY("sigmc_y", l, l, "");
902 Stat_t sumb = histb->Integral();
903 histb->Scale(1./sumb);
904
905 Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
906 delete hista;
907 delete histb;
908
909 // fill fHgON and fHgMC
910 for (Int_t k=1; k<=nbinssig; k++)
911 for (Int_t j=1; j<=nbinssig; j++)
912 {
913 Double_t a = fHga->GetBinContent(k,j);
914 fHga->SetBinContent (k, j, 0.0);
915 fHgON->SetBinContent(l, k, j, a);
916
917 Double_t b = fHgb->GetBinContent(k,j);
918 fHgb->SetBinContent (k, j, 0.0);
919 fHgMC->SetBinContent(l, k, j, b);
920 }
921
922
923 // fill target distribution fHSigmaTheta
924 // (only for plotting, it will not be used in the padding)
925 for (Int_t j=1; j<=nbinssig; j++)
926 {
927 Double_t a = histap->GetBinContent(j);
928 fHSigmaTheta->SetBinContent(l, j, a);
929 }
930
931 delete fHga;
932 delete fHgb;
933
934 delete histap;
935 }
936 //............ end of loop over Theta bins ....................
937
938
939 // Define the target distributions
940 // fHDiffPixTheta, fHSigmaPixTheta
941 // (they are calculated as
942 // averages of the ON and MCdistributions)
943 // fHBlindPixNTheta, fHBlindPixIdTheta
944 // (they are calculated as
945 // the sum of the ON and MC distributions)
946 // (fHDiffPixTheta will be used in the padding of MC events)
947
948 //............ start of new loop over Theta bins ....................
949 for (Int_t j=1; j<=nbinstheta; j++)
950 {
951 // fraction of ON/MC events to be padded : fracON, fracMC
952
953 Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, "");
954 Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, "");
955
956 Double_t normON;
957 Double_t normMC;
958
959 // set ranges for 2D-projections of 3D-histograms
960 TAxis *ax = diffpixthon->GetXaxis();
961 ax->SetRange(j, j);
962 ax = diffpixthmc->GetXaxis();
963 ax->SetRange(j, j);
964
965 TH2D *hist;
966 TH2D *histMC;
967
968 TH1D *hist1;
969 TH1D *hist1MC;
970
971 // weights for ON and MC distrubtions when
972 // calculating the weighted average
973 Double_t facON = 0.5 * (1. - fracON + fracMC);
974 Double_t facMC = 0.5 * (1. + fracON - fracMC);
975
976 *fLog << all << "Theta bin j : fracON, fracMC, facON, facMC = "
977 << j << ", " << fracON << ", " << fracMC << ", "
978 << facON << ", " << facMC << endl;
979
980 //------------------------------------------------------------------
981 // define target distribution 'sigma^2-sigmavar^2 vs. pixel number'
982 ay = diffpixthon->GetYaxis();
983 Int_t nbinspixel = ay->GetNbins();
984
985 TAxis *az = diffpixthon->GetZaxis();
986 Int_t nbinsdiff = az->GetNbins();
987
988 hist = (TH2D*)diffpixthon->Project3D("zy");
989 hist->SetName("dummy");
990 histMC = (TH2D*)diffpixthmc->Project3D("zy");
991
992 normON = hist->Integral();
993 normMC = histMC->Integral();
994 if (normON != 0.0) hist->Scale(1.0/normON);
995 if (normMC != 0.0) histMC->Scale(1.0/normMC);
996
997 // weighted average of ON and MC distributions
998 hist->Scale(facON);
999 hist->Add(histMC, facMC);
1000
1001 for (Int_t k=1; k<=nbinspixel; k++)
1002 for (Int_t l=1; l<=nbinsdiff; l++)
1003 {
1004 Double_t cont = hist->GetBinContent(k,l);
1005 fHDiffPixTheta->SetBinContent(j, k, l, cont);
1006 }
1007
1008 delete hist;
1009 delete histMC;
1010
1011 //------------------------------------------------------------------
1012 // define target distribution 'sigma vs. pixel number'
1013 ay = sigmapixthon->GetYaxis();
1014 nbinspixel = ay->GetNbins();
1015
1016 az = sigmapixthon->GetZaxis();
1017 Int_t nbinssigma = az->GetNbins();
1018
1019 hist = (TH2D*)sigmapixthon->Project3D("zy");
1020 hist->SetName("dummy");
1021 histMC = (TH2D*)sigmapixthmc->Project3D("zy");
1022
1023 normON = hist->Integral();
1024 normMC = histMC->Integral();
1025 if (normON != 0.0) hist->Scale(1.0/normON);
1026 if (normMC != 0.0) histMC->Scale(1.0/normMC);
1027
1028 // weighted average of ON and MC distributions
1029 hist->Scale(facON);
1030 hist->Add(histMC, facMC);
1031
1032 for (Int_t k=1; k<=nbinspixel; k++)
1033 for (Int_t l=1; l<=nbinssigma; l++)
1034 {
1035 Double_t cont = hist->GetBinContent(k,l);
1036 fHSigmaPixTheta->SetBinContent(j, k, l, cont);
1037 }
1038
1039 delete hist;
1040 delete histMC;
1041
1042
1043 //------------------------------------------------------------------
1044 // define target distribution 'number of blind pixels per event'
1045 ay = blindnthon->GetYaxis();
1046 Int_t nbinsn = ay->GetNbins();
1047
1048 hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, "");
1049 hist1->SetName("dummy");
1050 hist1MC = fHBlindPixNThetaMC->ProjectionY("", j, j, "");
1051
1052 normON = hist1->Integral();
1053 normMC = hist1MC->Integral();
1054 if (normON != 0.0) hist1->Scale(1.0/normON);
1055 if (normMC != 0.0) hist1MC->Scale(1.0/normMC);
1056
1057 // sum of ON and MC distributions
1058 hist1->Add(hist1MC, 1.0);
1059 hist1->Scale( 1.0/hist1->Integral() );
1060
1061 for (Int_t k=1; k<=nbinsn; k++)
1062 {
1063 Double_t cont = hist1->GetBinContent(k);
1064 fHBlindPixNTheta->SetBinContent(j, k, cont);
1065 }
1066
1067 delete hist1;
1068 delete hist1MC;
1069
1070 //------------------------------------------------------------------
1071 // define target distribution 'id of blind pixel'
1072 ay = blindidthon->GetYaxis();
1073 Int_t nbinsid = ay->GetNbins();
1074
1075 hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
1076 hist1->SetName("dummy");
1077 hist1MC = fHBlindPixIdThetaMC->ProjectionY("", j, j, "");
1078
1079 // divide by the number of events (from fHBlindPixNTheta)
1080 if (normON != 0.0) hist1->Scale(1.0/normON);
1081 if (normMC != 0.0) hist1MC->Scale(1.0/normMC);
1082
1083 // sum of ON and MC distributions
1084 hist1->Add(hist1MC, 1.0);
1085
1086 for (Int_t k=1; k<=nbinsid; k++)
1087 {
1088 Double_t cont = hist1->GetBinContent(k);
1089 fHBlindPixIdTheta->SetBinContent(j, k, cont);
1090 }
1091
1092 delete hist1;
1093 delete hist1MC;
1094 }
1095 //............ end of new loop over Theta bins ....................
1096
1097 *fLog << all
1098 << "MPad::MergeONMC(); The target distributions for the padding have been created"
1099 << endl;
1100 *fLog << all << "----------------------------------------------------------"
1101 << endl;
1102 //--------------------------------------------
1103
1104 fHSigmaThetaMC->SetDirectory(NULL);
1105 fHSigmaThetaON->SetDirectory(NULL);
1106
1107 fHDiffPixThetaMC->SetDirectory(NULL);
1108 fHDiffPixThetaON->SetDirectory(NULL);
1109
1110 fHSigmaPixThetaMC->SetDirectory(NULL);
1111 fHSigmaPixThetaON->SetDirectory(NULL);
1112
1113 fHBlindPixIdThetaMC->SetDirectory(NULL);
1114 fHBlindPixIdThetaON->SetDirectory(NULL);
1115
1116 fHBlindPixNThetaMC->SetDirectory(NULL);
1117 fHBlindPixNThetaON->SetDirectory(NULL);
1118
1119 fHgMC->SetDirectory(NULL);
1120 fHgON->SetDirectory(NULL);
1121
1122 return kTRUE;
1123}
1124
1125
1126// --------------------------------------------------------------------------
1127//
1128// Read padding distributions from a file
1129//
1130//
1131Bool_t MPad::ReadPaddingDist(const char* namefilein)
1132{
1133 *fLog << all << "MPad : Read padding histograms from file "
1134 << namefilein << endl;
1135
1136 fInfile = new TFile(namefilein);
1137 fInfile->ls();
1138
1139 //------------------------------------
1140 fHSigmaThetaMC =
1141 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC");
1142 if (!fHSigmaThetaMC)
1143 {
1144 *fLog << all
1145 << "MPad : Object '2D-ThetaSigmabarMC' not found on root file"
1146 << endl;
1147 return kFALSE;
1148 }
1149 *fLog << all
1150 << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl;
1151
1152 fHSigmaThetaON =
1153 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON");
1154 if (!fHSigmaThetaON)
1155 {
1156 *fLog << all
1157 << "MPad : Object '2D-ThetaSigmabarON' not found on root file"
1158 << endl;
1159 return kFALSE;
1160 }
1161 *fLog << all
1162 << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl;
1163
1164 fHSigmaThetaOFF =
1165 (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF");
1166 if (!fHSigmaThetaOFF)
1167 {
1168 *fLog << all
1169 << "MPad : Object '2D-ThetaSigmabarOFF' not found on root file"
1170 << endl;
1171 return kFALSE;
1172 }
1173 *fLog << all
1174 << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl;
1175
1176 //------------------------------------
1177 fHDiffPixTheta =
1178 (TH3D*) gROOT->FindObject("3D-ThetaPixDiff");
1179 if (!fHDiffPixTheta)
1180 {
1181 *fLog << all
1182 << "MPad : Object '3D-ThetaPixDiff' not found on root file"
1183 << endl;
1184 return kFALSE;
1185 }
1186 *fLog << all
1187 << "MPad : Object '3D-ThetaPixDiff' was read in" << endl;
1188
1189 fHDiffPixThetaMC =
1190 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC");
1191 if (!fHDiffPixThetaMC)
1192 {
1193 *fLog << all
1194 << "MPad : Object '3D-ThetaPixDiffMC' not found on root file"
1195 << endl;
1196 return kFALSE;
1197 }
1198 *fLog << all
1199 << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl;
1200
1201 fHDiffPixThetaON =
1202 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON");
1203 if (!fHDiffPixThetaON)
1204 {
1205 *fLog << all
1206 << "MPad : Object '3D-ThetaPixDiffON' not found on root file"
1207 << endl;
1208 return kFALSE;
1209 }
1210 *fLog << all
1211 << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl;
1212
1213 fHDiffPixThetaOFF =
1214 (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF");
1215 if (!fHDiffPixThetaOFF)
1216 {
1217 *fLog << all
1218 << "MPad : Object '3D-ThetaPixDiffOFF' not found on root file"
1219 << endl;
1220 return kFALSE;
1221 }
1222 *fLog << all
1223 << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl;
1224
1225
1226 //------------------------------------
1227 fHSigmaPixTheta =
1228 (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
1229 if (!fHSigmaPixTheta)
1230 {
1231 *fLog << all
1232 << "MPad : Object '3D-ThetaPixSigma' not found on root file"
1233 << endl;
1234 return kFALSE;
1235 }
1236 *fLog << all
1237 << "MPad : Object '3D-ThetaPixSigma' was read in" << endl;
1238
1239 fHSigmaPixThetaMC =
1240 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC");
1241 if (!fHSigmaPixThetaMC)
1242 {
1243 *fLog << all
1244 << "MPad : Object '3D-ThetaPixSigmaMC' not found on root file"
1245 << endl;
1246 return kFALSE;
1247 }
1248 *fLog << all
1249 << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl;
1250
1251 fHSigmaPixThetaON =
1252 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON");
1253 if (!fHSigmaPixThetaON)
1254 {
1255 *fLog << all
1256 << "MPad : Object '3D-ThetaPixSigmaON' not found on root file"
1257 << endl;
1258 return kFALSE;
1259 }
1260 *fLog << all
1261 << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl;
1262
1263 fHSigmaPixThetaOFF =
1264 (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF");
1265 if (!fHSigmaPixThetaOFF)
1266 {
1267 *fLog << all
1268 << "MPad : Object '3D-ThetaPixSigmaOFF' not found on root file"
1269 << endl;
1270 return kFALSE;
1271 }
1272 *fLog << all
1273 << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl;
1274
1275 //------------------------------------
1276 fHgMC =
1277 (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC");
1278 if (!fHgMC)
1279 {
1280 *fLog << all
1281 << "MPad : Object '3D-PaddingMatrixMC' not found on root file"
1282 << endl;
1283 return kFALSE;
1284 }
1285 *fLog << all
1286 << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl;
1287
1288 fHgON =
1289 (TH3D*) gROOT->FindObject("3D-PaddingMatrixON");
1290 if (!fHgON)
1291 {
1292 *fLog << all
1293 << "MPad : Object '3D-PaddingMatrixON' not found on root file"
1294 << endl;
1295 return kFALSE;
1296 }
1297 *fLog << all
1298 << "MPad : Object '3D-PaddingMatrixON' was read in" << endl;
1299
1300 fHgOFF =
1301 (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF");
1302 if (!fHgOFF)
1303 {
1304 *fLog << all
1305 << "MPad : Object '3D-PaddingMatrixOFF' not found on root file"
1306 << endl;
1307 return kFALSE;
1308 }
1309 *fLog << all
1310 << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl;
1311
1312 //------------------------------------
1313 fHBlindPixIdThetaMC =
1314 (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC");
1315 if (!fHBlindPixIdThetaMC)
1316 {
1317 *fLog << all
1318 << "MPad : Object '2D-ThetaBlindIdMC' not found on root file"
1319 << endl;
1320 return kFALSE;
1321 }
1322 *fLog << all
1323 << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl;
1324
1325 fHBlindPixIdThetaON =
1326 (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON");
1327 if (!fHBlindPixIdThetaON)
1328 {
1329 *fLog << all
1330 << "MPad : Object '2D-ThetaBlindIdON' not found on root file"
1331 << endl;
1332 return kFALSE;
1333 }
1334 *fLog << all
1335 << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl;
1336
1337 fHBlindPixIdThetaOFF =
1338 (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF");
1339 if (!fHBlindPixIdThetaOFF)
1340 {
1341 *fLog << all
1342 << "MPad : Object '2D-ThetaBlindIdOFF' not found on root file"
1343 << endl;
1344 return kFALSE;
1345 }
1346 *fLog << all
1347 << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl;
1348
1349 //------------------------------------
1350 fHBlindPixNThetaMC =
1351 (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC");
1352 if (!fHBlindPixNThetaMC)
1353 {
1354 *fLog << all
1355 << "MPad : Object '2D-ThetaBlindNMC' not found on root file"
1356 << endl;
1357 return kFALSE;
1358 }
1359 *fLog << all
1360 << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl;
1361
1362 fHBlindPixNThetaON =
1363 (TH2D*) gROOT->FindObject("2D-ThetaBlindNON");
1364 if (!fHBlindPixNThetaON)
1365 {
1366 *fLog << all
1367 << "MPad : Object '2D-ThetaBlindNON' not found on root file"
1368 << endl;
1369 return kFALSE;
1370 }
1371 *fLog << all
1372 << "MPad : Object '2D-ThetaBlindNON' was read in" << endl;
1373
1374 fHBlindPixNThetaOFF =
1375 (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF");
1376 if (!fHBlindPixNThetaOFF)
1377 {
1378 *fLog << all
1379 << "MPad : Object '2D-ThetaBlindNOFF' not found on root file"
1380 << endl;
1381 return kFALSE;
1382 }
1383 *fLog << all
1384 << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl;
1385
1386 //------------------------------------
1387
1388 return kTRUE;
1389}
1390
1391
1392// --------------------------------------------------------------------------
1393//
1394// Write padding distributions onto a file
1395//
1396//
1397Bool_t MPad::WritePaddingDist(const char* namefileout)
1398{
1399 *fLog << all << "MPad : Write padding histograms onto file "
1400 << namefileout << endl;
1401
1402 TFile outfile(namefileout, "RECREATE");
1403
1404 fHBlindPixNThetaMC->Write();
1405 fHBlindPixNThetaON->Write();
1406 fHBlindPixNThetaOFF->Write();
1407
1408 fHBlindPixIdThetaMC->Write();
1409 fHBlindPixIdThetaON->Write();
1410 fHBlindPixIdThetaOFF->Write();
1411
1412 fHSigmaThetaMC->Write();
1413 fHSigmaThetaON->Write();
1414 fHSigmaThetaOFF->Write();
1415
1416 fHDiffPixTheta->Write();
1417 fHDiffPixThetaMC->Write();
1418 fHDiffPixThetaON->Write();
1419 fHDiffPixThetaOFF->Write();
1420
1421 fHSigmaPixTheta->Write();
1422 fHSigmaPixThetaMC->Write();
1423 fHSigmaPixThetaON->Write();
1424 fHSigmaPixThetaOFF->Write();
1425
1426 fHgMC->Write();
1427 fHgON->Write();
1428 fHgOFF->Write();
1429
1430 *fLog << all
1431 << "MPad::WriteTargetDist(); padding histograms written onto file "
1432 << namefileout << endl;
1433
1434 return kTRUE;
1435}
1436
1437// --------------------------------------------------------------------------
1438//
1439// Set type of data to be padded
1440//
1441// this is not necessary if the type of the data can be recognized
1442// directly from the input
1443//
1444//
1445void MPad::SetDataType(const char* type)
1446{
1447 fType = type;
1448 *fLog << all << "MPad::SetDataType(); type of data to be padded : "
1449 << fType << endl;
1450
1451 return;
1452}
1453
1454// --------------------------------------------------------------------------
1455//
1456// Set the pointers and prepare the histograms
1457//
1458Int_t MPad::PreProcess(MParList *pList)
1459{
1460 if ( !fHSigmaThetaMC || !fHSigmaThetaON || !fHSigmaThetaOFF ||
1461 !fHDiffPixThetaMC || !fHDiffPixThetaON || !fHDiffPixThetaOFF ||
1462 !fHBlindPixIdThetaMC|| !fHBlindPixIdThetaON|| !fHBlindPixIdThetaOFF ||
1463 !fHBlindPixNThetaMC || !fHBlindPixNThetaON || !fHBlindPixNThetaOFF ||
1464 !fHgMC || !fHgON || !fHgOFF )
1465 {
1466 *fLog << err
1467 << "MPad : At least one of the histograms needed for the padding is not defined ... aborting."
1468 << endl;
1469 return kFALSE;
1470 }
1471
1472 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
1473 if (!fMcEvt)
1474 {
1475 *fLog << err << dbginf << "MPad : MMcEvt not found... aborting."
1476 << endl;
1477 return kFALSE;
1478 }
1479
1480 fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam");
1481 if (!fPed)
1482 {
1483 *fLog << err << "MPad : MPedPhotCam not found... aborting."
1484 << endl;
1485 return kFALSE;
1486 }
1487
1488 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
1489 if (!fCam)
1490 {
1491 *fLog << err
1492 << "MPad : MGeomCam not found (no geometry information available)... aborting."
1493 << endl;
1494 return kFALSE;
1495 }
1496
1497 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
1498 if (!fEvt)
1499 {
1500 *fLog << err << "MPad : MCerPhotEvt not found... aborting."
1501 << endl;
1502 return kFALSE;
1503 }
1504
1505 fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
1506 if (!fSigmabar)
1507 {
1508 *fLog << err << "MPad : MSigmabar not found... aborting." << endl;
1509 return kFALSE;
1510 }
1511
1512 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
1513 if (!fBlinds)
1514 {
1515 *fLog << err << "MPad : MBlindPixels not found... aborting."
1516 << endl;
1517 return kFALSE;
1518 }
1519
1520 if (fType !="ON" && fType !="OFF" && fType !="MC")
1521 {
1522 *fLog << err
1523 << "MPad : Type of data to be padded not defined... aborting."
1524 << endl;
1525 return kFALSE;
1526 }
1527
1528
1529 //--------------------------------------------------------------------
1530 // histograms for checking the padding
1531 //
1532 // plot pedestal sigmas
1533 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
1534 100, 0.0, 5.0, 100, 0.0, 5.0);
1535 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
1536 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
1537
1538 // plot no.of photons (before vs. after padding)
1539 fHPhotons = new TH2D("Photons","Photons: after vs.before padding",
1540 100, -10.0, 90.0, 100, -10, 90);
1541 fHPhotons->SetXTitle("no.of photons before padding");
1542 fHPhotons->SetYTitle("no.of photons after padding");
1543
1544 // plot of added NSB
1545 fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
1546 fHNSB->SetXTitle("no.of added NSB photons");
1547 fHNSB->SetYTitle("no.of pixels");
1548
1549
1550 //--------------------------------------------------------------------
1551
1552 fIter = 20;
1553
1554 memset(fErrors, 0, sizeof(fErrors));
1555 memset(fWarnings, 0, sizeof(fWarnings));
1556
1557 return kTRUE;
1558}
1559
1560// --------------------------------------------------------------------------
1561//
1562// Do the Padding (noise adjustment)
1563//
1564// input for the padding :
1565// - the matrices fHgON, fHgOFF, fHgMC
1566// - the original distributions fHBlindPixNTheta, fHBlindPixIdTheta
1567// fHSigmaTheta, fHDiffPixTheta
1568//
1569
1570Int_t MPad::Process()
1571{
1572 *fLog << all << "Entry MPad::Process();" << endl;
1573
1574 // this is the conversion factor from the number of photons
1575 // to the number of photo electrons
1576 // later to be taken from MCalibration
1577 // fPEperPhoton = PW * LG * QE * 1D
1578 Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9;
1579
1580 //------------------------------------------------
1581 // add pixels to MCerPhotEvt which are not yet in;
1582 // set their number of photons equal to zero
1583
1584 UInt_t imaxnumpix = fCam->GetNumPixels();
1585
1586 for (UInt_t i=0; i<imaxnumpix; i++)
1587 {
1588 Bool_t alreadythere = kFALSE;
1589 UInt_t npix = fEvt->GetNumPixels();
1590 for (UInt_t j=0; j<npix; j++)
1591 {
1592 MCerPhotPix &pix = (*fEvt)[j];
1593 UInt_t id = pix.GetPixId();
1594 if (i==id)
1595 {
1596 alreadythere = kTRUE;
1597 break;
1598 }
1599 }
1600 if (alreadythere)
1601 continue;
1602
1603 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms());
1604 }
1605
1606
1607
1608 //-----------------------------------------
1609 Int_t rc=0;
1610 Int_t rw=0;
1611
1612 const UInt_t npix = fEvt->GetNumPixels();
1613
1614 //-------------------------------------------
1615 // Calculate average sigma of the event
1616 //
1617 Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt);
1618 *fLog << all << "MPad::Process(); before padding : " << endl;
1619 fSigmabar->Print("");
1620 Double_t sigbarold = sigbarold_phot * fPEperPhoton;
1621 Double_t sigbarold2 = sigbarold*sigbarold;
1622
1623 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
1624 *fLog << all << "theta = " << theta << endl;
1625
1626 Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
1627 if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() )
1628 {
1629 *fLog << warn
1630 << "MPad::Process(); binNumber out of range : theta, binTheta = "
1631 << theta << ", " << binTheta << "; Skip event " << endl;
1632 // event cannot be padded; skip event
1633
1634 rc = 2;
1635 fErrors[rc]++;
1636 return kCONTINUE;
1637 }
1638
1639 //-------------------------------------------
1640 // get number of events in this theta bin
1641 // and number of events in this sigbarold_phot bin
1642
1643 Double_t nTheta;
1644 Double_t nSigma;
1645
1646 TH2D *st = NULL;
1647 if (fType == "MC") st = fHSigmaThetaMC;
1648 else if (fType == "ON") st = fHSigmaThetaON;
1649 else if (fType == "OFF") st = fHSigmaThetaOFF;
1650 else
1651 {
1652 *fLog << err << "MPad : illegal data type '" << fType
1653 << "'" << endl;
1654 return kERROR;
1655 }
1656
1657 TH1D *hn;
1658 hn = st->ProjectionY("", binTheta, binTheta, "");
1659 nTheta = hn->Integral();
1660 Int_t binSigma = hn->FindBin(sigbarold_phot);
1661 nSigma = hn->GetBinContent(binSigma);
1662
1663 *fLog << all
1664 << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1665 << theta << ", " << sigbarold_phot << ", " << binTheta << ", "
1666 << binSigma << ", " << nTheta << ", " << nSigma << endl;
1667
1668 delete hn;
1669
1670
1671 //-------------------------------------------
1672 // for the current theta :
1673 // generate blind pixels according to the histograms
1674 // fHBlindPixNTheta and fHBlindPixIDTheta
1675 //
1676 // ON : add the blind pixels from the OFF data
1677 // OFF : add the blind pixels from the ON data
1678 // MC : add the blind pixels from the ON and OFF data
1679
1680 //-----------------------------------
1681 if (fType == "ON" || fType == "MC")
1682 {
1683 // numBlind is the number of blind pixels in an event
1684 TH1D *nblind;
1685 UInt_t numBlind;
1686
1687 nblind = fHBlindPixNThetaOFF->ProjectionY("", binTheta, binTheta, "");
1688 if ( nblind->Integral() == 0.0 )
1689 {
1690 *fLog << warn << "MPad::Process(); projection for Theta bin "
1691 << binTheta << " has no entries; Skip event " << endl;
1692 // event cannot be padded; skip event
1693 delete nblind;
1694
1695 rc = 7;
1696 fErrors[rc]++;
1697 return kCONTINUE;
1698 }
1699 else
1700 {
1701 numBlind = (Int_t) (nblind->GetRandom()+0.5);
1702 }
1703 delete nblind;
1704
1705
1706 // throw the Id of numBlind different pixels in this event
1707 if ( numBlind > 0)
1708 {
1709 TH1D *hblind;
1710 UInt_t idBlind;
1711 UInt_t listId[npix];
1712 UInt_t nlist = 0;
1713 Bool_t equal;
1714
1715 hblind = fHBlindPixIdThetaOFF->ProjectionY("", binTheta, binTheta, "");
1716 if ( hblind->Integral() > 0.0 )
1717 for (UInt_t i=0; i<numBlind; i++)
1718 {
1719 while(1)
1720 {
1721 idBlind = (Int_t) (hblind->GetRandom()+0.5);
1722 equal = kFALSE;
1723 for (UInt_t j=0; j<nlist; j++)
1724 if (idBlind == listId[j])
1725 {
1726 equal = kTRUE;
1727 break;
1728 }
1729 if (!equal) break;
1730 }
1731 listId[nlist] = idBlind;
1732 nlist++;
1733
1734 fBlinds->SetPixelBlind(idBlind);
1735 *fLog << all << "idBlind = " << idBlind << endl;
1736 }
1737 fBlinds->SetReadyToSave();
1738
1739 delete hblind;
1740 }
1741 }
1742
1743 //------------------------------------
1744 if (fType == "OFF" || fType == "MC")
1745 {
1746 // throw numBlind;
1747 // numBlind is the number of blind pixels in an event
1748 TH1D *nblind;
1749 UInt_t numBlind;
1750
1751 nblind = fHBlindPixNThetaON->ProjectionY("", binTheta, binTheta, "");
1752 if ( nblind->Integral() == 0.0 )
1753 {
1754 *fLog << warn << "MPad::Process(); projection for Theta bin "
1755 << binTheta << " has no entries; Skip event " << endl;
1756 // event cannot be padded; skip event
1757 delete nblind;
1758
1759 rc = 7;
1760 fErrors[rc]++;
1761 return kCONTINUE;
1762 }
1763 else
1764 {
1765 numBlind = (Int_t) (nblind->GetRandom()+0.5);
1766 }
1767 delete nblind;
1768
1769
1770 // throw the Id of numBlind different pixels in this event
1771 if ( numBlind > 0)
1772 {
1773 TH1D *hblind;
1774 UInt_t idBlind;
1775 UInt_t listId[npix];
1776 UInt_t nlist = 0;
1777 Bool_t equal;
1778
1779 hblind = fHBlindPixIdThetaON->ProjectionY("", binTheta, binTheta, "");
1780 if ( hblind->Integral() > 0.0 )
1781 for (UInt_t i=0; i<numBlind; i++)
1782 {
1783 while(1)
1784 {
1785 idBlind = (Int_t) (hblind->GetRandom()+0.5);
1786 equal = kFALSE;
1787 for (UInt_t j=0; j<nlist; j++)
1788 if (idBlind == listId[j])
1789 {
1790 equal = kTRUE;
1791 break;
1792 }
1793 if (!equal) break;
1794 }
1795 listId[nlist] = idBlind;
1796 nlist++;
1797
1798 fBlinds->SetPixelBlind(idBlind);
1799 *fLog << all << "idBlind = " << idBlind << endl;
1800 }
1801 fBlinds->SetReadyToSave();
1802
1803 delete hblind;
1804 }
1805 }
1806
1807 //******************************************************************
1808 // has the event to be padded ?
1809 // if yes : to which sigmabar should it be padded ?
1810 //
1811
1812 Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot);
1813 *fLog << all << "binSig, sigbarold_phot = " << binSig << ", "
1814 << sigbarold_phot << endl;
1815
1816 Double_t prob;
1817 TH1D *hpad = NULL;
1818
1819 TH3D *Hg = NULL;
1820 if (fType == "MC") Hg = fHgMC;
1821 else if (fType == "ON") Hg = fHgON;
1822 else if (fType == "OFF") Hg = fHgOFF;
1823 else
1824 {
1825 *fLog << err << "MPad : illegal type of data '" << fType << "'"
1826 << endl;
1827 return kERROR;
1828 }
1829
1830 hpad = Hg->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");
1831
1832 //Int_t nb = hpad->GetNbinsX();
1833 //for (Int_t i=1; i<=nb; i++)
1834 //{
1835 // if (hpad->GetBinContent(i) != 0.0)
1836 // *fLog << all << "i, fHgON = " << i << ", "
1837 // << hpad->GetBinContent(i) << endl;
1838 //}
1839
1840 if (nSigma != 0.0)
1841 prob = hpad->Integral() * nTheta/nSigma;
1842 else
1843 prob = 0.0;
1844
1845 *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma
1846 << ", " << prob << endl;
1847
1848 if ( prob <= 0.0 || gRandom->Uniform() > prob )
1849 {
1850 delete hpad;
1851 // event should not be padded
1852 *fLog << all << "event is not padded" << endl;
1853
1854 rc = 8;
1855 fErrors[rc]++;
1856 return kTRUE;
1857 }
1858 // event should be padded
1859 *fLog << all << "event will be padded" << endl;
1860
1861
1862 //-------------------------------------------
1863 // for the current theta, generate a sigmabar :
1864 // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF
1865 //
1866 Double_t sigmabar_phot = 0;
1867 Double_t sigmabar = 0;
1868
1869
1870 //Int_t nbinsx = hpad->GetNbinsX();
1871 //for (Int_t i=1; i<=nbinsx; i++)
1872 //{
1873 // if (hpad->GetBinContent(i) != 0.0)
1874 // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i)
1875 // << endl;
1876 //}
1877
1878 sigmabar_phot = hpad->GetRandom();
1879 sigmabar = sigmabar_phot * fPEperPhoton;
1880
1881 *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl;
1882
1883 delete hpad;
1884
1885 const Double_t sigmabar2 = sigmabar*sigmabar;
1886
1887 //-------------------------------------------
1888
1889 *fLog << all << "MPad::Process(); sigbarold, sigmabar = "
1890 << sigbarold << ", "<< sigmabar << endl;
1891
1892 // Skip event if target sigmabar is <= sigbarold
1893 if (sigmabar <= sigbarold)
1894 {
1895 *fLog << all << "MPad::Process(); target sigmabar is less than sigbarold : "
1896 << sigmabar << ", " << sigbarold << ", Skip event" << endl;
1897
1898 rc = 4;
1899 fErrors[rc]++;
1900 return kCONTINUE;
1901 }
1902
1903
1904 //-------------------------------------------
1905 //
1906 // Calculate average number of NSB photo electrons to be added (lambdabar)
1907 // from the value of sigmabar,
1908 // - making assumptions about the average electronic noise (elNoise2) and
1909 // - using a fixed value (F2excess) for the excess noise factor
1910
1911 Double_t elNoise2; // [photo electrons]
1912 Double_t F2excess = 1.3;
1913 Double_t lambdabar; // [photo electrons]
1914
1915
1916 Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta);
1917 if (binTheta != bincheck)
1918 {
1919 *fLog << err
1920 << "MPad::Process(); The binnings of the 2 histograms are not identical; aborting"
1921 << endl;
1922 return kERROR;
1923 }
1924
1925 // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin.
1926 // The average electronic noise (to be added) has to be well above this RMS,
1927 // otherwise the electronic noise of an individual pixel (elNoise2Pix)
1928 // may become negative
1929
1930 TH1D *hnoise;
1931 if (fType == "MC")
1932 hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1933 else if (fType == "ON")
1934 hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1935 else if (fType == "OFF")
1936 hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, "");
1937 else
1938 {
1939 *fLog << err << "MPad::Process(); illegal data type... aborting"
1940 << endl;
1941 return kERROR;
1942 }
1943
1944 Double_t RMS_phot = hnoise->GetRMS(1);
1945 Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
1946 delete hnoise;
1947
1948 elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2);
1949 *fLog << all << "elNoise2 = " << elNoise2 << endl;
1950
1951 lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess;
1952
1953 // This value of lambdabar is the same for all pixels;
1954 // note that lambdabar is normalized to the area of pixel 0
1955
1956 //---------- start loop over pixels ---------------------------------
1957 // do the padding for each pixel
1958 //
1959 // pad only pixels - which are used (before image cleaning)
1960 //
1961
1962 Double_t sigma2 = 0;
1963
1964 Double_t diff_phot = 0;
1965 Double_t diff = 0;
1966
1967 Double_t addSig2_phot= 0;
1968 Double_t addSig2 = 0;
1969
1970 Double_t elNoise2Pix = 0;
1971
1972
1973 for (UInt_t i=0; i<npix; i++)
1974 {
1975 MCerPhotPix &pix = (*fEvt)[i];
1976 if ( !pix.IsPixelUsed() )
1977 continue;
1978
1979 //if ( pix.GetNumPhotons() == 0.0)
1980 //{
1981 // *fLog << warn
1982 // << "MPad::Process(); no.of photons is 0 for used pixel"
1983 // << endl;
1984 // continue;
1985 //}
1986
1987 Int_t j = pix.GetPixId();
1988
1989 // GetPixRatio returns (area of pixel 0 / area of current pixel)
1990 Double_t ratioArea = 1.0 / fCam->GetPixRatio(j);
1991
1992 MPedPhotPix &ppix = (*fPed)[j];
1993 Double_t oldsigma_phot = ppix.GetRms();
1994 Double_t oldsigma = oldsigma_phot * fPEperPhoton;
1995 Double_t oldsigma2 = oldsigma*oldsigma;
1996
1997 //---------------------------------
1998 // throw the Sigma for this pixel
1999 //
2000 Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );
2001
2002 Int_t count;
2003 Bool_t ok;
2004
2005 TH1D *hdiff;
2006
2007 // throw the Sigma for this pixel from the distribution fHDiffPixTheta
2008 // MC : from fHDiffPixTheta
2009 // ON : from fHDiffPixThetaOFF
2010 // OFF : from fHDiffPixThetaON
2011
2012 TH3D *sp = NULL;
2013
2014 if (fType == "MC") sp = fHDiffPixTheta;
2015 else if (fType == "ON") sp = fHDiffPixThetaOFF;
2016 else if (fType == "OFF") sp = fHDiffPixThetaON;
2017 else
2018 {
2019 *fLog << err << "MPad::Process(); illegal data type... aborting"
2020 << endl;
2021 return kERROR;
2022 }
2023
2024 hdiff = sp->ProjectionZ("", binTheta, binTheta,
2025 binPixel, binPixel, "");
2026
2027 if ( hdiff->GetEntries() == 0 )
2028 {
2029 *fLog << warn << "MPad::Process(); projection for Theta bin "
2030 << binTheta << " and pixel bin " << binPixel
2031 << " has no entries; aborting " << endl;
2032 delete hdiff;
2033
2034 rc = 5;
2035 fErrors[rc]++;
2036 return kCONTINUE;
2037 }
2038
2039 count = 0;
2040 ok = kFALSE;
2041 for (Int_t m=0; m<fIter; m++)
2042 {
2043 count += 1;
2044 diff_phot = hdiff->GetRandom();
2045 diff = diff_phot * fPEperPhoton * fPEperPhoton;
2046
2047 // the following condition ensures that elNoise2Pix > 0.0
2048 if ( (diff + sigmabar2 - oldsigma2/ratioArea
2049 - lambdabar*F2excess) > 0.0 )
2050 {
2051 ok = kTRUE;
2052 break;
2053 }
2054 }
2055
2056 if (!ok)
2057 {
2058 *fLog << all << "theta, j, count, sigmabar, diff = " << theta
2059 << ", "
2060 << j << ", " << count << ", " << sigmabar << ", "
2061 << diff << endl;
2062 diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2;
2063
2064 rw = 1;
2065 fWarnings[rw]++;
2066 }
2067 else
2068 {
2069 rw = 0;
2070 fWarnings[rw]++;
2071 }
2072
2073 delete hdiff;
2074 sigma2 = diff + sigmabar2;
2075
2076
2077 //---------------------------------
2078 // get the additional sigma^2 for this pixel (due to the padding)
2079
2080 addSig2 = sigma2*ratioArea - oldsigma2;
2081 addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
2082
2083 //---------------------------------
2084 // get the additional electronic noise for this pixel
2085
2086 elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;
2087
2088
2089 //---------------------------------
2090 // throw actual number of additional NSB photo electrons (NSB)
2091 // and its RMS (sigmaNSB)
2092
2093 Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
2094 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
2095 Double_t sigmaNSB0;
2096
2097 if (arg >= 0.0)
2098 {
2099 sigmaNSB0 = sqrt( arg );
2100 }
2101 else
2102 {
2103 sigmaNSB0 = 0.0000001;
2104 if (arg < -1.e-10)
2105 {
2106 *fLog << warn << "MPad::Process(); argument of sqrt < 0 : "
2107 << arg << endl;
2108 }
2109 }
2110
2111
2112 //---------------------------------
2113 // smear NSB0 according to sigmaNSB0
2114 // and subtract lambdabar because of AC coupling
2115
2116 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
2117 Double_t NSB_phot = NSB / fPEperPhoton;
2118
2119 //---------------------------------
2120
2121 // add additional NSB to the number of photons
2122 Double_t oldphotons_phot = pix.GetNumPhotons();
2123 Double_t oldphotons = oldphotons_phot * fPEperPhoton;
2124 Double_t newphotons = oldphotons + NSB;
2125 Double_t newphotons_phot = newphotons / fPEperPhoton;
2126 pix.SetNumPhotons( newphotons_phot );
2127
2128
2129 fHNSB->Fill( NSB_phot/ratioArea );
2130 fHPhotons->Fill( oldphotons_phot/ratioArea,
2131 newphotons_phot/ratioArea );
2132
2133
2134 // error: add sigma of padded noise quadratically
2135 Double_t olderror_phot = pix.GetErrorPhot();
2136 Double_t newerror_phot =
2137 sqrt( olderror_phot*olderror_phot + addSig2_phot );
2138 pix.SetErrorPhot( newerror_phot );
2139
2140
2141 Double_t newsigma = sqrt( oldsigma2 + addSig2 );
2142 Double_t newsigma_phot = newsigma / fPEperPhoton;
2143 ppix.SetRms( newsigma_phot );
2144
2145 fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot );
2146 }
2147 //---------- end of loop over pixels ---------------------------------
2148
2149 // Calculate sigmabar again and crosscheck
2150
2151 fSigmabar->Calc(*fCam, *fPed, *fEvt);
2152 *fLog << all << "MPad::Process(); after padding : " << endl;
2153 fSigmabar->Print("");
2154
2155 *fLog << all << "Exit MPad::Process();" << endl;
2156
2157 rc = 0;
2158 fErrors[rc]++;
2159 return kTRUE;
2160 //******************************************************************
2161}
2162
2163// --------------------------------------------------------------------------
2164//
2165//
2166Int_t MPad::PostProcess()
2167{
2168 if (GetNumExecutions() != 0)
2169 {
2170
2171 *fLog << inf << endl;
2172 *fLog << GetDescriptor() << " execution statistics:" << endl;
2173 *fLog << dec << setfill(' ');
2174
2175 int temp;
2176 if (fWarnings[0] != 0.0)
2177 temp = (int)(fWarnings[1]*100/fWarnings[0]);
2178 else
2179 temp = 100;
2180
2181 *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3)
2182 << temp
2183 << "%) Warning : iteration to find acceptable sigma failed"
2184 << ", fIter = " << fIter << endl;
2185
2186 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3)
2187 << (int)(fErrors[1]*100/GetNumExecutions())
2188 << "%) Evts skipped due to: Sigmabar_old > 0" << endl;
2189
2190 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
2191 << (int)(fErrors[2]*100/GetNumExecutions())
2192 << "%) Evts skipped due to: Zenith angle out of range" << endl;
2193
2194 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
2195 << (int)(fErrors[3]*100/GetNumExecutions())
2196 << "%) Evts skipped due to: No data for generating Sigmabar" << endl;
2197
2198 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
2199 << (int)(fErrors[4]*100/GetNumExecutions())
2200 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;
2201
2202 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3)
2203 << (int)(fErrors[5]*100/GetNumExecutions())
2204 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl;
2205
2206 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3)
2207 << (int)(fErrors[6]*100/GetNumExecutions())
2208 << "%) Evts skipped due to: No data for generating Sigma" << endl;
2209
2210 *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3)
2211 << (int)(fErrors[7]*100/GetNumExecutions())
2212 << "%) Evts skipped due to: No data for generating Blind pixels" << endl;
2213
2214 *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3)
2215 << (int)(fErrors[8]*100/GetNumExecutions())
2216 << "%) Evts didn't have to be padded" << endl;
2217
2218 *fLog << " " << fErrors[0] << " ("
2219 << (int)(fErrors[0]*100/GetNumExecutions())
2220 << "%) Evts were successfully padded!" << endl;
2221 *fLog << endl;
2222
2223 }
2224
2225 //---------------------------------------------------------------
2226 TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1500));
2227 c.Divide(3, 5);
2228 gROOT->SetSelectedPad(NULL);
2229
2230 c.cd(1);
2231 fHNSB->SetDirectory(NULL);
2232 fHNSB->DrawCopy();
2233 fHNSB->SetBit(kCanDelete);
2234
2235 c.cd(2);
2236 fHSigmaPedestal->SetDirectory(NULL);
2237 fHSigmaPedestal->DrawCopy();
2238 fHSigmaPedestal->SetBit(kCanDelete);
2239
2240 c.cd(3);
2241 fHPhotons->SetDirectory(NULL);
2242 fHPhotons->DrawCopy();
2243 fHPhotons->SetBit(kCanDelete);
2244
2245 //--------------------------------------------------------------------
2246
2247
2248 c.cd(4);
2249 fHSigmaTheta->SetDirectory(NULL);
2250 fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta");
2251 fHSigmaTheta->DrawCopy();
2252 fHSigmaTheta->SetBit(kCanDelete);
2253
2254
2255 //--------------------------------------------------------------------
2256
2257
2258 c.cd(7);
2259 fHBlindPixNTheta->SetDirectory(NULL);
2260 fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta");
2261 fHBlindPixNTheta->DrawCopy();
2262 fHBlindPixNTheta->SetBit(kCanDelete);
2263
2264 //--------------------------------------------------------------------
2265
2266
2267 c.cd(10);
2268 fHBlindPixIdTheta->SetDirectory(NULL);
2269 fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta");
2270 fHBlindPixIdTheta->DrawCopy();
2271 fHBlindPixIdTheta->SetBit(kCanDelete);
2272
2273
2274 //--------------------------------------------------------------------
2275 // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2
2276
2277 c.cd(5);
2278 TH2D *l1 = NULL;
2279 l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
2280 l1->SetDirectory(NULL);
2281 l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
2282 l1->SetXTitle("\\Theta [\\circ]");
2283 l1->SetYTitle("Sigma^2-Sigmabar^2");
2284
2285 l1->DrawCopy("box");
2286 l1->SetBit(kCanDelete);;
2287
2288 c.cd(8);
2289 TH2D *l2 = NULL;
2290 l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
2291 l2->SetDirectory(NULL);
2292 l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
2293 l2->SetXTitle("pixel");
2294 l2->SetYTitle("Sigma^2-Sigmabar^2");
2295
2296 l2->DrawCopy("box");
2297 l2->SetBit(kCanDelete);;
2298
2299 //--------------------------------------------------------------------
2300 // draw the 3D histogram (target): Theta, pixel, Sigma
2301
2302 c.cd(6);
2303 TH2D *k1 = NULL;
2304 k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
2305 k1->SetDirectory(NULL);
2306 k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)");
2307 k1->SetXTitle("\\Theta [\\circ]");
2308 k1->SetYTitle("Sigma");
2309
2310 k1->DrawCopy("box");
2311 k1->SetBit(kCanDelete);
2312
2313 c.cd(9);
2314 TH2D *k2 = NULL;
2315 k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
2316 k2->SetDirectory(NULL);
2317 k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)");
2318 k2->SetXTitle("pixel");
2319 k2->SetYTitle("Sigma");
2320
2321 k2->DrawCopy("box");
2322 k2->SetBit(kCanDelete);;
2323
2324
2325 //--------------------------------------------------------------------
2326
2327 c.cd(13);
2328 TH2D *m1;
2329 m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
2330 m1->SetDirectory(NULL);
2331 m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
2332 m1->SetXTitle("Sigmabar-old");
2333 m1->SetYTitle("Sigmabar-new");
2334
2335 m1->DrawCopy("box");
2336 m1->SetBit(kCanDelete);;
2337
2338 c.cd(14);
2339 TH2D *m2;
2340 m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
2341 m2->SetDirectory(NULL);
2342 m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
2343 m2->SetXTitle("Sigmabar-old");
2344 m2->SetYTitle("Sigmabar-new");
2345
2346 m2->DrawCopy("box");
2347 m2->SetBit(kCanDelete);;
2348
2349 c.cd(15);
2350 TH2D *m3;
2351 m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy");
2352 m3->SetDirectory(NULL);
2353 m3->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)");
2354 m3->SetXTitle("Sigmabar-old");
2355 m3->SetYTitle("Sigmabar-new");
2356
2357 m3->DrawCopy("box");
2358 m3->SetBit(kCanDelete);;
2359
2360
2361 return kTRUE;
2362}
2363
2364// --------------------------------------------------------------------------
2365
2366
2367
2368
2369
2370
Note: See TracBrowser for help on using the repository browser.