source: branches/Mars_use_drstimefiles/manalysis/MPad.cc

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