Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1748) +++ trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1748) @@ -0,0 +1,363 @@ +/* ======================================================================== *\ +! +! * +! * 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 1/2003 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHSigmaTheta (extension of Robert's MHSigmabarTheta) // +// // +// calculates - the 2D-histogram sigmabar vs. Theta, and // +// - the 3D-histogram sigma, pixel no., Theta // +// - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHSigmaTheta.h" + +#include + +#include "MTime.h" +#include "MMcEvt.hxx" + +#include "MBinning.h" +#include "MParList.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MCerPhotEvt.h" +#include "MCerPhotPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MSigmabar.h" + +ClassImp(MHSigmaTheta); + + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHSigmaTheta::MHSigmaTheta(const char *name, const char *title) + : fSigmaTheta(), fSigmaPixTheta() +{ + fName = name ? name : "MHSigmaTheta"; + fTitle = title ? title : "2D histogram sigmabar vs. Theta"; + + fSigmaTheta.SetDirectory(NULL); + fSigmaTheta.SetName("2D-ThetaSigmabar"); + fSigmaTheta.SetTitle("2D : Sigmabar, \\Theta"); + fSigmaTheta.SetXTitle("\\Theta [\\circ]"); + fSigmaTheta.SetYTitle("Sigmabar"); + + fSigmaPixTheta.SetDirectory(NULL); + fSigmaPixTheta.SetName("3D-ThetaPixSigma"); + fSigmaPixTheta.SetTitle("3D : \\Theta, pixel no., Sigma"); + fSigmaPixTheta.SetXTitle("\\Theta [\\circ]"); + fSigmaPixTheta.SetYTitle("pixel number"); + fSigmaPixTheta.SetZTitle("Sigma"); + + fDiffPixTheta.SetDirectory(NULL); + fDiffPixTheta.SetName("3D-ThetaPixDiff"); + fDiffPixTheta.SetTitle("3D : \\Theta, pixel, Sigma^2-Sigmabar^2"); + fDiffPixTheta.SetXTitle("\\Theta [\\circ]"); + fDiffPixTheta.SetYTitle("pixel number"); + fDiffPixTheta.SetZTitle("Sigma^2-sigmabar^2"); +} + +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histogram +// +Bool_t MHSigmaTheta::SetupFill(const MParList *plist) +{ + + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar"); + if (!fSigmabar) + { + *fLog << dbginf << "MHSigmaTheta::SetupFill : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get Theta Binning + MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binstheta) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningTheta not found... aborting." << endl; + return kFALSE; + } + + // Get Sigmabar binning + MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar"); + if (!binssigma) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningSigmabar not found... aborting." << endl; + return kFALSE; + } + + // Get binning for (sigma^2-sigmabar^2) + MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2"); + if (!binsdiff) + { + *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningDiffsigma2 not found... aborting." << endl; + return kFALSE; + } + + + // Set binnings in histograms + SetBinning(&fSigmaTheta, binstheta, binssigma); + + // Get binning for pixel number + UInt_t npix = fPed->GetSize(); + MBinning binspix("BinningPixel"); + MBinning* binspixel = &binspix; + binspixel->SetEdges(npix, -0.5, ((float)npix)-0.5 ); + + SetBinning(&fSigmaPixTheta, binstheta, binspixel, binssigma); + SetBinning(&fDiffPixTheta, binstheta, binspixel, binsdiff); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms +// +Bool_t MHSigmaTheta::Fill(const MParContainer *par) +{ + //*fLog << "entry Fill" << endl; + + Double_t Theta = fMcEvt->GetTelescopeTheta()*kRad2Deg; + Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt); + fSigmaTheta.Fill(Theta, mySig); + + //*fLog << "Theta, mySig = " << Theta << ", " << mySig << endl; + + const UInt_t npix = fEvt->GetNumPixels(); + for (Int_t i=0; ioperator[](i); + if (!cerpix.IsPixelUsed()) + continue; + + if (cerpix.GetNumPhotons() == 0) + continue; + + Int_t j = cerpix.GetPixId(); + const MPedestalPix pix = fPed->operator[](j); + + Double_t Sigma = pix.GetMeanRms(); + Double_t Area = fCam->GetPixRatio(j); + + fSigmaPixTheta.Fill(Theta, (Double_t)j, Sigma); + + Double_t Diff = Sigma*Sigma/Area - mySig*mySig; + fDiffPixTheta.Fill (Theta, (Double_t)j, Diff); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Plot the results +// +Bool_t MHSigmaTheta::Finalize() +{ + DrawClone(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw a copy of the histogram +// +TObject *MHSigmaTheta::DrawClone(Option_t *opt) +{ + TCanvas &c = *MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", + 900, 900); + c.Divide(3, 3); + + gROOT->SetSelectedPad(NULL); + + //-------------------------------------------------------------------- + // draw the 2D histogram Sigmabar versus Theta + TH1D *h; + + c.cd(1); + h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + gPad->SetLogy(); + + c.cd(4); + h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); + h->SetTitle("Distribution of Sigmabar"); + h->SetXTitle("Sigmabar"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + + c.cd(7); + ((TH2*)&fSigmaTheta)->DrawCopy(opt); + + //-------------------------------------------------------------------- + // draw the 3D histogram : Theta, pixel, Sigma^2-Sigmabar^2 + + TH2D *l; + + c.cd(2); + l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx"); + l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l->SetXTitle("\\Theta [\\circ]"); + l->SetYTitle("Sigma^2-Sigmabar^2"); + + l->Draw("box"); + l->SetBit(kCanDelete);; + + c.cd(5); + l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy"); + l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l->SetXTitle("pixel"); + l->SetYTitle("Sigma^2-Sigmabar^2"); + + l->Draw("box"); + l->SetBit(kCanDelete);; + + c.cd(8); + ((TH2*)&fDiffPixTheta)->DrawCopy(opt); + + + //-------------------------------------------------------------------- + // draw the 3D histogram : Theta, pixel, Sigma + + TH2D *k; + + c.cd(3); + k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx"); + k->SetTitle("Sigma vs. \\Theta (all pixels)"); + k->SetXTitle("\\Theta [\\circ]"); + k->SetYTitle("Sigma"); + + k->Draw("box"); + k->SetBit(kCanDelete);; + + c.cd(6); + k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy"); + k->SetTitle("Sigma vs. pixel number (all \\Theta)"); + k->SetXTitle("pixel"); + k->SetYTitle("Sigma"); + + k->Draw("box"); + k->SetBit(kCanDelete);; + + c.cd(9); + ((TH2*)&fSigmaPixTheta)->DrawCopy(opt); + + //-------------------------------------------------------------------- + c.Modified(); + c.Update(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHSigmaTheta::Draw(Option_t *opt) +{ + if (!gPad) + MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", + 600, 600); + + TH1D *h; + + gPad->Divide(2,2); + + gPad->cd(1); + h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + gPad->SetLogy(); + + gPad->cd(2); + h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); + h->SetTitle("Distribution of Sigmabar"); + h->SetXTitle("Sigmabar"); + h->SetYTitle("No.of events"); + + h->Draw(opt); + h->SetBit(kCanDelete);; + + gPad->cd(3); + fSigmaTheta.DrawCopy(opt); + + + + gPad->Modified(); + gPad->Update(); +} +// -------------------------------------------------------------------------- + + + + Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 1748) +++ trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 1748) @@ -0,0 +1,73 @@ +#ifndef MARS_MHSigmaTheta +#define MARS_MHSigmaTheta + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +#ifndef ROOT_TH3 +#include "TH3.h" +#endif + +#ifndef ROOT_TProfile2D +#include "TProfile2D.h" +#endif + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MPedestalCam; +class MSigmabar; +class MParList; + + +class MHSigmaTheta : public MH +{ +private: + const MGeomCam *fCam; + MPedestalCam *fPed; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + + TH2D fSigmaTheta; // 2D-distribution sigmabar versus Theta; + // sigmabar is the average pedestasl sigma + // in an event + TH3D fSigmaPixTheta;// 3D-distr.:Theta, pixel, pedestal sigma + TH3D fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 + + +public: + MHSigmaTheta(const char *name=NULL, const char *title=NULL); + + virtual Bool_t SetupFill(const MParList *plist); + virtual Bool_t Fill(const MParContainer *par); + virtual Bool_t Finalize(); + + const TH2D *GetSigmaTheta() { return &fSigmaTheta; } + const TH2D *GetSigmaTheta() const { return &fSigmaTheta; } + TH2D *GetSigmaThetaByName(const TString name) { return &fSigmaTheta; } + + const TH3D *GetSigmaPixTheta() { return &fSigmaPixTheta; } + const TH3D *GetSigmaPixTheta() const { return &fSigmaPixTheta; } + TH3D *GetSigmaPixThetaByName(const TString name) { return &fSigmaPixTheta; } + + const TH3D *GetDiffPixTheta() { return &fDiffPixTheta; } + const TH3D *GetDiffPixTheta() const { return &fDiffPixTheta; } + TH3D *GetDiffPixThetaByName(const TString name) { return &fDiffPixTheta; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option=""); + + ClassDef(MHSigmaTheta, 0) //2D-histogram sigmabar vs. Theta +}; + +#endif + + +