Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 5291) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 5431) @@ -51,4 +51,5 @@ #pragma link C++ class MPad+; #pragma link C++ class MMcCalibrationUpdate+; +#pragma link C++ class MMakePadHistograms+; #endif Index: trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc (revision 5431) +++ trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.cc (revision 5431) @@ -0,0 +1,350 @@ +/* ======================================================================== *\ +! +! * +! * 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, 09/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MMakePadHistograms +// +// reads calibrated data of a particular type (ON, OFF or MC) +// and produces the histograms for the padding +// +// these histograms are further treated in MPad +// +///////////////////////////////////////////////////////////////////////////// +#include "MMakePadHistograms.h" + +#include +#include + +#include +#include +#include +#include + +#include "MBadPixelsCalc.h" +#include "MBadPixelsTreat.h" +#include "MContinue.h" +#include "MEvtLoop.h" +#include "MFillH.h" +#include "MFSelBasic.h" +#include "MGeomApply.h" +#include "MHBadPixels.h" +#include "MHSigmaTheta.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MProgressBar.h" +#include "MReadMarsFile.h" +#include "MReadReports.h" +#include "MTaskList.h" +//#include "MExtrapolatePointingPos.h" +#include "MPointingPosCalc.h" + +ClassImp(MMakePadHistograms); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MMakePadHistograms::MMakePadHistograms(const char *name, const char *title) +{ + fName = name ? name : "MMakePadHistograms"; + fTitle = title ? title : "Class to make Pad histograms"; + + fHSigmaTheta = NULL; + fHSigmaThetaOuter = NULL; + fHDiffPixTheta = NULL; + fHSigmaPixTheta = NULL; + + fHBadIdTheta = NULL; + fHBadNTheta = NULL; + + // default values + fType = ""; + fMaxEvents = -1; + fNamePedPhotCam = "MPedPhotCamFromData"; + fNameInputFile = ""; + fNameOutputFile = ""; + fPedestalLevel = 3.0; + fUseInterpolation = kTRUE; + fProcessPedestal = kTRUE; + fProcessTime = kFALSE; +} + + +// -------------------------------------------------------------------------- +// +// SetDataType. +// +// the type may be ON, OFF or MC +// +void MMakePadHistograms::SetDataType(const char *type) +{ + fType = type; +} + +// -------------------------------------------------------------------------- +// +// SetMaxEvents. +// +// set the maximum number of events to be read +// +void MMakePadHistograms::SetMaxEvents(Int_t maxev) +{ + fMaxEvents = maxev; +} + +// -------------------------------------------------------------------------- +// +// SetNameInputFile. +// +// the input file contains the calibrated data +// +void MMakePadHistograms::SetNameInputFile(const char *input) +{ + fNameInputFile = input; +} + +// -------------------------------------------------------------------------- +// +// SetNameOutputFile. +// +// the output file contains the histgrams for the padding +// +void MMakePadHistograms::SetNameOutputFile(const char *output) +{ + fNameOutputFile = output; +} + +// -------------------------------------------------------------------------- +// +// SetNamePedPhotCam. +// +// set the name of the MPedPhotCam container; +// it is used in : MBadPixelsCalc, MBadPixelsTreat, MHSigmaTheta +// +void MMakePadHistograms::SetNamePedPhotCam(const char *name) +{ + fNamePedPhotCam = name; +} + +// -------------------------------------------------------------------------- +// +// SetPedestalLevel. +// +// set the pedestal level; +// it is used in : MBadPixelsCalc +// +void MMakePadHistograms::SetPedestalLevel(Float_t pedlevel) +{ + fPedestalLevel = pedlevel; +} + +// -------------------------------------------------------------------------- +// +// SetUseInterpolation. +// +// set the option "UseInterpolation"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetUseInterpolation(Bool_t useinter) +{ + fUseInterpolation = useinter; +} + +// -------------------------------------------------------------------------- +// +// SetProcessPedestal. +// +// set the option "ProcessPedestal"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetProcessPedestal(Bool_t procped) +{ + fProcessPedestal = procped; +} + +// -------------------------------------------------------------------------- +// +// SetProcessTime. +// +// set the option "ProcessTime"; +// it is used in : MBadPixelsTreat +// +void MMakePadHistograms::SetProcessTime(Bool_t proctim) +{ + fProcessTime = proctim; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MMakePadHistograms::MakeHistograms() +{ + if (fType == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; type of input data is not defined" + << endl; + return kFALSE; + } + + if (fNameInputFile == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; name of input file is not defined" + << endl; + return kFALSE; + } + + if (fNameOutputFile == "") + { + *fLog << err << "MMakePadHistograms::MakeHistograms; name of output file is not defined" + << endl; + return kFALSE; + } + + *fLog << "------------" << endl; + *fLog << fType << " events :" << endl; + *fLog << "------------" << endl; + *fLog << "MMakePadHistograms : calibrated " << fType + << " data are read from file " << fNameInputFile << endl; + + + MTaskList tlist; + MParList plist; + + MReadReports read; + read.AddTree("Drive"); + read.AddTree("Events","MTime.",kTRUE); + static_cast(read).AddFile(fNameInputFile); + + read.AddToBranchList("MReportDrive.*"); + read.AddToBranchList("MRawEvtHeader.*"); + +// MReadMarsFile read("Events", fNameInputFile); +// read.DisableAutoScheme(); + + + MGeomApply apply; + +// MExtrapolatePointingPos extraPointing(fNameInputFile); + MPointingPosCalc pointposcalc; + + MBadPixelsCalc badcalc; + badcalc.SetNamePedPhotContainer(fNamePedPhotCam); + badcalc.SetPedestalLevel(fPedestalLevel); + + MBadPixelsTreat badtreat; + badtreat.SetNamePedPhotCam(fNamePedPhotCam); + badtreat.SetUseInterpolation(fUseInterpolation); + badtreat.SetProcessPedestal(fProcessPedestal); + badtreat.SetProcessTimes(fProcessTime); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + MHBadPixels bad("BadPixels"); + bad.SetNamePedPhotCam(fNamePedPhotCam); + MFillH fillbad("BadPixels[MHBadPixels]", "MBadPixelsCam"); + fillbad.SetName("FillBad"); + + MHSigmaTheta sigth("SigmaTheta"); + sigth.SetNamePedPhotCam(fNamePedPhotCam); + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", ""); + fillsigtheta.SetName("FillSigTheta"); + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(&bad); + plist.AddToList(&sigth); + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + tlist.AddToList(&apply); + tlist.AddToList(&pointposcalc, "Events"); + +// tlist.AddToList(&badcalc); done in callisto +// tlist.AddToList(&badtreat); done in callisto + + tlist.AddToList(&contbasic, "Events"); + tlist.AddToList(&fillbad, "Events"); + tlist.AddToList(&fillsigtheta, "Events"); + //***************************** + + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetProgressBar(&bar); + + // if ( !evtloop.Eventloop(fMaxEvents) ) + // return kFALSE; + evtloop.Eventloop(fMaxEvents); + + tlist.PrintStatistics(0, kTRUE); + + bad.DrawClone(); + sigth.DrawClone(); + + // save the histograms for the padding + fHSigmaTheta = (TH2D*)sigth.GetSigmaTheta(); + fHSigmaThetaOuter= (TH2D*)sigth.GetSigmaThetaOuter(); + fHSigmaPixTheta = (TH3D*)sigth.GetSigmaPixTheta(); + fHDiffPixTheta = (TH3D*)sigth.GetDiffPixTheta(); + + fHBadIdTheta = (TH2D*)bad.GetBadId(); + fHBadNTheta = (TH2D*)bad.GetBadN(); + + TFile filewrite(fNameOutputFile, "RECREATE", ""); + fHSigmaTheta->Write(); + fHSigmaThetaOuter->Write(); + //fHSigmaPixTheta->Write(); + fHDiffPixTheta->Write(); + //fHBadIdTheta->Write(); + //fHBadNTheta->Write(); + filewrite.Close(); + + *fLog << "" << endl; + *fLog << "MMakePadHistograms : padding plots for " << fType + << " data written onto file " << fNameOutputFile << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h (revision 5431) +++ trunk/MagicSoft/Mars/manalysis/MMakePadHistograms.h (revision 5431) @@ -0,0 +1,73 @@ +#ifndef MARS_MMakePadHistograms +#define MARS_MMakePadHistograms + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH2D; +class TH3D; + +class MMakePadHistograms : public MH +{ +private: + + TString fType; // type of data to be padded (ON, OFF, MC) + Int_t fMaxEvents; // maximum no.of events to be read + TString fNamePedPhotCam; // name of the 'MPedPhotCam' container + TString fNameInputFile; // name of input file (calibrated data) + TString fNameOutputFile; // name of output file (Pad histograms) + Float_t fPedestalLevel; // pedestal level for MBadPixelsCalc + Bool_t fUseInterpolation; // option for MBadPixelsTreat + Bool_t fProcessPedestal; // option for MBadPixelsTreat + Bool_t fProcessTime; // option for MBadPixelsTreat + + + //---------------------------------- + // plots needed for the padding + + TH2D *fHSigmaTheta; // 2D-histogram (Theta, sigmabar(inner)) + TH2D *fHSigmaThetaOuter; // 2D-histogram (Theta, sigmabar(outer)) + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, (sigma^2-sigmabar^2)/Area) + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + + TH2D *fHBadIdTheta; // 2D-histogram (Theta, Id of bad pixel) + TH2D *fHBadNTheta; // 2D-histogram (Theta, no.of bad pixels) + + +public: + MMakePadHistograms(const char *name=NULL, const char *title=NULL); + + void SetDataType(const char *); // type of data (ON, OFF, MC) + void SetMaxEvents(Int_t); // mx. no.of events to be read + void SetNameInputFile(const char *); // Calibrated data + void SetNameOutputFile(const char *); // Pad histograms + void SetNamePedPhotCam(const char *); // "MPedPhotCamFromData" + + // parameters for MBadPixelsCalc + void SetPedestalLevel(Float_t); + + // parameters for MBadPixelsTreat + void SetUseInterpolation(Bool_t); + void SetProcessPedestal(Bool_t); + void SetProcessTime(Bool_t); + + Bool_t MakeHistograms(); + + ClassDef(MMakePadHistograms, 0) // class for making the padding histograms +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 5291) +++ (revision ) @@ -1,2282 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Wolfgang Wittek, 06/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// 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; -} - -// -------------------------------------------------------------------------- - - - - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.h (revision 5291) +++ (revision ) @@ -1,121 +1,0 @@ -#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 MSigmabar; -class MParList; -class MBadPixelsCam; -class MRead; -class MFilterList; - - -class MPad : public MTask -{ -private: - MGeomCam *fCam; - MCerPhotEvt *fEvt; - MSigmabar *fSigmabar; - 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 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" - - TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaMC; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaON; // 2D-histogram (sigmabar vs. Theta) - TH2D *fHSigmaThetaOFF; // 2D-histogram (sigmabar vs. Theta) - - //--------------------- - TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaMC; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaON; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - TH3D *fHDiffPixThetaOFF; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) - - //--------------------- - TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaMC; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaON; // 3D-histogram (Theta, pixel, sigma) - TH3D *fHSigmaPixThetaOFF; // 3D-histogram (Theta, pixel, sigma) - - //--------------------- - TH3D *fHgMC; // matrix (Theta, sigbarold, sigbarnew) for MC data - TH3D *fHgON; // matrix (Theta, sigbarold, sigbarnew) for ON data - TH3D *fHgOFF; // matrix (Theta, sigbarold, sigbarnew) for OFF data - - //------------------------------- - // plots for checking the padding - TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after - // versus before padding - TH2D *fHPhotons; // 2D-histogram : no.of photons after - // versus before padding - TH1D *fHNSB; // 1D-histogram : additional NSB - - - Bool_t Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, - TH2D *fHga, TH2D *fHgb, Int_t nbinssig); - - 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 SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } - - Bool_t MergeONOFFMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, - TH2D& sigthoff, TH3D& diffpixthoff,TH3D& sigmapixthoff); - - Bool_t MergeONMC( - TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, - TH2D& sigthon, TH3D&diffpixthon, TH3D& sigmapixthon); - - Bool_t ReadPaddingDist(const char *filein); - Bool_t WritePaddingDist(const char *fileout); - - void SetDataType(const char *type); // type of data to be padded - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - ClassDef(MPad, 0) // task for the padding -}; - -#endif - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 5291) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 5431) @@ -23,5 +23,5 @@ -I../mhistmc -I../mfileio -I../mmain -I../mhcalib -I../mcalib \ -I../msignal -I../mpointing -I../mtools -I../mfbase \ - -I../mbadpixels -I../mastro -I../mpedestal + -I../mbadpixels -I../mastro -I../mpedestal -I../mtemp/mucm/classes # mcalib: MGeomApply (MCalibrationCam) @@ -55,5 +55,6 @@ MFiltercutsCalc.cc \ MPad.cc \ - MMcCalibrationUpdate.cc + MMcCalibrationUpdate.cc \ + MMakePadHistograms.cc ############################################################ @@ -66,2 +67,3 @@ mrproper: clean rmbak +