Index: trunk/MagicSoft/Mars/mhist/MHGamma.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHGamma.cc (revision 1413) +++ trunk/MagicSoft/Mars/mhist/MHGamma.cc (revision 1413) @@ -0,0 +1,271 @@ +/* ======================================================================== *\ +! +! * +! * 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 time saving 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 6/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHGamma // +// // +// manipulation of alpha distributions // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHGamma.h" + +#include +#include + +#include + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHGamma); + + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHGamma::MHGamma() +{ +} + +// -------------------------------------------------------------------------- +// +// Dummy Fill +// +Bool_t MHGamma::Fill(const MParContainer *par) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the histogram as the difference of two histograms : +// fHist(gamma) = h1(source) - h2(antisource) +// +TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2, + const char *name, const char *title, + Bool_t Draw) +{ + TH3D *fHist; + fHist = new TH3D(); + fHist->SetName(name); + fHist->SetTitle(title); + + fHist->SetDirectory(NULL); + + SetBinning((TH3D*)fHist, (TH3D*)h1); + + TString strg1 = (((TH1*)h1)->GetXaxis())->GetTitle(); + TString strg2 = (((TH1*)h1)->GetYaxis())->GetTitle(); + TString strg3 = (((TH1*)h1)->GetZaxis())->GetTitle(); + fHist->SetXTitle(strg1); + fHist->SetYTitle(strg2); + fHist->SetZTitle(strg3); + + fHist->Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME! + + //........................................................... + // draw histogram + if (Draw == kTRUE) + { + TString strg7 = "3D-plot of "; + strg7 += strg1; + strg7 += ","; + strg7 += strg2; + strg7 += ","; + strg7 += strg3; + strg7 += " for \'gamma\' sample"; + + TCanvas &c = *MakeDefCanvas("Alpha", strg7); + + c.Divide(2, 2); + + gROOT->SetSelectedPad(NULL); + + TH1 *h; + + c.cd(1); + h = ((TH3D*)(fHist))->Project3D("ex"); + + TString strg0 = "SRC-ASRC : "; + TString strg4 = strg0 + strg1; + h->SetTitle(strg4); + h->SetXTitle(strg1); + h->SetYTitle("Counts"); + + h->Draw(); + h->SetBit(kCanDelete); + + c.cd(2); + h = ((TH3D*)(fHist))->Project3D("ey"); + + TString strg5 = strg0 + strg2; + h->SetTitle(strg5); + h->SetXTitle(strg2); + h->SetYTitle("Counts"); + + h->Draw(); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + c.cd(3); + h = ((TH3D*)(fHist))->Project3D("ez"); + + TString strg6 = strg0 + strg3; + h->SetTitle(strg6); + h->SetXTitle(strg3); + h->SetYTitle("Counts"); + + h->Draw(); + h->SetBit(kCanDelete); + + c.cd(4); + ((TH3D*)fHist)->DrawCopy(); + + c.Modified(); + c.Update(); + } + + return fHist; +} + +// -------------------------------------------------------------------------- +// +// Integrate fHist(gamma) in the alpha range (lo, up) +// +TH2D *MHGamma::GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up, + Bool_t Drawp) +{ + if (up < lo) + { + *fLog << err << fName << "MHGamma : Alpha projection not possible: lo=" << lo << " up=" << up << endl; + return NULL; + } + + TAxis &axe = *fHist->GetXaxis(); + + Int_t ilo = axe.FindFixBin(lo); + Int_t iup = axe.FindFixBin(up); + + const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo); + const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo; + + const Double_t epsup1 = up-axe.GetBinLowEdge(iup); + const Double_t epsup2 = axe.GetBinUpEdge(iup)-up; + + const Double_t epslo = epslo1epslo2) + ilo++; + + if (epsup10.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup)) + { + *fLog << err << fName << "MHGamma : binning is not adequate for the requested projection:" << endl; + *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl; + *fLog << " epslo = " << epslo << endl; + *fLog << " epsup = " << epsup << endl; + *fLog << " dwl = " << axe.GetBinWidth(ilo) << endl; + *fLog << " dwu = " << axe.GetBinWidth(iup) << endl; + return NULL; + } + + axe.SetRange(ilo, iup); + + TH2D &h2D = *(TH2D *)fHist->Project3D("ezy"); + + TString strg0 = "2D-plot of "; + TString strg1 = (fHist->GetYaxis())->GetTitle(); + TString strg2 = (fHist->GetZaxis())->GetTitle(); + strg0 += strg2; + strg0 += " vs. "; + strg0 += strg1; + h2D.SetTitle(strg0); + h2D.SetXTitle(strg1); + h2D.SetYTitle(strg2); + + + //........................................................... + // draw histogram + if (Drawp == kTRUE) + { + char txt[100]; + TString strg3 = "No.of Gammas vs. "; + strg3 += strg1; + strg3 += " and "; + strg3 += strg2; + sprintf(txt, " (%.1f < alpha < %.1f deg)", lo, up); + strg3 += txt; + + TCanvas &c = *MakeDefCanvas("Gamma", strg3); + + c.Divide(2, 2); + + gROOT->SetSelectedPad(NULL); + + TH1 *h; + + c.cd(1); + h = h2D.ProjectionX("xpro", -1, 9999, "E"); + TString strg0 = "No.of gammas : "; + TString strg7 = strg0 + strg1; + h->SetTitle(strg7); + h->SetXTitle(strg1); + h->SetYTitle("Counts"); + + h->Draw(); + h->SetBit(kCanDelete); + gPad->SetLogx(); + + c.cd(2); + h = h2D.ProjectionY("ypro", -1, 9999, "E"); + TString strg8 = strg0 + strg2; + h->SetTitle(strg8); + h->SetXTitle(strg2); + h->SetYTitle("Counts"); + + h->Draw(); + h->SetBit(kCanDelete); + + c.cd(3); + + h2D.DrawCopy(); + gPad->SetLogx(); + + c.Modified(); + c.Update(); + } + //........................................................... + + return &h2D; +} Index: trunk/MagicSoft/Mars/mhist/MHGamma.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHGamma.h (revision 1413) +++ trunk/MagicSoft/Mars/mhist/MHGamma.h (revision 1413) @@ -0,0 +1,42 @@ +#ifndef MARS_MHGamma +#define MARS_MHGamma + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH3 +#include "TH3.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +class TH2D; +class TH3D; + +class MHGamma : public MH +{ +public: + MHGamma(); + + Bool_t Fill(const MParContainer *par); + + TH3D *Subtract(const TH3D *h1, const TH3D *h2, + const char *name, const char *title, Bool_t Draw=kFALSE); + + TH2D *GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up, + Bool_t Drawp=kFALSE); + + + ClassDef(MHGamma, 1) // manipulation of alpha distributions +}; + +#endif + + + + + +