/* ======================================================================== *\ ! ! * ! * 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 (UInt_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(); } // --------------------------------------------------------------------------