Index: trunk/MagicSoft/Mars/mreflector/MHReflector.cc =================================================================== --- trunk/MagicSoft/Mars/mreflector/MHReflector.cc (revision 8315) +++ trunk/MagicSoft/Mars/mreflector/MHReflector.cc (revision 8315) @@ -0,0 +1,326 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 2/2007 +! +! Copyright: MAGIC Software Development, 2000-2007 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHReflector +// +///////////////////////////////////////////////////////////////////////////// +#include "MHReflector.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MBinning.h" + +#include "MGeomCam.h" + +#include "MRflEvtData.h" +#include "MRflEvtHeader.h" + + +ClassImp(MHReflector); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. Creates and initializes the histograms +// +MHReflector::MHReflector(const char *name, const char *title) +{ + fName = name ? name : "MHReflector"; + fTitle = title ? title : "Histogram for the reflected photons"; + + fHistXY.SetDirectory(NULL); + fHistXY.SetName("ReflXY"); + fHistXY.SetTitle("Histogram vs X/Y and Energy"); + fHistXY.SetXTitle("X [\\circ]"); + fHistXY.SetYTitle("Y [\\circ]"); + fHistXY.SetZTitle("E [GeV]"); + fHistXY.Sumw2(); + + fHistRad.SetDirectory(NULL); + fHistRad.SetName("ReflRad"); + fHistRad.SetTitle("Histogram vs Radius and Energy"); + fHistRad.SetXTitle("E [GeV]"); + fHistRad.SetYTitle("R [\\circ]"); + fHistRad.SetZTitle("Cnts/deg^{2}"); + fHistRad.Sumw2(); + + fHistSize.SetDirectory(NULL); + fHistSize.SetName("ReflSize"); + fHistSize.SetTitle("Histogram vs Size and Energy"); + fHistSize.SetXTitle("E [GeV]"); + fHistSize.SetYTitle("N\\gamma"); + fHistSize.SetZTitle("Cnts"); + fHistSize.Sumw2(); + + MBinning binse, binsd, binsr; + binse.SetEdgesLog(2*5, 10, 1000000); + binsd.SetEdges(50, -2.5, 2.5); + binsr.SetEdges(15, 0, 2.5); + + SetBinning(&fHistXY, &binsd, &binsd, &binse); + SetBinning(&fHistRad, &binse, &binsr); + SetBinning(&fHistSize, &binse, &binse); +} + +// -------------------------------------------------------------------------- +// +// Search for MRflEvtData, MRflEvtHeader and MGeomCam +// +Bool_t MHReflector::SetupFill(const MParList *pList) +{ + fHeader = (MRflEvtHeader*)pList->FindObject("MRflEvtHeader"); + if (!fHeader) + { + *fLog << err << "MRflEvtHeader not found... abort." << endl; + return kFALSE; + } + fData = (MRflEvtData*)pList->FindObject("MRflEvtData"); + if (!fData) + { + *fLog << err << "MRflEvtData not found... abort." << endl; + return kFALSE; + } + MGeomCam *geom = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!geom) + { + *fLog << err << "MGeomCam not found... abort." << endl; + return kFALSE; + } + + fMm2Deg = geom->GetConvMm2Deg(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill position and radius versus energy +// +#include "MRflSinglePhoton.h" +Bool_t MHReflector::Fill(const MParContainer *par, const Stat_t weight) +{ + const Double_t energy = fHeader->GetEnergy(); + + const Int_t n = fData->GetNumPhotons(); + + fHistSize.Fill(energy, n); + + for (int i=0; iGetPhoton(i); + + const Double_t x = ph.GetX()*fMm2Deg; + const Double_t y = ph.GetY()*fMm2Deg; + const Double_t r = TMath::Hypot(x, y); + const Double_t U = r*TMath::TwoPi(); + + if (U==0) + continue; + + fHistRad.Fill(energy, r, weight/U); + //fHistXY.Fill(x, y, energy); + } + + //fData->FillRad(fHistRad, energy/*x*/, fMm2Deg); + //fData->Fill(fHistXY, energy/*z*/, fMm2Deg); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Scale radial histogram with area +// +Bool_t MHReflector::Finalize() +{ + /* + const TAxis &axey = *fHistRad.GetYaxis(); + for (int y=1; y<=fHistRad.GetNbinsY(); y++) + { + const Double_t lo = axey.GetBinLowEdge(y); + const Double_t hi = axey.GetBinLowEdge(y+1); + + const Double_t A = (hi*hi-lo*lo)*TMath::Pi()*TMath::Pi(); + + for (int x=1; x<=fHistRad.GetNbinsX(); x++) + fHistRad.SetBinContent(x, y, fHistRad.GetBinContent(x, y)/A); + } + */ + return kTRUE; +} + +#include "../mmc/MMcEvt.hxx" +void MHReflector::Draw(Option_t *o) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(); + + pad->Divide(3,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetLogx(); + fHistRad.Draw("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + TH1 *h = fHistRad.ProjectionY("ProfRad", -1, 9999, "e"); + h->SetTitle("RadialProfile"); + h->SetDirectory(NULL); + h->SetBit(kCanDelete); + h->SetLineWidth(2); + h->SetLineColor(kBlue); + h->SetStats(kFALSE); + h->Draw(""); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + h = fHistSize.ProjectionY("ProfSize", -1, 9999, "e"); + h->SetTitle("Size distribution"); + h->SetDirectory(NULL); + h->SetBit(kCanDelete); + h->SetLineColor(kBlue); + h->SetLineWidth(2); + h->SetStats(kFALSE); + h->Draw(""); + + pad->cd(5); + gPad->SetBorderMode(0); + + TLegend *leg = new TLegend(0.01, 0.01, 0.99, 0.99, "Energy Range"); + for (int i=1; i<=fHistRad.GetNbinsX(); i++) + { + h = fHistRad.ProjectionY(Form("ProjRad%d", i), i, i, "e"); + h->SetDirectory(NULL); + h->SetBit(kCanDelete); + h->SetLineWidth(2); + h->SetStats(kFALSE); + h->Draw(i==1?"":"same"); + if (i==1) + h->SetTitle("Radial Profile"); + + const Float_t xlo = fHistRad.GetXaxis()->GetBinLowEdge(i); + const Float_t xhi = fHistRad.GetXaxis()->GetBinUpEdge(i); + + const TString lo = MMcEvt::GetEnergyStr(xlo); + const TString hi = MMcEvt::GetEnergyStr(xhi); + + leg->AddEntry(h, Form("%s - %s", lo.Data(), hi.Data())); + } + + pad->cd(4); + gPad->SetBorderMode(0); + leg->SetBit(kCanDelete); + leg->Draw(); + + pad->cd(6); + gPad->SetBorderMode(0); + + gPad->SetLogx(); + gPad->SetLogy(); + + for (int i=1; i<=fHistSize.GetNbinsX(); i++) + { + h = fHistSize.ProjectionY(Form("ProjSize%d", i), i, i, "e"); + h->SetDirectory(NULL); + h->SetBit(kCanDelete); + h->SetLineWidth(2); + h->SetStats(kFALSE); + h->Draw(i==1?"":"same"); + if (i==1) + h->SetTitle("Size distribution"); + } +} + +void MHReflector::Paint(Option_t *opt) +{ + TVirtualPad *pad = gPad; + + pad->cd(1); + SetPalette("pretty"); + + TH1 *h=0; + + pad->cd(2); + gPad->SetBorderMode(0); + if (gPad->FindObject("ProfRad")) + { + h = fHistRad.ProjectionY("ProfRad", -1, 9999, "e"); + h->Scale(1./h->Integral()); + } + + pad->cd(3); + gPad->SetBorderMode(0); + if (gPad->FindObject("ProfSize")) + { + h = fHistSize.ProjectionY("ProfSize", -1, 9999, "e"); + h->Scale(1./h->Integral()); + } + + pad->cd(5); + + for (int i=1; i<=fHistRad.GetNbinsX(); i++) + { + const TString name = Form("ProjRad%d", i); + if (!gPad->FindObject(name)) + continue; + + h = fHistRad.ProjectionY(name, i, i, "e"); + h->SetLineColor(i); + h->Scale(1./fHistRad.Integral()); + if (i==1) + h->SetMaximum(fHistRad.GetMaximum()/fHistRad.Integral()*1.05); + } + + pad->cd(6); + + for (int i=1; i<=fHistRad.GetNbinsX(); i++) + { + const TString name = Form("ProjSize%d", i); + if (!gPad->FindObject(name)) + continue; + + h = fHistSize.ProjectionY(name, i, i, "e"); + h->SetLineColor(i); + h->Scale(1./fHistSize.Integral()); + if (i==1) + { + h->SetMaximum(TMath::Power(fHistSize.GetMaximum()/fHistSize.Integral(), 1.05)); + h->SetMinimum(TMath::Power(fHistSize.GetMinimum(0)/fHistSize.Integral(), 0.95)); + } + } + +} Index: trunk/MagicSoft/Mars/mreflector/MHReflector.h =================================================================== --- trunk/MagicSoft/Mars/mreflector/MHReflector.h (revision 8315) +++ trunk/MagicSoft/Mars/mreflector/MHReflector.h (revision 8315) @@ -0,0 +1,47 @@ +#ifndef MARS_MHReflector +#define MARS_MHReflector + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include +#endif + +#ifndef ROOT_TH3 +#include +#endif + +class MRflEvtData; +class MRflEvtHeader; + +class MHReflector : public MH +{ +private: + MRflEvtData *fData; //! + MRflEvtHeader *fHeader; //! + + Double_t fMm2Deg; + + TH3D fHistXY; + TH2D fHistRad; + TH2D fHistSize; + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t weight=1); + Bool_t Finalize(); + +public: + MHReflector(const char *name=0, const char *title=0); + + TH2D *GetHistRad() { return &fHistRad; } + TH3D *GetHistXY() { return &fHistXY; } + + void Draw(Option_t *o=""); + void Paint(Option_t *o=""); + + ClassDef(MHReflector, 1) +}; + +#endif Index: trunk/MagicSoft/Mars/mreflector/Makefile =================================================================== --- trunk/MagicSoft/Mars/mreflector/Makefile (revision 8290) +++ trunk/MagicSoft/Mars/mreflector/Makefile (revision 8315) @@ -19,10 +19,11 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase -I../mgui -I../mgeom +INCLUDES = -I. -I../mbase -I../mgui -I../mgeom -I../mhbase SRCFILES = MRflEvtData.cc \ MRflEvtHeader.cc \ MRflRunHeader.cc \ - MRflSinglePhoton.cc + MRflSinglePhoton.cc \ + MHReflector.cc ############################################################ Index: trunk/MagicSoft/Mars/mreflector/ReflectorLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mreflector/ReflectorLinkDef.h (revision 8290) +++ trunk/MagicSoft/Mars/mreflector/ReflectorLinkDef.h (revision 8315) @@ -10,3 +10,5 @@ #pragma link C++ class MRflSinglePhoton+; +#pragma link C++ class MHReflector+; + #endif