/* ======================================================================== *\ ! ! * ! * 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 "MPedPhotCam.h" #include "MPedPhotPix.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 = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); if (!fPed) { *fLog << err << "MPad : MPedPhotCam 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].GetRms()); } //----------------------------------------- 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); MPedPhotPix &ppix = (*fPed)[j]; Double_t oldsigma_phot = ppix.GetRms(); 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.SetRms( 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; } // --------------------------------------------------------------------------