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

Last change on this file since 5436 was 5436, checked in by wittek, 20 years ago
*** empty log message ***
File size: 73.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Wolfgang Wittek, 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 return kTRUE;
916}
917
918
919// --------------------------------------------------------------------------
920//
921// Merge the 2D-histograms (Theta, sigmabar)
922// hA, hB and hC
923//
924// - put the merged distribution into hM
925// - and define the matrices hgA, hgB and hgC
926//
927// These matrices tell which fraction of events should be padded
928// from bin k of sigmabar to bin j
929//
930
931Bool_t MPad::MergeABC(TString tA, TString tB, TString tC,
932 TH2D *hA, TH2D *hB, TH2D *hC, TH2D *hM,
933 TH3D *hgA, TH3D *hgB, TH3D *hgC, TString canv)
934{
935 *fLog << all << "MPad::MergeABC(); Entry" << endl;
936
937 TAxis *ax = hM->GetXaxis();
938 Int_t nbinstheta = ax->GetNbins();
939
940 TAxis *ay = hM->GetYaxis();
941 Int_t nbinssig = ay->GetNbins();
942
943 TArrayD edgesy;
944 edgesy.Set(nbinssig+1);
945 for (Int_t i=0; i<=nbinssig; i++)
946 {
947 edgesy[i] = ay->GetBinLowEdge(i+1);
948 //*fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i]
949 // << endl;
950 }
951 MBinning binsg;
952 binsg.SetEdges(edgesy);
953
954
955 //............ loop over Theta bins ...........................
956
957 //*fLog << all << "MPad::MergeABC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
958 for (Int_t l=1; l<=nbinstheta; l++)
959 {
960 //-------------------------------------------
961 // merge A (ON) and B (OFF) distributions
962 // input : hista is the normalized 1D distr. of sigmabar for A (ON) data
963 // histb is the normalized 1D distr. of sigmabar for B (OFF) data
964 // output : histap will be the merged distribution (AB)
965 // fHga(k,j) will tell which fraction of A (ON) events should be
966 // padded from bin k to bin j
967 // fHgb(k,j) will tell which fraction of B (OFF) events should be
968 // padded from bin k to bin j
969
970 TH2D * fHga = new TH2D;
971 MH::SetBinning(fHga, &binsg, &binsg);
972 TH2D * fHgb = new TH2D;
973 MH::SetBinning(fHgb, &binsg, &binsg);
974
975 TH1D *hista = hA->ProjectionY("sigon_y", l, l, "");
976 TString tita(canv);
977 tita += "-A (orig)-";
978 tita += l;
979 hista->SetNameTitle(tita, tita);
980 Stat_t suma = hista->Integral();
981 if (suma != 0.0) hista->Scale(1./suma);
982
983 TH1D *histap = new TH1D( (const TH1D&)*hista );
984 TString titab(canv);
985 titab += "-AB (merged)-";
986 titab += l;
987 histap->SetNameTitle(titab, titab);
988
989 TH1D *histb = hB->ProjectionY("sigoff_y", l, l, "");
990 TString titb(canv);
991 titb += "-B (orig)-";
992 titb += l;
993 histb->SetNameTitle(titb, titb);
994 Stat_t sumb = histb->Integral();
995 if (sumb != 0.0) histb->Scale(1./sumb);
996
997 if (suma == 0.0 || sumb == 0.0)
998 {
999 delete hista;
1000 delete histb;
1001 delete fHga;
1002 delete fHgb;
1003 delete histap;
1004
1005 *fLog << err
1006 << "cannot call Merge2Distributions(A=" << tA << ", B="
1007 << tB << ") for theta bin l = "
1008 << l << "; NevA, NevB = " << suma << ", " << sumb << endl;
1009
1010 // clear the corresponding area of fHSigmaTheta
1011 for (Int_t j=1; j<=nbinssig; j++)
1012 {
1013 Double_t a = 0.0;
1014 hM->SetBinContent(l, j, a);
1015
1016 //*fLog << "MPad::MergeABC; l, j, a = " << l << ", " << j << ", "
1017 // << a << endl;
1018 }
1019
1020 // go to next theta bin (l)
1021 continue;
1022 }
1023
1024 //*fLog << "call Merge2Distributions(A=" << tA << ", B="
1025 // << tB << ") for theta bin l = "
1026 // << l << endl;
1027
1028 TString canvab(canv);
1029 canvab += "AB-";
1030 canvab += l;
1031 Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig, canvab);
1032
1033 // fill fHgA and fHgB
1034 for (Int_t k=1; k<=nbinssig; k++)
1035 for (Int_t j=1; j<=nbinssig; j++)
1036 {
1037 Double_t a = fHga->GetBinContent(k,j);
1038 hgA->SetBinContent(l, k, j, a);
1039
1040 Double_t b = fHgb->GetBinContent(k,j);
1041 hgB->SetBinContent(l, k, j, b);
1042 }
1043
1044
1045 //-------------------------------------------------
1046 // if there is no further distribution to be merged
1047 // fill target distribution fHSigmaTheta
1048 //
1049 if (!hC)
1050 {
1051 for (Int_t j=1; j<=nbinssig; j++)
1052 {
1053 Double_t a = histap->GetBinContent(j);
1054 hM->SetBinContent(l, j, a);
1055
1056 //*fLog << "MPad::MergeABC; l, j, a = " << l << ", " << j << ", "
1057 // << a << endl;
1058 }
1059
1060 delete fHga;
1061 delete fHgb;
1062 delete histap;
1063 }
1064
1065 //-------------- next merge (start) ----------------------------------
1066 if (hC)
1067 {
1068
1069 // now merge AB (ON-OFF) and C (MC) distributions
1070 // input : histe is the result of the merging of A (ON) and B (OFF)
1071 // distributions
1072 // histf is the normalized 1D distr. of sigmabar for C (MC) data
1073 // output : histep will be the merged distribution (target distribution)
1074 // fHge(k,j) will tell which fraction of A (ON), B (OFF) events
1075 // should be padded from bin k to bin j
1076 // fHgf(k,j) will tell which fraction of C (MC) events should be
1077 // padded from bin k to bin j
1078
1079 TH2D * fHge = new TH2D;
1080 MH::SetBinning(fHge, &binsg, &binsg);
1081 TH2D * fHgf = new TH2D;
1082 MH::SetBinning(fHgf, &binsg, &binsg);
1083
1084 TH1D *histe = new TH1D( (const TH1D&)*histap);
1085 TString titabm(canv);
1086 titabm += "-AB (merged)-";
1087 titabm += l;
1088 histe->SetNameTitle(titabm, titabm);
1089
1090 TH1D *histep = new TH1D( (const TH1D&)*histe);
1091 TString titabcm(canv);
1092 titabcm +="-ABC (merged)-";
1093 titabcm += l;
1094 histep->SetNameTitle(titabcm, titabcm);
1095
1096 TH1D *histf = hC->ProjectionY("sigmc_y", l, l, "");
1097 TString titc(canv);
1098 titc += "-C (orig)-";
1099 titc += l;
1100 histf->SetNameTitle(titc, titc);
1101 Double_t sumf = histf->Integral();
1102 if (sumf != 0.0) histf->Scale(1./sumf);
1103
1104 if (sumf == 0.0)
1105 {
1106 delete hista;
1107 delete histb;
1108 delete fHga;
1109 delete fHgb;
1110 delete histap;
1111
1112 delete histe;
1113 delete histf;
1114 delete fHge;
1115 delete fHgf;
1116 delete histep;
1117
1118 *fLog << err
1119 << "cannot call Merge2Distributions(AB=" << tA << "-" << tB
1120 << ", C=" << tC << ") for theta bin l = "
1121 << l << "; NevC = " << sumf << endl;
1122
1123 // go to next theta bin (l)
1124 continue;
1125 }
1126
1127 //*fLog << "call Merge2Distributions(AB=" << tA << "-" << tB
1128 // << ", C=" << tC << ") for theta bin l = "
1129 // << l << endl;
1130
1131 TString canvabc(canv);
1132 canvabc += "ABC-";
1133 canvabc += l;
1134 Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig, canvabc);
1135
1136 // update fHgA and fHgB
1137 UpdateHg(fHga, histap, fHge, hgA, nbinssig, l);
1138 UpdateHg(fHgb, histap, fHge, hgB, nbinssig, l);
1139
1140 // fill fHgC
1141 for (Int_t k=1; k<=nbinssig; k++)
1142 for (Int_t j=1; j<=nbinssig; j++)
1143 {
1144 Double_t f = fHgf->GetBinContent(k,j);
1145 hgC->SetBinContent(l, k, j, f);
1146 }
1147
1148 // fill target distribution fHSigmaTheta
1149 //
1150 for (Int_t j=1; j<=nbinssig; j++)
1151 {
1152 Double_t ep = histep->GetBinContent(j);
1153 hM->SetBinContent(l, j, ep);
1154
1155 //*fLog << "MPad::MergeABC; l, j, ep = " << l << ", " << j << ", "
1156 // << ep << endl;
1157 }
1158
1159
1160 //-------------------------------------------
1161
1162 delete hista;
1163 delete histb;
1164 delete fHga;
1165 delete fHgb;
1166 delete histap;
1167
1168 delete histe;
1169 delete histf;
1170 delete fHge;
1171 delete fHgf;
1172 delete histep;
1173
1174 }
1175 //-------------- next merge (end) ----------------------------------
1176
1177 }
1178 //............ end of loop over Theta bins ....................
1179
1180 *fLog << all << "MPad::MergeABC(); Return" << endl;
1181
1182 return kTRUE;
1183}
1184// --------------------------------------------------------------------------
1185//
1186// Merge the sigmabar distributions of 2 samples (samples A and B)
1187//
1188// input : the original distributions for samples A and B (hista, histb)
1189//
1190// output : the prescription which fraction of events should be padded
1191// from the sigmabar bin k to bin j (fHgMC, fHgON, fHgOFF)
1192//
1193
1194Bool_t MPad::Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap,
1195 TH2D *fHga, TH2D *fHgb, Int_t nbinssig,
1196 TString canv )
1197{
1198 *fLog << all << "MPad::Merge2Distributions(); Entry" << endl;
1199
1200
1201 // hista is the normalized 1D histogram of sigmabar for sample A
1202 // histb is the normalized 1D histogram of sigmabar for sample B
1203 // histc is the difference A-B
1204
1205 // at the beginning, histap is a copy of hista
1206 // at the end, it will be the 1D histogram for sample A after padding
1207
1208 // at the beginning, histbp is a copy of histb
1209 // at the end, it will be the 1D histogram for sample B after padding
1210
1211 // at the beginning, histcp is a copy of histc
1212 // at the end, it should be the difference histap-histbp,
1213 // which should be zero
1214
1215 // fHga[k][j] tells which fraction of events from sample A should be padded
1216 // from sigma_k to sigma_j
1217
1218
1219 // fHgb[k][j] tells which fraction of events from sample B should be padded
1220 // from sigma_k to sigma_j
1221
1222 //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl;
1223
1224 Double_t eps = 1.e-10;
1225
1226 TH1D *histbp = new TH1D( (const TH1D&)*histb );
1227
1228 TH1D *histc = new TH1D( (const TH1D&)*hista );
1229 histc->Add(histb, -1.0);
1230
1231 TH1D *histcp = new TH1D( (const TH1D&)*histc );
1232
1233
1234 //-------- start j loop ------------------------------------------------
1235 // loop over bins in histc,
1236 // starting from the end (i.e. at the largest sigmabar)
1237 Double_t v, s, w, t, x, u, a, b, arg;
1238
1239 //*fLog << "Content of histcp : " << endl;
1240
1241 for (Int_t j=nbinssig; j >= 1; j--)
1242 {
1243 //*fLog << "j, hista, histb , histap : " << j << ", "
1244 // << hista->GetBinContent(j) << ", "
1245 // << histb->GetBinContent(j) << ", "
1246 // << histap->GetBinContent(j) << endl;
1247
1248 v = histcp->GetBinContent(j);
1249 //*fLog << "cp : j, v = " << j << ", " << v << endl;
1250
1251 if ( fabs(v) < eps ) continue;
1252 if (v >= 0.0)
1253 s = 1.0;
1254 else
1255 s = -1.0;
1256
1257 //................ start k loop ......................................
1258 // look for a bin k which may compensate the content of bin j
1259 for (Int_t k=j-1; k >= 1; k--)
1260 {
1261 w = histcp->GetBinContent(k);
1262 if ( fabs(w) < eps ) continue;
1263 if (w >= 0.0)
1264 t = 1.0;
1265 else
1266 t = -1.0;
1267
1268 // if s==t bin k cannot compensate, go to next k bin
1269 if (s == t) continue;
1270
1271 x = v + w;
1272 if (x >= 0.0)
1273 u = 1.0;
1274 else
1275 u = -1.0;
1276
1277 // if u==s bin k will partly compensate : pad the whole fraction
1278 // w from bin k to bin j
1279 if (u == s)
1280 arg = -w;
1281
1282 // if u!=s bin k would overcompensate : pad only the fraction
1283 // v from bin k to bin j
1284 else
1285 arg = v;
1286
1287 if (arg <=0.0)
1288 {
1289 fHga->SetBinContent(k, j, -arg);
1290 fHgb->SetBinContent(k, j, 0.0);
1291 }
1292 else
1293 {
1294 fHga->SetBinContent(k, j, 0.0);
1295 fHgb->SetBinContent(k, j, arg);
1296 }
1297
1298 //*fLog << all << "k, j, arg = " << k << ", " << j
1299 // << ", " << arg << endl;
1300
1301 //......................................
1302 // this is for checking the procedure
1303 if (arg < 0.0)
1304 {
1305 a = histap->GetBinContent(k);
1306 histap->SetBinContent(k, a+arg);
1307 a = histap->GetBinContent(j);
1308 histap->SetBinContent(j, a-arg);
1309 }
1310 else
1311 {
1312 b = histbp->GetBinContent(k);
1313 histbp->SetBinContent(k, b-arg);
1314 b = histbp->GetBinContent(j);
1315 histbp->SetBinContent(j, b+arg);
1316 }
1317 //......................................
1318
1319 if (u == s)
1320 {
1321 histcp->SetBinContent(k, 0.0);
1322 histcp->SetBinContent(j, x);
1323
1324 v = histcp->GetBinContent(j);
1325 if ( fabs(v) < eps ) break;
1326
1327 // bin j was only partly compensated : go to next k bin
1328 continue;
1329 }
1330 else
1331 {
1332 histcp->SetBinContent(k, x);
1333 histcp->SetBinContent(j, 0.0);
1334
1335 // bin j was completely compensated : go to next j bin
1336 break;
1337 }
1338
1339 }
1340 //................ end k loop ......................................
1341 }
1342 //-------- end j loop ------------------------------------------------
1343
1344 // check results
1345
1346 //*fLog << "Content of histap, histbp, histcp : " << endl;
1347
1348 for (Int_t j=1; j<=nbinssig; j++)
1349 {
1350 Double_t a = histap->GetBinContent(j);
1351 Double_t b = histbp->GetBinContent(j);
1352 Double_t c = histcp->GetBinContent(j);
1353
1354 //*fLog << "j, a, b, c = " << j << ": " << a << ", " << b << ", "
1355 // << c << endl;
1356
1357 if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps )
1358 *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = "
1359 << j << ", " << a << ", " << b << ", " << c << endl;
1360 }
1361
1362 //---------------------------------------------------------------
1363 TCanvas *pad = MH::MakeDefCanvas(canv, canv, 600, 600);
1364 gROOT->SetSelectedPad(NULL);
1365
1366 pad->Divide(2, 2);
1367
1368 pad->cd(1);
1369 gPad->SetBorderMode(0);
1370 hista->SetDirectory(NULL);
1371 hista->UseCurrentStyle();
1372 hista->DrawClone();
1373 hista->SetBit(kCanDelete);
1374
1375 pad->cd(2);
1376 gPad->SetBorderMode(0);
1377 histb->SetDirectory(NULL);
1378 histb->UseCurrentStyle();
1379 histb->DrawClone();
1380 histb->SetBit(kCanDelete);
1381
1382 pad->cd(3);
1383 gPad->SetBorderMode(0);
1384 histap->SetDirectory(NULL);
1385 histap->UseCurrentStyle();
1386 histap->DrawClone();
1387 histap->SetBit(kCanDelete);
1388
1389 pad->Draw();
1390
1391 //--------------------------------------------------------------------
1392
1393 delete histc;
1394 delete histbp;
1395 delete histcp;
1396
1397 *fLog << all << "MPad::Merge2Distributions(); Return" << endl;
1398
1399 return kTRUE;
1400}
1401
1402
1403
1404// --------------------------------------------------------------------------
1405//
1406// Update the matrix fHgA
1407// which contains the final padding prescriptions
1408//
1409
1410Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA,
1411 Int_t nbinssig, Int_t l)
1412{
1413 // histap target distribution after ON-OFF merging
1414 // fHga padding prescription for ON (or OFF) data to get to histap
1415 // fHge padding prescription to get from histap to the target
1416 // distribution after the ON-OFF-MC merging
1417 // fHgA updated padding prescription for ON (or OFF) data to get
1418 // from the original ON (or OFF) histogram to the target
1419 // distribution after the ON-OFF-MC merging
1420 // l Theta bin
1421
1422 Double_t eps = 1.e-10;
1423
1424 for (Int_t k=1; k<=nbinssig; k++)
1425 {
1426 Double_t na = fHga->Integral(1, nbinssig, k, k, " ");
1427 Double_t ne = fHge->Integral(k, k, 1, nbinssig, " ");
1428 Double_t nap = histap->GetBinContent(k);
1429
1430 if (ne <= eps)
1431 {
1432 // go to next k
1433 continue;
1434 }
1435
1436 Double_t r = nap - na;
1437
1438 if (r >= ne-eps)
1439 {
1440 for (Int_t j=k+1; j<=nbinssig; j++)
1441 {
1442 Double_t e = fHge->GetBinContent(k,j);
1443 Double_t A = fHgA->GetBinContent(l,k,j);
1444 fHgA->SetBinContent(l, k, j, A + e);
1445 }
1446 // go to next k
1447 continue;
1448 }
1449
1450 for (Int_t j=k+1; j<=nbinssig; j++)
1451 {
1452 Double_t e = fHge->GetBinContent(k,j);
1453 Double_t A = fHgA->GetBinContent(l,k,j);
1454 fHgA->SetBinContent(l, k, j, A + r*e/ne);
1455 }
1456
1457 if (na <= eps)
1458 {
1459 // go to next k
1460 continue;
1461 }
1462
1463 for (Int_t i=1; i<k; i++)
1464 {
1465 Double_t a = fHga->GetBinContent(i,k);
1466 Double_t A = fHgA->GetBinContent(l,i,k);
1467 fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na);
1468 }
1469
1470 for (Int_t i=1; i<=nbinssig; i++)
1471 {
1472 if (i == k) continue;
1473 for (Int_t j=i+1; j<=nbinssig; j++)
1474 {
1475 if (j == k) continue;
1476 Double_t a = fHga->GetBinContent(i,k);
1477 Double_t e = fHge->GetBinContent(k,j);
1478 Double_t A = fHgA->GetBinContent(l,i,j);
1479 fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne) );
1480 }
1481 }
1482 }
1483
1484 return kTRUE;
1485}
1486
1487// --------------------------------------------------------------------------
1488//
1489// Write padding distributions onto a file
1490// these are the histograms produced in the member function
1491// MergeONOFFMC
1492// plus the distributions produced by MMakePadHistograms
1493//
1494Bool_t MPad::WritePaddingDist(TString namefileout)
1495{
1496 *fLog << all << "MPad::WritePaddingDist(); Padding distributions for ";
1497
1498 TFile outfile(namefileout, "RECREATE");
1499
1500 fHSigmaTheta->Write();
1501 fHSigmaThetaOuter->Write();
1502
1503 if (fHSigmaThetaMC)
1504 {
1505 *fLog << all << " MC ";
1506 fHSigmaThetaMC->Write();
1507 fHSigmaThetaOuterMC->Write();
1508 fHDiffPixThetaMC->Write();
1509 fHDiffPixThetaTargetMC->Write();
1510 fHgMC->Write();
1511 fHgOuterMC->Write();
1512 }
1513
1514 if (fHSigmaThetaON)
1515 {
1516 *fLog << all << " ON ";
1517 fHSigmaThetaON->Write();
1518 fHSigmaThetaOuterON->Write();
1519 fHDiffPixThetaON->Write();
1520 fHDiffPixThetaTargetON->Write();
1521 fHgON->Write();
1522 fHgOuterON->Write();
1523 }
1524
1525 if (fHSigmaThetaOFF)
1526 {
1527 *fLog << all << " OFF ";
1528 fHSigmaThetaOFF->Write();
1529 fHSigmaThetaOuterOFF->Write();
1530 fHDiffPixThetaOFF->Write();
1531 fHDiffPixThetaTargetOFF->Write();
1532 fHgOFF->Write();
1533 fHgOuterOFF->Write();
1534 }
1535
1536 *fLog << all << " data were written onto file "
1537 << namefileout << endl;
1538
1539 return kTRUE;
1540}
1541
1542// --------------------------------------------------------------------------
1543//
1544// Read padding distributions from a file
1545// thes are the distributions which were written out by
1546// the member function WritePaddingDist
1547//
1548Bool_t MPad::ReadPaddingDist(TString namefilein)
1549{
1550 *fLog << all << "MPad : Read padding histograms from file "
1551 << namefilein << endl;
1552
1553 Int_t OK = 0;
1554
1555 fInfile = new TFile(namefilein);
1556 fInfile->ls();
1557
1558 //------------------------------------
1559
1560 fHSigmaTheta = new TH2D;
1561 OK = fHSigmaTheta->Read("2D-ThetaSigmabar");
1562 if (OK)
1563 {
1564 *fLog << all
1565 << "MPad : Object '2D-ThetaSigmabar' was read in" << endl;
1566 }
1567
1568 fHSigmaThetaOuter = new TH2D;
1569 OK =fHSigmaThetaOuter->Read("2D-ThetaSigmabarOuter");
1570 if (OK)
1571 {
1572 *fLog << all
1573 << "MPad : Object '2D-ThetaSigmabarOuter' was read in" << endl;
1574 }
1575
1576
1577 //------------------------------------
1578
1579 fHSigmaThetaMC = new TH2D;
1580 OK = fHSigmaThetaMC->Read("2D-ThetaSigmabarMC");
1581 if (OK)
1582 {
1583 *fLog << all
1584 << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl;
1585 }
1586
1587 fHSigmaThetaOuterMC = new TH2D;
1588 OK = fHSigmaThetaOuterMC->Read("2D-ThetaSigmabarOuterMC");
1589 if (OK)
1590 {
1591 *fLog << all
1592 << "MPad : Object '2D-ThetaSigmabarOuterMC' was read in" << endl;
1593 }
1594
1595 fHDiffPixThetaMC = new TH3D;
1596 OK = fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC");
1597 if (OK)
1598 {
1599 *fLog << all
1600 << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl;
1601 }
1602
1603 fHDiffPixThetaTargetMC = new TH3D;
1604 OK = fHDiffPixThetaTargetMC->Read("3D-ThetaPixDiffTargetMC");
1605 if (OK)
1606 {
1607 *fLog << all
1608 << "MPad : Object '3D-ThetaPixDiffTargetMC' was read in" << endl;
1609 }
1610
1611 fHgMC = new TH3D;
1612 OK = fHgMC->Read("3D-PaddingMatrixMC");
1613 if (OK)
1614 {
1615 *fLog << all
1616 << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl;
1617 }
1618
1619 fHgOuterMC = new TH3D;
1620 OK = fHgOuterMC->Read("3D-PaddingMatrixOuterMC");
1621 if (OK)
1622 {
1623 *fLog << all
1624 << "MPad : Object '3D-PaddingMatrixOuterMC' was read in" << endl;
1625 }
1626
1627 //------------------------------------
1628
1629 fHSigmaThetaON = new TH2D;
1630 OK =fHSigmaThetaON->Read("2D-ThetaSigmabarON");
1631 if (OK)
1632 {
1633 *fLog << all
1634 << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl;
1635 }
1636
1637 fHSigmaThetaOuterON = new TH2D;
1638 OK = fHSigmaThetaOuterON->Read("2D-ThetaSigmabarOuterON");
1639 if (OK)
1640 {
1641 *fLog << all
1642 << "MPad : Object '2D-ThetaSigmabarOuterON' was read in" << endl;
1643 }
1644
1645 fHDiffPixThetaON = new TH3D;
1646 OK = fHDiffPixThetaON->Read("3D-ThetaPixDiffON");
1647 if (OK)
1648 {
1649 *fLog << all
1650 << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl;
1651 }
1652
1653 fHDiffPixThetaTargetON = new TH3D;
1654 OK = fHDiffPixThetaTargetON->Read("3D-ThetaPixDiffTargetON");
1655 if (OK)
1656 {
1657 *fLog << all
1658 << "MPad : Object '3D-ThetaPixDiffTargetON' was read in" << endl;
1659 }
1660
1661 fHgON = new TH3D;
1662 OK = fHgON->Read("3D-PaddingMatrixON");
1663 if (OK)
1664 {
1665 *fLog << all
1666 << "MPad : Object '3D-PaddingMatrixON' was read in" << endl;
1667 }
1668
1669 fHgOuterON = new TH3D;
1670 OK = fHgOuterON->Read("3D-PaddingMatrixOuterON");
1671 if (OK)
1672 {
1673 *fLog << all
1674 << "MPad : Object '3D-PaddingMatrixOuterON' was read in" << endl;
1675 }
1676
1677 //------------------------------------
1678
1679 fHSigmaThetaOFF = new TH2D;
1680 OK = fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF");
1681 if (OK)
1682 {
1683 *fLog << all
1684 << "MPad : Object '2D-ThetaSigmabarOFF' was read in" << endl;
1685 }
1686
1687 fHSigmaThetaOuterOFF = new TH2D;
1688 OK = fHSigmaThetaOuterOFF->Read("2D-ThetaSigmabarOuterOFF");
1689 if (OK)
1690 {
1691 *fLog << all
1692 << "MPad : Object '2D-ThetaSigmabarOuterOFF' was read in" << endl;
1693 }
1694
1695 fHDiffPixThetaOFF = new TH3D;
1696 OK = fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF");
1697 if (OK)
1698 {
1699 *fLog << all
1700 << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl;
1701 }
1702
1703 fHDiffPixThetaTargetOFF = new TH3D;
1704 OK = fHDiffPixThetaTargetOFF->Read("3D-ThetaPixDiffTargetOFF");
1705 if (OK)
1706 {
1707 *fLog << all
1708 << "MPad : Object '3D-ThetaPixDiffTargetOFF' was read in" << endl;
1709 }
1710
1711 fHgOFF = new TH3D;
1712 OK = fHgOFF->Read("3D-PaddingMatrixOFF");
1713 if (OK)
1714 {
1715 *fLog << all
1716 << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl;
1717 }
1718
1719 fHgOuterOFF = new TH3D;
1720 OK = fHgOuterOFF->Read("3D-PaddingMatrixOuterOFF");
1721 if (OK)
1722 {
1723 *fLog << all
1724 << "MPad : Object '3D-PaddingMatrixOuterOFF' was read in" << endl;
1725 }
1726
1727 //------------------------------------
1728
1729 return kTRUE;
1730}
1731
1732
1733
1734
1735// --------------------------------------------------------------------------
1736//
1737// Set the pointers and prepare the histograms
1738//
1739Int_t MPad::PreProcess(MParList *pList)
1740{
1741 if ( !fHSigmaThetaMC && !fHSigmaThetaON && !fHSigmaThetaOFF ||
1742 !fHDiffPixThetaMC && !fHDiffPixThetaON && !fHDiffPixThetaOFF ||
1743 !fHgMC && !fHgON && !fHgOFF )
1744 {
1745 *fLog << err
1746 << "MPad : There are no input histograms for the padding ... aborting."
1747 << endl;
1748 return kFALSE;
1749 }
1750
1751 fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
1752 if (!fPointPos)
1753 {
1754 *fLog << err << dbginf << "MPad : MPointingPos not found... aborting."
1755 << endl;
1756 return kFALSE;
1757 }
1758
1759 fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
1760 if (!fPed)
1761 {
1762 *fLog << err << AddSerialNumber(fNamePedPhotCam)
1763 << "[MPedPhotCam] not found... aborting." << endl;
1764 return kFALSE;
1765 }
1766 *fLog << all << "MPad::PreProcess(); name of MPedPhotCam container = "
1767 << fNamePedPhotCam << endl;
1768
1769
1770 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
1771 if (!fCam)
1772 {
1773 *fLog << err
1774 << "MPad : MGeomCam not found (no geometry information available)... aborting."
1775 << endl;
1776 return kFALSE;
1777 }
1778
1779 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
1780 if (!fEvt)
1781 {
1782 *fLog << err << "MPad : MCerPhotEvt not found... aborting."
1783 << endl;
1784 return kFALSE;
1785 }
1786
1787
1788 fBad = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam");
1789 if (!fBad)
1790 {
1791 *fLog << inf
1792 << "MPad : MBadPixelsCam container not present... continue."
1793 << endl;
1794 }
1795
1796
1797 if (fType !="ON" && fType !="OFF" && fType !="MC")
1798 {
1799 *fLog << err
1800 << "MPad : Type of data to be padded not defined... aborting."
1801 << endl;
1802 return kFALSE;
1803 }
1804
1805
1806 //--------------------------------------------------------------------
1807 // histograms for checking the padding
1808 //
1809 // plot pedestal sigmas
1810 fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding",
1811 100, 0.0, 75.0, 100, 0.0, 75.0);
1812 fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
1813 fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");
1814 fHSigmaPedestal->SetDirectory(NULL);
1815 fHSigmaPedestal->UseCurrentStyle();
1816 fHSigmaPedestal->GetYaxis()->SetTitleOffset(1.25);
1817
1818 // plot no.of photons (before vs. after padding)
1819 fHPhotons = new TH2D("Photons/area",
1820 "Photons/area: after vs.before padding",
1821 100, -100.0, 300.0, 100, -100, 300);
1822 fHPhotons->SetXTitle("no.of photons/area before padding");
1823 fHPhotons->SetYTitle("no.of photons/area after padding");
1824 fHPhotons->SetDirectory(NULL);
1825 fHPhotons->UseCurrentStyle();
1826 fHPhotons->GetYaxis()->SetTitleOffset(1.25);
1827
1828 // plot of added NSB
1829 fHNSB = new TH1D("NSB/area","Distribution of added NSB/area",
1830 100, -100.0, 200.0);
1831 fHNSB->SetXTitle("no.of added NSB photons/area");
1832 fHNSB->SetYTitle("no.of pixels");
1833 fHNSB->SetDirectory(NULL);
1834 fHNSB->UseCurrentStyle();
1835 fHNSB->GetYaxis()->SetTitleOffset(1.25);
1836
1837 //--------------------------------------------------------------------
1838
1839 *fLog << inf << "Type of data to be padded : " << fType << endl;
1840 *fLog << inf << "Name of MPedPhotCam container : " << fNamePedPhotCam
1841 << endl;
1842
1843 fIter = 10;
1844
1845 memset(fInf, 0, sizeof(fInf));
1846 memset(fErrors, 0, sizeof(fErrors));
1847 memset(fWarnings, 0, sizeof(fWarnings));
1848
1849 return kTRUE;
1850}
1851
1852// --------------------------------------------------------------------------
1853//
1854// Do the Padding (noise adjustment)
1855//
1856// input for the padding :
1857// - the matrices fHgON, fHgOFF, fHgMC
1858// - the original distributions fHSigmaTheta, fHDiffPixTheta
1859//
1860
1861Int_t MPad::Process()
1862{
1863 //*fLog << all << "Entry MPad::Process();" << endl;
1864
1865 // this is the conversion factor from the number of photons
1866 // to the number of photo electrons
1867 // later to be taken from MCalibration
1868 // fPEperPhoton = PW * LG * QE * 1D
1869 Double_t fPEperPhoton = 0.92 * 0.94 * 0.23 * 0.9; // (= 0.179)
1870
1871 //------------------------------------------------
1872 // Add pixels to MCerPhotEvt which are not yet in;
1873 // set their number of photons equal to zero.
1874 // Purpose : by the padding the number of photons is changed
1875 // so that cleaning and cuts may be affected.
1876 // However, no big effect is expectec unless the padding is strong
1877 // (strong increase of the noise level)
1878 // At present comment out this code
1879
1880 /*
1881 UInt_t imaxnumpix = fCam->GetNumPixels();
1882
1883 for (UInt_t i=0; i<imaxnumpix; i++)
1884 {
1885 Bool_t alreadythere = kFALSE;
1886 UInt_t npix = fEvt->GetNumPixels();
1887 for (UInt_t j=0; j<npix; j++)
1888 {
1889 MCerPhotPix &pix = (*fEvt)[j];
1890 UInt_t id = pix.GetPixId();
1891 if (i==id)
1892 {
1893 alreadythere = kTRUE;
1894 break;
1895 }
1896 }
1897 if (alreadythere)
1898 continue;
1899
1900 fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms());
1901 }
1902 */
1903
1904
1905 //-----------------------------------------
1906 Int_t rc=0;
1907 Int_t rw=0;
1908
1909 const UInt_t npix = fEvt->GetNumPixels();
1910
1911 //*fLog << all << "MPad::Process(); before padding : " << endl;
1912
1913 //-------------------------------------------
1914 // Calculate average sigma of the event
1915 //
1916 fPed->ReCalc(*fCam, fBad);
1917 //*fLog << "pedestalRMS, inner and outer = " << (fPed->GetArea(0)).GetRms()
1918 // << ", " << (fPed->GetArea(1)).GetRms() << endl;
1919
1920
1921 // inner sigmabar/sqrt(area)
1922 Double_t ratioinn = fCam->GetPixRatio(0);
1923 Double_t sigbarInnerold_phot = (fPed->GetArea(0)).GetRms();
1924 Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton * sqrt(ratioinn);
1925 Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold;
1926
1927 // outer sigmabar/sqrt(area)
1928 Double_t ratioout = fCam->GetPixRatio(500);
1929 Double_t sigbarOuterold_phot = (fPed->GetArea(1)).GetRms();
1930 Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton * sqrt(ratioout);
1931 Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold;
1932
1933 const Double_t theta = fPointPos->GetZd();
1934
1935 //*fLog << all << "theta = " << theta << endl;
1936
1937 Int_t binTheta = fHSigmaTheta->GetXaxis()->FindBin(theta);
1938 if ( binTheta < 1 || binTheta > fHSigmaTheta->GetNbinsX() )
1939 {
1940 *fLog << warn
1941 << "MPad::Process(); binNumber out of range : theta, binTheta = "
1942 << theta << ", " << binTheta << endl;
1943
1944 rc = 2;
1945 fErrors[rc]++;
1946 return kCONTINUE;
1947 }
1948
1949 //*fLog << all << "binTheta = " << binTheta << endl;
1950
1951
1952 //-------------------------------------------
1953 // get number of events in this theta bin (nTheta)
1954 // and number of events in this sigbarInnerold_phot bin (nSigma)
1955
1956 Double_t nTheta;
1957 Double_t nSigma;
1958
1959 TH2D *st = NULL;
1960 if (fType == "MC") st = fHSigmaThetaMC;
1961 else if (fType == "ON") st = fHSigmaThetaON;
1962 else if (fType == "OFF") st = fHSigmaThetaOFF;
1963 else
1964 {
1965 *fLog << err << "MPad : illegal data type '" << fType
1966 << "'" << endl;
1967 return kERROR;
1968 }
1969
1970 TH1D *hn;
1971 hn = st->ProjectionY("", binTheta, binTheta, "");
1972
1973 //*fLog << "Title of histogram : " << st->GetTitle() << endl;
1974 //for (Int_t i=1; i<=hn->GetNbinsX(); i++)
1975 //{
1976 // *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i)
1977 // << endl;
1978 //}
1979
1980 nTheta = hn->Integral();
1981 Int_t binSigma = hn->FindBin(sigbarInnerold_phot);
1982 nSigma = hn->GetBinContent(binSigma);
1983
1984 //*fLog << all
1985 // << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
1986 // << theta << ", " << sigbarInnerold_phot << ", " << binTheta
1987 // << ", "
1988 // << binSigma << ", " << nTheta << ", " << nSigma << endl;
1989
1990 delete hn;
1991
1992 //-------------------------------------------
1993
1994 //******************************************************************
1995 // has the event to be padded ?
1996 // if yes : to which sigmabar should it be padded ?
1997 //
1998
1999 TH3D *Hg = NULL;
2000 if (fType == "MC") Hg = fHgMC;
2001 else if (fType == "ON") Hg = fHgON;
2002 else if (fType == "OFF") Hg = fHgOFF;
2003 else
2004 {
2005 *fLog << err << "MPad : illegal type of data '" << fType << "'"
2006 << endl;
2007 return kERROR;
2008 }
2009
2010 Int_t binSig = Hg->GetYaxis()->FindBin(sigbarInnerold_phot);
2011 //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", "
2012 // << sigbarInnerold_phot << endl;
2013
2014 Double_t prob;
2015 TH1D *hpad = NULL;
2016
2017
2018 hpad = Hg->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");
2019
2020 //Int_t nb = hpad->GetNbinsX();
2021 //for (Int_t i=1; i<=nb; i++)
2022 //{
2023 // if (hpad->GetBinContent(i) != 0.0)
2024 // *fLog << all << "i, Hg = " << i << ", "
2025 // << hpad->GetBinContent(i) << endl;
2026 //}
2027
2028 if (nSigma != 0.0)
2029 prob = hpad->Integral() * nTheta/nSigma;
2030 else
2031 prob = 0.0;
2032
2033 //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma
2034 // << ", " << prob << endl;
2035
2036 if ( prob <= 0.0 || gRandom->Uniform() > prob )
2037 {
2038 delete hpad;
2039 // event should not be padded
2040 //*fLog << all << "event is not padded" << endl;
2041
2042 rc = 0;
2043 fInf[rc]++;
2044 return kTRUE;
2045 }
2046 // event should be padded
2047 //*fLog << all << "event will be padded" << endl;
2048
2049
2050 //-------------------------------------------
2051 // for the current theta, generate a sigmabar_inner :
2052 // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF
2053 //
2054 Double_t sigmabarInner_phot = 0;
2055 Double_t sigmabarInner = 0;
2056
2057
2058 //Int_t nbinsx = hpad->GetNbinsX();
2059 //for (Int_t i=1; i<=nbinsx; i++)
2060 //{
2061 // if (hpad->GetBinContent(i) != 0.0)
2062 // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i)
2063 // << endl;
2064 //}
2065
2066 sigmabarInner_phot = hpad->GetRandom();
2067 sigmabarInner = sigmabarInner_phot * fPEperPhoton * sqrt(ratioinn);
2068
2069 //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl;
2070
2071 delete hpad;
2072
2073 // new inner sigmabar2/area
2074 const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner;
2075
2076 //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = "
2077 // << sigbarInnerold << ", "<< sigmabarInner << endl;
2078
2079 // Stop if target sigmabar is <= sigbarold
2080 if (sigmabarInner <= sigbarInnerold)
2081 {
2082 *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : "
2083 << sigmabarInner << ", " << sigbarInnerold << ", aborting"
2084 << endl;
2085
2086 rc = 4;
2087 fErrors[rc]++;
2088 return kCONTINUE;
2089 }
2090
2091 //-------------------------------------------
2092 // estimate a sigmabar_outer from sigmabar_inner :
2093 // using the target distributions fHSigmaTheta and fHSigmaThetaOuter
2094 // for sigmabar(inner) and sigmabar(outer)
2095
2096 Double_t innerMeantarget = 0.0;
2097 Double_t outerMeantarget = 0.0;
2098 Double_t innerRMStarget = 0.0;
2099 Double_t outerRMStarget = 0.0;
2100
2101 // projection doesn't work
2102 //TH1D *hi = fHSigmaTheta->ProjectionY("proinner", binTheta, binTheta, "e");
2103 //TH1D *ho = fHSigmaThetaOuter->ProjectionY("proouter", binTheta, binTheta, "e");
2104 //sigmabarOuter2 = sigmabarInner2 + fPEperPhoton*fPEperPhoton *
2105 // ( ho->GetMean()*ho->GetMean()*ratioout
2106 // - hi->GetMean()*hi->GetMean()*ratioinn);
2107
2108 //innerMeantarget = hi->GetMean()*sqrt(ratioinn)*fPEperPhoton;
2109 //outerMeantarget = ho->GetMean()*sqrt(ratioout)*fPEperPhoton;
2110 //innerRMStarget = hi->GetRMS(1)*sqrt(ratioinn)*fPEperPhoton;
2111 //outerRMStarget = ho->GetRMS(1)*sqrt(ratioout)*fPEperPhoton;
2112
2113 Double_t y = 0.0;
2114 Double_t ybar = 0.0;
2115 Double_t y2bar = 0.0;
2116 Double_t w = 0.0;
2117
2118 Double_t isum = 0.0;
2119 Double_t isumy = 0.0;
2120 Double_t isumy2 = 0.0;
2121 for (Int_t i=1; i<=fHSigmaTheta->GetNbinsY(); i++)
2122 {
2123 w = fHSigmaTheta->GetBinContent(binTheta, i);
2124 y = fHSigmaTheta->GetYaxis()->GetBinCenter(i);
2125 isum += w;
2126 isumy += w * y;
2127 isumy2 += w * y*y;
2128 }
2129 if (isum == 0.0)
2130 {
2131 innerMeantarget = 0.0;
2132 innerRMStarget = 0.0;
2133 }
2134 else
2135 {
2136 ybar = isumy /isum;
2137 y2bar = isumy2/isum;
2138 innerMeantarget = ybar * sqrt(ratioinn)*fPEperPhoton;
2139 innerRMStarget = sqrt( y2bar - ybar*ybar ) * sqrt(ratioinn)*fPEperPhoton;
2140 }
2141
2142 Double_t osum = 0.0;
2143 Double_t osumy = 0.0;
2144 Double_t osumy2 = 0.0;
2145 for (Int_t i=1; i<=fHSigmaThetaOuter->GetNbinsY(); i++)
2146 {
2147 w = fHSigmaThetaOuter->GetBinContent(binTheta, i);
2148 y = fHSigmaThetaOuter->GetYaxis()->GetBinCenter(i);
2149 osum += w;
2150 osumy += w * y;
2151 osumy2 += w * y*y;
2152 }
2153 if (osum == 0.0)
2154 {
2155 outerMeantarget = 0.0;
2156 outerRMStarget = 0.0;
2157 }
2158 else
2159 {
2160 ybar = osumy /osum;
2161 y2bar = osumy2/osum;
2162 outerMeantarget = ybar * sqrt(ratioout)*fPEperPhoton;
2163 outerRMStarget = sqrt( y2bar - ybar*ybar ) * sqrt(ratioout)*fPEperPhoton;
2164 }
2165
2166
2167 // new outer sigmabar2/area
2168 Double_t sigmabarOuter2;
2169
2170 Double_t scal = ( innerMeantarget*innerRMStarget == 0.0 ||
2171 outerMeantarget*outerRMStarget == 0.0 ) ? 1.0 :
2172 (outerMeantarget*outerRMStarget)/(innerMeantarget*innerRMStarget);
2173 sigmabarOuter2 = outerMeantarget*outerMeantarget +
2174 scal * (sigmabarInner2 - innerMeantarget*innerMeantarget);
2175
2176 //*fLog << "innerMeantarget, innerRMStarget = " << innerMeantarget
2177 // << ", " << innerRMStarget << endl;
2178
2179 //*fLog << "outerMeantarget, outerRMStarget = " << outerMeantarget
2180 // << ", " << outerRMStarget << endl;
2181
2182 //*fLog << "sigmabarInner2, sigmabarOuter2, scal = " << sigmabarInner2
2183 // << ", " << sigmabarOuter2 << ", " << scal << endl;
2184
2185 //delete hi;
2186 //delete ho;
2187
2188 //-------------------------------------------
2189 //
2190 // Calculate average number of NSB photo electrons to be added (lambdabar)
2191 // from the value of sigmabar,
2192 // - using a fixed value (F2excess) for the excess noise factor
2193
2194 Double_t elNoise2; // [photo electrons]
2195 Double_t F2excess = 1.3;
2196 Double_t lambdabar; // [photo electrons]
2197
2198 //----------------
2199 TH3D *sp = NULL;
2200 if (fType == "MC") sp = fHDiffPixThetaTargetMC;
2201 else if (fType == "ON") sp = fHDiffPixThetaTargetON;
2202 else if (fType == "OFF") sp = fHDiffPixThetaTargetOFF;
2203
2204 //----------------
2205
2206 Int_t bincheck = sp->GetXaxis()->FindBin(theta);
2207 if (binTheta != bincheck)
2208 {
2209 *fLog << err
2210 << "MPad::Process(); The binnings of the 2 histograms are not identical; aborting"
2211 << endl;
2212 return kERROR;
2213 }
2214
2215 // In this Theta bin, get RMS of (Sigma^2-sigmabar^2)/area for inner pixels.
2216 // The average electronic noise (to be added) has to be in the order of
2217 // this RMS, otherwise the electronic noise of an individual pixel
2218 // (elNoise2Pix) may become negative
2219
2220 //----------------
2221
2222
2223 // Attention : maximum ID of inner pixels hard coded !!!
2224 Int_t idmaxpixinner = 396;
2225 Int_t binmaxpixinner =
2226 sp->GetYaxis()->FindBin( (Double_t)idmaxpixinner );
2227
2228 TH1D *hnoise = NULL;
2229 hnoise = sp->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
2230
2231 Double_t RMS_phot = hnoise->GetRMS(1);
2232 Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
2233 delete hnoise;
2234
2235 elNoise2 = TMath::Min(2.0*RMS, sigmabarInner2 - sigbarInnerold2);
2236 //*fLog << all << "RMS_phot, elNoise2 = " << RMS_phot << ", "
2237 // << elNoise2 << endl;
2238
2239 lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess;
2240
2241 if (lambdabar <= 0.0)
2242 {
2243 rc = 3;
2244 fErrors[rc]++;
2245 }
2246
2247 //*fLog << "lambdabar = " << lambdabar << endl;
2248
2249 // This value of lambdabar is the same for all pixels;
2250 // note that lambdabar is the NSB p.e. density
2251
2252 //---------- start loop over pixels ---------------------------------
2253 // do the padding for each pixel
2254 //
2255 // pad only pixels - which are used
2256 //
2257
2258 Double_t sigma2 = 0;
2259
2260 Double_t diff_phot = 0;
2261 Double_t diff = 0;
2262
2263 Double_t addSig2_phot= 0;
2264 Double_t addSig2 = 0;
2265
2266 Double_t elNoise2Pix = 0;
2267
2268
2269 // throw the Sigma for the pixels from the distribution fHDiffPixTheta
2270 // MC : from fHDiffPixThetaTargetMC
2271 // ON : from fHDiffPixThetaTaregtON
2272 // OFF : from fHDiffPixThetaTargetOFF
2273
2274
2275 Double_t sigbarold2;
2276 Double_t sigmabar2;
2277
2278 for (UInt_t i=0; i<npix; i++)
2279 {
2280
2281 MCerPhotPix &pix = (*fEvt)[i];
2282 if ( !pix.IsPixelUsed() )
2283 continue;
2284
2285 //if ( pix.GetNumPhotons() == 0.0)
2286 //{
2287 // *fLog << warn
2288 // << "MPad::Process(); no.of photons is 0 for used pixel"
2289 // << endl;
2290 // continue;
2291 //}
2292
2293 Int_t j = pix.GetPixId();
2294
2295 // GetPixRatio returns (area of pixel 0 / area of current pixel)
2296 Double_t ratio = fCam->GetPixRatio(j);
2297
2298 if (ratio > 0.5)
2299 {
2300 sigbarold2 = sigbarInnerold2;
2301 sigmabar2 = sigmabarInner2;
2302 }
2303 else
2304 {
2305 sigbarold2 = sigbarOuterold2;
2306 sigmabar2 = sigmabarOuter2;
2307 }
2308
2309 MPedPhotPix &ppix = (*fPed)[j];
2310
2311 // count number of pixels treated
2312 fWarnings[0]++;
2313
2314
2315 Double_t oldsigma_phot = ppix.GetRms();
2316 Double_t oldsigma = oldsigma_phot * fPEperPhoton * sqrt(ratio);
2317 Double_t oldsigma2 = oldsigma*oldsigma;
2318
2319 //---------------------------------
2320 // throw the Sigma for this pixel
2321 //
2322 Int_t binPixel = sp->GetYaxis()->FindBin( (Double_t)j );
2323
2324 Int_t count;
2325 Bool_t ok;
2326
2327 TH1D *hdiff = NULL;
2328
2329 hdiff = sp->ProjectionZ("", binTheta, binTheta,
2330 binPixel, binPixel, "");
2331 Double_t integral = hdiff->Integral();
2332 // if there are no entries in hdiff, diff cannot be thrown
2333 // in this case diff will be set to the old difference
2334 if ( integral == 0 )
2335 {
2336 //*fLog << warn << "MPad::Process(); fType = " << fType
2337 // << ", projection of '"
2338 // << sp->GetName() << "' for Theta bin "
2339 // << binTheta << " and pixel " << j
2340 // << " has no entries; set diff equal to the old difference "
2341 // << endl;
2342
2343 diff = TMath::Max(oldsigma2 - sigbarold2,
2344 lambdabar*F2excess + oldsigma2 - sigmabar2);
2345
2346 rc = 2;
2347 fWarnings[rc]++;
2348 }
2349
2350 // start of else -------------------
2351 else
2352 {
2353 count = 0;
2354 ok = kFALSE;
2355 for (Int_t m=0; m<fIter; m++)
2356 {
2357 count += 1;
2358 diff_phot = hdiff->GetRandom();
2359
2360 //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : "
2361 // << m << ", " << diff_phot << endl;
2362
2363 diff = diff_phot * fPEperPhoton * fPEperPhoton;
2364
2365 // the following condition ensures that elNoise2Pix > 0.0
2366 if ( (diff + sigmabar2 - oldsigma2
2367 - lambdabar*F2excess) > 0.0 )
2368 {
2369 ok = kTRUE;
2370 break;
2371 }
2372 }
2373
2374 if (!ok)
2375 {
2376 //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratio, lambdabar = "
2377 // << theta << ", "
2378 // << j << ", " << count << ", " << sigmabar2 << ", "
2379 // << diff << ", " << oldsigma2 << ", " << ratio << ", "
2380 // << lambdabar << endl;
2381 diff = lambdabar*F2excess + oldsigma2 - sigmabar2;
2382
2383 rw = 1;
2384 fWarnings[rw]++;
2385 }
2386 }
2387 // end of else --------------------
2388
2389 delete hdiff;
2390 sigma2 = diff + sigmabar2;
2391
2392
2393 //---------------------------------
2394 // get the additional sigma^2 for this pixel (due to the padding)
2395
2396 addSig2 = (sigma2 - oldsigma2) / ratio;
2397 addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);
2398
2399 //---------------------------------
2400 // get the additional electronic noise for this pixel
2401
2402 elNoise2Pix = addSig2 - lambdabar*F2excess/ratio;
2403
2404
2405 //---------------------------------
2406 // throw actual number of additional NSB photo electrons (NSB)
2407 // and its RMS (sigmaNSB)
2408
2409 Double_t NSB0 = gRandom->Poisson(lambdabar/ratio);
2410 Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix;
2411 Double_t sigmaNSB0;
2412
2413 if (arg >= 0.0)
2414 {
2415 sigmaNSB0 = sqrt( arg );
2416 }
2417 else
2418 {
2419 sigmaNSB0 = 0.0000001;
2420 if (arg < -1.e-10)
2421 {
2422 *fLog << warn << "MPad::Process(); argument of sqrt < 0 : "
2423 << arg << endl;
2424 }
2425 }
2426
2427
2428 //---------------------------------
2429 // smear NSB0 according to sigmaNSB0
2430 // and subtract lambdabar because of AC coupling
2431
2432 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar/ratio;
2433 Double_t NSB_phot = NSB / fPEperPhoton;
2434
2435 //---------------------------------
2436
2437 // add additional NSB to the number of photons
2438 Double_t oldphotons_phot = pix.GetNumPhotons();
2439 Double_t oldphotons = oldphotons_phot * fPEperPhoton;
2440 Double_t newphotons = oldphotons + NSB;
2441 Double_t newphotons_phot = newphotons / fPEperPhoton;
2442 pix.SetNumPhotons( newphotons_phot );
2443
2444
2445 fHNSB->Fill( NSB_phot*ratio );
2446 fHPhotons->Fill( oldphotons_phot*ratio,
2447 newphotons_phot*ratio );
2448
2449
2450 // error: add sigma of padded noise quadratically
2451 Double_t olderror_phot = pix.GetErrorPhot();
2452 Double_t newerror_phot =
2453 sqrt( olderror_phot*olderror_phot + addSig2_phot );
2454 pix.SetErrorPhot( newerror_phot );
2455
2456
2457 Double_t newsigma = sqrt( sigma2 / ratio );
2458 Double_t newsigma_phot = newsigma / fPEperPhoton;
2459 ppix.SetRms( newsigma_phot );
2460
2461 fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot );
2462 }
2463 //---------- end of loop over pixels ---------------------------------
2464
2465
2466 //*fLog << all << "MPad::Process(); after padding : " << endl;
2467
2468 // Calculate sigmabar again and crosscheck
2469 fPed->ReCalc(*fCam, fBad);
2470 //*fLog << "pedestalRMS, inner and outer = " << (fPed->GetArea(0)).GetRms()
2471 // << ", " << (fPed->GetArea(1)).GetRms() << endl;
2472
2473 //*fLog << all << "Exit MPad::Process();" << endl;
2474
2475 rc = 0;
2476 fErrors[rc]++;
2477
2478 return kTRUE;
2479 //******************************************************************
2480}
2481
2482// --------------------------------------------------------------------------
2483//
2484//
2485Int_t MPad::PostProcess()
2486{
2487 if (GetNumExecutions() != 0)
2488 {
2489
2490 *fLog << inf << endl;
2491 *fLog << GetDescriptor() << " execution statistics:" << endl;
2492 *fLog << dec << setfill(' ');
2493
2494 if (fWarnings[0] == 0 ) fWarnings[0] = 1;
2495
2496 *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3)
2497 << (int)(fWarnings[1]*100/fWarnings[0])
2498 << "%) Pixel: iteration to find acceptable sigma failed"
2499 << ", fIter = " << fIter << endl;
2500
2501 *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3)
2502 << (int)(fWarnings[2]*100/fWarnings[0])
2503 << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << endl;
2504
2505
2506 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3)
2507 << (int)(fErrors[2]*100/GetNumExecutions())
2508 << "%) Evts skipped due to: Zenith angle out of range" << endl;
2509
2510 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3)
2511 << (int)(fErrors[4]*100/GetNumExecutions())
2512 << "%) Evts skipped due to: new sigma <= old sigma" << endl;
2513
2514 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3)
2515 << (int)(fErrors[3]*100/GetNumExecutions())
2516 << "%) lambdabar = 0" << endl;
2517
2518 *fLog << " " << setw(7) << fInf[0] << " (" << setw(3)
2519 << (int)(fInf[0]*100/GetNumExecutions())
2520 << "%) Evts didn't have to be padded" << endl;
2521
2522 *fLog << " " << fErrors[0] << " ("
2523 << (int)(fErrors[0]*100/GetNumExecutions())
2524 << "%) Evts were successfully padded!" << endl;
2525 *fLog << endl;
2526
2527 }
2528
2529 //---------------------------------------------------------------
2530 TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1200));
2531 c.Divide(3, 4);
2532 gROOT->SetSelectedPad(NULL);
2533
2534 c.cd(1);
2535 fHNSB->DrawCopy();
2536 fHNSB->SetBit(kCanDelete);
2537
2538 c.cd(2);
2539 fHSigmaPedestal->DrawCopy();
2540 fHSigmaPedestal->SetBit(kCanDelete);
2541
2542 c.cd(3);
2543 fHPhotons->DrawCopy();
2544 fHPhotons->SetBit(kCanDelete);
2545
2546 //--------------------------------------------------------------------
2547
2548 if (fHgON)
2549 {
2550 c.cd(4);
2551 TH2D *m1;
2552 m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
2553 m1->SetDirectory(NULL);
2554 m1->UseCurrentStyle();
2555 m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
2556 m1->SetXTitle("Sigmabar-old");
2557 m1->SetYTitle("Sigmabar-new");
2558
2559 m1->DrawCopy("box");
2560 m1->SetBit(kCanDelete);;
2561 }
2562
2563 if (fHgOFF)
2564 {
2565 c.cd(5);
2566 TH2D *m2;
2567 m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
2568 m2->SetDirectory(NULL);
2569 m2->UseCurrentStyle();
2570 m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
2571 m2->SetXTitle("Sigmabar-old");
2572 m2->SetYTitle("Sigmabar-new");
2573
2574 m2->DrawCopy("box");
2575 m2->SetBit(kCanDelete);;
2576 }
2577
2578 if (fHgMC)
2579 {
2580 c.cd(6);
2581 TH2D *m3;
2582 m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy");
2583 m3->SetDirectory(NULL);
2584 m3->UseCurrentStyle();
2585 m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)");
2586 m3->SetXTitle("Sigmabar-old");
2587 m3->SetYTitle("Sigmabar-new");
2588
2589 m3->DrawCopy("box");
2590 m3->SetBit(kCanDelete);;
2591 }
2592
2593 //--------------------------------------------------------------------
2594
2595 if (fHgOuterON)
2596 {
2597 c.cd(7);
2598 TH2D *m1;
2599 m1 = (TH2D*) ((TH3*)fHgOuterON)->Project3D("zy");
2600 m1->SetDirectory(NULL);
2601 m1->UseCurrentStyle();
2602 m1->SetTitle("(fHgOuterON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)");
2603 m1->SetXTitle("Sigmabar-old");
2604 m1->SetYTitle("Sigmabar-new");
2605
2606 m1->DrawCopy("box");
2607 m1->SetBit(kCanDelete);;
2608 }
2609
2610 if (fHgOuterOFF)
2611 {
2612 c.cd(8);
2613 TH2D *m2;
2614 m2 = (TH2D*) ((TH3*)fHgOuterOFF)->Project3D("zy");
2615 m2->SetDirectory(NULL);
2616 m2->UseCurrentStyle();
2617 m2->SetTitle("(fHgOuterOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)");
2618 m2->SetXTitle("Sigmabar-old");
2619 m2->SetYTitle("Sigmabar-new");
2620
2621 m2->DrawCopy("box");
2622 m2->SetBit(kCanDelete);;
2623 }
2624
2625 if (fHgOuterMC)
2626 {
2627 c.cd(9);
2628 TH2D *m3;
2629 m3 = (TH2D*) ((TH3*)fHgOuterMC)->Project3D("zy");
2630 m3->SetDirectory(NULL);
2631 m3->UseCurrentStyle();
2632 m3->SetTitle("(fHgOuterMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)");
2633 m3->SetXTitle("Sigmabar-old");
2634 m3->SetYTitle("Sigmabar-new");
2635
2636 m3->DrawCopy("box");
2637 m3->SetBit(kCanDelete);;
2638 }
2639
2640 //--------------------------------------------------------------------
2641
2642 c.cd(10);
2643 fHSigmaTheta->SetDirectory(NULL);
2644 fHSigmaTheta->UseCurrentStyle();
2645 fHSigmaTheta->DrawCopy();
2646 fHSigmaTheta->SetBit(kCanDelete);
2647
2648 c.cd(11);
2649 fHSigmaThetaOuter->SetDirectory(NULL);
2650 fHSigmaThetaOuter->UseCurrentStyle();
2651 fHSigmaThetaOuter->DrawCopy();
2652 fHSigmaThetaOuter->SetBit(kCanDelete);
2653
2654
2655 return kTRUE;
2656}
2657
2658// --------------------------------------------------------------------------
2659
2660
2661
2662
2663
2664
2665
2666
2667
2668
Note: See TracBrowser for help on using the repository browser.