/* ======================================================================== *\ ! ! * ! * 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 #include "MBinning.h" #include "MPointingPos.h" #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MGeomCam.h" #include "MCerPhotPix.h" #include "MCerPhotEvt.h" #include "MH.h" #include "MPedPhotCam.h" #include "MPedPhotPix.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MRead.h" #include "MFilterList.h" #include "MTaskList.h" #include "MHBadPixels.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; fHSigmaPedestal = NULL; fHPhotons = NULL; fHNSB = NULL; fNamePedPhotCam = "MPedPhotCamFromData"; } // -------------------------------------------------------------------------- // // Destructor. // MPad::~MPad() { delete fHgMC; delete fHgON; delete fHgOFF; delete fHSigmaTheta; delete fHSigmaPixTheta; delete fHDiffPixTheta; 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 //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl; 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; //*fLog << "Content of histcp : " << endl; for (Int_t j=nbinssig; j >= 1; j--) { //*fLog << "j, hista, histb , histap : " << j << ", " // << hista->GetBinContent(j) << ", " // << histb->GetBinContent(j) << ", " // << histap->GetBinContent(j) << endl; v = histcp->GetBinContent(j); //*fLog << "cp : j, v = " << j << ", " << v << endl; 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 //*fLog << "Content of histap, histbp, histcp : " << endl; 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); //*fLog << "j, a, b, c = " << j << ": " << a << ", " << b << ", " // << c << endl; 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 *pad = MH::MakeDefCanvas("Merging", "", 600, 600); gROOT->SetSelectedPad(NULL); pad->Divide(2, 2); pad->cd(1); gPad->SetBorderMode(0); hista->SetDirectory(NULL); hista->Draw(); hista->SetBit(kCanDelete); pad->cd(2); gPad->SetBorderMode(0); histb->SetDirectory(NULL); histb->Draw(); histb->SetBit(kCanDelete); pad->cd(3); gPad->SetBorderMode(0); histap->SetDirectory(NULL); histap->Draw(); histap->SetBit(kCanDelete); pad->DrawClone(); //-------------------------------------------------------------------- delete histc; delete histbp; delete histcp; //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl; 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& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff) { //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ // for testing /* TAxis *xa = sigthon.GetXaxis(); Int_t nbitheta = xa->GetNbins(); TAxis *ya = sigthon.GetYaxis(); Int_t nbisig = ya->GetNbins(); for (Int_t i=1; i<=nbitheta; i++) { for (Int_t j=1; j<=nbisig; j++) { if (i>0) { sigthmc.SetBinContent( i, j, (Float_t) (625000.0+(nbisig*nbisig*nbisig-j*j*j)) ); sigthon.SetBinContent( i, j, (Float_t)(1.0) ); sigthoff.SetBinContent( i, j, (Float_t)( (0.5*nbisig-j)*(0.5*nbisig-j) ) ); } else { sigthmc.SetBinContent( i, j, 0.0); sigthon.SetBinContent( i, j, 0.0); sigthoff.SetBinContent(i, j, 0.0); } } } */ //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ *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( (const TH2D&) sigthon ); fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); *fLog << "after booking fHSigmaTheta" << endl; fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon ); fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); *fLog << "after booking fHDiffPixTheta" << endl; fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon ); fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); *fLog << "after booking fHSigmaPixTheta" << endl; *fLog << all << "Histograms for the merged padding plots were booked" << endl; //-------------------------- 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.)"); //*fLog << all << "addresses of SigmaTheta padding plots were copied" // << endl; //-------------------------- 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.)"); //*fLog << all << "addresses of DiffPixTheta padding plots were copied" // << endl; //-------------------------- 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.)"); //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" // << endl; //-------------------------- // 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"); //*fLog << all << "fHg histograms were booked" << endl; //............ 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++) { //------------------------------------------- // 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 TH2D * fHga = new TH2D; MH::SetBinning(fHga, &binsg, &binsg); TH2D * fHgb = new TH2D; MH::SetBinning(fHgb, &binsg, &binsg); TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); hista->SetNameTitle("ON-orig", "ON (orig)"); Stat_t suma = hista->Integral(); if (suma != 0.0) hista->Scale(1./suma); TH1D *histap = new TH1D( (const TH1D&)*hista ); histap->SetNameTitle("ON-OFF-merged)", "ON-OFF (merged)"); TH1D *histb = sigthoff.ProjectionY("sigoff_y", l, l, ""); histb->SetNameTitle("OFF-orig", "OFF (orig)"); Stat_t sumb = histb->Integral(); if (sumb != 0.0) histb->Scale(1./sumb); if (suma == 0.0 || sumb == 0.0) { delete hista; delete histb; delete fHga; delete fHgb; delete histap; *fLog << err << "cannot call Merge2Distributions(ON, OFF) for theta bin l = " << l << "; NevON, NevOFF = " << suma << ", " << sumb << endl; // go to next theta bin (l) continue; } //*fLog << "call Merge2Distributions(ON, OFF) for theta bin l = " // << l << endl; Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); // 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); fHgON->SetBinContent(l, k, j, a); Double_t b = fHgb->GetBinContent(k,j); fHgOFF->SetBinContent(l, k, j, b); } //------------------------------------------- // now merge ON-OFF and MC distributions // input : histe is the result of the merging of ON and OFF distributions // histf is the normalized 1D distr. of sigmabar for MC data // output : histep will be the merged distribution (target distribution) // fHge(k,j) will tell which fraction of ON, OFF events should be // padded from bin k to bin j // fHgf(k,j) will tell which fraction of MC events should be // padded from bin k to bin j TH2D * fHge = new TH2D; MH::SetBinning(fHge, &binsg, &binsg); TH2D * fHgf = new TH2D; MH::SetBinning(fHgf, &binsg, &binsg); TH1D *histe = new TH1D( (const TH1D&)*histap); histe->SetNameTitle("ON-OFF-merged", "ON-OFF (merged)"); TH1D *histep = new TH1D( (const TH1D&)*histe); histep->SetNameTitle("ON-OFF-MC-merged)", "ON-OFF-MC (merged)"); TH1D *histf = sigthmc.ProjectionY("sigmc_y", l, l, ""); histf->SetNameTitle("MC-orig", "MC (orig)"); Double_t sumf = histf->Integral(); if (sumf != 0.0) histf->Scale(1./sumf); if (sumf == 0.0) { delete hista; delete histb; delete fHga; delete fHgb; delete histap; delete histe; delete histf; delete fHge; delete fHgf; delete histep; *fLog << err << "cannot call Merge2Distributions(ON-OFF,MC) for theta bin l = " << l << "; NevMC = " << sumf << endl; // go to next theta bin (l) continue; } //*fLog << "call Merge2Distributions(ON-OFF, MC) for theta bin l = " // << l << endl; Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig); // update fHgON and fHgOFF UpdateHg(fHga, histap, fHge, fHgON, nbinssig, l); UpdateHg(fHgb, histap, fHge, fHgOFF, nbinssig, l); // fill fHgMC for (Int_t k=1; k<=nbinssig; k++) for (Int_t j=1; j<=nbinssig; j++) { Double_t f = fHgf->GetBinContent(k,j); fHgMC->SetBinContent(l, k, j, f); } // 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 ep = histep->GetBinContent(j); fHSigmaTheta->SetBinContent(l, j, ep); } //------------------------------------------- delete hista; delete histb; delete fHga; delete fHgb; delete histap; delete histe; delete histf; delete fHge; delete fHgf; delete histep; } //............ end of loop over Theta bins .................... // Define the target distributions // fHDiffPixTheta, fHSigmaPixTheta // (they are calculated as // averages 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; // 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, fracMC, 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; } //............ 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); fHgMC->SetDirectory(NULL); fHgON->SetDirectory(NULL); fHgOFF->SetDirectory(NULL); return kTRUE; } // -------------------------------------------------------------------------- // // Update the matrix fHgA // which contains the final padding prescriptions // Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, Int_t nbinssig, Int_t l) { // histap target distribution after ON-OFF merging // fHga padding prescription for ON (or OFF) data to get to histap // fHge padding prescription to get from histap to the target // distribution after the ON-OFF-MC merging // fHgA updated padding prescription for ON (or OFF) data to get // from the original ON (or OFF) histogram to the target // distribution after the ON-OFF-MC merging // l Theta bin Double_t eps = 1.e-10; for (Int_t k=1; k<=nbinssig; k++) { Double_t na = fHga->Integral(1, nbinssig, k, k, " "); Double_t ne = fHge->Integral(k, k, 1, nbinssig, " "); Double_t nap = histap->GetBinContent(k); if (ne <= eps) { // go to next k continue; } Double_t r = nap - na; if (r >= ne-eps) { for (Int_t j=k+1; j<=nbinssig; j++) { Double_t e = fHge->GetBinContent(k,j); Double_t A = fHgA->GetBinContent(l,k,j); fHgA->SetBinContent(l, k, j, A + e); } // go to next k continue; } for (Int_t j=k+1; j<=nbinssig; j++) { Double_t e = fHge->GetBinContent(k,j); Double_t A = fHgA->GetBinContent(l,k,j); fHgA->SetBinContent(l, k, j, A + r*e/ne); } if (na <= eps) { // go to next k continue; } for (Int_t i=1; iGetBinContent(i,k); Double_t A = fHgA->GetBinContent(l,i,k); fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na); } for (Int_t i=1; i<=nbinssig; i++) { if (i == k) continue; for (Int_t j=i+1; j<=nbinssig; j++) { if (j == k) continue; Double_t a = fHga->GetBinContent(i,k); Double_t e = fHge->GetBinContent(k,j); Double_t A = fHgA->GetBinContent(l,i,j); fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne) ); } } } 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& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon) { *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)"); //-------------------------- fHSigmaThetaMC = &sigthmc; fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); fHSigmaThetaON = &sigthon; fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (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(); if (suma != 0.0) hista->Scale(1./suma); TH1D *histap = new TH1D( (const TH1D&)*hista ); TH1D *histb = sigthmc.ProjectionY("sigmc_y", l, l, ""); Stat_t sumb = histb->Integral(); if (sumb != 0.0) 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) // (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; // 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; //------------------------------------------------------------------ } //............ 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); 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; } */ fHSigmaThetaMC = new TH2D; fHSigmaThetaMC->Read("2D-ThetaSigmabarMC"); *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; } */ fHSigmaThetaON = new TH2D; fHSigmaThetaON->Read("2D-ThetaSigmabarON"); *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; } */ fHSigmaThetaOFF = new TH2D; fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF"); *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; } */ fHDiffPixTheta = new TH3D; fHDiffPixTheta->Read("3D-ThetaPixDiff"); *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; } */ fHDiffPixThetaMC = new TH3D; fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC"); *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; } */ fHDiffPixThetaON = new TH3D; fHDiffPixThetaON->Read("3D-ThetaPixDiffON"); *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; } */ fHDiffPixThetaOFF = new TH3D; fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF"); *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; } */ fHSigmaPixTheta = new TH3D; fHSigmaPixTheta->Read("3D-ThetaPixSigma"); *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; } */ fHSigmaPixThetaMC = new TH3D; fHSigmaPixThetaMC->Read("3D-ThetaPixSigmaMC"); *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; } */ fHSigmaPixThetaON = new TH3D; fHSigmaPixThetaON->Read("3D-ThetaPixSigmaON"); *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; } */ fHSigmaPixThetaOFF = new TH3D; fHSigmaPixThetaOFF->Read("3D-ThetaPixSigmaOFF"); *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; } */ fHgMC = new TH3D; fHgMC->Read("3D-PaddingMatrixMC"); *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; } */ fHgON = new TH3D; fHgON->Read("3D-PaddingMatrixON"); *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; } */ fHgOFF = new TH3D; fHgOFF->Read("3D-PaddingMatrixOFF"); *fLog << all << "MPad : Object '3D-PaddingMatrixOFF' 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"); 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 || !fHgMC || !fHgON || !fHgOFF ) { *fLog << err << "MPad : At least one of the histograms needed for the padding is not defined ... aborting." << endl; return kFALSE; } fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); if (!fPointPos) { *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." << endl; return kFALSE; } fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); if (!fPed) { *fLog << err << AddSerialNumber(fNamePedPhotCam) << "[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; } fBad = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); if (!fBad) { *fLog << inf << "MPad : MBadPixelsCam container not present... continue." << endl; } 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, 75.0, 100, 0.0, 75.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, -100.0, 300.0, 100, -100, 300); 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, -100.0, 200.0); fHNSB->SetXTitle("no.of added NSB photons"); fHNSB->SetYTitle("no.of pixels"); //-------------------------------------------------------------------- fIter = 10; 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 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.92 * 0.94 * 0.23 * 0.9; //------------------------------------------------ // Add pixels to MCerPhotEvt which are not yet in; // set their number of photons equal to zero. // Purpose : by the padding the number of photons is changed // so that cleaning and cuts may be affected. // However, no big effect is expectec unless the padding is strong // (strong increase of the noise level) // At present comment out this code /* 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 // //*fLog << all << "MPad::Process(); before padding : " << endl; // inner sigmabar/sqrt(area) Double_t ratioinn = fCam->GetPixRatio(0); Double_t sigbarInnerold_phot = (fPed->GetArea(0)).GetRms() * sqrt(ratioinn); Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton; Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold; // outer sigmabar/sqrt(area) Double_t ratioout = fCam->GetPixRatio(500); Double_t sigbarOuterold_phot = (fPed->GetArea(1)).GetRms() * sqrt(ratioout); Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton; Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold; const Double_t theta = fPointPos->GetZd(); //*fLog << all << "theta = " << theta << endl; Int_t binTheta = fHSigmaThetaON->GetXaxis()->FindBin(theta); if ( binTheta < 1 || binTheta > fHSigmaThetaON->GetNbinsX() ) { *fLog << warn << "MPad::Process(); binNumber out of range : theta, binTheta = " << theta << ", " << binTheta << endl; rc = 2; fErrors[rc]++; return kCONTINUE; } //*fLog << all << "binTheta = " << binTheta << endl; //------------------------------------------- // get number of events in this theta bin // and number of events in this sigbarInnerold_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, ""); //*fLog << "Title of histogram : " << st->GetTitle() << endl; //for (Int_t i=1; i<=hn->GetNbinsX(); i++) //{ // *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i) // << endl; //} nTheta = hn->Integral(); Int_t binSigma = hn->FindBin(sigbarInnerold_phot); nSigma = hn->GetBinContent(binSigma); //*fLog << all // << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " // << theta << ", " << sigbarInnerold_phot << ", " << binTheta // << ", " // << binSigma << ", " << nTheta << ", " << nSigma << endl; delete hn; //------------------------------------------- //****************************************************************** // has the event to be padded ? // if yes : to which sigmabar should it be padded ? // Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot); //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", " // << sigbarInnerold_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 sigmabarInner_phot = 0; Double_t sigmabarInner = 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; //} sigmabarInner_phot = hpad->GetRandom(); sigmabarInner = sigmabarInner_phot * fPEperPhoton; //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl; delete hpad; // new inner sigmabar2/area const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner; //------------------------------------------- //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " // << sigbarInnerold << ", "<< sigmabarInner << endl; // Stop if target sigmabar is <= sigbarold if (sigmabarInner <= sigbarInnerold) { *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : " << sigmabarInner << ", " << sigbarInnerold << ", aborting" << 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; } // In this Theta bin, get RMS of (Sigma^2-sigmabar^2) for inner pixels // The average electronic noise (to be added) has to be in the order of // this RMS, otherwise the electronic noise of an individual pixel // (elNoise2Pix) may become negative // Attention : maximum ID of inner pixels hard coded !!! Int_t idmaxpixinner = 396; Int_t binmaxpixinner = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)idmaxpixinner ); TH1D *hnoise; if (fType == "MC") hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "ON") hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "OFF") hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); 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(2.0*RMS, sigmabarInner2 - sigbarInnerold2); //*fLog << all << "elNoise2 = " << elNoise2 << endl; lambdabar = (sigmabarInner2 - sigbarInnerold2 - 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 and not bad (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; // throw the Sigma for the pixels 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; } Double_t sigbarold2; Double_t sigmabar2; Double_t sigmabarOuter2; for (UInt_t i=0; iGetPixRatio(j); if (ratioArea < 1.5) { sigbarold2 = sigbarInnerold2; sigmabar2 = sigmabarInner2; } else { sigbarold2 = sigbarOuterold2; //sigbarOuter2 = sigmabarInner2 * sigbarOuterold2 / sigbarInnerold2; sigmabarOuter2 = sigmabarInner2 + (sigbarOuterold2 - sigbarInnerold2); sigmabar2 = sigmabarOuter2; } MPedPhotPix &ppix = (*fPed)[j]; if ( fBad != NULL && (*fBad)[j].IsUnsuitable() ) { // this should never occur, because bad pixels should have // been set unused *fLog << all << "MPad::Process; bad pixel found which is used, j = " << j << "... go to next pixel." << endl; continue; } // count number of pixels treated fWarnings[0]++; Double_t oldsigma_phot = ppix.GetRms(); Double_t oldsigma = oldsigma_phot * fPEperPhoton; Double_t oldsigma2 = oldsigma*oldsigma / ratioArea; //--------------------------------- // throw the Sigma for this pixel // Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); Int_t count; Bool_t ok; TH1D *hdiff; hdiff = sp->ProjectionZ("", binTheta, binTheta, binPixel, binPixel, ""); Double_t integral = hdiff->Integral(); // if there are no entries in hdiff, diff cannot be thrown // in this case diff will be set to zero if ( integral == 0 ) { //*fLog << warn << "MPad::Process(); fType = " << fType // << ", projection of '" // << sp->GetName() << "' for Theta bin " // << binTheta << " and pixel " << j // << " has no entries; set diff equal to the old difference " // << endl; diff = TMath::Max(oldsigma2 - sigbarold2, lambdabar*F2excess + oldsigma2 - sigmabar2); rc = 2; fWarnings[rc]++; } // start of else ------------------- else { count = 0; ok = kFALSE; for (Int_t m=0; mGetRandom(); //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : " // << m << ", " << diff_phot << endl; diff = diff_phot * fPEperPhoton * fPEperPhoton; // the following condition ensures that elNoise2Pix > 0.0 if ( (diff + sigmabar2 - oldsigma2 - lambdabar*F2excess) > 0.0 ) { ok = kTRUE; break; } } if (!ok) { //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratioArea, lambdabar = " // << theta << ", " // << j << ", " << count << ", " << sigmabar2 << ", " // << diff << ", " << oldsigma2 << ", " << ratioArea << ", " // << lambdabar << endl; diff = lambdabar*F2excess + oldsigma2 - sigmabar2; rw = 1; fWarnings[rw]++; } } // end of else -------------------- delete hdiff; sigma2 = diff + sigmabar2; //--------------------------------- // get the additional sigma^2 for this pixel (due to the padding) addSig2 = (sigma2 - oldsigma2) * ratioArea; 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( sigma2 * ratioArea ); 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 //*fLog << all << "MPad::Process(); after padding : " << endl; //*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(' '); if (fWarnings[0] == 0 ) fWarnings[0] = 1; *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) << (int)(fWarnings[1]*100/fWarnings[0]) << "%) Pixel: iteration to find acceptable sigma failed" << ", fIter = " << fIter << endl; *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3) << (int)(fWarnings[2]*100/fWarnings[0]) << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << 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[4] << " (" << setw(3) << (int)(fErrors[4]*100/GetNumExecutions()) << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << 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, 600)); c.Divide(3, 2); 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); */ //-------------------------------------------------------------------- // 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(4); TH2D *m1; m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); m1->SetDirectory(NULL); m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); m1->SetXTitle("Sigmabar-old"); m1->SetYTitle("Sigmabar-new"); m1->DrawCopy("box"); m1->SetBit(kCanDelete);; c.cd(5); TH2D *m2; m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); m2->SetDirectory(NULL); m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); m2->SetXTitle("Sigmabar-old"); m2->SetYTitle("Sigmabar-new"); m2->DrawCopy("box"); m2->SetBit(kCanDelete);; c.cd(6); TH2D *m3; m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); m3->SetDirectory(NULL); m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); m3->SetXTitle("Sigmabar-old"); m3->SetYTitle("Sigmabar-new"); m3->DrawCopy("box"); m3->SetBit(kCanDelete);; return kTRUE; } // --------------------------------------------------------------------------