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

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