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

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