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

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