/* ======================================================================== *\ ! ! * ! * 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 ! Author(s): Thomas Bretz 4/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // 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 "MSigmabar.h" #include "MGeomCam.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MCerPhotEvt.h" #include "MCerPhotPix.h" #include "MLog.h" #include "MLogManip.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 Id, Sigma"); fSigmaPixTheta.SetXTitle("\\Theta [\\circ]"); fSigmaPixTheta.SetYTitle("Pixel Id"); fSigmaPixTheta.SetZTitle("Sigma"); fDiffPixTheta.SetDirectory(NULL); fDiffPixTheta.SetName("3D-ThetaPixDiff"); fDiffPixTheta.SetTitle("3D : \\Theta, Pixel Id, Sigma^2-Sigmabar^2"); fDiffPixTheta.SetXTitle("\\Theta [\\circ]"); fDiffPixTheta.SetYTitle("Pixel Id"); fDiffPixTheta.SetZTitle("Sigma^2-Sigmabar^2"); } // -------------------------------------------------------------------------- // // Default Destructor // MHSigmaTheta::~MHSigmaTheta() { } // -------------------------------------------------------------------------- // // 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 << "MMcEvt not found... aborting." << endl; return kFALSE; } fPed = (MPedestalCam*)plist->FindObject("MPedestalCam"); if (!fPed) { *fLog << err << "MPedestalCam not found... aborting." << endl; return kFALSE; } fCam = (MGeomCam*)plist->FindObject("MGeomCam"); if (!fCam) { *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; return kFALSE; } fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt"); if (!fEvt) { *fLog << err << "MCerPhotEvt not found... aborting." << endl; return kFALSE; } fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar"); if (!fSigmabar) { *fLog << err << "MSigmabar not found... aborting." << endl; return kFALSE; } // Get Theta Binning MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); if (!binstheta) { *fLog << err << "BinningTheta [MBinning] not found... aborting." << endl; return kFALSE; } // Get Sigmabar binning MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar"); if (!binssigma) { *fLog << err << "BinningSigmabar [MBinning] not found... aborting." << endl; return kFALSE; } // Get binning for (sigma^2-sigmabar^2) MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2"); if (!binsdiff) { *fLog << err << "BinningDiffsigma2 [MBinning] not found... aborting." << endl; return kFALSE; } // Set binnings in histograms SetBinning(&fSigmaTheta, binstheta, binssigma); // Get binning for pixel number const UInt_t npix = fPed->GetSize(); MBinning binspix("BinningPixel"); binspix.SetEdges(npix+1, -0.5, 0.5+npix ); SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma); SetBinning(&fDiffPixTheta, binstheta, &binspix, binsdiff); return kTRUE; } // -------------------------------------------------------------------------- // // Fill the histograms // Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w) { Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg; Double_t mysig = fSigmabar->Calc(*fCam, *fPed, *fEvt); fSigmaTheta.Fill(theta, mysig); const UInt_t npix = fEvt->GetNumPixels(); for (UInt_t i=0; iGetPixRatio(id); fSigmaPixTheta.Fill(theta, (Double_t)id, sigma); const Double_t diff = sigma*sigma/area - mysig*mysig; fDiffPixTheta.Fill(theta, (Double_t)id, diff); } return kTRUE; } // -------------------------------------------------------------------------- // // Draw the histogram // void MHSigmaTheta::Draw(Option_t *opt) { TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this, 900, 900); pad->SetBorderMode(0); AppendPad(""); pad->Divide(3, 2); // draw the 2D histogram Sigmabar versus Theta TH1 *h; pad->cd(1); gPad->SetLogy(); h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E"); h->SetDirectory(NULL); h->SetTitle("Distribution of \\Theta"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("No.of events"); h->Draw(opt); h->SetBit(kCanDelete); pad->cd(2); h = fDiffPixTheta.Project3D("zx"); h->SetDirectory(NULL); h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("\\sigma_{ped}^2-\\bar{\\sigma}_{ped}^{2}"); h->Draw("box"); h->SetBit(kCanDelete); pad->cd(3); h = fSigmaPixTheta.Project3D("zx"); h->SetDirectory(NULL); h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("\\sigma_{ped}"); h->Draw("box"); h->SetBit(kCanDelete); pad->cd(4); h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E"); h->SetDirectory(NULL); h->SetTitle("Distribution of \\bar{\\sigma}_{ped}"); h->SetXTitle("\\bar{\\sigma}_{ped}"); h->SetYTitle("No.of events"); h->Draw(opt); h->SetBit(kCanDelete); pad->cd(5); h = fDiffPixTheta.Project3D("zy"); h->SetDirectory(NULL); h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)"); h->SetXTitle("Id"); h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}"); h->Draw("box"); h->SetBit(kCanDelete); pad->cd(6); h = fSigmaPixTheta.Project3D("zy"); h->SetDirectory(NULL); h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)"); h->SetXTitle("Id"); h->SetYTitle("\\sigma_{ped}"); h->Draw("box"); h->SetBit(kCanDelete); //pad->cd(7); //fSigmaTheta.Draw(opt); //pad->cd(8); //fDiffPixTheta.Draw(opt); //pad->cd(9); //fSigmaPixTheta.Draw(opt); }