Index: trunk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 5432) +++ trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 5432) @@ -0,0 +1,2658 @@ +/* ======================================================================== *\ +! +! * +! * 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, 10/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// 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; + +// -------------------------------------------------------------------------- +// +// Constructor. +// +MPad::MPad(const char *name, const char *title) +{ + fName = name ? name : "MPad"; + fTitle = title ? title : "Task for the padding"; + + fType = ""; + + fHSigmaTheta = NULL; + fHSigmaThetaOuter = NULL; + + fHSigmaThetaMC = NULL; + fHSigmaThetaOuterMC = NULL; + fHDiffPixThetaMC = NULL; + fHDiffPixThetaTargetMC = NULL; + + fHSigmaThetaON = NULL; + fHSigmaThetaOuterON = NULL; + fHDiffPixThetaON = NULL; + fHDiffPixThetaTargetON = NULL; + + fHSigmaThetaOFF = NULL; + fHSigmaThetaOuterOFF = NULL; + fHDiffPixThetaOFF = NULL; + fHDiffPixThetaTargetOFF = NULL; + + fHgMC = NULL; + fHgON = NULL; + fHgOFF = NULL; + + fHgOuterMC = NULL; + fHgOuterON = NULL; + fHgOuterOFF = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; + + fNamePedPhotCam = "MPedPhotCamFromData"; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPad::~MPad() +{ + delete fHSigmaTheta; + delete fHSigmaThetaOuter; + + if (fHSigmaThetaMC) + { + delete fHSigmaThetaMC; + delete fHSigmaThetaOuterMC; + delete fHDiffPixThetaMC; + delete fHDiffPixThetaTargetMC; + delete fHgMC; + delete fHgOuterMC; + } + + if (fHSigmaThetaON) + { + delete fHSigmaThetaON; + delete fHSigmaThetaOuterON; + delete fHDiffPixThetaON; + delete fHDiffPixThetaTargetON; + delete fHgON; + delete fHgOuterON; + } + + if (fHSigmaThetaOFF) + { + delete fHSigmaThetaOFF; + delete fHSigmaThetaOuterOFF; + delete fHDiffPixThetaOFF; + delete fHDiffPixThetaTargetOFF; + delete fHgOFF; + delete fHgOuterOFF; + } + + delete fHSigmaPedestal; + delete fHPhotons; + delete fHNSB; + + delete fInfile; +} + +// -------------------------------------------------------------------------- +// +// SetNamePedPhotCam +// +// set the name of the MPedPhotCam container + +void MPad::SetNamePedPhotCam(const char *name) +{ + fNamePedPhotCam = name; +} + +// -------------------------------------------------------------------------- +// +// 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; + + return; +} + +// -------------------------------------------------------------------------- +// +// Read pad histograms from a file +// these are the histograms generated by MMakePadHistograms +// +Bool_t MPad::ReadPadHistograms(TString type, TString namepad) +{ + + //------------------------------------ + Int_t OK = 0; + + if (type == "ON") + { + *fLog << inf << "MPad : Read pad histograms for " << type + << " data from file " << namepad << endl; + + TH2D *sigth = new TH2D; + TH2D *sigthout = new TH2D; + TH3D *diff = new TH3D; + + TFile *fInfile = new TFile(namepad); + fInfile->ls(); + + OK = sigth->Read("2D-ThetaSigmabar(Inner)"); + fHSigmaThetaON = new TH2D( (const TH2D&)(*sigth) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", + "2D-ThetaSigmabarON (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' was read in" + << endl; + + + OK = sigthout->Read("2D-ThetaSigmabar(Outer)"); + fHSigmaThetaOuterON = new TH2D( (const TH2D&)(*sigthout) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaOuterON->SetNameTitle("2D-ThetaSigmabarOuterON", + "2D-ThetaSigmabarOuterON (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' was read in" + << endl; + + + OK = diff->Read("3D-ThetaPixDiff"); + fHDiffPixThetaON = new TH3D( (const TH3D&)(*diff) ); + if (!OK) + { + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' not found" + << endl; + return kFALSE; + } + + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", + "3D-ThetaPixDiffON (orig.)"); + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' was read in" << endl; + + + return kTRUE; + } + + //------------------------------------ + if (type == "OFF") + { + *fLog << inf << "MPad : Read pad histograms for " << type + << " data from file " << namepad << endl; + + + TH2D *sigth = new TH2D; + TH2D *sigthout = new TH2D; + TH3D *diff = new TH3D; + + TFile *fInfile = new TFile(namepad); + fInfile->ls(); + + OK = sigth->Read("2D-ThetaSigmabar(Inner)"); + fHSigmaThetaOFF = new TH2D( (const TH2D&)(*sigth) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", + "2D-ThetaSigmabarOFF (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' was read in" + << endl; + + + OK = sigthout->Read("2D-ThetaSigmabar(Outer)"); + fHSigmaThetaOuterOFF = new TH2D( (const TH2D&)(*sigthout) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaOuterOFF->SetNameTitle("2D-ThetaSigmabarOuterOFF", + "2D-ThetaSigmabarOuterOFF (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' was read in" + << endl; + + + OK = diff->Read("3D-ThetaPixDiff"); + fHDiffPixThetaOFF = new TH3D( (const TH3D&)(*diff) ); + if (!OK) + { + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' not found" + << endl; + return kFALSE; + } + + fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", + "3D-ThetaPixDiffOFF (orig.)"); + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' was read in" << endl; + + + + return kTRUE; + } + + //------------------------------------ + if (type == "MC") + { + + TH2D *sigth = new TH2D; + TH2D *sigthout = new TH2D; + TH3D *diff = new TH3D; + + *fLog << inf << "MPad : Read pad histograms for " << type + << " data from file " << namepad << endl; + + TFile *fInfile = new TFile(namepad); + fInfile->ls(); + + OK = sigth->Read("2D-ThetaSigmabar(Inner)"); + fHSigmaThetaMC = new TH2D( (const TH2D&)(*sigth) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", + "2D-ThetaSigmabarMC (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Inner)' was read in" + << endl; + + OK = sigthout->Read("2D-ThetaSigmabar(Outer)"); + fHSigmaThetaOuterMC = new TH2D( (const TH2D&)(*sigthout) ); + if (!OK) + { + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' not found" + << endl; + return kFALSE; + } + + fHSigmaThetaOuterMC->SetNameTitle("2D-ThetaSigmabarOuterMC", + "2D-ThetaSigmabarOuterMC (orig.)"); + *fLog << all + << " MPad : Object '2D-ThetaSigmabar(Outer)' was read in" + << endl; + + + + OK = diff->Read("3D-ThetaPixDiff"); + fHDiffPixThetaMC = new TH3D( (const TH3D&)(*diff) ); + if (!OK) + { + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' not found" + << endl; + return kFALSE; + } + + fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", + "3D-ThetaPixDiffMC (orig.)"); + + *fLog << all + << " MPad : Object '3D-ThetaPixDiff' was read in" << endl; + + return kTRUE; + } + + 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(TString nameon, TString nameoff, + TString namemc, TString fileout) +{ + // read the histograms produced by MMakePadHistograms + TH2D *hist2patternin = NULL; + TH2D *hist2patternout = NULL; + TH3D *hist3pattern = NULL; + + Bool_t OK = 0; + + if (namemc != "") + { + OK = ReadPadHistograms("MC", namemc); + if (!OK) + { + *fLog << err << "MPad::MergeONOFFMC : Pad histograms for MC not available" + << endl; + return kFALSE; + } + + hist2patternin = fHSigmaThetaMC; + hist2patternout = fHSigmaThetaOuterMC; + hist3pattern = fHDiffPixThetaMC; + *fLog << inf << "" << endl; + *fLog << inf << "MPad::MergeONOFFMC : Pad histograms for MC data read from file " + << namemc << endl; + } + + if (nameoff != "") + { + OK = ReadPadHistograms("OFF", nameoff); + if (!OK) + { + *fLog << err << "MPad::MergeONOFFMC : Pad histograms for OFF not available" + << endl; + return kFALSE; + } + + hist2patternin = fHSigmaThetaOFF; + hist2patternout = fHSigmaThetaOuterOFF; + hist3pattern = fHDiffPixThetaOFF; + *fLog << inf << "" << endl; + *fLog << inf << "MPad::MergeONOFFMC : Pad histograms for OFF data read from file " + << nameoff << endl; + } + + if (nameon != "") + { + OK = ReadPadHistograms("ON", nameon); + if (!OK) + { + *fLog << err << "MPad::MergeONOFFMC : Pad histograms for ON not available" + << endl; + return kFALSE; + } + + hist2patternin = fHSigmaThetaON; + hist2patternout = fHSigmaThetaOuterON; + hist3pattern = fHDiffPixThetaON; + *fLog << inf << "" << endl; + *fLog << inf << "MPad::MergeONOFFMC : Pad histograms for ON data read from file " + << nameon << endl; + } + + *fLog << "hist2patternin = " << hist2patternin << endl; + *fLog << "hist2patternout = " << hist2patternout << endl; + *fLog << "hist3pattern = " << hist3pattern << endl; + + + TAxis *axs = hist2patternin->GetXaxis(); + Int_t nbinsthetas = axs->GetNbins(); + *fLog << "vor fHSigmaTheta : nbinsthetas = " << nbinsthetas << endl; + + TArrayD edgess; + edgess.Set(nbinsthetas+1); + for (Int_t i=0; i<=nbinsthetas; i++) + { + edgess[i] = axs->GetBinLowEdge(i+1); + *fLog << all << "i, theta low edge = " << i << ", " << edgess[i] + << endl; + } + + + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + // 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; + hist2patternin->Copy(*fHSigmaTheta); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + //------------------------- + fHSigmaThetaOuter = new TH2D; + hist2patternout->Copy(*fHSigmaThetaOuter); + fHSigmaThetaOuter->SetNameTitle("2D-ThetaSigmabarOuter", "2D-ThetaSigmabarOuter (target)"); + + //-------------------------- + fHDiffPixThetaTargetMC = new TH3D; + hist3pattern->Copy(*fHDiffPixThetaTargetMC); + fHDiffPixThetaTargetMC->SetNameTitle("3D-ThetaPixDiffTargetMC", "3D-ThetaPixDiffMC (target)"); + + //------------------------- + fHDiffPixThetaTargetON = new TH3D; + hist3pattern->Copy(*fHDiffPixThetaTargetON); + fHDiffPixThetaTargetON->SetNameTitle("3D-ThetaPixDiffTargetON", "3D-ThetaPixDiffON (target)"); + + //------------------------- + fHDiffPixThetaTargetOFF = new TH3D; + hist3pattern->Copy(*fHDiffPixThetaTargetOFF); + fHDiffPixThetaTargetOFF->SetNameTitle("3D-ThetaPixDiffTargetOFF", "3D-ThetaPixDiffOFF (target)"); + + *fLog << all << "MPad::MergeONOFFMC(); Histograms for the merged padding plots were booked" + << endl; + + + //-------------------------- + + // get binning for fHgON, fHgOFF and fHgMC from fHSigmaTheta + // binning in Theta + TAxis *ax = hist2patternin->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + + //*fLog << "nbinstheta = " << nbinstheta << endl; + + 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 = hist2patternin->GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + + //*fLog << "nbinssig = " << nbinssig << endl; + + 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; + + //----------------------- + // get binning for fHgOuterON, fHgOuterOFF and fHgOuterMC + // from fHSigmaThetaOuter + // binning in sigmabar + TAxis *ayout = hist2patternout->GetYaxis(); + Int_t nbinssigout = ayout->GetNbins(); + TArrayD edgesyout; + edgesyout.Set(nbinssigout+1); + for (Int_t i=0; i<=nbinssigout; i++) + { + edgesyout[i] = ayout->GetBinLowEdge(i+1); + //*fLog << all << "i, sigmabar low edge = " << i << ", " << edgesyout[i] + // << endl; + } + MBinning binsgout; + binsgout.SetEdges(edgesyout); + + + fHgOuterMC = new TH3D; + MH::SetBinning(fHgOuterMC, &binth, &binsgout, &binsgout); + fHgOuterMC->SetNameTitle("3D-PaddingMatrixOuterMC", "3D-PadMatrixThetaOuterMC"); + + fHgOuterON = new TH3D; + MH::SetBinning(fHgOuterON, &binth, &binsgout, &binsgout); + fHgOuterON->SetNameTitle("3D-PaddingMatrixOuterON", "3D-PadMatrixThetaOuterON"); + + fHgOuterOFF = new TH3D; + MH::SetBinning(fHgOuterOFF, &binth, &binsgout, &binsgout); + fHgOuterOFF->SetNameTitle("3D-PaddingMatrixOuterOFF", "3D-PadMatrixThetaOuterOFF"); + + *fLog << all << "fHgOuter histograms were booked" << endl; + + + //-------------------------------------------------------------------- + // define the first (A), second (B) and third (C) data set + // (A, B, C) may be (ON, OFF, MC) + // or (ON, OFF, "") + // or (ON, MC, "") + // or (OFF,MC, "") + // or the same sets in any other order + + TString tagA = ""; + TString tagB = ""; + TString tagC = ""; + + TH2D *histA = NULL; + TH2D *histB = NULL; + TH2D *histC = NULL; + + TH2D *histOuterA = NULL; + TH2D *histOuterB = NULL; + TH2D *histOuterC = NULL; + + TH3D *histDiffA = NULL; + TH3D *histDiffB = NULL; + TH3D *histDiffC = NULL; + + TH3D *fHgA = NULL; + TH3D *fHgB = NULL; + TH3D *fHgC = NULL; + + TH3D *fHgOuterA = NULL; + TH3D *fHgOuterB = NULL; + TH3D *fHgOuterC = NULL; + + if (fHSigmaThetaON) + { + tagA = "ON"; + histA = fHSigmaThetaON; + histOuterA = fHSigmaThetaOuterON; + histDiffA = fHDiffPixThetaON; + fHgA = fHgON; + fHgOuterA = fHgOuterON; + + if (fHSigmaThetaOFF) + { + tagB = "OFF"; + histB = fHSigmaThetaOFF; + histOuterB = fHSigmaThetaOuterOFF; + histDiffB = fHDiffPixThetaOFF; + fHgB = fHgOFF; + fHgOuterB = fHgOuterOFF; + + if (fHSigmaThetaMC) + { + tagC = "MC"; + histC = fHSigmaThetaMC; + histOuterC = fHSigmaThetaOuterMC; + histDiffC = fHDiffPixThetaMC; + fHgC = fHgMC; + fHgOuterC = fHgOuterMC; + } + } + + else if (fHSigmaThetaMC) + { + tagB = "MC"; + histB = fHSigmaThetaMC; + histOuterB = fHSigmaThetaOuterMC; + histDiffB = fHDiffPixThetaMC; + fHgB = fHgMC; + fHgOuterB = fHgOuterMC; + } + } + + else if (fHSigmaThetaOFF) + { + tagA = "OFF"; + histA = fHSigmaThetaOFF; + histOuterA = fHSigmaThetaOuterOFF; + histDiffA = fHDiffPixThetaOFF; + fHgA = fHgOFF; + fHgOuterA = fHgOuterOFF; + + if (fHSigmaThetaMC) + { + tagB = "MC"; + histB = fHSigmaThetaMC; + histOuterB = fHSigmaThetaOuterMC; + histDiffB = fHDiffPixThetaMC; + fHgB = fHgMC; + fHgOuterB = fHgOuterMC; + } + + else + { + *fLog << err << "MPad::MergeONOFFMC; there are no data sets to be merged : tagA, tagB, tagC = " + << tagA << ", " << tagB << ", " << tagC << endl; + return kFALSE; + } + } + + else + { + *fLog << err << "MPad::MergeONOFFMC; there are no data sets to be merged" + << endl; + return kFALSE; + } + + *fLog << inf << "MPad::MergeONOFFMC; the following data sets will be merged : " + << tagA << " " << tagB << " " << tagC << endl; + + + //------------------------------------------- + // merge the 2D histograms (theta, Sigmabar) + // from the data sets A, B and C + // for the inner pixels and for the outer pixels + + TString canv("Inner"); + MergeABC(tagA, tagB, tagC, histA, histB, histC, + fHSigmaTheta, fHgA, fHgB, fHgC, canv); + + TString canvout("Outer"); + MergeABC(tagA, tagB, tagC, histOuterA, histOuterB, histOuterC, + fHSigmaThetaOuter, fHgOuterA, fHgOuterB, fHgOuterC, canvout); + + + //------------------------------------------- + // Define the target distributions fHDiffPixThetaTargetON, ..OFF, ..MC + // (fHDiffPixThetaTarget.. will be used in the padding) + + + //............ start of new loop over Theta bins .................... + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of A/B/C events to be padded : fracA, fracB, fracC + Double_t fracA = 0.0; + Double_t fracB = 0.0; + Double_t fracC = 0.0; + + fracA = fHgA->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + fracB = fHgB->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + if (tagC != "") fracC = fHgC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + *fLog << all << "Theta bin j : tagA, fracA, tagB, fracB, tagC, fracC = " + << j << ". " << tagA << ": " << fracA << "; " << tagB << ": " + << fracB << "; " << tagC << ": " << fracC << endl; + + + //------------------------------------------------------------------ + // for the current theta bin : + // define target distribution 'sigma^2-sigmabar^2 vs. pixel number' + // take the distribution from that sample which has the lowest 'frac.' + // At present, for all samples (ON, OFF, MC) the same target distribution + // is used + + Double_t fracmin = fracA; + TString tag = tagA; + + if (tagB != "") + { + if (fracB < fracmin) + { + fracmin = fracB; + tag = tagB; + } + if (tagC != "") + { + if (fracC < fracmin) + { + fracmin = fracC; + tag = tagC; + } + } + } + + *fLog << all << "MPad::MergeONOFFMC; thetabin " << j + << " : sample with smallest frac = " << tag << endl; + + TH3D *hist3 = NULL; + if (tag == "ON") + hist3 = fHDiffPixThetaON; + else if (tag == "OFF") + hist3 = fHDiffPixThetaOFF; + else if (tag == "MC") + hist3 = fHDiffPixThetaMC; + + ay = hist3->GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + TAxis *az = hist3->GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + for (Int_t k=1; k<=nbinspixel; k++) + { + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist3->GetBinContent(j, k,l); + fHDiffPixThetaTargetON-> SetBinContent(j, k, l, cont); + fHDiffPixThetaTargetOFF->SetBinContent(j, k, l, cont); + fHDiffPixThetaTargetMC-> SetBinContent(j, k, l, cont); + } + } + } + //............ end of new loop over Theta bins .................... + + *fLog << all + << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaTheta->SetDirectory(NULL); + fHSigmaThetaOuter->SetDirectory(NULL); + + if (tagA == "MC" || tagB == "MC" || tagC == "MC") + { + fHSigmaThetaMC->SetDirectory(NULL); + fHSigmaThetaOuterMC->SetDirectory(NULL); + fHDiffPixThetaMC->SetDirectory(NULL); + fHgMC->SetDirectory(NULL); + fHgOuterMC->SetDirectory(NULL); + } + + if (tagA == "ON" || tagB == "ON" || tagC == "ON") + { + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOuterON->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHgON->SetDirectory(NULL); + fHgOuterON->SetDirectory(NULL); + } + + if (tagA == "OFF" || tagB == "OFF" || tagC == "OFF") + { + fHSigmaThetaOFF->SetDirectory(NULL); + fHSigmaThetaOuterOFF->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + fHgOuterOFF->SetDirectory(NULL); + } + + // write the target padding histograms onto a file --------- + WritePaddingDist(fileout); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Merge the 2D-histograms (Theta, sigmabar) +// hA, hB and hC +// +// - put the merged distribution into hM +// - and define the matrices hgA, hgB and hgC +// +// These matrices tell which fraction of events should be padded +// from bin k of sigmabar to bin j +// + +Bool_t MPad::MergeABC(TString tA, TString tB, TString tC, + TH2D *hA, TH2D *hB, TH2D *hC, TH2D *hM, + TH3D *hgA, TH3D *hgB, TH3D *hgC, TString canv) +{ + *fLog << all << "MPad::MergeABC(); Entry" << endl; + + TAxis *ax = hM->GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + + TAxis *ay = hM->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); + + + //............ loop over Theta bins ........................... + + //*fLog << all << "MPad::MergeABC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + //------------------------------------------- + // merge A (ON) and B (OFF) distributions + // input : hista is the normalized 1D distr. of sigmabar for A (ON) data + // histb is the normalized 1D distr. of sigmabar for B (OFF) data + // output : histap will be the merged distribution (AB) + // fHga(k,j) will tell which fraction of A (ON) events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of B (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 = hA->ProjectionY("sigon_y", l, l, ""); + TString tita(canv); + tita += "-A (orig)-"; + tita += l; + hista->SetNameTitle(tita, tita); + Stat_t suma = hista->Integral(); + if (suma != 0.0) hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + TString titab(canv); + titab += "-AB (merged)-"; + titab += l; + histap->SetNameTitle(titab, titab); + + TH1D *histb = hB->ProjectionY("sigoff_y", l, l, ""); + TString titb(canv); + titb += "-B (orig)-"; + titb += l; + histb->SetNameTitle(titb, titb); + 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(A=" << tA << ", B=" + << tB << ") for theta bin l = " + << l << "; NevA, NevB = " << suma << ", " << sumb << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(A=" << tA << ", B=" + // << tB << ") for theta bin l = " + // << l << endl; + + TString canvab(canv); + canvab += "AB-"; + canvab += l; + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig, canvab); + + // fill fHgA and fHgB + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent(k,j); + hgA->SetBinContent(l, k, j, a); + + Double_t b = fHgb->GetBinContent(k,j); + hgB->SetBinContent(l, k, j, b); + } + + + //------------------------------------------------- + // if there is no further distribution to be merged + // fill target distribution fHSigmaTheta + // + if (!hC) + { + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = histap->GetBinContent(j); + hM->SetBinContent(l, j, a); + + //*fLog << "MPad::MergeABC; l, j, a = " << l << ", " << j << ", " + // << a << endl; + } + + delete fHga; + delete fHgb; + delete histap; + } + + //-------------- next merge (start) ---------------------------------- + if (hC) + { + + // now merge AB (ON-OFF) and C (MC) distributions + // input : histe is the result of the merging of A (ON) and B (OFF) + // distributions + // histf is the normalized 1D distr. of sigmabar for C (MC) data + // output : histep will be the merged distribution (target distribution) + // fHge(k,j) will tell which fraction of A (ON), B (OFF) events + // should be padded from bin k to bin j + // fHgf(k,j) will tell which fraction of C (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); + TString titabm(canv); + titabm += "-AB (merged)-"; + titabm += l; + histe->SetNameTitle(titabm, titabm); + + TH1D *histep = new TH1D( (const TH1D&)*histe); + TString titabcm(canv); + titabcm +="-ABC (merged)-"; + titabcm += l; + histep->SetNameTitle(titabcm, titabcm); + + TH1D *histf = hC->ProjectionY("sigmc_y", l, l, ""); + TString titc(canv); + titc += "-C (orig)-"; + titc += l; + histf->SetNameTitle(titc, titc); + 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(AB=" << tA << "-" << tB + << ", C=" << tC << ") for theta bin l = " + << l << "; NevC = " << sumf << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(AB=" << tA << "-" << tB + // << ", C=" << tC << ") for theta bin l = " + // << l << endl; + + TString canvabc(canv); + canvabc += "ABC-"; + canvabc += l; + Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig, canvabc); + + // update fHgA and fHgB + UpdateHg(fHga, histap, fHge, hgA, nbinssig, l); + UpdateHg(fHgb, histap, fHge, hgB, nbinssig, l); + + // fill fHgC + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t f = fHgf->GetBinContent(k,j); + hgC->SetBinContent(l, k, j, f); + } + + // fill target distribution fHSigmaTheta + // + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t ep = histep->GetBinContent(j); + hM->SetBinContent(l, j, ep); + + //*fLog << "MPad::MergeABC; l, j, ep = " << l << ", " << j << ", " + // << ep << endl; + } + + + //------------------------------------------- + + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + delete histe; + delete histf; + delete fHge; + delete fHgf; + delete histep; + + } + //-------------- next merge (end) ---------------------------------- + + } + //............ end of loop over Theta bins .................... + + *fLog << all << "MPad::MergeABC(); Return" << endl; + + return kTRUE; +} +// -------------------------------------------------------------------------- +// +// 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, + TString canv ) +{ + *fLog << all << "MPad::Merge2Distributions(); Entry" << endl; + + + // 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(canv, canv, 600, 600); + gROOT->SetSelectedPad(NULL); + + pad->Divide(2, 2); + + pad->cd(1); + gPad->SetBorderMode(0); + hista->SetDirectory(NULL); + hista->UseCurrentStyle(); + hista->DrawClone(); + hista->SetBit(kCanDelete); + + pad->cd(2); + gPad->SetBorderMode(0); + histb->SetDirectory(NULL); + histb->UseCurrentStyle(); + histb->DrawClone(); + histb->SetBit(kCanDelete); + + pad->cd(3); + gPad->SetBorderMode(0); + histap->SetDirectory(NULL); + histap->UseCurrentStyle(); + histap->DrawClone(); + histap->SetBit(kCanDelete); + + pad->Draw(); + + //-------------------------------------------------------------------- + + delete histc; + delete histbp; + delete histcp; + + *fLog << all << "MPad::Merge2Distributions(); Return" << endl; + + 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; +} + +// -------------------------------------------------------------------------- +// +// Write padding distributions onto a file +// these are the histograms produced in the member function +// MergeONOFFMC +// plus the distributions produced by MMakePadHistograms +// +Bool_t MPad::WritePaddingDist(TString namefileout) +{ + *fLog << all << "MPad::WritePaddingDist(); Padding distributions for "; + + TFile outfile(namefileout, "RECREATE"); + + fHSigmaTheta->Write(); + fHSigmaThetaOuter->Write(); + + if (fHSigmaThetaMC) + { + *fLog << all << " MC "; + fHSigmaThetaMC->Write(); + fHSigmaThetaOuterMC->Write(); + fHDiffPixThetaMC->Write(); + fHDiffPixThetaTargetMC->Write(); + fHgMC->Write(); + fHgOuterMC->Write(); + } + + if (fHSigmaThetaON) + { + *fLog << all << " ON "; + fHSigmaThetaON->Write(); + fHSigmaThetaOuterON->Write(); + fHDiffPixThetaON->Write(); + fHDiffPixThetaTargetON->Write(); + fHgON->Write(); + fHgOuterON->Write(); + } + + if (fHSigmaThetaOFF) + { + *fLog << all << " OFF "; + fHSigmaThetaOFF->Write(); + fHSigmaThetaOuterOFF->Write(); + fHDiffPixThetaOFF->Write(); + fHDiffPixThetaTargetOFF->Write(); + fHgOFF->Write(); + fHgOuterOFF->Write(); + } + + *fLog << all << " data were written onto file " + << namefileout << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read padding distributions from a file +// thes are the distributions which were written out by +// the member function WritePaddingDist +// +Bool_t MPad::ReadPaddingDist(TString namefilein) +{ + *fLog << all << "MPad : Read padding histograms from file " + << namefilein << endl; + + Int_t OK = 0; + + fInfile = new TFile(namefilein); + fInfile->ls(); + + //------------------------------------ + + fHSigmaTheta = new TH2D; + OK = fHSigmaTheta->Read("2D-ThetaSigmabar"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabar' was read in" << endl; + } + + fHSigmaThetaOuter = new TH2D; + OK =fHSigmaThetaOuter->Read("2D-ThetaSigmabarOuter"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOuter' was read in" << endl; + } + + + //------------------------------------ + + fHSigmaThetaMC = new TH2D; + OK = fHSigmaThetaMC->Read("2D-ThetaSigmabarMC"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl; + } + + fHSigmaThetaOuterMC = new TH2D; + OK = fHSigmaThetaOuterMC->Read("2D-ThetaSigmabarOuterMC"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOuterMC' was read in" << endl; + } + + fHDiffPixThetaMC = new TH3D; + OK = fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl; + } + + fHDiffPixThetaTargetMC = new TH3D; + OK = fHDiffPixThetaTargetMC->Read("3D-ThetaPixDiffTargetMC"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffTargetMC' was read in" << endl; + } + + fHgMC = new TH3D; + OK = fHgMC->Read("3D-PaddingMatrixMC"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl; + } + + fHgOuterMC = new TH3D; + OK = fHgOuterMC->Read("3D-PaddingMatrixOuterMC"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixOuterMC' was read in" << endl; + } + + //------------------------------------ + + fHSigmaThetaON = new TH2D; + OK =fHSigmaThetaON->Read("2D-ThetaSigmabarON"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl; + } + + fHSigmaThetaOuterON = new TH2D; + OK = fHSigmaThetaOuterON->Read("2D-ThetaSigmabarOuterON"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOuterON' was read in" << endl; + } + + fHDiffPixThetaON = new TH3D; + OK = fHDiffPixThetaON->Read("3D-ThetaPixDiffON"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl; + } + + fHDiffPixThetaTargetON = new TH3D; + OK = fHDiffPixThetaTargetON->Read("3D-ThetaPixDiffTargetON"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffTargetON' was read in" << endl; + } + + fHgON = new TH3D; + OK = fHgON->Read("3D-PaddingMatrixON"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixON' was read in" << endl; + } + + fHgOuterON = new TH3D; + OK = fHgOuterON->Read("3D-PaddingMatrixOuterON"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixOuterON' was read in" << endl; + } + + //------------------------------------ + + fHSigmaThetaOFF = new TH2D; + OK = fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOFF' was read in" << endl; + } + + fHSigmaThetaOuterOFF = new TH2D; + OK = fHSigmaThetaOuterOFF->Read("2D-ThetaSigmabarOuterOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '2D-ThetaSigmabarOuterOFF' was read in" << endl; + } + + fHDiffPixThetaOFF = new TH3D; + OK = fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl; + } + + fHDiffPixThetaTargetOFF = new TH3D; + OK = fHDiffPixThetaTargetOFF->Read("3D-ThetaPixDiffTargetOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-ThetaPixDiffTargetOFF' was read in" << endl; + } + + fHgOFF = new TH3D; + OK = fHgOFF->Read("3D-PaddingMatrixOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl; + } + + fHgOuterOFF = new TH3D; + OK = fHgOuterOFF->Read("3D-PaddingMatrixOuterOFF"); + if (OK) + { + *fLog << all + << "MPad : Object '3D-PaddingMatrixOuterOFF' was read in" << endl; + } + + //------------------------------------ + + return kTRUE; +} + + + + +// -------------------------------------------------------------------------- +// +// 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 : There are no input histograms for the padding ... 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; + } + *fLog << all << "MPad::PreProcess(); name of MPedPhotCam container = " + << fNamePedPhotCam << endl; + + + 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"); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->UseCurrentStyle(); + fHSigmaPedestal->GetYaxis()->SetTitleOffset(1.25); + + // plot no.of photons (before vs. after padding) + fHPhotons = new TH2D("Photons/area", + "Photons/area: after vs.before padding", + 100, -100.0, 300.0, 100, -100, 300); + fHPhotons->SetXTitle("no.of photons/area before padding"); + fHPhotons->SetYTitle("no.of photons/area after padding"); + fHPhotons->SetDirectory(NULL); + fHPhotons->UseCurrentStyle(); + fHPhotons->GetYaxis()->SetTitleOffset(1.25); + + // plot of added NSB + fHNSB = new TH1D("NSB/area","Distribution of added NSB/area", + 100, -100.0, 200.0); + fHNSB->SetXTitle("no.of added NSB photons/area"); + fHNSB->SetYTitle("no.of pixels"); + fHNSB->SetDirectory(NULL); + fHNSB->UseCurrentStyle(); + fHNSB->GetYaxis()->SetTitleOffset(1.25); + + //-------------------------------------------------------------------- + + *fLog << inf << "Type of data to be padded : " << fType << endl; + *fLog << inf << "Name of MPedPhotCam container : " << fNamePedPhotCam + << endl; + + fIter = 10; + + memset(fInf, 0, sizeof(fInf)); + 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; // (= 0.179) + + //------------------------------------------------ + // 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(); + + //*fLog << all << "MPad::Process(); before padding : " << endl; + + //------------------------------------------- + // Calculate average sigma of the event + // + fPed->ReCalc(*fCam, fBad); + //*fLog << "pedestalRMS, inner and outer = " << (fPed->GetArea(0)).GetRms() + // << ", " << (fPed->GetArea(1)).GetRms() << endl; + + + // inner sigmabar/sqrt(area) + Double_t ratioinn = fCam->GetPixRatio(0); + Double_t sigbarInnerold_phot = (fPed->GetArea(0)).GetRms(); + Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton * sqrt(ratioinn); + Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold; + + // outer sigmabar/sqrt(area) + Double_t ratioout = fCam->GetPixRatio(500); + Double_t sigbarOuterold_phot = (fPed->GetArea(1)).GetRms(); + Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton * sqrt(ratioout); + Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold; + + const Double_t theta = fPointPos->GetZd(); + + //*fLog << all << "theta = " << theta << endl; + + Int_t binTheta = fHSigmaTheta->GetXaxis()->FindBin(theta); + if ( binTheta < 1 || binTheta > fHSigmaTheta->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 (nTheta) + // and number of events in this sigbarInnerold_phot bin (nSigma) + + 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 ? + // + + 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; + } + + Int_t binSig = Hg->GetYaxis()->FindBin(sigbarInnerold_phot); + //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", " + // << sigbarInnerold_phot << endl; + + Double_t prob; + TH1D *hpad = NULL; + + + 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, Hg = " << 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 = 0; + fInf[rc]++; + return kTRUE; + } + // event should be padded + //*fLog << all << "event will be padded" << endl; + + + //------------------------------------------- + // for the current theta, generate a sigmabar_inner : + // 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 * sqrt(ratioinn); + + //*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; + } + + //------------------------------------------- + // estimate a sigmabar_outer from sigmabar_inner : + // using the target distributions fHSigmaTheta and fHSigmaThetaOuter + // for sigmabar(inner) and sigmabar(outer) + + Double_t innerMeantarget = 0.0; + Double_t outerMeantarget = 0.0; + Double_t innerRMStarget = 0.0; + Double_t outerRMStarget = 0.0; + + // projection doesn't work + //TH1D *hi = fHSigmaTheta->ProjectionY("proinner", binTheta, binTheta, "e"); + //TH1D *ho = fHSigmaThetaOuter->ProjectionY("proouter", binTheta, binTheta, "e"); + //sigmabarOuter2 = sigmabarInner2 + fPEperPhoton*fPEperPhoton * + // ( ho->GetMean()*ho->GetMean()*ratioout + // - hi->GetMean()*hi->GetMean()*ratioinn); + + //innerMeantarget = hi->GetMean()*sqrt(ratioinn)*fPEperPhoton; + //outerMeantarget = ho->GetMean()*sqrt(ratioout)*fPEperPhoton; + //innerRMStarget = hi->GetRMS(1)*sqrt(ratioinn)*fPEperPhoton; + //outerRMStarget = ho->GetRMS(1)*sqrt(ratioout)*fPEperPhoton; + + Double_t y = 0.0; + Double_t ybar = 0.0; + Double_t y2bar = 0.0; + Double_t w = 0.0; + + Double_t isum = 0.0; + Double_t isumy = 0.0; + Double_t isumy2 = 0.0; + for (Int_t i=1; i<=fHSigmaTheta->GetNbinsY(); i++) + { + w = fHSigmaTheta->GetBinContent(binTheta, i); + y = fHSigmaTheta->GetYaxis()->GetBinCenter(i); + isum += w; + isumy += w * y; + isumy2 += w * y*y; + } + if (isum == 0.0) + { + innerMeantarget = 0.0; + innerRMStarget = 0.0; + } + else + { + ybar = isumy /isum; + y2bar = isumy2/isum; + innerMeantarget = ybar * sqrt(ratioinn)*fPEperPhoton; + innerRMStarget = sqrt( y2bar - ybar*ybar ) * sqrt(ratioinn)*fPEperPhoton; + } + + Double_t osum = 0.0; + Double_t osumy = 0.0; + Double_t osumy2 = 0.0; + for (Int_t i=1; i<=fHSigmaThetaOuter->GetNbinsY(); i++) + { + w = fHSigmaThetaOuter->GetBinContent(binTheta, i); + y = fHSigmaThetaOuter->GetYaxis()->GetBinCenter(i); + osum += w; + osumy += w * y; + osumy2 += w * y*y; + } + if (osum == 0.0) + { + outerMeantarget = 0.0; + outerRMStarget = 0.0; + } + else + { + ybar = osumy /osum; + y2bar = osumy2/osum; + outerMeantarget = ybar * sqrt(ratioout)*fPEperPhoton; + outerRMStarget = sqrt( y2bar - ybar*ybar ) * sqrt(ratioout)*fPEperPhoton; + } + + + // new outer sigmabar2/area + Double_t sigmabarOuter2; + + Double_t scal = ( innerMeantarget*innerRMStarget == 0.0 || + outerMeantarget*outerRMStarget == 0.0 ) ? 1.0 : + (outerMeantarget*outerRMStarget)/(innerMeantarget*innerRMStarget); + sigmabarOuter2 = outerMeantarget*outerMeantarget + + scal * (sigmabarInner2 - innerMeantarget*innerMeantarget); + + //*fLog << "innerMeantarget, innerRMStarget = " << innerMeantarget + // << ", " << innerRMStarget << endl; + + //*fLog << "outerMeantarget, outerRMStarget = " << outerMeantarget + // << ", " << outerRMStarget << endl; + + //*fLog << "sigmabarInner2, sigmabarOuter2, scal = " << sigmabarInner2 + // << ", " << sigmabarOuter2 << ", " << scal << endl; + + //delete hi; + //delete ho; + + //------------------------------------------- + // + // Calculate average number of NSB photo electrons to be added (lambdabar) + // from the value of sigmabar, + // - 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] + + //---------------- + TH3D *sp = NULL; + if (fType == "MC") sp = fHDiffPixThetaTargetMC; + else if (fType == "ON") sp = fHDiffPixThetaTargetON; + else if (fType == "OFF") sp = fHDiffPixThetaTargetOFF; + + //---------------- + + Int_t bincheck = sp->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)/area 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 = + sp->GetYaxis()->FindBin( (Double_t)idmaxpixinner ); + + TH1D *hnoise = NULL; + hnoise = sp->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); + + 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 << "RMS_phot, elNoise2 = " << RMS_phot << ", " + // << elNoise2 << endl; + + lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess; + + if (lambdabar <= 0.0) + { + rc = 3; + fErrors[rc]++; + } + + //*fLog << "lambdabar = " << lambdabar << endl; + + // This value of lambdabar is the same for all pixels; + // note that lambdabar is the NSB p.e. density + + //---------- start loop over pixels --------------------------------- + // do the padding for each pixel + // + // pad only pixels - which are used + // + + 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 fHDiffPixThetaTargetMC + // ON : from fHDiffPixThetaTaregtON + // OFF : from fHDiffPixThetaTargetOFF + + + Double_t sigbarold2; + Double_t sigmabar2; + + for (UInt_t i=0; iGetPixRatio(j); + + if (ratio > 0.5) + { + sigbarold2 = sigbarInnerold2; + sigmabar2 = sigmabarInner2; + } + else + { + sigbarold2 = sigbarOuterold2; + sigmabar2 = sigmabarOuter2; + } + + MPedPhotPix &ppix = (*fPed)[j]; + + // count number of pixels treated + fWarnings[0]++; + + + Double_t oldsigma_phot = ppix.GetRms(); + Double_t oldsigma = oldsigma_phot * fPEperPhoton * sqrt(ratio); + Double_t oldsigma2 = oldsigma*oldsigma; + + //--------------------------------- + // throw the Sigma for this pixel + // + Int_t binPixel = sp->GetYaxis()->FindBin( (Double_t)j ); + + Int_t count; + Bool_t ok; + + TH1D *hdiff = NULL; + + 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 the old difference + 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, ratio, lambdabar = " + // << theta << ", " + // << j << ", " << count << ", " << sigmabar2 << ", " + // << diff << ", " << oldsigma2 << ", " << ratio << ", " + // << 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) / ratio; + addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess/ratio; + + + //--------------------------------- + // throw actual number of additional NSB photo electrons (NSB) + // and its RMS (sigmaNSB) + + Double_t NSB0 = gRandom->Poisson(lambdabar/ratio); + 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/ratio; + 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*ratio ); + fHPhotons->Fill( oldphotons_phot*ratio, + newphotons_phot*ratio ); + + + // 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 / ratio ); + Double_t newsigma_phot = newsigma / fPEperPhoton; + ppix.SetRms( newsigma_phot ); + + fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot ); + } + //---------- end of loop over pixels --------------------------------- + + + //*fLog << all << "MPad::Process(); after padding : " << endl; + + // Calculate sigmabar again and crosscheck + fPed->ReCalc(*fCam, fBad); + //*fLog << "pedestalRMS, inner and outer = " << (fPed->GetArea(0)).GetRms() + // << ", " << (fPed->GetArea(1)).GetRms() << 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: new sigma <= old sigma" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) + << (int)(fErrors[3]*100/GetNumExecutions()) + << "%) lambdabar = 0" << endl; + + *fLog << " " << setw(7) << fInf[0] << " (" << setw(3) + << (int)(fInf[0]*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, 1200)); + c.Divide(3, 4); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + fHNSB->DrawCopy(); + fHNSB->SetBit(kCanDelete); + + c.cd(2); + fHSigmaPedestal->DrawCopy(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(3); + fHPhotons->DrawCopy(); + fHPhotons->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + if (fHgON) + { + c.cd(4); + TH2D *m1; + m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); + m1->SetDirectory(NULL); + m1->UseCurrentStyle(); + 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);; + } + + if (fHgOFF) + { + c.cd(5); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->UseCurrentStyle(); + 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);; + } + + if (fHgMC) + { + c.cd(6); + TH2D *m3; + m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); + m3->SetDirectory(NULL); + m3->UseCurrentStyle(); + 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);; + } + + //-------------------------------------------------------------------- + + if (fHgOuterON) + { + c.cd(7); + TH2D *m1; + m1 = (TH2D*) ((TH3*)fHgOuterON)->Project3D("zy"); + m1->SetDirectory(NULL); + m1->UseCurrentStyle(); + m1->SetTitle("(fHgOuterON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetXTitle("Sigmabar-old"); + m1->SetYTitle("Sigmabar-new"); + + m1->DrawCopy("box"); + m1->SetBit(kCanDelete);; + } + + if (fHgOuterOFF) + { + c.cd(8); + TH2D *m2; + m2 = (TH2D*) ((TH3*)fHgOuterOFF)->Project3D("zy"); + m2->SetDirectory(NULL); + m2->UseCurrentStyle(); + m2->SetTitle("(fHgOuterOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetXTitle("Sigmabar-old"); + m2->SetYTitle("Sigmabar-new"); + + m2->DrawCopy("box"); + m2->SetBit(kCanDelete);; + } + + if (fHgOuterMC) + { + c.cd(9); + TH2D *m3; + m3 = (TH2D*) ((TH3*)fHgOuterMC)->Project3D("zy"); + m3->SetDirectory(NULL); + m3->UseCurrentStyle(); + m3->SetTitle("(fHgOuterMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); + m3->SetXTitle("Sigmabar-old"); + m3->SetYTitle("Sigmabar-new"); + + m3->DrawCopy("box"); + m3->SetBit(kCanDelete);; + } + + //-------------------------------------------------------------------- + + c.cd(10); + fHSigmaTheta->SetDirectory(NULL); + fHSigmaTheta->UseCurrentStyle(); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + c.cd(11); + fHSigmaThetaOuter->SetDirectory(NULL); + fHSigmaThetaOuter->UseCurrentStyle(); + fHSigmaThetaOuter->DrawCopy(); + fHSigmaThetaOuter->SetBit(kCanDelete); + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.h (revision 5432) +++ trunk/MagicSoft/Mars/manalysis/MPad.h (revision 5432) @@ -0,0 +1,131 @@ +#ifndef MARS_MPad +#define MARS_MPad + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class TH3D; + +class MGeomCam; +class MCerPhotEvt; +class MPedPhotCam; +class MPointingPos; +class MParList; +class MBadPixelsCam; +class MRead; +class MFilterList; + + +class MPad : public MTask +{ +private: + MGeomCam *fCam; + MCerPhotEvt *fEvt; + MPointingPos *fPointPos; + MPedPhotCam *fPed; + MBadPixelsCam *fBad; + + TString fNamePedPhotCam; // name of the 'MPedPhotCam' container + TString fType; // type of data to be padded + TFile *fInfile; // input file containing padding histograms + + Int_t fIter; + + Int_t fInf[9]; + Int_t fErrors[9]; + Int_t fWarnings[3]; + + //---------------------------------- + // plots used for the padding + // for all plots it is assumed that the pedestal RMS is given in units of "number of photons" + + // original distributions + TH2D *fHSigmaThetaMC; // 2D-histogram (sigmabar_inner vs. Theta) + TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar_inner vs. Theta) + TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar_inner vs. Theta) + + TH2D *fHSigmaThetaOuterMC; // 2D-histogram (sigmabar_outer vs. Theta) + TH2D *fHSigmaThetaOuterON; // 2D-histogram (sigmabar_outer vs. Theta) + TH2D *fHSigmaThetaOuterOFF; // 2D-histogram (sigmabar_outer vs. Theta) + + TH3D *fHDiffPixThetaMC; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + + //--------------------- + // target distributions + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar_inner vs. Theta) + TH2D *fHSigmaThetaOuter; // 2D-histogram (sigmabar_outer vs. Theta) + + TH3D *fHDiffPixThetaTargetMC; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + TH3D *fHDiffPixThetaTargetON; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + TH3D *fHDiffPixThetaTargetOFF; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/area ) + + //--------------------- + // matrices according to which the padding is performed + TH3D *fHgMC; // matrix (Theta, sigbarold, sigbarnew) for MC data + TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data + TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data + + TH3D *fHgOuterMC; // matrix (Theta, sigbarold, sigbarnew) for MC data + TH3D *fHgOuterON; // matrix (Theta, sigbarold, sigbarnew) for ON data + TH3D *fHgOuterOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data + + //------------------------------- + // plots for checking the padding + TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after + // versus before padding + TH2D *fHPhotons; // 2D-histogram : no.of photons/area after + // versus before padding + TH1D *fHNSB; // 1D-histogram : additional NSB/area + + //------------------------------- + Bool_t MergeABC(TString tA, TString tB, TString tC, + TH2D *hA, TH2D *hB, TH2D *hC, TH2D *hM, + TH3D *hgA, TH3D *hgB, TH3D *hgC, TString canv); + + Bool_t Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, + TH2D *fHga, TH2D *fHgb, Int_t nbinssig, + TString canv); + + Bool_t UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, + Int_t nbinssig, Int_t l); + +public: + MPad(const char *name=NULL, const char *title=NULL); + ~MPad(); + + void SetDataType(const char *type); // type of data to be padded + void SetNamePedPhotCam(const char *name); // name of MPedPhotCam container + + Bool_t ReadPadHistograms(TString type, TString filein); + + Bool_t MergeONOFFMC(TString nameon="", TString nameoff="", + TString namemc="", TString fileout=""); + + Bool_t WritePaddingDist(TString fileout); + Bool_t ReadPaddingDist(TString filein); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + ClassDef(MPad, 0) // task for the padding +}; + +#endif + + + + + + + +