Index: trunk/MagicSoft/Mars/macros/ONAnalysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONAnalysis.C (revision 2776) +++ trunk/MagicSoft/Mars/macros/ONAnalysis.C (revision 2777) @@ -761,14 +761,9 @@ contbasic.SetName("SelBasic"); - - // There are 2 options for Thomas Schweizer's padding - // fPadFlag = 1 get Sigmabar from fHSigmaTheta - // and Sigma from fHDiffPixTheta - // fPadFlag = 2 get Sigma from fHSigmaPixTheta - - MPadON padthomas("MPadON","Task for the padding"); - padthomas.SetHistograms(fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta, - fHIdBlindPixels, fHNBlindPixels); - padthomas.SetPadFlag(1); + MPad padthomas("MPad","Task for the padding"); + padthomas.MergeONMC(fHSigmaThetaMC, fHDiffPixThetaMC, fHSigmaPixThetaMC, + fHIdBlindPixelsMC, fHNBlindPixelsMC, + fHSigmaThetaON, fHDiffPixThetaON, fHSigmaPixThetaON, + fHIdBlindPixelsON, fHNBlindPixelsON); MFillH fillblind("MCBlindPixels[MHBlindPixels]", "MBlindPixels"); Index: trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 2776) +++ trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 2777) @@ -357,7 +357,6 @@ // and merge these histograms - MPadONOFF pad; - pad.SetName("MPadONOFF"); - pad.SetPadFlag(1); + MPad pad; + pad.SetName("MPad"); pad.SetDataType(typeInput); @@ -531,14 +530,16 @@ gLog << "=====================================================" << endl; - pad.MergeHistograms(sigthon, sigthoff, - sigpixthon, sigpixthoff, - diffpixthon, diffpixthoff, - blindidthon, blindidthoff, - blindnthon, blindnthoff); + pad.MergeONOFFMC(sigthmc, diffpixthmc, sigpixthmc, + blindidthmc, blindnthmc, + sigthon, diffpixthon, sigpixthon, + blindidthon, blindnthon, + sigthoff, diffpixthoff, sigpixthoff, + blindidthoff, blindnthoff); + if (WPad) { // write the padding histograms onto a file --------- - pad.WriteTargetDist(outNameSigTh); + pad.WritePaddingDist(outNameSigTh); } } @@ -547,5 +548,5 @@ if (RPad) { - pad.ReadTargetDist(outNameSigTh); + pad.ReadPaddingDist(outNameSigTh); } Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2776) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2777) @@ -48,6 +48,5 @@ #pragma link C++ class MCT1PadSchweizer+; #pragma link C++ class MCT1PadONOFF+; -#pragma link C++ class MPadON+; -#pragma link C++ class MPadONOFF+; +#pragma link C++ class MPad+; #pragma link C++ class MCT1PointingCorrCalc+; Index: trunk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 2777) +++ trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 2777) @@ -0,0 +1,2369 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Wolfgang Wittek, 06/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPad +// +// The aim of the padding is to make the noise level (NSB + electronic noise) +// equal for the MC, ON and OFF data samples +// +// The target noise level is determined by 'merging' the (sigmabar vs. Theta) +// distributions of MC, ON and OFF data. +// +// The prescription on which fraction of events has to padded from +// bin k of sigmabar to bin j is stored in the matrices +// fHgMC, fHgON and fHgOFF. +// +// By the padding the number of photons, its error and the pedestal sigmas +// are altered. On average, the number of photons added is zero. +// +// The formulas used can be found in Thomas Schweizer's Thesis, +// Section 2.2.1 +// +// The padding is done as follows : +// +// Increase the sigmabar according to the matrices fHg.... Then generate +// a pedestal sigma for each pixel using the 3D-histogram Theta, pixel no., +// Sigma^2-Sigmabar^2 (fHDiffPixTheta). +// +// The padding has to be done before the image cleaning because the +// image cleaning depends on the pedestal sigmas. +// +///////////////////////////////////////////////////////////////////////////// +#include "MPad.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MSigmabar.h" +#include "MMcEvt.hxx" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MBlindPixels.h" + +#include "MRead.h" +#include "MFilterList.h" +#include "MTaskList.h" +#include "MBlindPixelCalc.h" +#include "MHBlindPixels.h" +#include "MFillH.h" +#include "MHSigmaTheta.h" +#include "MEvtLoop.h" + +ClassImp(MPad); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPad::MPad(const char *name, const char *title) +{ + fName = name ? name : "MPad"; + fTitle = title ? title : "Task for the padding"; + + fType = ""; + + fHSigmaThetaMC = NULL; + fHSigmaThetaON = NULL; + fHSigmaThetaOFF = NULL; + + fHDiffPixThetaMC = NULL; + fHDiffPixThetaON = NULL; + fHDiffPixThetaOFF = NULL; + + fHSigmaPixThetaMC = NULL; + fHSigmaPixThetaON = NULL; + fHSigmaPixThetaOFF = NULL; + + fHgMC = NULL; + fHgON = NULL; + fHgOFF = NULL; + + fHBlindPixIdThetaMC = NULL; + fHBlindPixIdThetaON = NULL; + fHBlindPixIdThetaOFF = NULL; + + fHBlindPixNThetaMC = NULL; + fHBlindPixNThetaON = NULL; + fHBlindPixNThetaOFF = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPad::~MPad() +{ + delete fHgMC; + delete fHgON; + delete fHgOFF; + + delete fHSigmaTheta; + delete fHSigmaPixTheta; + delete fHDiffPixTheta; + delete fHBlindPixNTheta; + delete fHBlindPixIdTheta; + + delete fHSigmaPedestal; + delete fHPhotons; + delete fHNSB; + + delete fInfile; +} +// -------------------------------------------------------------------------- +// +// Merge the sigmabar distributions of 2 samples (samples A and B) +// +// input : the original distributions for samples A and B (hista, histb) +// +// output : the prescription which fraction of events should be padded +// from the sigmabar bin k to bin j (fHgMC, fHgON, fHgOFF) +// + +Bool_t MPad::Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, + TH2D *fHga, TH2D *fHgb, Int_t nbinssig ) +{ + // hista is the normalized 1D histogram of sigmabar for sample A + // histb is the normalized 1D histogram of sigmabar for sample B + // histc is the difference A-B + + // at the beginning, histap is a copy of hista + // at the end, it will be the 1D histogram for sample A after padding + + // at the beginning, histbp is a copy of histb + // at the end, it will be the 1D histogram for sample B after padding + + // at the beginning, histcp is a copy of histc + // at the end, it should be the difference histap-histbp, + // which should be zero + + // fHga[k][j] tells which fraction of events from sample A should be padded + // from sigma_k to sigma_j + + + // fHgb[k][j] tells which fraction of events from sample B should be padded + // from sigma_k to sigma_j + + Double_t eps = 1.e-10; + + TH1D *histbp = new TH1D( (const TH1D&)*histb ); + + TH1D *histc = new TH1D( (const TH1D&)*hista ); + histc->Add(histb, -1.0); + + TH1D *histcp = new TH1D( (const TH1D&)*histc ); + + + //-------- start j loop ------------------------------------------------ + // loop over bins in histc, + // starting from the end (i.e. at the largest sigmabar) + Double_t v, s, w, t, x, u, a, b, arg; + + for (Int_t j=nbinssig; j >= 1; j--) + { + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) continue; + if (v >= 0.0) + s = 1.0; + else + s = -1.0; + + //................ start k loop ...................................... + // look for a bin k which may compensate the content of bin j + for (Int_t k=j-1; k >= 1; k--) + { + w = histcp->GetBinContent(k); + if ( fabs(w) < eps ) continue; + if (w >= 0.0) + t = 1.0; + else + t = -1.0; + + // if s==t bin k cannot compensate, go to next k bin + if (s == t) continue; + + x = v + w; + if (x >= 0.0) + u = 1.0; + else + u = -1.0; + + // if u==s bin k will partly compensate : pad the whole fraction + // w from bin k to bin j + if (u == s) + arg = -w; + + // if u!=s bin k would overcompensate : pad only the fraction + // v from bin k to bin j + else + arg = v; + + if (arg <=0.0) + { + fHga->SetBinContent(k, j, -arg); + fHgb->SetBinContent(k, j, 0.0); + } + else + { + fHga->SetBinContent(k, j, 0.0); + fHgb->SetBinContent(k, j, arg); + } + + *fLog << all << "k, j, arg = " << k << ", " << j + << ", " << arg << endl; + + //...................................... + // this is for checking the procedure + if (arg < 0.0) + { + a = histap->GetBinContent(k); + histap->SetBinContent(k, a+arg); + a = histap->GetBinContent(j); + histap->SetBinContent(j, a-arg); + } + else + { + b = histbp->GetBinContent(k); + histbp->SetBinContent(k, b-arg); + b = histbp->GetBinContent(j); + histbp->SetBinContent(j, b+arg); + } + //...................................... + + if (u == s) + { + histcp->SetBinContent(k, 0.0); + histcp->SetBinContent(j, x); + + v = histcp->GetBinContent(j); + if ( fabs(v) < eps ) break; + + // bin j was only partly compensated : go to next k bin + continue; + } + else + { + histcp->SetBinContent(k, x); + histcp->SetBinContent(j, 0.0); + + // bin j was completely compensated : go to next j bin + break; + } + + } + //................ end k loop ...................................... + } + //-------- end j loop ------------------------------------------------ + + // check results + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = histap->GetBinContent(j); + Double_t b = histbp->GetBinContent(j); + Double_t c = histcp->GetBinContent(j); + + if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) + *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = " + << j << ", " << a << ", " << b << ", " << c << endl; + } + + //--------------------------------------------------------------- + TCanvas &c = *(MH::MakeDefCanvas("Merging", "", 600, 600)); + c.Divide(2, 2); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + hista->SetDirectory(NULL); + hista->DrawCopy(); + hista->SetBit(kCanDelete); + + c.cd(2); + histb->SetDirectory(NULL); + histb->DrawCopy(); + histb->SetBit(kCanDelete); + + c.cd(3); + histap->SetDirectory(NULL); + histap->DrawCopy(); + histap->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + delete histc; + delete histbp; + delete histcp; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Merge the distributions +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// +// ===> of ON, OFF and MC data <============== +// +// and define the matrices fHgMC, fHgON and fHgOFF +// +// These matrices tell which fraction of events should be padded +// from bin k of sigmabar to bin j +// + +Bool_t MPad::MergeONOFFMC( + TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, + TH2D *blindidthmc, TH2D *blindnthmc, + TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, + TH2D *blindidthon, TH2D *blindnthon, + TH2D *sigthoff, TH3D *diffpixthoff, TH3D *sigmapixthoff, + TH2D *blindidthoff, TH2D *blindnthoff) +{ + *fLog << all << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding" + << endl; + + fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + //-------------------------- + fHSigmaThetaMC = sigthmc; + fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); + fHSigmaThetaON = sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); + fHSigmaThetaOFF = sigthoff; + fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)"); + + //-------------------------- + fHBlindPixNThetaMC = blindnthmc; + fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); + fHBlindPixNThetaON = blindnthon; + fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); + fHBlindPixNThetaOFF = blindnthoff; + fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)"); + + //-------------------------- + fHBlindPixIdThetaMC = blindidthmc; + fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); + fHBlindPixIdThetaON = blindidthon; + fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); + fHBlindPixIdThetaOFF = blindidthoff; + fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)"); + + //-------------------------- + fHDiffPixThetaMC = diffpixthmc; + fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); + fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); + fHDiffPixThetaOFF = diffpixthoff; + fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)"); + + //-------------------------- + fHSigmaPixThetaMC = sigmapixthmc; + fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); + fHSigmaPixThetaON = sigmapixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); + fHSigmaPixThetaOFF = sigmapixthoff; + fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)"); + //-------------------------- + + + // get binning for fHgON, fHgOFF and fHgMC from sigthon + // binning in Theta + TAxis *ax = sigthon->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon->GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgMC = new TH3D; + MH::SetBinning(fHgMC, &binth, &binsg, &binsg); + fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); + fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); + + //............ loop over Theta bins ........................... + + + + *fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + TH2D * fHga = new TH2D; + MH::SetBinning(fHga, &binsg, &binsg); + TH2D * fHgb = new TH2D; + MH::SetBinning(fHgb, &binsg, &binsg); + + //------------------------------------------- + // merge ON and OFF distributions + // input : hista is the normalized 1D distr. of sigmabar for ON data + // histb is the normalized 1D distr. of sigmabar for OFF data + // output : histap will be the merged distribution (ON-OFF) + // fHga(k,j) will tell which fraction of ON events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of OFF events should be + // padded from bin k to bin j + + TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + Stat_t suma = hista->Integral(); + hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + + TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); + Stat_t sumb = histb->Integral(); + histb->Scale(1./sumb); + + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); + delete hista; + delete histb; + + // fill fHgON and fHgOFF + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent(k,j); + fHga->SetBinContent (k, j, 0.0); + fHgON->SetBinContent(l, k, j, a); + + Double_t b = fHgb->GetBinContent(k,j); + fHgb->SetBinContent (k, j, 0.0); + fHgOFF->SetBinContent(l, k, j, b); + } + + //------------------------------------------- + // now merge ON-OFF and MC distributions + // input : hista is the result of the merging of ON and OFF distributions + // histb is the normalized 1D distr. of sigmabar for MC data + // output : histap will be the merged distribution (target distribution) + // fHga(k,j) will tell which fraction of ON, OFF events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of MC events should be + // padded from bin k to bin j + + hista = new TH1D( (const TH1D&)*histap); + + histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); + sumb = histb->Integral(); + histb->Scale(1./sumb); + + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); + delete hista; + delete histb; + + // update fHgON and fHgOFF and fill fHgMC + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent (k,j); + + Double_t aON = fHgON->GetBinContent(l,k,j); + fHgON->SetBinContent(l, k, j, aON+a); + + Double_t aOFF = fHgOFF->GetBinContent(l,k,j); + fHgOFF->SetBinContent(l, k, j, aOFF+a); + + Double_t b = fHgb->GetBinContent(k,j); + fHgMC->SetBinContent(l, k, j, b); + } + + // fill target distribution fHSigmaTheta + // (only for plotting, it will not be used in the padding) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = histap->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, a); + } + + delete fHga; + delete fHgb; + + delete histap; + } + //............ end of loop over Theta bins .................... + + + // Define the target distributions + // fHDiffPixTheta, fHSigmaPixTheta + // (they are calculated as + // averages of the ON and OFF distributions) + // fHBlindPixNTheta, fHBlindPixIdTheta + // (they are calculated as + // the sum of the ON and OFF distributions) + // (fHDiffPixTheta will be used in the padding of MC events) + + //............ start of new loop over Theta bins .................... + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + Double_t normON; + Double_t normOFF; + + // set ranges for 2D-projections of 3D-histograms + TAxis *ax = diffpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthoff->GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histOFF; + + TH1D *hist1; + TH1D *hist1OFF; + + // weights for ON and OFF distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracOFF); + Double_t facOFF = 0.5 * (1. + fracON - fracOFF); + + *fLog << all << "Theta bin j : fracON, fracOFF, facON, facOFF = " + << j << ", " << fracON << ", " << fracOFF << ", " + << fracMC << ", " << facON << ", " << facOFF << endl; + + //------------------------------------------------------------------ + // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' + ay = diffpixthon->GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + TAxis *az = diffpixthon->GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)diffpixthoff->Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------------------------------ + // define target distribution 'sigma vs. pixel number' + ay = sigmapixthon->GetYaxis(); + nbinspixel = ay->GetNbins(); + + az = sigmapixthon->GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigmapixthon->Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)sigmapixthoff->Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + + //------------------------------------------------------------------ + // define target distribution 'number of blind pixels per event' + ay = blindnthon->GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, ""); + + normON = hist1->Integral(); + normOFF = hist1OFF->Integral(); + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + hist1->Scale( 1.0/hist1->Integral() ); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------------------------------ + // define target distribution 'id of blind pixel' + ay = blindidthon->GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, ""); + + // divide by the number of events (from fHBlindPixNTheta) + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + } + //............ end of new loop over Theta bins .................... + + *fLog << all + << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaThetaMC->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOFF->SetDirectory(NULL); + + fHDiffPixThetaMC->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + + fHSigmaPixThetaMC->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + fHSigmaPixThetaOFF->SetDirectory(NULL); + + fHBlindPixIdThetaMC->SetDirectory(NULL); + fHBlindPixIdThetaON->SetDirectory(NULL); + fHBlindPixIdThetaOFF->SetDirectory(NULL); + + fHBlindPixNThetaMC->SetDirectory(NULL); + fHBlindPixNThetaON->SetDirectory(NULL); + fHBlindPixNThetaOFF->SetDirectory(NULL); + + fHgMC->SetDirectory(NULL); + fHgON->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Merge the distributions +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// +// ===> of ON and MC data <============== +// +// and define the matrices fHgMC and fHgON +// +// These matrices tell which fraction of events should be padded +// from bin k of sigmabar to bin j +// + +Bool_t MPad::MergeONMC( + TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, + TH2D *blindidthmc, TH2D *blindnthmc, + TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, + TH2D *blindidthon, TH2D *blindnthon) +{ + *fLog << all << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" + << endl; + + fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + //-------------------------- + fHSigmaThetaMC = sigthmc; + fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); + fHSigmaThetaON = sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); + + //-------------------------- + fHBlindPixNThetaMC = blindnthmc; + fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); + fHBlindPixNThetaON = blindnthon; + fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); + + //-------------------------- + fHBlindPixIdThetaMC = blindidthmc; + fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); + fHBlindPixIdThetaON = blindidthon; + fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); + + //-------------------------- + fHDiffPixThetaMC = diffpixthmc; + fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); + fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); + + //-------------------------- + fHSigmaPixThetaMC = sigmapixthmc; + fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); + fHSigmaPixThetaON = sigmapixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); + + //-------------------------- + + + // get binning for fHgON and fHgMC from sigthon + // binning in Theta + TAxis *ax = sigthon->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon->GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgMC = new TH3D; + MH::SetBinning(fHgMC, &binth, &binsg, &binsg); + fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + + //............ loop over Theta bins ........................... + + + + *fLog << all << "MPad::MergeONMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + TH2D * fHga = new TH2D; + MH::SetBinning(fHga, &binsg, &binsg); + TH2D * fHgb = new TH2D; + MH::SetBinning(fHgb, &binsg, &binsg); + + //------------------------------------------- + // merge ON and MC distributions + // input : hista is the normalized 1D distr. of sigmabar for ON data + // histb is the normalized 1D distr. of sigmabar for MC data + // output : histap will be the merged distribution (ON-MC) + // fHga(k,j) will tell which fraction of ON events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of MC events should be + // padded from bin k to bin j + + TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + Stat_t suma = hista->Integral(); + hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + + TH1D *histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); + Stat_t sumb = histb->Integral(); + histb->Scale(1./sumb); + + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); + delete hista; + delete histb; + + // fill fHgON and fHgMC + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent(k,j); + fHga->SetBinContent (k, j, 0.0); + fHgON->SetBinContent(l, k, j, a); + + Double_t b = fHgb->GetBinContent(k,j); + fHgb->SetBinContent (k, j, 0.0); + fHgMC->SetBinContent(l, k, j, b); + } + + + // fill target distribution fHSigmaTheta + // (only for plotting, it will not be used in the padding) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = histap->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, a); + } + + delete fHga; + delete fHgb; + + delete histap; + } + //............ end of loop over Theta bins .................... + + + // Define the target distributions + // fHDiffPixTheta, fHSigmaPixTheta + // (they are calculated as + // averages of the ON and MCdistributions) + // fHBlindPixNTheta, fHBlindPixIdTheta + // (they are calculated as + // the sum of the ON and MC distributions) + // (fHDiffPixTheta will be used in the padding of MC events) + + //............ start of new loop over Theta bins .................... + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/MC events to be padded : fracON, fracMC + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + Double_t normON; + Double_t normMC; + + // set ranges for 2D-projections of 3D-histograms + TAxis *ax = diffpixthon->GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthmc->GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histMC; + + TH1D *hist1; + TH1D *hist1MC; + + // weights for ON and MC distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracMC); + Double_t facMC = 0.5 * (1. + fracON - fracMC); + + *fLog << all << "Theta bin j : fracON, fracMC, facON, facMC = " + << j << ", " << fracON << ", " << fracMC << ", " + << facON << ", " << facMC << endl; + + //------------------------------------------------------------------ + // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' + ay = diffpixthon->GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + TAxis *az = diffpixthon->GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon->Project3D("zy"); + hist->SetName("dummy"); + histMC = (TH2D*)diffpixthmc->Project3D("zy"); + + normON = hist->Integral(); + normMC = histMC->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normMC != 0.0) histMC->Scale(1.0/normMC); + + // weighted average of ON and MC distributions + hist->Scale(facON); + hist->Add(histMC, facMC); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histMC; + + //------------------------------------------------------------------ + // define target distribution 'sigma vs. pixel number' + ay = sigmapixthon->GetYaxis(); + nbinspixel = ay->GetNbins(); + + az = sigmapixthon->GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigmapixthon->Project3D("zy"); + hist->SetName("dummy"); + histMC = (TH2D*)sigmapixthmc->Project3D("zy"); + + normON = hist->Integral(); + normMC = histMC->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normMC != 0.0) histMC->Scale(1.0/normMC); + + // weighted average of ON and MC distributions + hist->Scale(facON); + hist->Add(histMC, facMC); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histMC; + + + //------------------------------------------------------------------ + // define target distribution 'number of blind pixels per event' + ay = blindnthon->GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1MC = fHBlindPixNThetaMC->ProjectionY("", j, j, ""); + + normON = hist1->Integral(); + normMC = hist1MC->Integral(); + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + + // sum of ON and MC distributions + hist1->Add(hist1MC, 1.0); + hist1->Scale( 1.0/hist1->Integral() ); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1MC; + + //------------------------------------------------------------------ + // define target distribution 'id of blind pixel' + ay = blindidthon->GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1MC = fHBlindPixIdThetaMC->ProjectionY("", j, j, ""); + + // divide by the number of events (from fHBlindPixNTheta) + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + + // sum of ON and MC distributions + hist1->Add(hist1MC, 1.0); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1MC; + } + //............ end of new loop over Theta bins .................... + + *fLog << all + << "MPad::MergeONMC(); The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaThetaMC->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + + fHDiffPixThetaMC->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + + fHSigmaPixThetaMC->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + + fHBlindPixIdThetaMC->SetDirectory(NULL); + fHBlindPixIdThetaON->SetDirectory(NULL); + + fHBlindPixNThetaMC->SetDirectory(NULL); + fHBlindPixNThetaON->SetDirectory(NULL); + + fHgMC->SetDirectory(NULL); + fHgON->SetDirectory(NULL); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Read padding distributions from a file +// +// +Bool_t MPad::ReadPaddingDist(const char* namefilein) +{ + *fLog << all << "MPad : Read padding histograms from file " + << namefilein << endl; + + fInfile = new TFile(namefilein); + fInfile->ls(); + + //------------------------------------ + fHSigmaThetaMC = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC"); + if (!fHSigmaThetaMC) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl; + + fHSigmaThetaON = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); + if (!fHSigmaThetaON) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl; + + fHSigmaThetaOFF = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); + if (!fHSigmaThetaOFF) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl; + + //------------------------------------ + fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiff' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixDiff' was read in" << endl; + + fHDiffPixThetaMC = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC"); + if (!fHDiffPixThetaMC) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl; + + fHDiffPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); + if (!fHDiffPixThetaON) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl; + + fHDiffPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); + if (!fHDiffPixThetaOFF) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl; + + + //------------------------------------ + fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + *fLog << all + << "MPad : Object '3D-ThetaPixSigma' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixSigma' was read in" << endl; + + fHSigmaPixThetaMC = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC"); + if (!fHSigmaPixThetaMC) + { + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl; + + fHSigmaPixThetaON = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); + if (!fHSigmaPixThetaON) + { + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl; + + fHSigmaPixThetaOFF = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); + if (!fHSigmaPixThetaOFF) + { + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl; + + //------------------------------------ + fHgMC = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC"); + if (!fHgMC) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl; + + fHgON = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); + if (!fHgON) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-PaddingMatrixON' was read in" << endl; + + fHgOFF = + (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); + if (!fHgOFF) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl; + + //------------------------------------ + fHBlindPixIdThetaMC = + (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC"); + if (!fHBlindPixIdThetaMC) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindIdMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + + fHBlindPixIdThetaON = + (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON"); + if (!fHBlindPixIdThetaON) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindIdON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl; + + fHBlindPixIdThetaOFF = + (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF"); + if (!fHBlindPixIdThetaOFF) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindIdOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + + //------------------------------------ + fHBlindPixNThetaMC = + (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC"); + if (!fHBlindPixNThetaMC) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindNMC' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl; + + fHBlindPixNThetaON = + (TH2D*) gROOT->FindObject("2D-ThetaBlindNON"); + if (!fHBlindPixNThetaON) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindNON' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindNON' was read in" << endl; + + fHBlindPixNThetaOFF = + (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF"); + if (!fHBlindPixNThetaOFF) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindNOFF' not found on root file" + << endl; + return kFALSE; + } + *fLog << all + << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl; + + //------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Write padding distributions onto a file +// +// +Bool_t MPad::WritePaddingDist(const char* namefileout) +{ + *fLog << all << "MPad : Write padding histograms onto file " + << namefileout << endl; + + TFile outfile(namefileout, "RECREATE"); + + fHBlindPixNThetaMC->Write(); + fHBlindPixNThetaON->Write(); + fHBlindPixNThetaOFF->Write(); + + fHBlindPixIdThetaMC->Write(); + fHBlindPixIdThetaON->Write(); + fHBlindPixIdThetaOFF->Write(); + + fHSigmaThetaMC->Write(); + fHSigmaThetaON->Write(); + fHSigmaThetaOFF->Write(); + + fHDiffPixTheta->Write(); + fHDiffPixThetaMC->Write(); + fHDiffPixThetaON->Write(); + fHDiffPixThetaOFF->Write(); + + fHSigmaPixTheta->Write(); + fHSigmaPixThetaMC->Write(); + fHSigmaPixThetaON->Write(); + fHSigmaPixThetaOFF->Write(); + + fHgMC->Write(); + fHgON->Write(); + fHgOFF->Write(); + + *fLog << all + << "MPad::WriteTargetDist(); padding histograms written onto file " + << namefileout << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set type of data to be padded +// +// this is not necessary if the type of the data can be recognized +// directly from the input +// +// +void MPad::SetDataType(const char* type) +{ + fType = type; + *fLog << all << "MPad::SetDataType(); type of data to be padded : " + << fType << endl; + + return; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers and prepare the histograms +// +Int_t MPad::PreProcess(MParList *pList) +{ + if ( !fHSigmaThetaMC || !fHSigmaThetaON || !fHSigmaThetaOFF || + !fHDiffPixThetaMC || !fHDiffPixThetaON || !fHDiffPixThetaOFF || + !fHBlindPixIdThetaMC|| !fHBlindPixIdThetaON|| !fHBlindPixIdThetaOFF || + !fHBlindPixNThetaMC || !fHBlindPixNThetaON || !fHBlindPixNThetaOFF || + !fHgMC || !fHgON || !fHgOFF ) + { + *fLog << err + << "MPad : At least one of the histograms needed for the padding is not defined ... aborting." + << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MPad : MMcEvt not found... aborting." + << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPad : MPedestalCam not found... aborting." + << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err + << "MPad : MGeomCam not found (no geometry information available)... aborting." + << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << err << "MPad : MCerPhotEvt not found... aborting." + << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << err << "MPad : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlinds) + { + *fLog << err << "MPad : MBlindPixels not found... aborting." + << endl; + return kFALSE; + } + + if (fType !="ON" && fType !="OFF" && fType !="MC") + { + *fLog << err + << "MPad : Type of data to be padded not defined... aborting." + << endl; + return kFALSE; + } + + + //-------------------------------------------------------------------- + // histograms for checking the padding + // + // plot pedestal sigmas + fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", + 100, 0.0, 5.0, 100, 0.0, 5.0); + fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); + fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); + + // plot no.of photons (before vs. after padding) + fHPhotons = new TH2D("Photons","Photons: after vs.before padding", + 100, -10.0, 90.0, 100, -10, 90); + fHPhotons->SetXTitle("no.of photons before padding"); + fHPhotons->SetYTitle("no.of photons after padding"); + + // plot of added NSB + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB->SetXTitle("no.of added NSB photons"); + fHNSB->SetYTitle("no.of pixels"); + + + //-------------------------------------------------------------------- + + fIter = 20; + + memset(fErrors, 0, sizeof(fErrors)); + memset(fWarnings, 0, sizeof(fWarnings)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the Padding (noise adjustment) +// +// input for the padding : +// - the matrices fHgON, fHgOFF, fHgMC +// - the original distributions fHBlindPixNTheta, fHBlindPixIdTheta +// fHSigmaTheta, fHDiffPixTheta +// + +Int_t MPad::Process() +{ + *fLog << all << "Entry MPad::Process();" << endl; + + // this is the conversion factor from the number of photons + // to the number of photo electrons + // later to be taken from MCalibration + // fPEperPhoton = PW * LG * QE * 1D + Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + + //------------------------------------------------ + // add pixels to MCerPhotEvt which are not yet in; + // set their number of photons equal to zero + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetNumPixels(); + for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); + } + + + + //----------------------------------------- + Int_t rc=0; + Int_t rw=0; + + const UInt_t npix = fEvt->GetNumPixels(); + + //------------------------------------------- + // Calculate average sigma of the event + // + Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "MPad::Process(); before padding : " << endl; + fSigmabar->Print(""); + Double_t sigbarold = sigbarold_phot * fPEperPhoton; + Double_t sigbarold2 = sigbarold*sigbarold; + + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + *fLog << all << "theta = " << theta << endl; + + Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) + { + *fLog << warn + << "MPad::Process(); binNumber out of range : theta, binTheta = " + << theta << ", " << binTheta << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + //------------------------------------------- + // get number of events in this theta bin + // and number of events in this sigbarold_phot bin + + Double_t nTheta; + Double_t nSigma; + + TH2D *st = NULL; + if (fType == "MC") st = fHSigmaThetaMC; + else if (fType == "ON") st = fHSigmaThetaON; + else if (fType == "OFF") st = fHSigmaThetaOFF; + else + { + *fLog << err << "MPad : illegal data type '" << fType + << "'" << endl; + return kERROR; + } + + TH1D *hn; + hn = st->ProjectionY("", binTheta, binTheta, ""); + nTheta = hn->Integral(); + Int_t binSigma = hn->FindBin(sigbarold_phot); + nSigma = hn->GetBinContent(binSigma); + + *fLog << all + << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + << theta << ", " << sigbarold_phot << ", " << binTheta << ", " + << binSigma << ", " << nTheta << ", " << nSigma << endl; + + delete hn; + + + //------------------------------------------- + // for the current theta : + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // + // ON : add the blind pixels from the OFF data + // OFF : add the blind pixels from the ON data + // MC : add the blind pixels from the ON and OFF data + + //----------------------------------- + if (fType == "ON" || fType == "MC") + { + // numBlind is the number of blind pixels in an event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNThetaOFF->ProjectionY("", binTheta, binTheta, ""); + if ( nblind->Integral() == 0.0 ) + { + *fLog << warn << "MPad::Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + } + delete nblind; + + + // throw the Id of numBlind different pixels in this event + if ( numBlind > 0) + { + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdThetaOFF->ProjectionY("", binTheta, binTheta, ""); + if ( hblind->Integral() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + *fLog << all << "idBlind = " << idBlind << endl; + } + fBlinds->SetReadyToSave(); + + delete hblind; + } + } + + //------------------------------------ + if (fType == "OFF" || fType == "MC") + { + // throw numBlind; + // numBlind is the number of blind pixels in an event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNThetaON->ProjectionY("", binTheta, binTheta, ""); + if ( nblind->Integral() == 0.0 ) + { + *fLog << warn << "MPad::Process(); projection for Theta bin " + << binTheta << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + } + delete nblind; + + + // throw the Id of numBlind different pixels in this event + if ( numBlind > 0) + { + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdThetaON->ProjectionY("", binTheta, binTheta, ""); + if ( hblind->Integral() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + *fLog << all << "idBlind = " << idBlind << endl; + } + fBlinds->SetReadyToSave(); + + delete hblind; + } + } + + //****************************************************************** + // has the event to be padded ? + // if yes : to which sigmabar should it be padded ? + // + + Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); + *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " + << sigbarold_phot << endl; + + Double_t prob; + TH1D *hpad = NULL; + + TH3D *Hg = NULL; + if (fType == "MC") Hg = fHgMC; + else if (fType == "ON") Hg = fHgON; + else if (fType == "OFF") Hg = fHgOFF; + else + { + *fLog << err << "MPad : illegal type of data '" << fType << "'" + << endl; + return kERROR; + } + + hpad = Hg->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); + + //Int_t nb = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nb; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHgON = " << i << ", " + // << hpad->GetBinContent(i) << endl; + //} + + if (nSigma != 0.0) + prob = hpad->Integral() * nTheta/nSigma; + else + prob = 0.0; + + *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + << ", " << prob << endl; + + if ( prob <= 0.0 || gRandom->Uniform() > prob ) + { + delete hpad; + // event should not be padded + *fLog << all << "event is not padded" << endl; + + rc = 8; + fErrors[rc]++; + return kTRUE; + } + // event should be padded + *fLog << all << "event will be padded" << endl; + + + //------------------------------------------- + // for the current theta, generate a sigmabar : + // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF + // + Double_t sigmabar_phot = 0; + Double_t sigmabar = 0; + + + //Int_t nbinsx = hpad->GetNbinsX(); + //for (Int_t i=1; i<=nbinsx; i++) + //{ + // if (hpad->GetBinContent(i) != 0.0) + // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i) + // << endl; + //} + + sigmabar_phot = hpad->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + + *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; + + delete hpad; + + const Double_t sigmabar2 = sigmabar*sigmabar; + + //------------------------------------------- + + *fLog << all << "MPad::Process(); sigbarold, sigmabar = " + << sigbarold << ", "<< sigmabar << endl; + + // Skip event if target sigmabar is <= sigbarold + if (sigmabar <= sigbarold) + { + *fLog << all << "MPad::Process(); target sigmabar is less than sigbarold : " + << sigmabar << ", " << sigbarold << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; + return kCONTINUE; + } + + + //------------------------------------------- + // + // Calculate average number of NSB photo electrons to be added (lambdabar) + // from the value of sigmabar, + // - making assumptions about the average electronic noise (elNoise2) and + // - using a fixed value (F2excess) for the excess noise factor + + Double_t elNoise2; // [photo electrons] + Double_t F2excess = 1.3; + Double_t lambdabar; // [photo electrons] + + + Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta); + if (binTheta != bincheck) + { + *fLog << err + << "MPad::Process(); The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. + // The average electronic noise (to be added) has to be well above this RMS, + // otherwise the electronic noise of an individual pixel (elNoise2Pix) + // may become negative + + TH1D *hnoise; + if (fType == "MC") + hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "ON") + hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else if (fType == "OFF") + hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + else + { + *fLog << err << "MPad::Process(); illegal data type... aborting" + << endl; + return kERROR; + } + + Double_t RMS_phot = hnoise->GetRMS(1); + Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; + delete hnoise; + + elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); + *fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; + + // This value of lambdabar is the same for all pixels; + // note that lambdabar is normalized to the area of pixel 0 + + //---------- start loop over pixels --------------------------------- + // do the padding for each pixel + // + // pad only pixels - which are used (before image cleaning) + // + + Double_t sigma2 = 0; + + Double_t diff_phot = 0; + Double_t diff = 0; + + Double_t addSig2_phot= 0; + Double_t addSig2 = 0; + + Double_t elNoise2Pix = 0; + + + for (UInt_t i=0; iGetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; + Double_t oldsigma_phot = ppix.GetPedestalRms(); + Double_t oldsigma = oldsigma_phot * fPEperPhoton; + Double_t oldsigma2 = oldsigma*oldsigma; + + //--------------------------------- + // throw the Sigma for this pixel + // + Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); + + Int_t count; + Bool_t ok; + + TH1D *hdiff; + + // throw the Sigma for this pixel from the distribution fHDiffPixTheta + // MC : from fHDiffPixTheta + // ON : from fHDiffPixThetaOFF + // OFF : from fHDiffPixThetaON + + TH3D *sp = NULL; + + if (fType == "MC") sp = fHDiffPixTheta; + else if (fType == "ON") sp = fHDiffPixThetaOFF; + else if (fType == "OFF") sp = fHDiffPixThetaON; + else + { + *fLog << err << "MPad::Process(); illegal data type... aborting" + << endl; + return kERROR; + } + + hdiff = sp->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + + if ( hdiff->GetEntries() == 0 ) + { + *fLog << warn << "MPad::Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hdiff; + + rc = 5; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + diff = diff_phot * fPEperPhoton * fPEperPhoton; + + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2/ratioArea + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, diff = " << theta + << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; + + rw = 1; + fWarnings[rw]++; + } + else + { + rw = 0; + fWarnings[rw]++; + } + + delete hdiff; + sigma2 = diff + sigmabar2; + + + //--------------------------------- + // get the additional sigma^2 for this pixel (due to the padding) + + addSig2 = sigma2*ratioArea - oldsigma2; + addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; + + + //--------------------------------- + // throw actual number of additional NSB photo electrons (NSB) + // and its RMS (sigmaNSB) + + Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); + Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; + Double_t sigmaNSB0; + + if (arg >= 0.0) + { + sigmaNSB0 = sqrt( arg ); + } + else + { + sigmaNSB0 = 0.0000001; + if (arg < -1.e-10) + { + *fLog << warn << "MPad::Process(); argument of sqrt < 0 : " + << arg << endl; + } + } + + + //--------------------------------- + // smear NSB0 according to sigmaNSB0 + // and subtract lambdabar because of AC coupling + + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; + Double_t NSB_phot = NSB / fPEperPhoton; + + //--------------------------------- + + // add additional NSB to the number of photons + Double_t oldphotons_phot = pix.GetNumPhotons(); + Double_t oldphotons = oldphotons_phot * fPEperPhoton; + Double_t newphotons = oldphotons + NSB; + Double_t newphotons_phot = newphotons / fPEperPhoton; + pix.SetNumPhotons( newphotons_phot ); + + + fHNSB->Fill( NSB_phot/ratioArea ); + fHPhotons->Fill( oldphotons_phot/ratioArea, + newphotons_phot/ratioArea ); + + + // error: add sigma of padded noise quadratically + Double_t olderror_phot = pix.GetErrorPhot(); + Double_t newerror_phot = + sqrt( olderror_phot*olderror_phot + addSig2_phot ); + pix.SetErrorPhot( newerror_phot ); + + + Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + Double_t newsigma_phot = newsigma / fPEperPhoton; + ppix.SetPedestalRms( newsigma_phot ); + + fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); + } + //---------- end of loop over pixels --------------------------------- + + // Calculate sigmabar again and crosscheck + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "MPad::Process(); after padding : " << endl; + fSigmabar->Print(""); + + *fLog << all << "Exit MPad::Process();" << endl; + + rc = 0; + fErrors[rc]++; + return kTRUE; + //****************************************************************** +} + +// -------------------------------------------------------------------------- +// +// +Int_t MPad::PostProcess() +{ + if (GetNumExecutions() != 0) + { + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + + int temp; + if (fWarnings[0] != 0.0) + temp = (int)(fWarnings[1]*100/fWarnings[0]); + else + temp = 100; + + *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) + << temp + << "%) Warning : iteration to find acceptable sigma failed" + << ", fIter = " << fIter << endl; + + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) + << (int)(fErrors[2]*100/GetNumExecutions()) + << "%) Evts skipped due to: Zenith angle out of range" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) + << (int)(fErrors[3]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigmabar" << endl; + + *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) + << (int)(fErrors[4]*100/GetNumExecutions()) + << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; + + *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) + << (int)(fErrors[5]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; + + *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) + << (int)(fErrors[6]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma" << endl; + + *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) + << (int)(fErrors[7]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Blind pixels" << endl; + + *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) + << (int)(fErrors[8]*100/GetNumExecutions()) + << "%) Evts didn't have to be padded" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts were successfully padded!" << endl; + *fLog << endl; + + } + + //--------------------------------------------------------------- + TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1500)); + c.Divide(3, 5); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + fHNSB->SetDirectory(NULL); + fHNSB->DrawCopy(); + fHNSB->SetBit(kCanDelete); + + c.cd(2); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->DrawCopy(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(3); + fHPhotons->SetDirectory(NULL); + fHPhotons->DrawCopy(); + fHPhotons->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(4); + fHSigmaTheta->SetDirectory(NULL); + fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta"); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 + + c.cd(5); + TH2D *l1 = NULL; + l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); + l1->SetDirectory(NULL); + l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l1->SetXTitle("\\Theta [\\circ]"); + l1->SetYTitle("Sigma^2-Sigmabar^2"); + + l1->DrawCopy("box"); + l1->SetBit(kCanDelete);; + + c.cd(8); + TH2D *l2 = NULL; + l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); + l2->SetDirectory(NULL); + l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l2->SetXTitle("pixel"); + l2->SetYTitle("Sigma^2-Sigmabar^2"); + + l2->DrawCopy("box"); + l2->SetBit(kCanDelete);; + + //-------------------------------------------------------------------- + // draw the 3D histogram (target): Theta, pixel, Sigma + + c.cd(6); + TH2D *k1 = NULL; + k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); + k1->SetDirectory(NULL); + k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)"); + k1->SetXTitle("\\Theta [\\circ]"); + k1->SetYTitle("Sigma"); + + k1->DrawCopy("box"); + k1->SetBit(kCanDelete); + + c.cd(9); + TH2D *k2 = NULL; + k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); + k2->SetDirectory(NULL); + k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)"); + k2->SetXTitle("pixel"); + k2->SetYTitle("Sigma"); + + k2->DrawCopy("box"); + k2->SetBit(kCanDelete);; + + + //-------------------------------------------------------------------- + + c.cd(13); + TH2D *m1; + m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); + m1->SetDirectory(NULL); + m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetXTitle("Sigmabar-old"); + m1->SetYTitle("Sigmabar-new"); + + m1->DrawCopy("box"); + m1->SetBit(kCanDelete);; + + c.cd(14); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetXTitle("Sigmabar-old"); + m2->SetYTitle("Sigmabar-new"); + + m2->DrawCopy("box"); + m2->SetBit(kCanDelete);; + + c.cd(15); + TH2D *m3; + m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); + m3->SetDirectory(NULL); + m3->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); + m3->SetXTitle("Sigmabar-old"); + m3->SetYTitle("Sigmabar-new"); + + m3->DrawCopy("box"); + m3->SetBit(kCanDelete);; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.h (revision 2777) +++ trunk/MagicSoft/Mars/manalysis/MPad.h (revision 2777) @@ -0,0 +1,131 @@ +#ifndef MARS_MPad +#define MARS_MPad + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class TH3D; + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MSigmabar; +class MParList; +class MBlindPixels; +class MRead; +class MFilterList; + + +class MPad : public MTask +{ +private: + MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + MBlindPixels *fBlinds; + + TString fType; // type of data to be padded + TFile *fInfile; // input file containing padding histograms + + Int_t fIter; + + Int_t fErrors[9]; + Int_t fWarnings[2]; + + //---------------------------------- + // plots used for the padding + // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" + TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixIdThetaMC; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixIdThetaON; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixIdThetaOFF; // 2D-histogram (blind pixel Id vs. Theta) + + //--------------------- + TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) + TH2D *fHBlindPixNThetaMC; // 2D-histogram (no.of blind pixels vs. Theta) + TH2D *fHBlindPixNThetaON; // 2D-histogram (no.of blind pixels vs. Theta) + TH2D *fHBlindPixNThetaOFF; // 2D-histogram (no.of blind pixels vs. Theta) + + //--------------------- + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaMC; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) + TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) + + //--------------------- + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaMC; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + + //--------------------- + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaMC; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) + + //--------------------- + TH3D *fHgMC; // matrix (Theta, sigbarold, sigbarnew) for MC data + TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data + TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data + + //------------------------------- + // plots for checking the padding + TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after + // versus before padding + TH2D *fHPhotons; // 2D-histogram : no.of photons after + // versus before padding + TH1D *fHNSB; // 1D-histogram : additional NSB + + + Bool_t Merge2Distributions(TH1D *hista, TH1D * histb, TH1D * histap, + TH2D *fHga, TH2D * fHgb, Int_t nbinssig); + + +public: + MPad(const char *name=NULL, const char *title=NULL); + ~MPad(); + + Bool_t MergeONOFFMC( + TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, + TH2D *blindidthmc, TH2D *blindnthmc, + TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, + TH2D *blindidthon, TH2D *blindnthon, + TH2D *sigthoff=NULL, TH3D *diffpixthoff=NULL,TH3D *sigmapixthoff=NULL, + TH2D *blindidthoff=NULL, TH2D *blindnthoff=NULL); + + Bool_t MergeONMC( + TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, + TH2D *blindidthmc, TH2D *blindnthmc, + TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, + TH2D *blindidthon, TH2D *blindnthon); + + Bool_t ReadPaddingDist(const char *filein); + Bool_t WritePaddingDist(const char *fileout); + + void SetDataType(const char *type); // type of data to be padded + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + ClassDef(MPad, 0) // task for the padding +}; + +#endif + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadON.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadON.cc (revision 2776) +++ (revision ) @@ -1,940 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Wolfgang Wittek, 12/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MPadON -// -// This task applies padding such that for a given pixel and for a given -// Theta bin the resulting distribution of the pedestal sigma is identical -// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta. -// -// The number of photons, its error and the pedestal sigmas are altered. -// On average, the number of photons added is zero. -// -// The formulas used can be found in Thomas Schweizer's Thesis, -// Section 2.2.1 -// -// There are 2 options for the padding : -// -// 1) fPadFlag = 1 : -// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta -// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using -// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 -// (fHDiffPixTheta). -// -// This is the preferred option as it takes into account the -// correlations between the Sigma of a pixel and Sigmabar. -// -// THE FOLLOWING OPTION HAS NOT BEEN TESTED : -// 2) fPadFlag = 2 : -// Generate a pedestal sigma for each pixel using the 3D-histogram -// Theta, pixel no., Sigma (fHSigmaPixTheta). -// -// -// The padding has to be done before the image cleaning because the -// image cleaning depends on the pedestal sigmas. -// -// For random numbers gRandom is used. -// -///////////////////////////////////////////////////////////////////////////// -#include "MPadON.h" - -#include - -#include -#include -#include -#include -#include - -#include "MBinning.h" -#include "MSigmabar.h" -#include "MMcEvt.hxx" -#include "MLog.h" -#include "MLogManip.h" -#include "MParList.h" -#include "MGeomCam.h" - -#include "MCerPhotPix.h" -#include "MCerPhotEvt.h" - -#include "MPedestalCam.h" -#include "MPedestalPix.h" -#include "MBlindPixels.h" - -ClassImp(MPadON); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPadON::MPadON(const char *name, const char *title) -{ - fName = name ? name : "MPadON"; - fTitle = title ? title : "Task for the padding (Schweizer)"; - - fHSigmaTheta = NULL; - fHSigmaPixTheta = NULL; - fHDiffPixTheta = NULL; - fHBlindPixIdTheta = NULL; - fHBlindPixNTheta = NULL; - - fHSigmaPedestal = NULL; - fHPhotons = NULL; - fHNSB = NULL; -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MPadON::~MPadON() -{ - if (fHSigmaPedestal != NULL) delete fHSigmaPedestal; - if (fHPhotons != NULL) delete fHPhotons; - if (fHNSB != NULL) delete fHNSB; -} - -// -------------------------------------------------------------------------- -// -// Set the references to the histograms to be used in the padding -// -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) -// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) -// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) -// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) -// -// -void MPadON::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, - TH2D *hist2Pix, TH2D *hist2PixN) -{ - fHSigmaTheta = hist2; - fHSigmaPixTheta = hist3; - fHDiffPixTheta = hist3Diff; - fHBlindPixIdTheta = hist2Pix; - fHBlindPixNTheta = hist2PixN; - - fHSigmaTheta->SetDirectory(NULL); - fHSigmaPixTheta->SetDirectory(NULL); - fHDiffPixTheta->SetDirectory(NULL); - fHBlindPixIdTheta->SetDirectory(NULL); - fHBlindPixNTheta->SetDirectory(NULL); - - Print(); -} - -// -------------------------------------------------------------------------- -// -// Set the option for the padding -// -// There are 2 options for the padding : -// -// 1) fPadFlag = 1 : -// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta -// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using -// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 -// (fHDiffPixTheta). -// -// This is the preferred option as it takes into account the -// correlations between the Sigma of a pixel and Sigmabar. -// -// NOT YET TESTED : -// 2) fPadFlag = 2 : -// Generate a pedestal sigma for each pixel using the 3D-histogram -// Theta, pixel no., Sigma (fHSigmaPixTheta). -// -void MPadON::SetPadFlag(Int_t padflag) -{ - fPadFlag = padflag; - *fLog << inf << "MPadON::SetPadFlag(); choose option " << fPadFlag << endl; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers and prepare the histograms -// -Int_t MPadON::PreProcess(MParList *pList) -{ - if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta || - !fHBlindPixIdTheta || !fHBlindPixNTheta) - { - *fLog << err - << "MPadON : At least one of the histograms needed for the padding is not defined ... aborting." - << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MPadON : MMcEvt not found... aborting." - << endl; - return kFALSE; - } - - fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); - if (!fPed) - { - *fLog << err << "MPadON : MPedestalCam not found... aborting." << endl; - return kFALSE; - } - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << err - << "MPadON : MGeomCam not found (no geometry information available)... aborting." - << endl; - return kFALSE; - } - - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << err << "MPadON : MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); - if (!fSigmabar) - { - *fLog << err << "MPadON : MSigmabar not found... aborting." << endl; - return kFALSE; - } - - fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); - if (!fBlinds) - { - *fLog << err << "MPadON : MBlindPixels not found... aborting." << endl; - return kFALSE; - } - - - //-------------------------------------------------------------------- - // histograms for checking the padding - // - // plot pedestal sigmas - fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", - 100, 0.0, 5.0, 100, 0.0, 5.0); - fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); - fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); - - // plot no.of photons (before vs. after padding) - fHPhotons = new TH2D("Photons","Photons: after vs.before padding", - 100, -10.0, 90.0, 100, -10, 90); - fHPhotons->SetXTitle("no.of photons before padding"); - fHPhotons->SetYTitle("no.of photons after padding"); - - // plot of added NSB - fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); - fHNSB->SetXTitle("no.of added NSB photons"); - fHNSB->SetYTitle("no.of pixels"); - - - //-------------------------------------------------------------------- - - fIter = 20; - - memset(fErrors, 0, sizeof(fErrors)); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Do the Padding -// idealy the events to be padded should have been generated without NSB -// -Int_t MPadON::Process() -{ - *fLog << all << "Entry MPadON::Process();" << endl; - - // this is the conversion factor from the number of photons - // to the number of photo electrons - // later to be taken from MCalibration - // fPEperPhoton = PW * LG * QE * 1D - Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; - - //------------------------------------------------ - // add pixels to MCerPhotEvt which are not yet in; - // set their numnber of photons equal to zero - - UInt_t imaxnumpix = fCam->GetNumPixels(); - - for (UInt_t i=0; iGetNumPixels(); - for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); - } - - - - //----------------------------------------- - Int_t rc=0; - - const UInt_t npix = fEvt->GetNumPixels(); - - //fSigmabar->Calc(*fCam, *fPed, *fEvt); - //*fLog << all << "before padding : " << endl; - //fSigmabar->Print(""); - - - //$$$$$$$$$$$$$$$$$$$$$$$$$$ - // to simulate the situation that before the padding the NSB and - // electronic noise are zero : set Sigma = 0 for all pixels - //for (UInt_t i=0; ioperator[](i); - // Int_t j = pix.GetPixId(); - - // MPedestalPix &ppix = fPed->operator[](j); - // ppix.SetMeanRms(0.0); - //} - //$$$$$$$$$$$$$$$$$$$$$$$$$$ - - //------------------------------------------- - // Calculate average sigma of the event - // - Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "before padding : " << endl; - fSigmabar->Print(""); - Double_t sigbarold = sigbarold_phot * fPEperPhoton; - - Double_t sigbarold2 = sigbarold*sigbarold; - - - // for MC data : expect sigmabar to be zero before the padding - if (sigbarold_phot > 0) - { - *fLog << err - << "MPadON::Process(); sigmabar of event to be padded is > 0 : " - << sigbarold_phot << ". Stop event loop " << endl; - // input data should have sigmabar = 0; stop event loop - - rc = 1; - fErrors[rc]++; - return kERROR; - } - - const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - *fLog << all << "theta = " << theta << endl; - - - //------------------------------------------- - // for the current theta, - // generate blind pixels according to the histograms - // fHBlindPixNTheta and fHBlindPixIDTheta - // - - - Int_t binPix = fHBlindPixNTheta->GetXaxis()->FindBin(theta); - - if ( binPix < 1 || binPix > fHBlindPixNTheta->GetNbinsX() ) - { - *fLog << warn - << "MPadON::Process(); binNumber out of range : theta, binPix = " - << theta << ", " << binPix << "; Skip event " << endl; - // event cannot be padded; skip event - - rc = 2; - fErrors[rc]++; - return kCONTINUE; - } - - // numBlind is the number of blind pixels in this event - TH1D *nblind; - UInt_t numBlind; - - nblind = fHBlindPixNTheta->ProjectionY("", binPix, binPix, ""); - if ( nblind->GetEntries() == 0.0 ) - { - *fLog << warn << "MPadON::Process(); projection for Theta bin " - << binPix << " has no entries; Skip event " << endl; - // event cannot be padded; skip event - delete nblind; - - rc = 7; - fErrors[rc]++; - return kCONTINUE; - } - else - { - numBlind = (Int_t) (nblind->GetRandom()+0.5); - } - delete nblind; - - - // throw the Id of numBlind different pixels in this event - TH1D *hblind; - UInt_t idBlind; - UInt_t listId[npix]; - UInt_t nlist = 0; - Bool_t equal; - - hblind = fHBlindPixIdTheta->ProjectionY("", binPix, binPix, ""); - if ( hblind->GetEntries() > 0.0 ) - for (UInt_t i=0; iGetRandom()+0.5); - equal = kFALSE; - for (UInt_t j=0; jSetPixelBlind(idBlind); - *fLog << "idBlind = " << idBlind << endl; - } - fBlinds->SetReadyToSave(); - - delete hblind; - - - //------------------------------------------- - // for the current theta, - // generate a sigmabar according to the histogram fHSigmaTheta - // - Double_t sigmabar_phot = 0; - Double_t sigmabar = 0; - Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta); - - if (binPix != binNumber) - { - *fLog << err - << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting" - << endl; - return kERROR; - } - - TH1D *hsigma; - - hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); - if ( hsigma->GetEntries() == 0.0 ) - { - *fLog << warn << "MPadON::Process(); projection for Theta bin " - << binNumber << " has no entries; Skip event " << endl; - // event cannot be padded; skip event - delete hsigma; - - rc = 3; - fErrors[rc]++; - return kCONTINUE; - } - else - { - sigmabar_phot = hsigma->GetRandom(); - sigmabar = sigmabar_phot * fPEperPhoton; - *fLog << all << "Theta, bin number = " << theta << ", " << binNumber - << ", sigmabar_phot = " << sigmabar_phot << endl; - } - delete hsigma; - - const Double_t sigmabar2 = sigmabar*sigmabar; - - //------------------------------------------- - - *fLog << all << "MPadON::Process(); sigbarold, sigmabar = " - << sigbarold << ", "<< sigmabar << endl; - - // Skip event if target sigmabar is <= sigbarold - if (sigmabar <= sigbarold) - { - *fLog << all - << "MPadON::Process(); target sigmabar is less than sigbarold : " - << sigmabar << ", " << sigbarold << ", Skip event" << endl; - - rc = 4; - fErrors[rc]++; - return kCONTINUE; - } - - - //------------------------------------------- - // - // Calculate average number of NSB photons to be added (lambdabar) - // from the value of sigmabar, - // - making assumptions about the average electronic noise (elNoise2) and - // - using a fixed value (F2excess) for the excess noise factor - - Double_t elNoise2; // [photo electrons] - Double_t F2excess = 1.3; - Double_t lambdabar; // [photo electrons] - - - - Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta); - if (binTheta != binNumber) - { - *fLog << err - << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting" - << endl; - return kERROR; - } - - // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. - // The average electronic noise (to be added) has to be well above this RMS, - // otherwise the electronic noise of an individual pixel (elNoise2Pix) - // may become negative - - TH1D *hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, - 0, 9999, ""); - Double_t RMS_phot = hnoise->GetRMS(1); - Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; - delete hnoise; - - elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); - *fLog << all << "elNoise2 = " << elNoise2 << endl; - - lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; - - // This value of lambdabar is the same for all pixels; - // note that lambdabar is normalized to the area of pixel 0 - - //---------- start loop over pixels --------------------------------- - // do the padding for each pixel - // - // pad only pixels - which are used (before image cleaning) - // - Double_t sig_phot = 0; - Double_t sig = 0; - - Double_t sigma2 = 0; - - Double_t diff_phot = 0; - Double_t diff = 0; - - Double_t addSig2_phot= 0; - Double_t addSig2 = 0; - - Double_t elNoise2Pix = 0; - - - for (UInt_t i=0; iGetPixRatio(j); - - MPedestalPix &ppix = (*fPed)[j]; - Double_t oldsigma_phot = ppix.GetPedestalRms(); - Double_t oldsigma = oldsigma_phot * fPEperPhoton; - Double_t oldsigma2 = oldsigma*oldsigma; - - //--------------------------------- - // throw the Sigma for this pixel - // - Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); - - Int_t count; - Bool_t ok; - - TH1D *hdiff; - TH1D *hsig; - - switch (fPadFlag) - { - case 1 : - // throw the Sigma for this pixel from the distribution fHDiffPixTheta - - hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - if ( hdiff->GetEntries() == 0 ) - { - *fLog << warn << "MPadON::Process(); projection for Theta bin " - << binTheta << " and pixel bin " << binPixel - << " has no entries; aborting " << endl; - delete hdiff; - - rc = 5; - fErrors[rc]++; - return kCONTINUE; - } - - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - diff = diff_phot * fPEperPhoton * fPEperPhoton; - - // the following condition ensures that elNoise2Pix > 0.0 - if ( (diff + sigmabar2 - oldsigma2/ratioArea - - lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - if (!ok) - { - *fLog << all << "theta, j, count, sigmabar, diff = " << theta << ", " - << j << ", " << count << ", " << sigmabar << ", " - << diff << endl; - diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; - } - delete hdiff; - sigma2 = diff + sigmabar2; - break; - - case 2 : - // throw the Sigma for this pixel from the distribution fHSigmaPixTheta - - hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - if ( hsig->GetEntries() == 0 ) - { - *fLog << err << "MPadON::Process(); projection for Theta bin " - << binTheta << " and pixel bin " << binPixel - << " has no entries; aborting " << endl; - delete hsig; - - rc = 6; - fErrors[rc]++; - return kCONTINUE; - } - - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - sig = sig_phot * fPEperPhoton; - - sigma2 = sig*sig/ratioArea; - // the following condition ensures that elNoise2Pix > 0.0 - - if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - if (!ok) - { - - *fLog << "theta, j, count, sigmabar, sig = " << theta << ", " - << j << ", " << count << ", " << sigmabar << ", " - << sig << endl; - sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; - } - delete hsig; - break; - } - - //--------------------------------- - // get the additional sigma^2 for this pixel (due to the padding) - - addSig2 = sigma2*ratioArea - oldsigma2; - addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); - - //--------------------------------- - // get the additional electronic noise for this pixel - - elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; - - - //--------------------------------- - // throw actual number of additional NSB photons (NSB) - // and its RMS (sigmaNSB) - - Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); - Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; - Double_t sigmaNSB0; - - if (arg >= 0) - { - sigmaNSB0 = sqrt( arg ); - } - else - { - sigmaNSB0 = 0.0000001; - if (arg < -1.e-10) - { - *fLog << warn << "MPadON::Process(); argument of sqrt < 0 : " - << arg << endl; - } - } - - - //--------------------------------- - // smear NSB0 according to sigmaNSB0 - // and subtract lambdabar because of AC coupling - - Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; - Double_t NSB_phot = NSB / fPEperPhoton; - - //--------------------------------- - - // add additional NSB to the number of photons - Double_t oldphotons_phot = pix.GetNumPhotons(); - Double_t oldphotons = oldphotons_phot * fPEperPhoton; - Double_t newphotons = oldphotons + NSB; - Double_t newphotons_phot = newphotons / fPEperPhoton; - pix.SetNumPhotons(newphotons_phot); - - - fHNSB->Fill( NSB_phot/sqrt(ratioArea) ); - fHPhotons->Fill( oldphotons_phot/ratioArea, - newphotons_phot/ratioArea ); - - - // error: add sigma of padded noise quadratically - Double_t olderror_phot = pix.GetErrorPhot(); - Double_t newerror_phot = - sqrt(olderror_phot*olderror_phot + addSig2_phot); - pix.SetErrorPhot(newerror_phot); - - - Double_t newsigma = sqrt(oldsigma2 + addSig2); - Double_t newsigma_phot = newsigma / fPEperPhoton; - ppix.SetPedestalRms(newsigma_phot); - - fHSigmaPedestal->Fill(oldsigma_phot, newsigma_phot); - } - //---------- end of loop over pixels --------------------------------- - - // Calculate sigmabar again and crosscheck - - fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPadON::Process(); after padding : " << endl; - fSigmabar->Print(""); - - - *fLog << all << "Exit MPadON::Process();" << endl; - - rc = 0; - fErrors[rc]++; - - return kTRUE; - -} - -// -------------------------------------------------------------------------- -// -// -Int_t MPadON::PostProcess() -{ - if (GetNumExecutions() != 0) - { - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) - << (int)(fErrors[1]*100/GetNumExecutions()) - << "%) Evts skipped due to: Sigmabar_old > 0" << endl; - - *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) - << (int)(fErrors[2]*100/GetNumExecutions()) - << "%) Evts skipped due to: Zenith angle out of range" << endl; - - *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) - << (int)(fErrors[3]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigmabar" << endl; - - *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) - << (int)(fErrors[4]*100/GetNumExecutions()) - << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; - - *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) - << (int)(fErrors[5]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; - - *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) - << (int)(fErrors[6]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma" << endl; - - *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) - << (int)(fErrors[7]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Blind pixels" << endl; - - *fLog << " " << fErrors[0] << " (" - << (int)(fErrors[0]*100/GetNumExecutions()) - << "%) Evts survived the padding!" << endl; - *fLog << endl; - - } - - //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("PadON", "", 900, 1200)); - c.Divide(3, 4); - gROOT->SetSelectedPad(NULL); - - c.cd(1); - fHNSB->SetDirectory(NULL); - fHNSB->DrawCopy(); - fHNSB->SetBit(kCanDelete); - - c.cd(2); - fHSigmaPedestal->SetDirectory(NULL); - fHSigmaPedestal->DrawCopy(); - fHSigmaPedestal->SetBit(kCanDelete); - - c.cd(3); - fHPhotons->SetDirectory(NULL); - fHPhotons->DrawCopy(); - fHPhotons->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - - c.cd(4); - fHSigmaTheta->SetDirectory(NULL); - fHSigmaTheta->SetTitle("(Input) 2D : Sigmabar, \\Theta"); - fHSigmaTheta->DrawCopy(); - fHSigmaTheta->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - - c.cd(7); - fHBlindPixNTheta->SetDirectory(NULL); - fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta"); - fHBlindPixNTheta->DrawCopy(); - fHBlindPixNTheta->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - - c.cd(10); - fHBlindPixIdTheta->SetDirectory(NULL); - fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta"); - fHBlindPixIdTheta->DrawCopy(); - fHBlindPixIdTheta->SetBit(kCanDelete); - - - //-------------------------------------------------------------------- - // draw the 3D histogram (input): Theta, pixel, Sigma^2-Sigmabar^2 - - c.cd(5); - TH2D *l1; - l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); - l1->SetDirectory(NULL); - l1->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); - l1->SetXTitle("\\Theta [\\circ]"); - l1->SetYTitle("Sigma^2-Sigmabar^2"); - - l1->DrawCopy("box"); - l1->SetBit(kCanDelete);; - - c.cd(8); - TH2D *l2; - l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); - l2->SetDirectory(NULL); - l2->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); - l2->SetXTitle("pixel"); - l2->SetYTitle("Sigma^2-Sigmabar^2"); - - l2->DrawCopy("box"); - l2->SetBit(kCanDelete);; - - //-------------------------------------------------------------------- - // draw the 3D histogram (input): Theta, pixel, Sigma - - c.cd(6); - TH2D *k1; - k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); - k1->SetDirectory(NULL); - k1->SetTitle("(Input) Sigma vs. \\Theta (all pixels)"); - k1->SetXTitle("\\Theta [\\circ]"); - k1->SetYTitle("Sigma"); - - k1->DrawCopy("box"); - k1->SetBit(kCanDelete); - - c.cd(9); - TH2D *k2; - k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); - k2->SetDirectory(NULL); - k2->SetTitle("(Input) Sigma vs. pixel number (all \\Theta)"); - k2->SetXTitle("pixel"); - k2->SetYTitle("Sigma"); - - k2->DrawCopy("box"); - k2->SetBit(kCanDelete);; - - - //-------------------------------------------------------------------- - - - return kTRUE; -} - -// -------------------------------------------------------------------------- - - - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MPadON.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadON.h (revision 2776) +++ (revision ) @@ -1,82 +1,0 @@ -#ifndef MARS_MPadON -#define MARS_MPadON - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1D; -class TH2D; -class TH3D; - -class MGeomCam; -class MCerPhotEvt; -class MPedestalCam; -class MMcEvt; -class MSigmabar; -class MParList; -class MBlindPixels; - -class MPadON : public MTask -{ -private: - MGeomCam *fCam; - MCerPhotEvt *fEvt; - MSigmabar *fSigmabar; - MMcEvt *fMcEvt; - MPedestalCam *fPed; - MBlindPixels *fBlinds; - - Int_t fPadFlag; - Int_t fIter; - - Int_t fRunType; - Int_t fGroup; - - Int_t fErrors[8]; - - // plots used for the padding - TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) - TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) - TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) - TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - - // plots for checking the padding - TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after - // versus before padding - TH2D *fHPhotons; // 2D-histogram : no.of photons after - // versus before padding - TH1D *fHNSB; // 1D-histogram : additional NSB - - -public: - MPadON(const char *name=NULL, const char *title=NULL); - ~MPadON(); - - void SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, - TH2D *hist2Pix, TH2D *hist2PixN); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - void SetPadFlag(Int_t padflag); - - ClassDef(MPadON, 0) // task for the padding (Schweizer) -}; - -#endif - - - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2776) +++ (revision ) @@ -1,1922 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Wolfgang Wittek, 06/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MPadONOFF -// -// This task applies padding such that for a given pixel and for a given -// Theta bin the resulting distribution of the pedestal sigma is identical -// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta. -// -// The number of photons, its error and the pedestal sigmas are altered. -// On average, the number of photons added is zero. -// -// The formulas used can be found in Thomas Schweizer's Thesis, -// Section 2.2.1 -// -// There are 2 options for the padding : -// -// 1) fPadFlag = 1 : -// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta -// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using -// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 -// (fHDiffPixTheta). -// -// This is the preferred option as it takes into account the -// correlations between the Sigma of a pixel and Sigmabar. -// -// THE FOLLOWING OPTION HAS NOT BEEN TESTED : -// 2) fPadFlag = 2 : -// Generate a pedestal sigma for each pixel using the 3D-histogram -// Theta, pixel no., Sigma (fHSigmaPixTheta). -// -// -// The padding has to be done before the image cleaning because the -// image cleaning depends on the pedestal sigmas. -// -// For random numbers gRandom is used. -// -///////////////////////////////////////////////////////////////////////////// -#include "MPadONOFF.h" - -#include -#include - -#include -#include -#include -#include -#include -#include - -#include "MBinning.h" -#include "MSigmabar.h" -#include "MMcEvt.hxx" -#include "MLog.h" -#include "MLogManip.h" -#include "MParList.h" -#include "MGeomCam.h" - -#include "MCerPhotPix.h" -#include "MCerPhotEvt.h" - -#include "MPedestalCam.h" -#include "MPedestalPix.h" -#include "MBlindPixels.h" - -#include "MRead.h" -#include "MFilterList.h" -#include "MTaskList.h" -#include "MBlindPixelCalc.h" -#include "MHBlindPixels.h" -#include "MFillH.h" -#include "MHSigmaTheta.h" -#include "MEvtLoop.h" - -ClassImp(MPadONOFF); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MPadONOFF::MPadONOFF(const char *name, const char *title) -{ - fName = name ? name : "MPadONOFF"; - fTitle = title ? title : "Task for the ON-OFF padding"; - - fPadFlag = 1; - *fLog << all << "MadONOFF : fPadFlag = " << fPadFlag << endl; - - fType = ""; - - fHSigmaTheta = NULL; - fHSigmaThetaON = NULL; - fHSigmaThetaOFF = NULL; - - fHSigmaPixTheta = NULL; - fHSigmaPixThetaON = NULL; - fHSigmaPixThetaOFF = NULL; - - fHDiffPixTheta = NULL; - fHDiffPixThetaON = NULL; - fHDiffPixThetaOFF = NULL; - - fHgON = NULL; - fHgOFF = NULL; - - fHBlindPixIdTheta = NULL; - fHBlindPixNTheta = NULL; - - fHSigmaPedestal = NULL; - fHPhotons = NULL; - fHNSB = NULL; -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MPadONOFF::~MPadONOFF() -{ - if (fHBlindPixIdTheta != NULL) delete fHBlindPixIdTheta; - if (fHBlindPixNTheta != NULL) delete fHBlindPixNTheta; - - if (fHSigmaTheta != NULL) delete fHSigmaTheta; - if (fHSigmaPixTheta != NULL) delete fHSigmaPixTheta; - if (fHDiffPixTheta != NULL) delete fHDiffPixTheta; - - if (fHSigmaPedestal != NULL) delete fHSigmaPedestal; - if (fHPhotons != NULL) delete fHPhotons; - if (fHNSB != NULL) delete fHNSB; - - if (fInfile != NULL) delete fInfile; -} - -// -------------------------------------------------------------------------- -// -// Merge the distributions of ON and OFF data -// -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) -// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) -// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) -// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) -// -// and define the target distributions for the padding -// -Bool_t MPadONOFF::MergeHistograms(TH2D *sigthon, TH2D *sigthoff, - TH3D *sigpixthon, TH3D *sigpixthoff, - TH3D *diffpixthon, TH3D *diffpixthoff, - TH2D *blindidthon, TH2D *blindidthoff, - TH2D *blindnthon, TH2D *blindnthoff) -{ - *fLog << all << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "MPadONOFF::MergeHistograms(); Merge the ON and OFF histograms to obtain the target distributions for the padding" - << endl; - - //------------------------------------------------------------- - // merge the distributions of ON and OFF events - // to obtain the target distribution fHSigmaTheta - // - - Double_t eps = 1.e-10; - - fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); - fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - - // get binning for fHgON and fHgOFF from sigthon - // binning in Theta - TAxis *ax = sigthon->GetXaxis(); - Int_t nbinstheta = ax->GetNbins(); - TArrayD edgesx; - edgesx.Set(nbinstheta+1); - for (Int_t i=0; i<=nbinstheta; i++) - { - edgesx[i] = ax->GetBinLowEdge(i+1); - *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] - << endl; - } - MBinning binth; - binth.SetEdges(edgesx); - - // binning in sigmabar - TAxis *ay = sigthon->GetYaxis(); - Int_t nbinssig = ay->GetNbins(); - TArrayD edgesy; - edgesy.Set(nbinssig+1); - for (Int_t i=0; i<=nbinssig; i++) - { - edgesy[i] = ay->GetBinLowEdge(i+1); - *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] - << endl; - } - MBinning binsg; - binsg.SetEdges(edgesy); - - - fHgON = new TH3D; - MH::SetBinning(fHgON, &binth, &binsg, &binsg); - fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - - fHgOFF = new TH3D; - MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); - fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); - - //............ loop over Theta bins ........................... - - // hista is the normalized 1D histogram of sigmabar for ON data - // histb is the normalized 1D histogram of sigmabar for OFF data - // histc is the difference ON-OFF - - // at the beginning, histap is a copy of hista - // at the end, it will be the 1D histogram for ON data after padding - - // at the beginning, histbp is a copy of histb - // at the end, it will be the 1D histogram for OFF data after padding - - // at the beginning, histcp is a copy of histc - // at the end, it should be zero - - *fLog << all << "MPadONOFF::MergeHistograms(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; - for (Int_t l=1; l<=nbinstheta; l++) - { - TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); - Stat_t suma = hista->Integral(); - hista->Scale(1./suma); - - TH1D *histap = new TH1D( (const TH1D&)*hista ); - - TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); - Stat_t sumb = histb->Integral(); - histb->Scale(1./sumb); - - TH1D *histbp = new TH1D( (const TH1D&)*histb ); - - TH1D *histc = new TH1D( (const TH1D&)*hista ); - histc->Add(histb, -1.0); - - TH1D *histcp = new TH1D( (const TH1D&)*histc ); - - - // calculate matrix g - - // fHgON[k][j] tells which fraction of ON events should be padded - // from sigma_k to sigma_j - - - // fHgOFF[k][j] tells which fraction of OFF events should be padded - // from sigma_k to sigma_j - - //-------- start j loop ------------------------------------------------ - // loop over bins in histc, starting from the end - Double_t v, s, w, t, x, u, a, b, c, arg; - - for (Int_t j=nbinssig; j >= 1; j--) - { - v = histcp->GetBinContent(j); - if ( fabs(v) < eps ) continue; - if (v >= 0.0) - s = 1.0; - else - s = -1.0; - - //................ start k loop ...................................... - // look for a bin k which may compensate the content of bin j - for (Int_t k=j-1; k >= 1; k--) - { - w = histcp->GetBinContent(k); - if ( fabs(w) < eps ) continue; - if (w >= 0.0) - t = 1.0; - else - t = -1.0; - - if (s==t) continue; - - x = v + w; - if (x >= 0.0) - u = 1.0; - else - u = -1.0; - - if (u == s) - { - arg = -w; - - if (arg <=0.0) - { - fHgON->SetBinContent (l, k, j, -arg); - fHgOFF->SetBinContent(l, k, j, 0.0); - } - else - { - fHgON->SetBinContent (l, k, j, 0.0); - fHgOFF->SetBinContent(l, k, j, arg); - } - - - *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j - << ", " << arg << endl; - - // g(k-1, j-1) = arg; - //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " - // << arg << endl; - - //...................................... - // this is for checking the procedure - if (arg < 0.0) - { - a = histap->GetBinContent(k); - histap->SetBinContent(k, a+arg); - a = histap->GetBinContent(j); - histap->SetBinContent(j, a-arg); - } - else - { - b = histbp->GetBinContent(k); - histbp->SetBinContent(k, b-arg); - b = histbp->GetBinContent(j); - histbp->SetBinContent(j, b+arg); - } - //...................................... - - histcp->SetBinContent(k, 0.0); - histcp->SetBinContent(j, x); - - //...................................... - // redefine v - v = histcp->GetBinContent(j); - if ( fabs(v) < eps ) break; - if (v >= 0.0) - s = 1.0; - else - s = -1.0; - //...................................... - - continue; - } - - arg = v; - - if (arg <=0.0) - { - fHgON->SetBinContent (l, k, j, -arg); - fHgOFF->SetBinContent(l, k, j, 0.0); - } - else - { - fHgON->SetBinContent (l, k, j, 0.0); - fHgOFF->SetBinContent(l, k, j, arg); - } - - *fLog << "l, k, j, arg = " << l << ", " << k << ", " << j - << ", " << arg << endl; - - //g(k-1, j-1) = arg; - //cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " - // << arg << endl; - - //...................................... - // this is for checking the procedure - if (arg < 0.0) - { - a = histap->GetBinContent(k); - histap->SetBinContent(k, a+arg); - a = histap->GetBinContent(j); - histap->SetBinContent(j, a-arg); - } - else - { - b = histbp->GetBinContent(k); - - histbp->SetBinContent(k, b-arg); - b = histbp->GetBinContent(j); - histbp->SetBinContent(j, b+arg); - } - //...................................... - - histcp->SetBinContent(k, x); - histcp->SetBinContent(j, 0.0); - - break; - } - //................ end k loop ...................................... - } - //-------- end j loop ------------------------------------------------ - - // check results for this Theta bin - for (Int_t j=1; j<=nbinssig; j++) - { - a = histap->GetBinContent(j); - b = histbp->GetBinContent(j); - c = histcp->GetBinContent(j); - - if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) - *fLog << err << "MPadONOFF::MergeHistograms(); inconsistency in results; a, b, c = " - << a << ", " << b << ", " << c << endl; - } - - - // fill target distribution SigmaTheta - // for this Theta bin - // - for (Int_t j=1; j<=nbinssig; j++) - { - a = histap->GetBinContent(j); - fHSigmaTheta->SetBinContent(l, j, a); - } - - delete hista; - delete histb; - delete histc; - - delete histap; - delete histbp; - delete histcp; - } - //............ end of loop over Theta bins .................... - - - // target distributions for MC - // SigmaPixTheta and DiffPixTheta - // BlindPixIdTheta and BlindPixNTheta - // are calculated as averages of the ON and OFF distributions - - // Note : in all histograms the pedestal RMS is assumed to be given in units of "number of photons" - - fHSigmaThetaON = sigthon; - fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (target)"); - - fHSigmaThetaOFF = sigthoff; - fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (target)"); - - - fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); - fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); - - fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); - fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); - - fHSigmaPixTheta = new TH3D( (TH3D&)*sigpixthon ); - fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - fHSigmaPixThetaON = sigpixthon; - fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (target)"); - - fHSigmaPixThetaOFF = sigpixthoff; - fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (target)"); - - fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); - fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - - fHDiffPixThetaON = diffpixthon; - fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (target)"); - - fHDiffPixThetaOFF = diffpixthoff; - fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (target)"); - - - for (Int_t j=1; j<=nbinstheta; j++) - { - // fraction of ON/OFF events to be padded to a sigmabar - // of the OFF/ON events : fracON, fracOFF - - Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - - TAxis *ax; - TAxis *ay; - TAxis *az; - - Double_t entON; - Double_t entOFF; - - Double_t normON; - Double_t normOFF; - - // set ranges for 2D-projections of 3D-histograms - ax = sigpixthon->GetXaxis(); - ax->SetRange(j, j); - ax = sigpixthoff->GetXaxis(); - ax->SetRange(j, j); - - ax = diffpixthon->GetXaxis(); - ax->SetRange(j, j); - ax = diffpixthoff->GetXaxis(); - ax->SetRange(j, j); - - TH2D *hist; - TH2D *histOFF; - - TH1D *hist1; - TH1D *hist1OFF; - - // weights for ON and OFF distrubtions when - // calculating the weighted average - Double_t facON = 0.5 * (1. - fracON + fracOFF); - Double_t facOFF = 0.5 * (1. + fracON - fracOFF); - - *fLog << all << fracON << ", " << fracOFF << ", " - << facON << ", " << facOFF << endl; - //------------------------------------------- - ay = blindnthon->GetYaxis(); - Int_t nbinsn = ay->GetNbins(); - - hist1 = (TH1D*)blindnthon->ProjectionY ("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = (TH1D*)blindnthoff->ProjectionY("", j, j, ""); - - // number of events in this theta bin - entON = hist1->Integral(); - entOFF = hist1OFF->Integral(); - - *fLog << all << "BlindPixNTheta : j, entON, entOFF = " << j << ", " - << entON << ", " << entOFF << endl; - - normON = entON<=0.0 ? 0.0 : 1.0/entON; - normOFF = entOFF<=0.0 ? 0.0 : 1.0/entOFF; - - hist1->Scale(normON); - hist1OFF->Scale(normOFF); - - // weighted average of ON and OFF distributions - hist1->Scale(facON); - hist1->Add(hist1OFF, facOFF); - - for (Int_t k=1; k<=nbinsn; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixNTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - - //------------------------------------------- - ay = blindidthon->GetYaxis(); - Int_t nbinsid = ay->GetNbins(); - - hist1 = (TH1D*)blindidthon->ProjectionY ("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = (TH1D*)blindidthoff->ProjectionY("", j, j, ""); - - hist1->Scale(normON); - hist1OFF->Scale(normOFF); - - // weighted average of ON and OFF distributions - hist1->Scale(facON); - hist1->Add(hist1OFF, facOFF); - - for (Int_t k=1; k<=nbinsid; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixIdTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - - //------------------------------------------- - ay = sigpixthon->GetYaxis(); - Int_t nbinspix = ay->GetNbins(); - - az = sigpixthon->GetZaxis(); - Int_t nbinssigma = az->GetNbins(); - - hist = (TH2D*)sigpixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)sigpixthoff->Project3D("zy"); - - hist->Scale(normON); - histOFF->Scale(normOFF); - - // weighted average of ON and OFF distributions - - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspix; k++) - for (Int_t l=1; l<=nbinssigma; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHSigmaPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - //------------------------------------------- - ay = diffpixthon->GetYaxis(); - Int_t nbinspixel = ay->GetNbins(); - - az = diffpixthon->GetZaxis(); - Int_t nbinsdiff = az->GetNbins(); - - hist = (TH2D*)diffpixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)diffpixthoff->Project3D("zy"); - - hist->Scale(normON); - histOFF->Scale(normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinsdiff; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHDiffPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - //------------------------------------------- - } - - - *fLog << all - << "MPadONOFF::MergeHistograms(); The target distributions for the padding have been created" - << endl; - *fLog << all << "----------------------------------------------------------" - << endl; - //-------------------------------------------- - - fHSigmaTheta->SetDirectory(NULL); - fHSigmaThetaON->SetDirectory(NULL); - fHSigmaThetaOFF->SetDirectory(NULL); - - fHSigmaPixTheta->SetDirectory(NULL); - fHSigmaPixThetaON->SetDirectory(NULL); - fHSigmaPixThetaOFF->SetDirectory(NULL); - - fHDiffPixTheta->SetDirectory(NULL); - fHDiffPixThetaON->SetDirectory(NULL); - fHDiffPixThetaOFF->SetDirectory(NULL); - - fHBlindPixIdTheta->SetDirectory(NULL); - fHBlindPixNTheta->SetDirectory(NULL); - - - fHgON->SetDirectory(NULL); - fHgOFF->SetDirectory(NULL); - - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Read target distributions from a file -// -// -Bool_t MPadONOFF::ReadTargetDist(const char* namefilein) -{ - *fLog << all << "MPadONOFF : Read padding histograms from file " - << namefilein << endl; - - fInfile = new TFile(namefilein); - fInfile->ls(); - - //------------------------------------ - - fHSigmaTheta = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); - if (!fHSigmaTheta) - { - *fLog << all - << "MPadONOFF : Object '2D-ThetaSigmabar' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '2D-ThetaSigmabar' was read in" << endl; - - fHSigmaThetaON = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); - if (!fHSigmaThetaON) - { - *fLog << all - << "MPadONOFF : Object '2D-ThetaSigmabarON' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '2D-ThetaSigmabarON' was read in" << endl; - - fHSigmaThetaOFF = - (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); - if (!fHSigmaThetaOFF) - { - *fLog << all - << "MPadONOFF : Object '2D-ThetaSigmabarOFF' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF :Object '2D-ThetaSigmabarOFF' was read in" << endl; - - fHSigmaPixTheta = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); - if (!fHSigmaPixTheta) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigma' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigma' was read in" << endl; - - fHSigmaPixThetaON = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); - if (!fHSigmaPixThetaON) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigmaON' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigmaON' was read in" << endl; - - fHSigmaPixThetaOFF = - (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); - if (!fHSigmaPixThetaOFF) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' was read in" << endl; - - fHDiffPixTheta = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); - if (!fHDiffPixTheta) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiff' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiff' was read in" << endl; - - fHDiffPixThetaON = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); - if (!fHDiffPixThetaON) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiffON' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiffON' was read in" << endl; - - fHDiffPixThetaOFF = - (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); - if (!fHDiffPixThetaOFF) - { - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiffOFF' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-ThetaPixDiffOFF' was read in" << endl; - - fHgON = - (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); - if (!fHgON) - { - *fLog << all - << "MPadONOFF : Object '3D-PaddingMatrixON' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-PaddingMatrixON' was read in" << endl; - - fHgOFF = - (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); - if (!fHgOFF) - { - *fLog << all - << "MPadONOFF : Object '3D-PaddingMatrixOFF' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '3D-PaddingMatrixOFF' was read in" << endl; - - - fHBlindPixIdTheta = - (TH2D*) gROOT->FindObject("2D-ThetaBlindId"); - if (!fHBlindPixIdTheta) - { - *fLog << all - << "MPadONOFF : Object '2D-ThetaBlindId' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '2D-ThetaBlindId' was read in" << endl; - - fHBlindPixNTheta = - (TH2D*) gROOT->FindObject("2D-ThetaBlindN"); - if (!fHBlindPixNTheta) - { - *fLog << all - << "MPadONOFF : Object '2D-ThetaBlindN' not found on root file" - << endl; - return kFALSE; - } - *fLog << all - << "MPadONOFF : Object '2D-ThetaBlindN' was read in" << endl; - - - //------------------------------------ - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Write target distributions onto a file -// -// -Bool_t MPadONOFF::WriteTargetDist(const char* namefileout) -{ - *fLog << all << "MPadONOFF : Write padding histograms onto file " - << namefileout << endl; - - TFile outfile(namefileout, "RECREATE"); - - fHBlindPixNTheta->Write(); - fHBlindPixIdTheta->Write(); - - fHSigmaTheta->Write(); - fHSigmaThetaON->Write(); - fHSigmaThetaOFF->Write(); - - fHSigmaPixTheta->Write(); - fHSigmaPixThetaON->Write(); - fHSigmaPixThetaOFF->Write(); - - fHDiffPixTheta->Write(); - fHDiffPixThetaON->Write(); - fHDiffPixThetaOFF->Write(); - - fHgON->Write(); - fHgOFF->Write(); - - *fLog << all - << "MPadONOFF::WriteTargetDist(); target histograms written onto file " - << namefileout << endl; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Set type of data to be padded -// -// this is not necessary if the type of the data can be recognized -// directly from the input -// -// -void MPadONOFF::SetDataType(const char* type) -{ - fType = type; - *fLog << all << "MPadONOFF::SetDataType(); type of data to be padded : " - << fType << endl; - - return; -} - - -// -------------------------------------------------------------------------- -// -// Set the option for the padding -// -// There are 2 options for the padding : -// -// 1) fPadFlag = 1 : -// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta -// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using -// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 -// (fHDiffPixTheta). -// -// This is the preferred option as it takes into account the -// correlations between the Sigma of a pixel and Sigmabar. -// -// NOT YET TESTED : -// 2) fPadFlag = 2 : -// Generate a pedestal sigma for each pixel using the 3D-histogram -// Theta, pixel no., Sigma (fHSigmaPixTheta). -// -void MPadONOFF::SetPadFlag(Int_t padflag) -{ - fPadFlag = padflag; - *fLog << all << "MPadONOFF::SetPadFlag(); choose option " << fPadFlag - << endl; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers and prepare the histograms -// -Int_t MPadONOFF::PreProcess(MParList *pList) -{ - if ( !fHSigmaTheta || !fHSigmaThetaON || !fHSigmaThetaOFF || - !fHSigmaPixTheta || !fHSigmaPixThetaON || !fHSigmaPixThetaOFF || - !fHDiffPixTheta || !fHDiffPixThetaON || !fHDiffPixThetaOFF || - !fHBlindPixIdTheta || !fHBlindPixNTheta || - !fHgON || !fHgOFF) - { - *fLog << err - << "MPadONOFF : At least one of the histograms needed for the padding is not defined ... aborting." - << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MPadONOFF : MMcEvt not found... aborting." - << endl; - return kFALSE; - } - - fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); - if (!fPed) - { - *fLog << err << "MPadONOFF : MPedestalCam not found... aborting." - << endl; - return kFALSE; - } - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << err - << "MPadONOFF : MGeomCam not found (no geometry information available)... aborting." - << endl; - return kFALSE; - } - - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << err << "MPadONOFF : MCerPhotEvt not found... aborting." - << endl; - return kFALSE; - } - - fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); - if (!fSigmabar) - { - *fLog << err << "MPadONOFF : MSigmabar not found... aborting." << endl; - return kFALSE; - } - - fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); - if (!fBlinds) - { - *fLog << err << "MPadONOFF : MBlindPixels not found... aborting." - << endl; - return kFALSE; - } - - if (fType !="ON" && fType !="OFF" && fType !="MC") - { - *fLog << err - << "MPadONOFF : Type of data to be padded not defined... aborting." - << endl; - return kFALSE; - } - - - //-------------------------------------------------------------------- - // histograms for checking the padding - // - // plot pedestal sigmas - fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", - 100, 0.0, 5.0, 100, 0.0, 5.0); - fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); - fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); - - // plot no.of photons (before vs. after padding) - fHPhotons = new TH2D("Photons","Photons: after vs.before padding", - 100, -10.0, 90.0, 100, -10, 90); - fHPhotons->SetXTitle("no.of photons before padding"); - fHPhotons->SetYTitle("no.of photons after padding"); - - // plot of added NSB - fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); - fHNSB->SetXTitle("no.of added NSB photons"); - fHNSB->SetYTitle("no.of pixels"); - - - //-------------------------------------------------------------------- - - fIter = 20; - - memset(fErrors, 0, sizeof(fErrors)); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Do the Padding -// idealy the events to be padded should have been generated without NSB -// -Int_t MPadONOFF::Process() -{ - *fLog << all << "Entry MPadONOFF::Process();" << endl; - - // this is the conversion factor from the number of photons - // to the number of photo electrons - // later to be taken from MCalibration - // fPEperPhoton = PW * LG * QE * 1D - Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; - - //------------------------------------------------ - // add pixels to MCerPhotEvt which are not yet in; - // set their number of photons equal to zero - - UInt_t imaxnumpix = fCam->GetNumPixels(); - - for (UInt_t i=0; iGetNumPixels(); - for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); - } - - - - //----------------------------------------- - Int_t rc=0; - Int_t rw=0; - - const UInt_t npix = fEvt->GetNumPixels(); - - //fSigmabar->Calc(*fCam, *fPed, *fEvt); - // *fLog << all << "before padding : " << endl; - //fSigmabar->Print(""); - - - //$$$$$$$$$$$$$$$$$$$$$$$$$$ - // to simulate the situation that before the padding the NSB and - // electronic noise are zero : set Sigma = 0 for all pixels - //for (UInt_t i=0; ioperator[](i); - // Int_t j = pix.GetPixId(); - - // MPedestalPix &ppix = fPed->operator[](j); - // ppix.SetMeanRms(0.0); - //} - //$$$$$$$$$$$$$$$$$$$$$$$$$$ - - //------------------------------------------- - // Calculate average sigma of the event - // - Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPadONOFF::Process(); before padding : " << endl; - fSigmabar->Print(""); - Double_t sigbarold = sigbarold_phot * fPEperPhoton; - - - Double_t sigbarold2 = sigbarold*sigbarold; - - - // for MC data : expect sigmabar to be zero before the padding - if (fType == "MC") - { - if (sigbarold_phot > 0) - { - *fLog << err - << "MPadONOFF::Process(); sigmabar of event to be padded is > 0 : " - << sigbarold_phot << ". Stop event loop " << endl; - // input data should have sigmabar = 0; stop event loop - - rc = 1; - fErrors[rc]++; - return kERROR; - } - } - - const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - *fLog << all << "theta = " << theta << endl; - - Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); - if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) - { - *fLog << warn - << "MPadONOFF::Process(); binNumber out of range : theta, binTheta = " - << theta << ", " << binTheta << "; Skip event " << endl; - // event cannot be padded; skip event - - rc = 2; - fErrors[rc]++; - return kCONTINUE; - } - - //------------------------------------------- - // get number of events in this theta bin - // and number of events in this sigbarold_phot bin - - Double_t nTheta; - Double_t nSigma; - if (fType == "ON") - { - TH1D *hn; - - hn = fHSigmaThetaON->ProjectionY("", binTheta, binTheta, ""); - nTheta = hn->Integral(); - Int_t binSigma = hn->FindBin(sigbarold_phot); - nSigma = hn->GetBinContent(binSigma); - - *fLog << all - << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " - << theta << ", " << sigbarold_phot << ", " << binTheta << ", " - << binSigma << ", " << nTheta << ", " << nSigma << endl; - - delete hn; - } - - else if (fType == "OFF") - { - TH1D *hn; - - hn = fHSigmaThetaOFF->ProjectionY("", binTheta, binTheta, ""); - nTheta = hn->Integral(); - Int_t binSigma = hn->FindBin(sigbarold_phot); - nSigma = hn->GetBinContent(binSigma); - - *fLog << all - << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " - << theta << ", " << sigbarold_phot << ", " << binTheta << ", " - << binSigma << ", " << nTheta << ", " << nSigma << endl; - - delete hn; - } - - else - { - nTheta = 0.0; - nSigma = 0.0; - } - - //------------------------------------------- - // for the current theta, - // generate blind pixels according to the histograms - // fHBlindPixNTheta and fHBlindPixIDTheta - // do this only for MC data - - - if (fType == "MC") - { - - // numBlind is the number of blind pixels in this event - TH1D *nblind; - UInt_t numBlind; - - nblind = fHBlindPixNTheta->ProjectionY("", binTheta, binTheta, ""); - if ( nblind->Integral() == 0.0 ) - { - *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " - << binTheta << " has no entries; Skip event " << endl; - // event cannot be padded; skip event - delete nblind; - - rc = 7; - fErrors[rc]++; - return kCONTINUE; - } - else - { - numBlind = (Int_t) (nblind->GetRandom()+0.5); - } - delete nblind; - - //warn Code commented out due no compilation errors on Alpha OSF (tgb) - - // throw the Id of numBlind different pixels in this event - TH1D *hblind; - UInt_t idBlind; - UInt_t listId[npix]; - UInt_t nlist = 0; - Bool_t equal; - - hblind = fHBlindPixIdTheta->ProjectionY("", binTheta, binTheta, ""); - if ( hblind->Integral() > 0.0 ) - for (UInt_t i=0; iGetRandom()+0.5); - equal = kFALSE; - for (UInt_t j=0; jSetPixelBlind(idBlind); - *fLog << all << "idBlind = " << idBlind << endl; - } - fBlinds->SetReadyToSave(); - - delete hblind; - - } - - //****************************************************************** - // has the event to be padded ? - // if yes : to which sigmabar should it be padded ? - // - - Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); - *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " - << sigbarold_phot << endl; - - Double_t prob; - TH1D *hpad = NULL; - if (fType == "ON") - { - hpad = fHgON->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, ""); - - //Int_t nb = hpad->GetNbinsX(); - //for (Int_t i=1; i<=nb; i++) - //{ - // if (hpad->GetBinContent(i) != 0.0) - // *fLog << all << "i, fHgON = " << i << ", " - // << hpad->GetBinContent(i) << endl; - //} - - if (nSigma != 0.0) - prob = hpad->Integral() * nTheta/nSigma; - else - prob = 0.0; - - *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma - << ", " << prob << endl; - } - else if (fType == "OFF") - { - hpad = fHgOFF->ProjectionZ("zOFF", binTheta, binTheta, binSig, binSig, ""); - if (nSigma != 0.0) - prob = hpad->Integral() * nTheta/nSigma; - else - prob = 0.0; - } - else - prob = 1.0; - - if ( fType != "MC" && - (prob <= 0.0 || gRandom->Uniform() > prob) ) - { - delete hpad; - // event should not be padded - *fLog << all << "event is not padded" << endl; - - rc = 8; - fErrors[rc]++; - return kTRUE; - } - // event should be padded - *fLog << all << "event is padded" << endl; - - - //------------------------------------------- - // for the current theta, generate a sigmabar - // - // for ON/OFF data according to the matrix fHgON/OFF - // for MC data according to the histogram fHSigmaTheta - // - Double_t sigmabar_phot = 0; - Double_t sigmabar = 0; - if (fType == "ON" || fType == "OFF") - { - //Int_t nbinsx = hpad->GetNbinsX(); - //for (Int_t i=1; i<=nbinsx; i++) - //{ - // if (hpad->GetBinContent(i) != 0.0) - // *fLog << all << "i, fHg = " << i << ", " << hpad->GetBinContent(i) - // << endl; - //} - - sigmabar_phot = hpad->GetRandom(); - sigmabar = sigmabar_phot * fPEperPhoton; - - *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; - - delete hpad; - } - - else if (fType == "MC") - { - Int_t bincheck = fHSigmaTheta->GetXaxis()->FindBin(theta); - - if (binTheta != bincheck) - { - cout << "The binnings of the 2 histograms are not identical; aborting" - << endl; - return kERROR; - } - - TH1D *hsigma; - - hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, ""); - if ( hsigma->GetEntries() == 0.0 ) - { - *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " - << binTheta << " has no entries; Skip event " << endl; - // event cannot be padded; skip event - delete hsigma; - - rc = 3; - fErrors[rc]++; - return kCONTINUE; - } - else - { - sigmabar_phot = hsigma->GetRandom(); - sigmabar = sigmabar_phot * fPEperPhoton; - - *fLog << all << "Theta, bin number = " << theta << ", " << binTheta - << ", sigmabar_phot = " << sigmabar_phot << endl; - } - delete hsigma; - } - - const Double_t sigmabar2 = sigmabar*sigmabar; - - //------------------------------------------- - - *fLog << all << "MPadONOFF::Process(); sigbarold, sigmabar = " - << sigbarold << ", "<< sigmabar << endl; - - // Skip event if target sigmabar is <= sigbarold - if (sigmabar <= sigbarold) - { - *fLog << all << "MPadONOFF::Process(); target sigmabar is less than sigbarold : " - << sigmabar << ", " << sigbarold << ", Skip event" << endl; - - rc = 4; - fErrors[rc]++; - return kCONTINUE; - } - - - //------------------------------------------- - // - // Calculate average number of NSB photo electrons to be added (lambdabar) - // from the value of sigmabar, - // - making assumptions about the average electronic noise (elNoise2) and - // - using a fixed value (F2excess) for the excess noise factor - - Double_t elNoise2; // [photo electrons] - Double_t F2excess = 1.3; - Double_t lambdabar; // [photo electrons] - - - - Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta); - if (binTheta != bincheck) - { - *fLog << err - << "MPadONOFF::Process(); The binnings of the 2 histograms are not identical; aborting" - << endl; - return kERROR; - } - - // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. - // The average electronic noise (to be added) has to be well above this RMS, - // otherwise the electronic noise of an individual pixel (elNoise2Pix) - // may become negative - - TH1D *hnoise; - if (fType == "MC") - hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); - else if (fType == "ON") - hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); - else if (fType == "OFF") - hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); - else - { - *fLog << err << "MPadONOFF::Process(); illegal data type... aborting" - << endl; - return kERROR; - } - - Double_t RMS_phot = hnoise->GetRMS(1); - Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; - delete hnoise; - - elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); - *fLog << all << "elNoise2 = " << elNoise2 << endl; - - lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; - - // This value of lambdabar is the same for all pixels; - // note that lambdabar is normalized to the area of pixel 0 - - //---------- start loop over pixels --------------------------------- - // do the padding for each pixel - // - // pad only pixels - which are used (before image cleaning) - // - Double_t sig_phot = 0; - Double_t sig = 0; - - Double_t sigma2 = 0; - - Double_t diff_phot = 0; - Double_t diff = 0; - - Double_t addSig2_phot= 0; - Double_t addSig2 = 0; - - Double_t elNoise2Pix = 0; - - - for (UInt_t i=0; iGetPixRatio(j); - - MPedestalPix &ppix = (*fPed)[j]; - Double_t oldsigma_phot = ppix.GetPedestalRms(); - Double_t oldsigma = oldsigma_phot * fPEperPhoton; - Double_t oldsigma2 = oldsigma*oldsigma; - - //--------------------------------- - // throw the Sigma for this pixel - // - Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); - - Int_t count; - Bool_t ok; - - TH1D *hdiff; - TH1D *hsig; - - switch (fPadFlag) - { - case 1 : - // throw the Sigma for this pixel from the distribution fHDiffPixTheta - - if (fType == "MC") - hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - else if (fType == "ON") - hdiff = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - else - hdiff = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - - if ( hdiff->GetEntries() == 0 ) - { - *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " - << binTheta << " and pixel bin " << binPixel - << " has no entries; aborting " << endl; - delete hdiff; - - rc = 5; - fErrors[rc]++; - return kCONTINUE; - } - - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - diff = diff_phot * fPEperPhoton * fPEperPhoton; - - // the following condition ensures that elNoise2Pix > 0.0 - if ( (diff + sigmabar2 - oldsigma2/ratioArea - - lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - if (!ok) - { - *fLog << all << "theta, j, count, sigmabar, diff = " << theta - << ", " - << j << ", " << count << ", " << sigmabar << ", " - << diff << endl; - diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; - - rw = 1; - fWarnings[rw]++; - } - else - { - rw = 0; - fWarnings[rw]++; - } - - delete hdiff; - sigma2 = diff + sigmabar2; - break; - - case 2 : - // throw the Sigma for this pixel from the distribution fHSigmaPixTheta - - if (fType == "MC") - hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - else if (fType == "ON") - hsig = fHSigmaPixThetaOFF->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - else - hsig = fHSigmaPixThetaON->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - - if ( hsig->GetEntries() == 0 ) - { - *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " - << binTheta << " and pixel bin " << binPixel - << " has no entries; aborting " << endl; - delete hsig; - - rc = 6; - fErrors[rc]++; - return kCONTINUE; - } - - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - sig = sig_phot * fPEperPhoton; - - sigma2 = sig*sig/ratioArea; - // the following condition ensures that elNoise2Pix > 0.0 - - if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - if (!ok) - { - *fLog << all << "theta, j, count, sigmabar, sig = " << theta - << ", " - << j << ", " << count << ", " << sigmabar << ", " - << sig << endl; - sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; - - rw = 1; - fWarnings[rw]++; - } - else - { - rw = 0; - fWarnings[rw]++; - } - - delete hsig; - break; - } - - //--------------------------------- - // get the additional sigma^2 for this pixel (due to the padding) - - addSig2 = sigma2*ratioArea - oldsigma2; - addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); - - //--------------------------------- - // get the additional electronic noise for this pixel - - elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; - - - //--------------------------------- - // throw actual number of additional NSB photo electrons (NSB) - // and its RMS (sigmaNSB) - - Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); - Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; - Double_t sigmaNSB0; - - if (arg >= 0.0) - { - sigmaNSB0 = sqrt( arg ); - } - else - { - sigmaNSB0 = 0.0000001; - if (arg < -1.e-10) - { - *fLog << warn << "MPadONOFF::Process(); argument of sqrt < 0 : " - << arg << endl; - } - } - - - //--------------------------------- - // smear NSB0 according to sigmaNSB0 - // and subtract lambdabar because of AC coupling - - Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; - Double_t NSB_phot = NSB / fPEperPhoton; - - //--------------------------------- - - // add additional NSB to the number of photons - Double_t oldphotons_phot = pix.GetNumPhotons(); - Double_t oldphotons = oldphotons_phot * fPEperPhoton; - Double_t newphotons = oldphotons + NSB; - Double_t newphotons_phot = newphotons / fPEperPhoton; - pix.SetNumPhotons( newphotons_phot ); - - - fHNSB->Fill( NSB_phot/ratioArea ); - fHPhotons->Fill( oldphotons_phot/ratioArea, - newphotons_phot/ratioArea ); - - - // error: add sigma of padded noise quadratically - Double_t olderror_phot = pix.GetErrorPhot(); - Double_t newerror_phot = - sqrt( olderror_phot*olderror_phot + addSig2_phot ); - pix.SetErrorPhot( newerror_phot ); - - - Double_t newsigma = sqrt( oldsigma2 + addSig2 ); - Double_t newsigma_phot = newsigma / fPEperPhoton; - ppix.SetPedestalRms( newsigma_phot ); - - fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); - } - //---------- end of loop over pixels --------------------------------- - - // Calculate sigmabar again and crosscheck - - fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPadONOFF::Process(); after padding : " << endl; - fSigmabar->Print(""); - - *fLog << all << "Exit MPadONOFF::Process();" << endl; - - rc = 0; - fErrors[rc]++; - return kTRUE; - //****************************************************************** -} - -// -------------------------------------------------------------------------- -// -// -Int_t MPadONOFF::PostProcess() -{ - if (GetNumExecutions() != 0) - { - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - - int temp; - if (fWarnings[0] != 0.0) - temp = (int)(fWarnings[1]*100/fWarnings[0]); - else - temp = 100; - - *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) - << temp - << "%) Warning : iteration to find acceptable sigma failed" - << ", fIter = " << fIter << endl; - - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) - << (int)(fErrors[1]*100/GetNumExecutions()) - << "%) Evts skipped due to: Sigmabar_old > 0" << endl; - - *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) - << (int)(fErrors[2]*100/GetNumExecutions()) - << "%) Evts skipped due to: Zenith angle out of range" << endl; - - *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) - << (int)(fErrors[3]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigmabar" << endl; - - *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) - << (int)(fErrors[4]*100/GetNumExecutions()) - << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; - - *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) - << (int)(fErrors[5]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; - - *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) - << (int)(fErrors[6]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma" << endl; - - *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) - << (int)(fErrors[7]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Blind pixels" << endl; - - *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) - << (int)(fErrors[8]*100/GetNumExecutions()) - << "%) Evts didn't have to be padded" << endl; - - *fLog << " " << fErrors[0] << " (" - << (int)(fErrors[0]*100/GetNumExecutions()) - << "%) Evts were successfully padded!" << endl; - *fLog << endl; - - } - - //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("PadONOFF", "", 900, 1200)); - c.Divide(3, 4); - gROOT->SetSelectedPad(NULL); - - c.cd(1); - fHNSB->SetDirectory(NULL); - fHNSB->DrawCopy(); - fHNSB->SetBit(kCanDelete); - - c.cd(2); - fHSigmaPedestal->SetDirectory(NULL); - fHSigmaPedestal->DrawCopy(); - fHSigmaPedestal->SetBit(kCanDelete); - - c.cd(3); - fHPhotons->SetDirectory(NULL); - fHPhotons->DrawCopy(); - fHPhotons->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - - c.cd(4); - fHSigmaTheta->SetDirectory(NULL); - fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta"); - fHSigmaTheta->DrawCopy(); - fHSigmaTheta->SetBit(kCanDelete); - - - //-------------------------------------------------------------------- - - - c.cd(7); - fHBlindPixNTheta->SetDirectory(NULL); - fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta"); - fHBlindPixNTheta->DrawCopy(); - fHBlindPixNTheta->SetBit(kCanDelete); - - //-------------------------------------------------------------------- - - - c.cd(10); - fHBlindPixIdTheta->SetDirectory(NULL); - fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta"); - fHBlindPixIdTheta->DrawCopy(); - fHBlindPixIdTheta->SetBit(kCanDelete); - - - //-------------------------------------------------------------------- - // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 - - c.cd(5); - TH2D *l1; - l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); - l1->SetDirectory(NULL); - l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); - l1->SetXTitle("\\Theta [\\circ]"); - l1->SetYTitle("Sigma^2-Sigmabar^2"); - - l1->DrawCopy("box"); - l1->SetBit(kCanDelete);; - - c.cd(8); - TH2D *l2; - l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); - l2->SetDirectory(NULL); - l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); - l2->SetXTitle("pixel"); - l2->SetYTitle("Sigma^2-Sigmabar^2"); - - l2->DrawCopy("box"); - l2->SetBit(kCanDelete);; - - //-------------------------------------------------------------------- - // draw the 3D histogram (target): Theta, pixel, Sigma - - c.cd(6); - TH2D *k1; - k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); - k1->SetDirectory(NULL); - k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)"); - k1->SetXTitle("\\Theta [\\circ]"); - k1->SetYTitle("Sigma"); - - k1->DrawCopy("box"); - k1->SetBit(kCanDelete); - - c.cd(9); - TH2D *k2; - k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); - k2->SetDirectory(NULL); - k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)"); - k2->SetXTitle("pixel"); - k2->SetYTitle("Sigma"); - - k2->DrawCopy("box"); - k2->SetBit(kCanDelete);; - - - //-------------------------------------------------------------------- - - c.cd(11); - TH2D *m1; - m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); - m1->SetDirectory(NULL); - m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); - m1->SetXTitle("Sigmabar-old"); - m1->SetYTitle("Sigmabar-new"); - - m1->DrawCopy("box"); - m1->SetBit(kCanDelete);; - - c.cd(12); - TH2D *m2; - m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); - m2->SetDirectory(NULL); - m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); - m2->SetXTitle("Sigmabar-old"); - m2->SetYTitle("Sigmabar-new"); - - m2->DrawCopy("box"); - m2->SetBit(kCanDelete);; - - - return kTRUE; -} - -// -------------------------------------------------------------------------- - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.h (revision 2776) +++ (revision ) @@ -1,105 +1,0 @@ -#ifndef MARS_MPadONOFF -#define MARS_MPadONOFF - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1D; -class TH2D; -class TH3D; - -class MGeomCam; -class MCerPhotEvt; -class MPedestalCam; -class MMcEvt; -class MSigmabar; -class MParList; -class MBlindPixels; -class MRead; -class MFilterList; - - -class MPadONOFF : public MTask -{ -private: - MGeomCam *fCam; - MCerPhotEvt *fEvt; - MSigmabar *fSigmabar; - MMcEvt *fMcEvt; - MPedestalCam *fPed; - MBlindPixels *fBlinds; - - TString fType; // type of data to be padded - TFile *fInfile; // input file containing padding histograms - - Int_t fPadFlag; - Int_t fIter; - - Int_t fErrors[9]; - Int_t fWarnings[2]; - - - // plots used for the padding - // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" - TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) - TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) - - TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) - - TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) - - TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - - TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data - TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data - - // plots for checking the padding - TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after - // versus before padding - TH2D *fHPhotons; // 2D-histogram : no.of photons after - // versus before padding - TH1D *fHNSB; // 1D-histogram : additional NSB - - -public: - MPadONOFF(const char *name=NULL, const char *title=NULL); - ~MPadONOFF(); - - Bool_t MergeHistograms(TH2D *sigthon, TH2D *sigthoff, - TH3D *sigpixthon, TH3D *sigpixthoff, - TH3D *diffpixthon, TH3D *diffpixthoff, - TH2D *blindidthon, TH2D *blindidthoff, - TH2D *blindnthon, TH2D *blindnthoff); - - Bool_t ReadTargetDist(const char *filein); - Bool_t WriteTargetDist(const char *fileout); - - void SetDataType(const char *type); // type of data to be padded - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - void SetPadFlag(Int_t padflag); - - ClassDef(MPadONOFF, 0) // task for the ON-OFF padding -}; - -#endif - - - - - - Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2776) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2777) @@ -74,6 +74,5 @@ MFiltercutsCalc.cc \ MCT1PadONOFF.cc \ - MPadON.cc \ - MPadONOFF.cc \ + MPad.cc \ MPedestalWorkaround.cc \ MExtractedSignalCam.cc \