Index: trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 6977) +++ trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 7001) @@ -107,11 +107,17 @@ TVirtualPad *save = gPad; gPad = 0; + /* #if ROOT_VERSION_CODE < ROOT_VERSION(3,01,06) SetPalette(1, 0); -#else +#endif + */ +/* +#if ROOT_VERSION_CODE < ROOT_VERSION(4,04,00) SetPrettyPalette(); +#elese // WORAROUND - FIXME: Calling it many times becomes slower and slower - //SetInvDeepBlueSeaPalette(); + SetInvDeepBlueSeaPalette(); #endif +*/ gPad = save; } @@ -121,5 +127,5 @@ // Default Constructor. To be used by the root system ONLY. // -MHCamera::MHCamera() : TH1D(), fGeomCam(NULL), fColors(kItemsLegend), fAbberation(0) +MHCamera::MHCamera() : TH1D(), fGeomCam(NULL), fAbberation(0) { Init(); @@ -133,5 +139,5 @@ // MHCamera::MHCamera(const MGeomCam &geom, const char *name, const char *title) -: fGeomCam(NULL), fColors(kItemsLegend), fAbberation(0) +: fGeomCam(NULL), fAbberation(0) { //fGeomCam = (MGeomCam*)geom.Clone(); @@ -463,4 +469,7 @@ // GeMinimum() ((val-min)/(max-min)) // 'proj' Display the y-projection of the histogram +// 'pal0' Use Pretty palette +// 'pal1' Use Deep Blue Sea palette +// 'pal2' Use Inverse Depp Blue Sea palette // 'same' Draw trandparent pixels on top of an existing pad. This // makes it possible to draw the camera image on top of an @@ -970,4 +979,12 @@ { opt.ReplaceAll("hist", ""); + opt.ReplaceAll("box", ""); + opt.ReplaceAll("pixelindex", ""); + opt.ReplaceAll("sectorindex", ""); + opt.ReplaceAll("content", ""); + opt.ReplaceAll("proj", ""); + opt.ReplaceAll("pal0", ""); + opt.ReplaceAll("pal1", ""); + opt.ReplaceAll("pal2", ""); TH1D::Paint(opt); return; @@ -1007,4 +1024,16 @@ } + const Bool_t pal1 = opt.Contains("pal1"); + const Bool_t pal2 = opt.Contains("pal2"); + + if (!pal1 && !pal2) + SetPrettyPalette(); + + if (pal1) + SetDeepBlueSeaPalette(); + + if (pal2) + SetInvDeepBlueSeaPalette(); + // Update Contents of the pixels and paint legend Update(gPad->GetLogy(), hasbox, hascol, hassame); @@ -1021,4 +1050,23 @@ } +void MHCamera::SetDrawOption(Option_t *option) +{ + // This is a workaround. For some reason MHCamera is + // stored in a TObjLink instead of a TObjOptLink + if (!option || !gPad) + return; + + TListIter next(gPad->GetListOfPrimitives()); + delete gPad->FindObject("Tframe"); + TObject *obj; + while ((obj = next())) + if (obj == this && (TString)next.GetOption()!=(TString)option) + { + gPad->GetListOfPrimitives()->Remove(this); + gPad->GetListOfPrimitives()->AddFirst(this, option); + return; + } +} + // ------------------------------------------------------------------------ // @@ -1053,8 +1101,4 @@ else gStyle->SetPalette(ncolors, colors); - - fColors.Set(kItemsLegend); - for (int i=0; iGetColorPalette(i); } @@ -1068,6 +1112,13 @@ void MHCamera::SetPrettyPalette() { - if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase)) + TString opt(GetDrawOption()); + + if (!opt.Contains("hist", TString::kIgnoreCase)) SetPalette(1, 0); + + opt.ReplaceAll("pal1", ""); + opt.ReplaceAll("pal2", ""); + + SetDrawOption(opt); } @@ -1080,6 +1131,14 @@ void MHCamera::SetDeepBlueSeaPalette() { - if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase)) + TString opt(GetDrawOption()); + + if (!opt.Contains("hist", TString::kIgnoreCase)) SetPalette(51, 0); + + opt.ReplaceAll("pal1", ""); + opt.ReplaceAll("pal2", ""); + opt += "pal1"; + + SetDrawOption(opt); } @@ -1092,6 +1151,14 @@ void MHCamera::SetInvDeepBlueSeaPalette() { - if (!TString(GetDrawOption()).Contains("hist", TString::kIgnoreCase)) + TString opt(GetDrawOption()); + + if (!opt.Contains("hist", TString::kIgnoreCase)) SetPalette(52, 0); + + opt.ReplaceAll("pal1", ""); + opt.ReplaceAll("pal2", ""); + opt += "pal2"; + + SetDrawOption(opt); } @@ -1515,8 +1582,8 @@ if (val >= max) - return fColors[maxcolidx]; + return gStyle->GetColorPalette(maxcolidx); if (val <= min) - return fColors[0]; + return gStyle->GetColorPalette(0); // @@ -1530,5 +1597,5 @@ const Int_t colidx = (Int_t)(ratio*maxcolidx + .5); - return fColors[colidx]; + return gStyle->GetColorPalette(colidx); } @@ -1612,5 +1679,5 @@ for (Int_t i=0; iGetColorPalette(i)); newbox.PaintBox(range, H*(i*h-1)-offset, range+w, H*((i+1)*h-1)-offset); } Index: trunk/MagicSoft/Mars/mhist/MHCamera.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.h (revision 6977) +++ trunk/MagicSoft/Mars/mhist/MHCamera.h (revision 7001) @@ -48,5 +48,4 @@ TArrayI fBinEntries; // number of entries per bin - TArrayI fColors; //! Color conversion table TList *fNotify; //! @@ -214,4 +213,5 @@ char *GetObjectInfo(Int_t px, Int_t py) const; void ExecuteEvent(Int_t event, Int_t px, Int_t py); + void SetDrawOption(Option_t *option); //*MENU* void SetPalette(Int_t ncolors, Int_t *colors); @@ -222,6 +222,4 @@ void SetAutoScale() { fMinimum = fMaximum = -1111; } // *MENU* - void DisplayAsHistogram() { SetDrawOption("histEP"); } // *MENU* - void DisplayAsCamera() { SetDrawOption(""); } // *MENU* void SetFreezed(Bool_t f=kTRUE) { f ? SetBit(kFreezed) : ResetBit(kFreezed); } // *TOGGLE* *GETTER=IsFreezed Index: trunk/MagicSoft/Mars/mhist/MHGamma.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHGamma.cc (revision 6977) +++ (revision ) @@ -1,286 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MHAlphaEnergyTheta.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHGamma); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MHGamma::MHGamma(const TString &name, const TString &title) - : fHist(NULL), fProject(NULL) -{ - fName = name.IsNull() ? (TString)"MHGamma" : name; - fTitle = title.IsNull() ? (TString)"3D-histogram of Alpha, E_est, Theta (Gamma sample)" : title; -} - -TH3D *MHGamma::Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2) -{ - return Subtract(h1.GetHist(), h2.GetHist()); -} - -TObject *MHGamma::DrawClone(Option_t *opt) const -{ - DrawClone1(); - DrawClone2(); - - return NULL; -} - -void MHGamma::DrawClone1() const -{ - if (!fHist) - return; - - // - // ------------- Part 1 --------------------- - // - TString titlex = fHist->GetXaxis()->GetTitle(); - TString titley = fHist->GetYaxis()->GetTitle(); - TString titlez = fHist->GetYaxis()->GetTitle(); - - TString canvtitle = "3D-plot "; //of "; - /* - canvtitle += titlex; - canvtitle += ", "; - canvtitle += titley; - canvtitle += ", "; - canvtitle += titlez+" "; - */ - canvtitle += "for 'gamma' sample"; - - TCanvas &c = *MakeDefCanvas("Alpha", canvtitle); - - c.Divide(2, 2); - - gROOT->SetSelectedPad(NULL); - - TH1 *h; - - c.cd(1); - h = ((TH3D*)(fHist))->Project3D(fName+"_ex"); - - TString title= "Source-Antisource: "; - h->SetTitle(title + titlex); - h->SetXTitle(titlex); - h->SetYTitle("Counts"); - - h->Draw(); - h->SetBit(kCanDelete); - - c.cd(2); - h = ((TH3D*)(fHist))->Project3D(fName+"_ey"); - - h->SetTitle(title + titley); - h->SetXTitle(titley); - h->SetYTitle("Counts"); - - h->Draw(); - h->SetBit(kCanDelete); - gPad->SetLogx(); - - c.cd(3); - h = ((TH3D*)(fHist))->Project3D(fName+"_ez"); - - h->SetTitle(title + titlez); - h->SetXTitle(titlez); - h->SetYTitle("Counts"); - - h->Draw(); - h->SetBit(kCanDelete); - - c.cd(4); - ((TH3D*)fHist)->DrawCopy(); - - c.Modified(); - c.Update(); -} - - -// -------------------------------------------------------------------------- -// -// Calculate the histogram as the difference of two histograms : -// fHist(gamma) = h1(source) - h2(antisource) -// -TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2) -{ - if (fHist) - delete fHist; - - fHist = new TH3D; - fHist->SetName(fName); - fHist->SetTitle(fTitle); - fHist->SetDirectory(NULL); - - SetBinning((TH1*)fHist, (TH1*)h1); - - fHist->SetXTitle((((TH1*)h1)->GetXaxis())->GetTitle()); - fHist->SetYTitle((((TH1*)h1)->GetYaxis())->GetTitle()); - fHist->SetZTitle((((TH1*)h1)->GetZaxis())->GetTitle()); - - fHist->Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME! - - return fHist; -} - -// -------------------------------------------------------------------------- -// -// Integrate fHist(gamma) in the alpha range (lo, up) -// -TH2D *MHGamma::GetAlphaProjection(Axis_t lo, Axis_t up) -{ - 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); - - fLo = lo; - fHi = up; - - if (fProject) - delete fProject; - fProject = (TH2D*)fHist->Project3D(fName+"_ezy"); - - const TString title = "2D-plot of "; - const TString titley = fHist->GetYaxis()->GetTitle(); - const TString titlez = fHist->GetZaxis()->GetTitle(); - - fProject->SetTitle(title + titley + " vs. " + titlez); - fProject->SetXTitle(titley); - fProject->SetYTitle(titlez); - - return fProject; -} - -void MHGamma::DrawClone2() const -{ - if (!fProject) - return; - - const TString titley = fHist->GetYaxis()->GetTitle(); - const TString titlez = fHist->GetZaxis()->GetTitle(); - - TString canvtitle = "No.of Gammas ";//vs. "; - /* - canvtitle += titley; - canvtitle += " and "; - canvtitle += titlez; - */ - canvtitle += Form("(%.1f < alpha < %.1f deg)", fLo, fHi); - - TCanvas &c = *MakeDefCanvas("Gamma", canvtitle); - - c.Divide(2, 2); - - gROOT->SetSelectedPad(NULL); - - TH1 *h; - - c.cd(1); - h = fProject->ProjectionX(fName+"_xpro", -1, 9999, "E"); - TString title = "No.of gammas: "; - h->SetTitle(title+titley); - h->SetXTitle(titley); - h->SetYTitle("Counts"); - - h->Draw(); - h->SetBit(kCanDelete); - gPad->SetLogx(); - - c.cd(2); - h = fProject->ProjectionY(fName+"_ypro", -1, 9999, "E"); - h->SetTitle(title+titlez); - h->SetXTitle(titlez); - h->SetYTitle("Counts"); - - h->Draw(); - h->SetBit(kCanDelete); - - c.cd(3); - - fProject->DrawCopy(); - gPad->SetLogx(); - - c.Modified(); - c.Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHGamma.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHGamma.h (revision 6977) +++ (revision ) @@ -1,46 +1,0 @@ -#ifndef MARS_MHGamma -#define MARS_MHGamma - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH2D; -class TH3D; - -class MHAlphaEnergyTheta; - -class MHGamma : public MH -{ -private: - TH3D *fHist; //! - TH2D *fProject; //! - - Axis_t fLo; //! - Axis_t fHi; //! - -public: - MHGamma(const TString &name="", const TString &title=""); - - TH3D *Subtract(const TH3D *h1, const TH3D *h2); - - TH3D *Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2); - - TH2D *GetAlphaProjection(Axis_t lo, Axis_t up); - - TObject *DrawClone(Option_t *opt="") const; - void DrawClone1() const; - void DrawClone2() const; - - const TH2D *GetProject() const { return fProject; } - - ClassDef(MHGamma, 0) // manipulation of alpha distributions -}; - -#endif - - - - - - Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 6977) +++ (revision ) @@ -1,1767 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 expressed -! * or implied warranty. -! * -! -! -! Author(s): Robert Wagner 9/2002 -! Author(s): Robert Wagner 3/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHOnSubtraction // -// // -// // -// extracts the gamma signal from a pure ON-signal given in an // -// ALPHA-Energy-Theta histogram. The class will take this histogram from // -// the parameter list and will provide result histograms in the // -// parameter list. // -// // -// This class still is work in progress. // -// // -// // -////////////////////////////////////////////////////////////////////////////// - -// This part of MARS is code still evolving. Please do not change the code -// without prior feedback by the author. - -#include "MHOnSubtraction.h" - -#include -#include -#include -#include -#include -#include -#include - -#include "MBinning.h" -#include "MParList.h" -#include "MHArray.h" -#include "MH3.h" - -#include "MHAlphaEnergyTheta.h" -#include "MHAlphaEnergyTime.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHOnSubtraction); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0), fMaxRedChiSq(0), fSlope(20.0) -{ - fName = name ? name : "MHOnSubtraction"; - fTitle = title ? title : "Extracts Gamma signal from pure ON data"; - fHistogramType = "Theta"; - fChiSquareHisto=0; - fSignificanceHisto=0; - fSummedAlphaPlots=0; - fThetaBin=0; - -} - - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MHOnSubtraction::~MHOnSubtraction() -{ - if (fChiSquareHisto) delete fChiSquareHisto; - if (fSignificanceHisto) delete fSignificanceHisto; - if (fSummedAlphaPlots) delete fSummedAlphaPlots; -} - -// ----------------------------------------------------------------------- -// -// Calculate Significance according to Li and Ma, ApJ 272 (1983) 317, eq -// (17). n_{On}, n_{Off} and the ratio of On-times for On and Off -// measurements have to be provided. -// -// This function underestimates the true significance for it does not take -// into account errors on the event numbers. A more accurate variation wil -// be provided soon -// -Double_t MHOnSubtraction::CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta) -{ - if (nOn<=0) - *fLog << warn << "Got " << nOn << " total events, " << flush; - if (nOff<=0) - *fLog << warn << "Got " << nOff << " background events, " << flush; - if (nOn<=0 || nOff<=0.0) { - *fLog << warn << "returning significance of 0.0" << endl; - return 0.0; - } - return sqrt(2*(nOn*log((1+theta)*nOn/(theta*(nOn+nOff))) - +nOff*log((1+theta)*(nOff/(nOn+nOff))))); -} - -// ----------------------------------------------------------------------- -// -// This function takes a first look at a given ALPHA distribution -// and determines if a fit is applicable. -// -// Fits a given TH1 containing an ALPHA distribution with a combined -// gaussian plus polynomial of third grade and returns the fitted -// function. By signalRegionFactor the width of the region in which -// signal extraction is to be performed can be specified in units of -// sigma of the gaussian which is fitted to the distribution (default: -// 3.5). Accordingly, FitHistogram returns the range in which it suggests -// signal extraction (lowerBin < ALPHA < upperBin). An estimated -// significance of the signal is also provided. draw specifies whether -// FitHistogram should draw the distribution. -// -Bool_t MHOnSubtraction::FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa, - Double_t &lowerBin, Double_t &upperBin, - Float_t signalRegionFactor, const Bool_t draw, - TString funcName) -{ - if (alphaHisto.GetEntries() == 0) { - *fLog << warn << "Histogram contains no entries. Returning. " << endl; - lowerBin=0; - upperBin=0; - sigLiMa=0; - if (draw) { - TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no events)"); - lab->SetFillStyle(0); - lab->SetBorderSize(0); - lab->Draw(); - } - return kFALSE; - } - - //cosmetics - alphaHisto.GetXaxis()->SetTitle("ALPHA"); - alphaHisto.GetYaxis()->SetTitle("events"); - - Double_t outerEvents = 0.0; - Double_t innerEvents = 0.0; - Int_t outerBins = 0; - Int_t innerBins = 0; - Int_t outerBinsZero = 0; - Int_t innerBinsZero = 0; - - for (Int_t alphaBin = 1; alphaBin < alphaHisto.GetNbinsX(); alphaBin++) { - if ((alphaHisto.GetBinLowEdge(alphaBin) >= -35.)&& //inner Region - (alphaHisto.GetBinLowEdge(alphaBin+1) <= 35.0)) { - innerEvents+=alphaHisto.GetBinContent(alphaBin); - innerBins++; - if (alphaHisto.GetBinContent(alphaBin)==0.0) - innerBinsZero++; - } else { - if ((alphaHisto.GetBinLowEdge(alphaBin) >= -89.5)&& //outer Region - (alphaHisto.GetBinLowEdge(alphaBin+1) <=89.5)) { - outerEvents+=alphaHisto.GetBinContent(alphaBin); - outerBins++; - if (alphaHisto.GetBinContent(alphaBin)==0.0) - outerBinsZero++; - } - } - } - - *fLog << dbg << "Plot contains " << - outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " << - innerEvents << " inner ev (" << innerEvents/innerBins << "/bin) " << endl; - - if ((outerBinsZero!=0) || (innerBinsZero != 0)) - *fLog << warn << "There are "; - if (outerBinsZero != 0) - *fLog << dbg << outerBinsZero << " empty outer bins "; - if (innerBinsZero != 0) - *fLog << dbg <SetLineWidth(2); - po->SetLineColor(2); - po->SetParNames("c"); - po->SetParameter(0,0.0); - alphaHisto.GetListOfFunctions()->Add(po); - alphaHisto.SetLineColor(97); - if (draw) { - alphaHisto.DrawCopy(); //rwagner - po->Draw("SAME"); - } - sigLiMa = 0; - lowerBin = 1; - upperBin = alphaHisto.GetNbinsX()-1; - signalRegionFactor = 0; - - return kFALSE; //No gaus fit applied - } - - if (outerEvents/outerBins < 2.5) { - *fLog << warn << "Only " << /*setprecision(2) <<*/ outerEvents/outerBins - << " events/bin in OFF region: " - << "Assuming this as background." << endl; - - TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5); - po->SetLineWidth(2); - po->SetLineColor(94); - po->SetParNames("c"); - po->SetParameter(0,outerEvents/outerBins); - alphaHisto.GetListOfFunctions()->Add(po); - if (draw) { - alphaHisto.DrawCopy(); //rwagner - po->Draw("SAME"); - } - - Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0); - Int_t maxBin = centerBin > alphaHisto.GetNbinsX() - centerBin ? - alphaHisto.GetNbinsX()- centerBin : centerBin; - - Int_t lowerSignalEdge = centerBin; - for (Int_t i=3; i < maxBin; i++) { - if ((Double_t)alphaHisto.GetBinContent(centerBin-i) - < outerEvents/outerBins) { - lowerSignalEdge = centerBin-i+1; - break; - } - } - if (centerBinupperSignalEdge) upperSignalEdge=centerBin; - - Double_t nOnInt = 0; - for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++) - nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins; - - Double_t nOffInt = (upperSignalEdge - lowerSignalEdge + 1) - * (outerEvents/outerBins); - - sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); - - if (sigLiMa < 3) - alphaHisto.SetLineColor(2); - - - *fLog << inf << "Estimated significance is " << sigLiMa - << " sigma " << endl; - *fLog << inf << "Signal region is " - << lowerBin << " < ALPHA < " << upperBin - << " (Most likely you wanna ignore this)" << endl; - - return kFALSE; //No gaus fit applied - } - - //fit combined gaus+pol3 to data - TF1 *gp = new TF1("gauspol3"+funcName,"gaus(0)+pol3(3)",-89.5,89.5); - gp->SetLineWidth(2); - gp->SetLineColor(2); - gp->SetParNames("Excess","A","sigma","a","b","c","d"); - gp->SetParLimits(0, 200,2000); - gp->SetParLimits(1, -4.,4.); - gp->SetParLimits(2, 2.,20.); - // gp->SetParameter(6, 0.0000); // include for slope(0)=0 constrain - TString gpDrawOptions = draw ? "RIQ" : "RIQ0"; - gpDrawOptions = "RIQ0"; // rwagner - alphaHisto.Fit("gauspol3"+funcName,gpDrawOptions); - alphaHisto.DrawCopy(); //rwagner - alphaHisto.SetDirectory(NULL); //rwagner - - Double_t gausMean = gp->GetParameter(1); - Double_t gausSigma = gp->GetParameter(2); - -// TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-signalRegionFactor*gausSigma+gausMean, -// signalRegionFactor*gausSigma+gausMean); - TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-100, - 100); - p3->SetParameters(gp->GetParameter(3),gp->GetParameter(4), - gp->GetParameter(5),gp->GetParameter(6)); - // p3->SetLineStyle(2); - p3->SetLineColor(4); - p3->SetLineWidth(1); - if (draw) p3->Draw("SAME"); - - TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40); - ga->SetParameters(gp->GetParameter(0),gp->GetParameter(1),gp->GetParameter(2)); - ga->SetLineColor(93); - ga->SetLineWidth(2); - // if (draw) ga->Draw("SAME"); - - // identify the signal region: signalRegionFactor*gausSigma - // this allows us to - // (1) determine n_{ON}, n_{OFF} - Double_t scalingFactor = - (alphaHisto.GetXaxis()->GetBinLowEdge(alphaHisto.GetNbinsX()+1) - - alphaHisto.GetXaxis()->GetBinLowEdge(1)) / alphaHisto.GetNbinsX(); - - Double_t nOnInt = (gp->Integral(-signalRegionFactor*gausSigma+gausMean, - signalRegionFactor*gausSigma+gausMean) / scalingFactor); - Double_t nOffInt = (p3->Integral(-signalRegionFactor*gausSigma+gausMean, - signalRegionFactor*gausSigma+gausMean) / scalingFactor); - - // (2) determine the signal region from fit in degrees - // we do it a bit more complicated: assuming that the binning in all - // histograms is the same, we want to be sure that summing up is always - // done over the same bins. - - lowerBin = alphaHisto.GetXaxis()->FindBin(-signalRegionFactor*gausSigma+gausMean); - upperBin = alphaHisto.GetXaxis()->FindBin( signalRegionFactor*gausSigma+gausMean); - - lowerBin = alphaHisto.GetBinLowEdge((Int_t)lowerBin); - upperBin = alphaHisto.GetBinLowEdge((Int_t)upperBin)+alphaHisto.GetBinWidth((Int_t)upperBin); - - sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); -// if (sigLiMa < 3) -// alphaHisto.SetLineColor(2); - - - *fLog << inf << "Fit estimates significance to be " - << sigLiMa << " sigma " << endl; - - *fLog << inf << "Fit yields signal region to be " - << lowerBin << " < ALPHA < " << upperBin - << " (Chisquare/dof=" << gp->GetChisquare()/gp->GetNDF() - << ", prob=" << gp->GetProb() << ")" << endl; - - return kTRUE; //returning gaus fit -} - - - - - -// ----------------------------------------------------------------------- -// -// Does the actual extraction of the gamma signal. For performance -// reasons, fits already done by MHOnSubtraction::FitHistogram are used -// and not redone. -// From it, a polynomial function for the background is evaluated and the -// gamma signal is extracted in the region given by lowerBin < ALPHA < -// upperBin. -// Significance of the signal is also provided. draw specifies whether -// FitHistogram should draw the distribution. -// -Bool_t MHOnSubtraction::ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa, - Double_t &lowerBin, Double_t &upperBin, - Double_t &gammaSignal, Double_t &errorGammaSignal, - Double_t &off, Double_t &errorOff, - Float_t signalRegionFactor, const Bool_t draw, TString funcName, - TPad *drawPad, Int_t drawBase) -{ - TF1 *gausPol = alphaHisto.GetFunction("gauspol3"+funcName); - TF1 *pol = alphaHisto.GetFunction("pol0"+funcName); - - if (!gausPol && !pol) { - *fLog << err << "Fatal: ALPHA histogram has no gauspol or pol " - << " fit attached to it." << endl; - TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no fit)"); - lab->SetFillStyle(3000); - lab->SetBorderSize(0); - lab->DrawClone(); - lab->SetBit(kCanDelete); - return kFALSE; - } - - TF1* po; - po = NULL; - - if (gausPol) { - Double_t gausMean = gausPol->GetParameter(1); - Double_t gausSigma = gausPol->GetParameter(2); - po = new TF1("po"+funcName,"pol3(0)", - -signalRegionFactor*gausSigma+gausMean, - signalRegionFactor*gausSigma+gausMean); - po->SetParameters(gausPol->GetParameter(3),gausPol->GetParameter(4), - gausPol->GetParameter(5),gausPol->GetParameter(6)); - - TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40); - ga->SetParameters(gausPol->GetParameter(0),gausPol->GetParameter(1), - gausPol->GetParameter(2)); - - if (draw) { - alphaHisto.Draw(); - gausPol->Draw("SAME"); //Maybe not even necessary? - - po->SetLineColor(4); - po->SetLineWidth(2); - po->Draw("SAME"); - - ga->SetLineColor(93); - ga->SetLineWidth(2); - ga->Draw("SAME"); - - char legendTitle[80]; - sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f", - lowerBin, upperBin); - - TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle); - - legend->SetBorderSize(0); - legend->SetFillColor(10); - legend->SetFillStyle(0); - legend->AddEntry(gausPol, "combined N_{on}","l"); - legend->AddEntry(po,"polynomial N_{bg} Signal region","l"); - legend->AddEntry(ga, "putative gaussian N_{S}","l"); - legend->Draw(); - } - } // gausPol - - - if (pol) { - po = pol; - - if (draw) { - alphaHisto.Draw(); - - po->SetLineColor(6); - po->SetLineWidth(2); - po->Draw("SAME"); - - char legendTitle[80]; - sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f", - lowerBin, upperBin); - - TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle); - - legend->SetBorderSize(0); - legend->SetFillColor(10); - legend->SetFillStyle(0); - legend->AddEntry(po,"polynomial N_{bg} Signal region","l"); - legend->Draw(); - } - } // pol - - Double_t nOn = 0; - Double_t eNOn = 0; - - Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin); - Int_t binNumberHigh = alphaHisto.GetXaxis()->FindBin(upperBin); - - for (Int_t bin=binNumberLow; binEval(x); - //cout << bin << ": " << offEvts << endl; - nOff += offEvts; - } //for bin - eNOff = sqrt(fabs(nOff)); - - if (nOn==0) // there should not be a negative number of signal events - nOff=0; - - if (nOff<0) { // there should not be a negative number of off events - nOff=0; - eNOff=0; - } - - *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", "; - - off = nOff; errorOff = eNOff; - gammaSignal = nOn-nOff; errorGammaSignal = sqrt(eNOn*eNOn + eNOff*eNOff); - - *fLog << inf << "nBack = " << nOff << "+-" << eNOff << ", "; - *fLog << inf << "nSig = " << gammaSignal << "+-" << errorGammaSignal << endl; - - sigLiMa = CalcSignificance(nOn, nOff, 1); - // Double_t sigFit=(ga->GetParameter(1))/(ga->GetParError(1)); //Mean / sigMean - - *fLog << inf << "Significance: "<AddText(tx); - sprintf(tx, "Off: %2.2f #pm %2.2f", nOff, eNOff); - pt->AddText(tx); - sprintf(tx, "Significance: %2.2f #sigma", sigLiMa); - pt->AddText(tx); - pt->SetFillStyle(0); - pt->SetBorderSize(0); - pt->SetTextAlign(12); - pt->Draw(""); - } - if (draw||drawPad) { - - // Plot significance vs alpha_cut - - Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0); - Int_t maxBin = centerBin > alphaHisto.GetNbinsX()- centerBin ? - alphaHisto.GetNbinsX()- centerBin : centerBin; - - const Int_t totalBins = centerBin; - Float_t alpha[totalBins]; - Float_t signi[totalBins]; - Float_t maxSigni = 0; - - for (Int_t i=0; i < maxBin; i++) { - Double_t nOn = 0; Double_t eNOn = 0; - Double_t nOff = 0; Double_t eNOff = 0; - for (Int_t bin=centerBin-i; binEval(x); - nOff += offEvts; - } //for bin - eNOn = sqrt(eNOn); - eNOff = sqrt(nOff); - alpha[i] = - (alphaHisto.GetXaxis()->GetBinLowEdge(centerBin+i+1)- - alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2; - signi[i] = CalcSignificance(nOn, nOff, 1); - maxSigni = maxSigni > signi[i] ? maxSigni : signi[i]; - } - - if (!drawPad) { - TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600); - c3->cd(); - } - - if (drawPad) - drawPad->cd(drawBase-1); - - TGraph* gr = new TGraph(totalBins-2,alpha,signi); - TString drawOpt = "L"; - - if (draw || (drawPad && fSigniPlotIndex == 0)) - drawOpt += "A"; - - gr->Draw(drawOpt); - if (drawPad && fSigniPlotIndex == 0) { - gr->GetXaxis()->SetTitle("|#alpha_{max}|"); - gr->GetYaxis()->SetTitle("Significance"); - } - gr->SetMarkerStyle(2); - gr->SetMarkerSize(1); - gr->SetMarkerColor(4+fSigniPlotIndex); - gr->SetLineColor(4+fSigniPlotIndex); - gr->SetLineWidth(1); - gr->SetTitle("Significance vs ALPHA integration range"); - gr->Draw("L"); - - fSigniPlotIndex++; - - } //draw - - return kTRUE; -} - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a MHAlphaEnergyTheta or -// MHAlphaEnergyTime object expected in the ParList. -// -Bool_t MHOnSubtraction::Calc(MParList *parlist, const Bool_t Draw) -{ - //Find source histograms - fHistogramType = "Theta"; - MHAlphaEnergyTheta* mhisttheta = (MHAlphaEnergyTheta*)parlist->FindObject("MHAlphaEnergyTheta"); - if (mhisttheta) return Calc((TH3D*)mhisttheta->GetHist(), parlist, Draw); - - fHistogramType = "Time"; - MHAlphaEnergyTime* mhisttime = (MHAlphaEnergyTime*)parlist->FindObject("MHAlphaEnergyTime"); - if (mhisttime) return Calc((TH3D*)mhisttime->GetHist(), parlist, Draw); - - *fLog << err << "No MHAlphaEnergyTheta type object found in the parameter list. Aborting." << endl; - return kFALSE; -} - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a TH3D histogram in -// ALPHA, Energy and Theta. -// -Bool_t MHOnSubtraction::CalcAET(TH3D *aetHisto, MParList *parlist, const Bool_t Draw) -{ - // Analyze aetHisto - // ------------------------------------- - Int_t alphaBins = aetHisto->GetNbinsX(); // # of alpha bins - Int_t energyBins = aetHisto->GetNbinsY(); - Int_t thetaBins = aetHisto->GetNbinsZ(); - - // We output an array of histograms to the parlist. - // Create a template for such a histogram, needed by MHArray - // ------------------------------------- - MBinning *binsfft = new MBinning("BinningMH3X"); - binsfft->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1), - aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1)); - parlist->AddToList(binsfft); - MH3 *energyHistogram = new MH3(1); // The template histogram in energy - // for a given theta value - energyHistogram->SetName("MHOnSubtractionEnergyHistogramTemplate"); - parlist->AddToList(energyHistogram); - - fThetaHistoArray = new MHArray("MHOnSubtractionEnergyHistogramTemplate", kTRUE, - "MHOnSubtractionGammaSignalArray", - "Array of histograms in energy bins"); - fThetaHistoArray->SetupFill(parlist); - parlist->AddToList(fThetaHistoArray); - - // Canvases---direct debug output for the time being - // ------------------------------------- - TCanvas *c3 = new TCanvas("c3", "Plots by MHOnSubtraction::ExtractSignal", 800,600); - cout << thetaBins << " x " << energyBins << endl; - c3->Divide(thetaBins,energyBins); - - TCanvas *c4a = new TCanvas("c4a", "Energy distributions for different ZA", 800,600); - - TH1D* histalphaon[energyBins*thetaBins]; // ALPHA histograms - - fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5); - fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5); - fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", - alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1), - aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - - // ------------------------------------- - - fThetaLegend = new TLegend(0.83, 0.07, 0.98, 0.42, "Energy Distributions"); - fThetaLegend->SetBorderSize(1); - - Double_t overallSigLiMa = 0; - - for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) { - - char hname[80]; - sprintf(hname, "Energy distribution for %s bin %d", fHistogramType.Data(), - thetaBin); - *fLog << all << "Calculating " << hname << endl; - - Double_t minLowerBin = -10; // minimum integration range - Double_t maxUpperBin = 10; // minimum integration range - Double_t maxAlpha = 70; // maximum integration range - - Double_t sumSigLiMa = 0; - - // This loop just fixes the integration range - // And alerts when no significant excess is found. - - for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { - - sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin); - histalphaon[(thetaBin-1)*(energyBins)+energyBin] = - new TH1D(hname, "ALPHA histogram", - alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1), - aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - histalphaon[(thetaBin-1)*(energyBins)+energyBin]->Sumw2(); - histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetTitle(hname); - - sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d", - energyBin, thetaBin); - *fLog << inf << hname << endl; - - // now we fill one histogram for one particular set - // of Energy and Theta... - - if (aetHisto->InheritsFrom("TH3D")||aetHisto->InheritsFrom("TH2D")) { - aetHisto->GetYaxis()->SetRange(energyBin,energyBin); // E - if (aetHisto->InheritsFrom("TH3D")) - aetHisto->GetZaxis()->SetRange(thetaBin,thetaBin); // theta - histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto->Project3D("xe"); - } else { - histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto; - } - - *fLog << inf << "This histogram has " - << histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetEntries() - << " entries" << endl; - - sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin); - histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetName(hname); - - // Test: Find signal region and significance in histogram - Double_t lowerBin, upperBin, sSigLiMa; - FitHistogram(*histalphaon[(thetaBin-1)*(energyBins)+energyBin], - sSigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE); - - Double_t redChiSq = histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetChisquare()/histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetNDF(); - - fChiSquareHisto->Fill(redChiSq); - fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq; - - // histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetProb() <Add(histalphaon[(thetaBin-1)*(energyBins)+energyBin], 1); - - if (sSigLiMa < 3) - *fLog << inf << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"< maxAlpha) ? maxAlpha : upperBin; - maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin; - - } //energyBin - - *fLog << inf << "Integration range for this " << fHistogramType - << " value will be " << minLowerBin << " < ALPHA < " - << maxUpperBin << endl; - - *fLog << inf << "Summed estd. significance for this " << fHistogramType - << " value is " << sumSigLiMa << endl; - - // Create Histogram in Energy for this Theta value - cout << "STARTIGN REAL CALC1 "<< endl; - fThetaHistoArray->AddHistogram(); - cout << "STARTIGN REAL CALC1a "<< endl; - fThetaHistoArray->SetIndex(thetaBin-1); - - MH3 &mThetaHisto = *(MH3*)(fThetaHistoArray->GetH()); - mThetaHisto.Print(); - - TH1D &thetaHisto = (TH1D&)(mThetaHisto.GetHist()); - sprintf(hname,"Energy distribution for theta bin %d", thetaBin); - thetaHisto.SetTitle(hname); - thetaHisto.SetEntries(0); - - // we have a rough idea of what is going on in the ALPHA plot... - // So now we can indeed extract the signal. - -// cout << "STARTIGN REAL CALC "<< endl; - - sumSigLiMa = 0; - for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { - - sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d", - energyBin, thetaBin); - *fLog << inf << hname << endl; - - Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa; - - c3->cd((thetaBin-1)*(energyBins)+energyBin); - - ExtractSignal(*histalphaon[(thetaBin-1)*(energyBins)+energyBin], - sigLiMa, minLowerBin, maxUpperBin, - gammaSignal, errorGammaSignal, - off, errorOff, (Float_t)3, kTRUE); - - fSignificanceHisto->Fill(sigLiMa); - fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif; - - thetaHisto.SetBinContent(energyBin, gammaSignal*TMath::Exp(7-energyBin)); - thetaHisto.SetBinError(energyBin, errorGammaSignal); - thetaHisto.SetEntries(thetaHisto.GetEntries()+gammaSignal); - - sumSigLiMa += sigLiMa; - - }//energyBin - - *fLog << inf << "Summed significance for this " << fHistogramType - << " value is " << sumSigLiMa << endl; - - //fit exponential function to data - TF1* e = new TF1("expF","expo",0,5); - e->SetLineWidth(3); - e->SetLineColor(thetaBin); - // e->SetParLimits(1, -0.5,3); //limits on slope - - if (fSlope!=20.0) { - e->SetParameter(1, fSlope); - *fLog << inf << "Fixing slope to " << e->GetParameter(1) << endl; - } - - TString eDrawOptions = Draw ? "RIQ" : "RIQ0"; - cout << "doing the fit" << endl; - thetaHisto.Fit("expF"); - - Double_t expoSlope = e->GetParameter(1); - - *fLog << inf << "Determined slope for energy distribution is " - << expoSlope << endl; - - Double_t integralEvts = - e->Integral(aetHisto->GetYaxis()->GetBinLowEdge(1), - aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1)); - - *fLog << inf << "Integral in E range [" << - aetHisto->GetYaxis()->GetBinLowEdge(1) << ";" << - aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1) << "] is " << - integralEvts << endl; - - // Plot Energy histogram - - c4a->cd(0); - thetaHisto.SetLineColor(thetaBin); - if (thetaBin==1) { - thetaHisto.Draw(); - } else { - thetaHisto.Draw("SAME"); - } - - sprintf(hname,"Theta bin %d",thetaBin); - fThetaLegend->AddEntry(&thetaHisto, hname, "l"); - - overallSigLiMa += sumSigLiMa; - - } //thetaBin - - fThetaLegend->Draw(); - - Double_t sigLiMa, lowerBin, upperBin; - FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin); - fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution"); - - *fLog << inf << "Summed overall significance is " << overallSigLiMa << endl; - - *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl; - // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1); - - // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl; - // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5); - - return kTRUE; -} - - - - - - - - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a TH3D histogram in -// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms -// in energy and Theta. -// -Bool_t MHOnSubtraction::Calc(TH3 *aetHisto, MParList *parlist, const Bool_t Draw) -{ - // Analyze aetHisto - // ------------------------------------- - Int_t energyBins = aetHisto->GetNbinsY(); - Int_t thetaBins = aetHisto->GetNbinsZ(); - - *fLog << "Total events: " << aetHisto->GetEntries() << endl; - *fLog << "Total energy bins: " << energyBins << endl; - *fLog << "Total theta bins: " << thetaBins << endl; - - // We output results in an array of histograms to the parameter list. - // Create a template for such a histogram, needed by MH3 - // ------------------------------------- - MBinning *binsmh3x = new MBinning("BinningMH3X"); - // dummy binning, real one follows some lines down - binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1), - aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1)); - parlist->AddToList(binsmh3x); - - MBinning *binsmh3y = new MBinning("BinningMH3Y"); - // dummy binning, real one follows some lines down - binsmh3y->SetEdges(thetaBins,aetHisto->GetZaxis()->GetBinLowEdge(1), - aetHisto->GetZaxis()->GetBinLowEdge(thetaBins+1)); - parlist->AddToList(binsmh3y); - - MH3 *signalHisto = new MH3(2); // Signal(E,t) - signalHisto->SetupFill(parlist); - parlist->AddToList(signalHisto); - signalHisto->SetName("MHOnSubtractionSignal"); - signalHisto->SetTitle("Gamma Events"); - - MH3 *offHisto = new MH3(2); // Off(E,t) - offHisto->SetupFill(parlist); - parlist->AddToList(offHisto); - offHisto->SetName("MHOnSubtractionOff"); - offHisto->SetTitle("Background Events"); - - MH3 *signifHisto = new MH3(2); // Significance(E,t) - signifHisto->SetupFill(parlist); - parlist->AddToList(signifHisto); - signifHisto->SetName("MHOnSubtractionSignif"); - signifHisto->SetTitle("Significance"); - -// if (!fChiSquareHisto) -// fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5); -// if (!fSignificanceHisto) -// fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41); -// if (!fSummedAlphaPlots) -// fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", -// alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1), -// aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - - TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600); - c4->Divide(energyBins+3,thetaBins,0,0,0); - - TH2D& signalTH2DHist = (TH2D&)signalHisto->GetHist(); - TH2D& offTH2DHist = (TH2D&)offHisto->GetHist(); - TH2D& signifTH2DHist = (TH2D&)signifHisto->GetHist(); - - signalTH2DHist.Reset(); - signalTH2DHist.SetTitle("Gamma Signal"); - signalTH2DHist.SetXTitle("E [GeV]"); - signalTH2DHist.SetYTitle("theta [deg]"); - signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); - signalTH2DHist.Sumw2(); - - offTH2DHist.Reset(); - offTH2DHist.SetTitle("Off Signal"); - offTH2DHist.SetXTitle("E [GeV]"); - offTH2DHist.SetYTitle("theta [deg]"); - offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); - offTH2DHist.Sumw2(); - - signifTH2DHist.Reset(); - signifTH2DHist.SetTitle("Significance"); - signifTH2DHist.SetXTitle("E [GeV]"); - signifTH2DHist.SetYTitle("theta [deg]"); - signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); - signifTH2DHist.Sumw2(); - - for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) { - - *fLog << dbg << "--------------------------------------" << endl; - *fLog << dbg << "Processing ALPHA plots for theta bin " << thetaBin << endl; - *fLog << dbg << "--------------------------------------" << endl; - aetHisto->GetZaxis()->SetRange(thetaBin, thetaBin); - TH2F* aeHisto = (TH2F*)aetHisto->Project3D("yxe"); - aeHisto->SetDirectory(NULL); - char hname[80]; - sprintf(hname, "%2.1f < #theta < %2.1f", - aetHisto->GetZaxis()->GetBinLowEdge(thetaBin), - aetHisto->GetZaxis()->GetBinLowEdge(thetaBin+1)); - *fLog << inf << "There are " << aeHisto->GetEntries() - << " entries in " << hname << endl; - aeHisto->SetTitle(hname); - sprintf(hname, "MHOnSubtractionThetaBin%d", thetaBin); - aeHisto->SetName(hname); - - c4->cd((energyBins+3)*(thetaBin-1)+1); - aeHisto->SetMarkerColor(3); - aeHisto->DrawCopy(); - - c4->Update(); - - // We hand over a nonstandard parameter list, which - // contails lower-dimensional result histograms - // signalHisto, offHisto and signifHisto - -// cout << fLog->GetDebugLevel() << endl; - fLog->SetDebugLevel(1); - - MParList *parlistth2 = new MParList(); - // parlistth2->SetNullOutput(); - parlistth2->AddToList(binsmh3x); - - MH3* signalHisto2 = new MH3(1); // Signal(E) - signalHisto2->SetupFill(parlistth2); - parlistth2->AddToList(signalHisto2); - signalHisto2->SetName("MHOnSubtractionSignal"); - signalHisto2->SetTitle("Gamma Events"); - - MH3* offHisto2 = new MH3(1); // Off(E) - offHisto2->SetupFill(parlistth2); - parlistth2->AddToList(offHisto2); - offHisto2->SetName("MHOnSubtractionOff"); - offHisto2->SetTitle("Background Events"); - - MH3* signifHisto2 = new MH3(1); // Significance(E) - signifHisto2->SetupFill(parlistth2); - parlistth2->AddToList(signifHisto2); - signifHisto2->SetName("MHOnSubtractionSignif"); - signifHisto2->SetTitle("Significance"); - - fLog->SetDebugLevel(-1); - - TH1D& signalTH1DHist = (TH1D&)signalHisto2->GetHist(); - TH1D& offTH1DHist = (TH1D&)offHisto2->GetHist(); - TH1D& signifTH1DHist = (TH1D&)signifHisto2->GetHist(); - - signalTH1DHist.Sumw2(); - offTH1DHist.Sumw2(); - signifTH1DHist.Sumw2(); - - TH2Calc(aeHisto, parlistth2, kFALSE, c4, - (energyBins+3)*(thetaBin-1)+4); - - - for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { - - // cout << "filling " << thetaBin << " " << energyBin << ": " -// << signalTH1DHist.GetBinContent(energyBin) << "+-" -// << signalTH1DHist.GetBinError(energyBin) << " " -// << offTH1DHist.GetBinContent(energyBin) << "+-" -// << offTH1DHist.GetBinError(energyBin) << endl; - - if (signalTH1DHist.GetBinContent(energyBin)>=0.0) { - signalTH2DHist.SetBinContent(energyBin, thetaBin, - signalTH1DHist.GetBinContent(energyBin)); - signalTH2DHist.SetBinError(energyBin, thetaBin, - signalTH1DHist.GetBinError(energyBin)); - } - - if (offTH1DHist.GetBinContent(energyBin)>=0.0) { - offTH2DHist.SetBinContent(energyBin, thetaBin, - offTH1DHist.GetBinContent(energyBin)); - offTH2DHist.SetBinError(energyBin, thetaBin, - offTH1DHist.GetBinError(energyBin)); - } - - signifTH2DHist.SetBinContent(energyBin, thetaBin, - signifTH1DHist.GetBinContent(energyBin)); - signifTH2DHist.SetBinError(energyBin, thetaBin, - signifTH1DHist.GetBinError(energyBin)); - } //energyBin - - - c4->cd((energyBins+3)*(thetaBin-1)+2); - sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2); - TPad* tp = (TPad*)(gROOT->FindObject(hname)); - tp->SetLogx(); - signalTH1DHist.SetLineColor(2); - signalTH1DHist.DrawCopy(); - offTH1DHist.SetLineColor(4); - offTH1DHist.DrawCopy("SAME"); - c4->Update(); - - signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries()); - offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries()); - signifTH2DHist.SetEntries(signifTH2DHist.GetEntries()+signifTH1DHist.GetEntries()); - - delete signifHisto2; - delete offHisto2; - delete signalHisto2; - delete parlistth2; - - } - - - c4->Update(); - - return kTRUE; -} - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/* - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a TH3D histogram in -// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms -// in energy and Theta. -// -Bool_t MHOnSubtraction::CalcLightCurve(TH3 *aetHisto, MParList *parlist, const Bool_t draw) -{ - // Analyze aetHisto - // ------------------------------------- - Int_t alphaBins = aetHisto->GetNbinsX(); - Int_t energyBins = aetHisto->GetNbinsY(); - Int_t timeBins = aetHisto->GetNbinsZ(); - - *fLog << "Total events: " << aetHisto->GetEntries() << endl; - *fLog << "Total energy bins: " << energyBins << endl; - *fLog << "Total time bins: " << timeBins << endl; - - // We output results in an array of histograms to the parameter list. - // Create a template for such a histogram, needed by MH3 - // ------------------------------------- - MBinning *binsmh3x = new MBinning("BinningMH3X"); - // dummy binning, real one follows some lines down - binsmh3x->SetEdges(timeBins,aetHisto->GetZaxis()->GetBinLowEdge(1), - aetHisto->GetZaxis()->GetBinLowEdge(timeBins+1)); - parlist->AddToList(binsmh3x); - - MH3 *signalHisto = new MH3(1); // Signal(t) - signalHisto->SetupFill(parlist); - parlist->AddToList(signalHisto); - signalHisto->SetName("MHOnSubtractionSignal"); - signalHisto->SetTitle("Gamma Events"); - - MH3 *offHisto = new MH3(1); // Off(t) - offHisto->SetupFill(parlist); - parlist->AddToList(offHisto); - offHisto->SetName("MHOnSubtractionOff"); - offHisto->SetTitle("Background Events"); - - MH3 *signifHisto = new MH3(1); // Significance(t) - signifHisto->SetupFill(parlist); - parlist->AddToList(signifHisto); - signifHisto->SetName("MHOnSubtractionSignif"); - signifHisto->SetTitle("Significance"); - - TH1D& signalTH1DHist = (TH1D&)signalHisto->GetHist(); - TH1D& offTH1DHist = (TH1D&)offHisto->GetHist(); - TH1D& signifTH1DHist = (TH1D&)signifHisto->GetHist(); - - signalTH1DHist.Reset(); - signalTH1DHist.SetTitle("Gamma Signal"); - signalTH1DHist.SetXTitle("Time [MJD]"); - signalTH1DHist.SetYTitle("Events"); - signalHisto->SetBinning(&signalTH1DHist, aetHisto->GetZaxis()); - signalTH1DHist.Sumw2(); - - offTH1DHist.Reset(); - offTH1DHist.SetTitle("Background Signal"); - offTH1DHist.SetXTitle("Time [MJD]"); - offTH1DHist.SetYTitle("Events"); - offHisto->SetBinning(&offTH1DHist, aetHisto->GetZaxis()); - offTH1DHist.Sumw2(); - - signifTH1DHist.Reset(); - signifTH1DHist.SetTitle("Significance"); - signifTH1DHist.SetXTitle("Time [MJD]"); - signifTH1DHist.SetYTitle("Significance #sigma"); - signifHisto->SetBinning(&signifTH1DHist, aetHisto->GetZaxis()); - signifTH1DHist.Sumw2(); - - //The following histogram is an additional histogram needed for - //the light curve - - TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600); - c4->Divide(8,7,0,0,0); - - - - // The following loop fixes a common integration region for each - // energy bin and alerts when no significant excess is found. - - Double_t minLowerBin = -10; // minimum integration range - Double_t maxUpperBin = 10; // minimum integration range - Double_t maxAlpha = 70; // maximum integration range - Double_t sumSigLiMa = 0; - - TH1F* alphaHisto[timeBins]; - - for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) { - - *fLog << dbg << "--------------------------------------" << endl; - *fLog << dbg << "Processing ALPHA plots for time bin " << timeBin << endl; - *fLog << dbg << "--------------------------------------" << endl; - - aetHisto->GetZaxis()->SetRange(timeBin, timeBin); - - char hname[80]; - sprintf(hname, "MHOnSubtractionTimeBin%d", timeBin); - alphaHisto[timeBin-1] = - new TH1F(hname, "ALPHA histogram", - alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1), - aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - alphaHisto[timeBin-1]->SetDirectory(NULL); - alphaHisto[timeBin-1]->Sumw2(); - alphaHisto[timeBin-1] = (TH1F*)aetHisto->Project3D("xe"); - alphaHisto[timeBin-1]->SetName(hname); - alphaHisto[timeBin-1]->SetDirectory(NULL); - - sprintf(hname, "%6.0f < t < %6.0f", - aetHisto->GetZaxis()->GetBinLowEdge(timeBin), - aetHisto->GetZaxis()->GetBinLowEdge(timeBin+1)); - *fLog << inf << "There are " << alphaHisto[timeBin-1]->GetEntries() - << " entries in " << hname << endl; - alphaHisto[timeBin-1]->SetTitle(hname); - - // Find signal region and significance in histogram - Double_t lowerBin, upperBin, sSigLiMa; - char funcName[20]; - sprintf(funcName, "%2d", timeBin); - - Bool_t drawFit = kTRUE; - - c4->cd(timeBin); -// alphaHisto[timeBin-1]->SetMarkerColor(3); - alphaHisto[timeBin-1]->DrawCopy(); - - c4->Update(); - - fSigniPlotColor = 0; - ; - Bool_t fit = FitHistogram(*alphaHisto[timeBin-1], sSigLiMa, - lowerBin, upperBin, (Float_t)3.0, drawFit, - funcName); - - if (fit) { - if (sSigLiMa < 3) - *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"< maxAlpha) ? maxAlpha : upperBin; - maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin; - } - - - - - - - } //timeBin - - *fLog << inf << "=> Integration range will be " << minLowerBin << " < ALPHA < " - << maxUpperBin << "," << endl; - *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl; - - // we have an idea of what is going on in the ALPHA plot... - // So now we can indeed extract the signal. - - sumSigLiMa = 0; - for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) { - *fLog << inf << "Processing ALPHA distribution for time bin " << timeBin << endl; - if (alphaHisto[timeBin-1]->GetEntries() == 0) { - *fLog << warn << "No attempt to extract a signal from 0 events." << endl; - if (draw) { - TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-"); - lab->SetBorderSize(0); - lab->Draw(); - } - } else { - char funcName[20]; - sprintf(funcName, "%2d", timeBin); - - Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa; - - Bool_t drawFit = kFALSE; - - ExtractSignal(*alphaHisto[timeBin-1], - sigLiMa, minLowerBin, maxUpperBin, - gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit, - funcName, NULL, 1); - - sumSigLiMa += sigLiMa; - - fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif; - - // Fill Signal - TH1D &h = (TH1D&)(signalHisto->GetHist()); - h.SetBinContent(timeBin, gammaSignal); - h.SetBinError(timeBin, errorGammaSignal); - h.SetEntries(h.GetEntries()+gammaSignal); - - // Fill Off Events - TH1D &ho = (TH1D&)(offHisto->GetHist()); - ho.SetBinContent(timeBin, off); - ho.SetBinError(timeBin, errorOff); - ho.SetEntries(ho.GetEntries()+off); - - // Fill Significance - TH1D &hg = (TH1D&)(signifHisto->GetHist()); - hg.SetBinContent(timeBin, sigLiMa); - hg.SetEntries(hg.GetEntries()+sigLiMa); - } - - }//timeBin - - *fLog << inf << "Summed significance is " << sumSigLiMa << endl; - - - if (draw) { - Double_t sigLiMa, lowerBin, upperBin; - FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin); - fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution"); - - *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl; - // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1); - - // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl; - // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5); - - fChiSquareHisto->DrawClone(); - fSignificanceHisto->DrawClone(); - fSummedAlphaPlots->DrawClone(); - } - - return kTRUE; -} - - - - - -*/ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a TH2 histogram in -// ALPHA and Energy. Output to parameter list is TH1 histogram in -// energy -// -Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase) -{ - - fSigniPlotColor = 0; - - // Analyze aeHisto - // ------------------------------------- - const Int_t alphaBins = aeHisto->GetNbinsX(); // # of alpha bins - const Int_t energyBins = aeHisto->GetNbinsY(); - - // Check availability of result histograms - // ------------------------------------- - - MH3* signalHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignal","MH3"); - MH3* offHisto = (MH3*)parlist->FindObject("MHOnSubtractionOff","MH3"); - MH3* signifHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignif","MH3"); - - if (signalHisto && offHisto && signifHisto) { -// *fLog << dbg << "Result histograms are present in parameter list " << -// "and will be used further." << endl; - } else { - - *fLog << warn << "Result histograms are not present in parameter " << - "list and thus are going to be created now." << endl; - - MBinning *binsmh3x = new MBinning("BinningMH3X"); - binsmh3x->SetEdgesLog(energyBins,aeHisto->GetYaxis()->GetBinLowEdge(1), - aeHisto->GetYaxis()->GetBinLowEdge(energyBins+1)); - parlist->AddToList(binsmh3x); - - signalHisto = new MH3(1); // Signal(E) - signalHisto->SetName("MHOnSubtractionSignal"); - signalHisto->SetTitle("Extracted Signal"); - parlist->AddToList(signalHisto); - signalHisto->SetBinning(&((TH1D&)signalHisto->GetHist()), - aeHisto->GetYaxis()); - - offHisto = new MH3(1); // Off(E) - offHisto->SetName("MHOnSubtractionOff"); - offHisto->SetTitle("Off Signal"); - parlist->AddToList(offHisto); - offHisto->SetBinning(&((TH1D&)offHisto->GetHist()), - aeHisto->GetYaxis()); - - signifHisto = new MH3(1); // Significance(E) - signifHisto->SetName("MHOnSubtractionSignif"); - signifHisto->SetTitle("Significance"); - parlist->AddToList(signifHisto); - signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()), - aeHisto->GetYaxis()); - - } - if (fChiSquareHisto==0x0) - fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5); - if (fSignificanceHisto==0x0) - fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5); - if (fSummedAlphaPlots==0x0) - fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", - alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1), - aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - - // The following loop fixes a common integration region for each - // energy bin and alerts when no significant excess is found. - - Double_t minLowerBin = -10; // minimum integration range - Double_t maxUpperBin = 10; // minimum integration range - Double_t maxAlpha = 70; // maximum integration range - Double_t sumSigLiMa = 0; - - TH1F* alphaHisto[energyBins]; - - fSigniPlotIndex = 0; // cf. ExtractSignal - - for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { - *fLog << inf << "Preprocessing ALPHA distribution for E bin " << energyBin << endl; - char hname[80]; - sprintf(hname, "MHOnSubtractionAlpha%d", energyBin); - alphaHisto[energyBin-1] = - new TH1F(hname, "ALPHA histogram", - alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1), - aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1)); - alphaHisto[energyBin-1]->SetDirectory(NULL); - alphaHisto[energyBin-1]->Sumw2(); - alphaHisto[energyBin-1] = (TH1F*)aeHisto->ProjectionX("xe", energyBin, energyBin); - alphaHisto[energyBin-1]->SetName(hname); - alphaHisto[energyBin-1]->SetDirectory(NULL); - - sprintf(hname, "%2.1f < E < %2.1f", - aeHisto->GetYaxis()->GetBinLowEdge(energyBin), - aeHisto->GetYaxis()->GetBinLowEdge(energyBin+1)); - alphaHisto[energyBin-1]->SetTitle(hname); - // alphaHisto[energyBin-1]->DrawCopy(); - alphaHisto[energyBin-1]->SetDirectory(NULL); - - // Find signal region and significance in histogram - Double_t lowerBin, upperBin, sSigLiMa; - char funcName[20]; - sprintf(funcName, "%2d", energyBin); - - Bool_t drawFit = kFALSE; - - if (drawPad) { - drawPad->cd(drawBase+energyBin-1); - drawFit = kTRUE; - } - - Bool_t fit = FitHistogram(*alphaHisto[energyBin-1], sSigLiMa, - lowerBin, upperBin, (Float_t)3.0, drawFit, - funcName); - - if (fit) { - Double_t redChiSq = alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetChisquare()/ - alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetNDF(); - fChiSquareHisto->Fill(redChiSq); - fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq; - fSummedAlphaPlots->Add(alphaHisto[energyBin-1], 1); - if (sSigLiMa < 3) - *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"< maxAlpha) ? maxAlpha : upperBin; - maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin; - } else { - fChiSquareHisto->Fill(0); - } - -// debugOut->Update(); - - } //energyBin - - *fLog << inf << "=> Integration range for this theta bin will be " << minLowerBin << " < ALPHA < " - << maxUpperBin << "," << endl; - *fLog << inf << " Summed estimated significance is " << sumSigLiMa << endl; - - // we have an idea of what is going on in the ALPHA plot... - // So now we can indeed extract the signal. - - sumSigLiMa = 0; - for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { - *fLog << inf << "Processing ALPHA distribution for E bin " << energyBin << endl; - if (alphaHisto[energyBin-1]->GetEntries() == 0) { - *fLog << warn << "No attempt to extract a signal from 0 events." << endl; - if (draw) { - TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-"); - lab->SetBorderSize(0); - lab->Draw(); - } - } else { - char funcName[20]; - sprintf(funcName, "%2d", energyBin); - - Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa; - - Bool_t drawFit = kFALSE; - -// if (drawPad) { -// cout << "Going to use pad " <cd(drawBase+energyBin-1); -// drawFit = kTRUE; -// } - - - ExtractSignal(*alphaHisto[energyBin-1], - sigLiMa, minLowerBin, maxUpperBin, - gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit, - funcName, drawPad, drawBase); - - sumSigLiMa += sigLiMa; - fSignificanceHisto->Fill(sigLiMa); - fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif; - - // Fill Signal - TH1D &h = (TH1D&)(signalHisto->GetHist()); - h.SetBinContent(energyBin, gammaSignal); - h.SetBinError(energyBin, errorGammaSignal); - h.SetEntries(h.GetEntries()+gammaSignal); - - // Fill Off Events - TH1D &ho = (TH1D&)(offHisto->GetHist()); - ho.SetBinContent(energyBin, off); - ho.SetBinError(energyBin, errorOff); - ho.SetEntries(ho.GetEntries()+off); - - // Fill Significance - TH1D &hg = (TH1D&)(signifHisto->GetHist()); - hg.SetBinContent(energyBin, sigLiMa); - hg.SetEntries(hg.GetEntries()+sigLiMa); - } - - }//energyBin - - *fLog << inf << "Summed significance is " << sumSigLiMa << endl; - - - if (draw) { - Double_t sigLiMa, lowerBin, upperBin; - FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin); - fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution"); - - *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl; - // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1); - - // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl; - // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5); - - fChiSquareHisto->DrawClone(); - fSignificanceHisto->DrawClone(); - fSummedAlphaPlots->DrawClone(); - } - - return kTRUE; -} - - - - - - - - - - - - - - - - -// ----------------------------------------------------------------------- -// -// Extraction of Gamma signal is performed on a TH1 histogram in ALPHA -// -// -Bool_t MHOnSubtraction::Calc(TH1 *aHisto, MParList *parlist, - const Bool_t draw, Float_t signalRegion) -{ - // Find signal region and estimate significance - Double_t lowerBin, upperBin, sigLiMa; - FitHistogram(*aHisto, sigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE); - - //if (sigLiMa < 3) - // *fLog << err << "No significant excess (sigma=" << sigLiMa <<")!"<SetSelectedPad(NULL); - gStyle->SetOptStat(0); - - c.cd(1); - fSummedAlphaPlots->DrawCopy(opt); - - c.cd(2); - fSignificanceHisto->DrawCopy(opt); - - c.cd(3); - TString drawopt=""; - for (fThetaHistoArray->SetIndex(0); - MH3* h=(MH3*)(fThetaHistoArray->GetH()); - fThetaHistoArray->IncIndex()) - { - // Get the current histogram - TH1D& hist = (TH1D&)h->GetHist(); - fSummedAlphaPlots->SetTitle("Energy distributions for different theta bins"); - hist.Draw(drawopt); - drawopt="SAME"; - } - fThetaLegend->Draw(); - - c.cd(4); - fChiSquareHisto->DrawCopy(opt); - - c.Modified(); - c.Update(); - - return &c; -} - -// ------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHOnSubtraction::Draw(Option_t *opt) -{ - if (!gPad) - MakeDefCanvas("OnSubtraction", "Results from OnSubtraction"); - - gPad->Divide(2,3); - -// Ok, at some point this will be implemented - - gPad->Modified(); - gPad->Update(); -} - - - - - Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h (revision 6977) +++ (revision ) @@ -1,75 +1,0 @@ -#ifndef MARS_MHOnSubtraction -#define MARS_MHOnSubtraction - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1D; -class TH3D; -class TPad; -class TGraph; -class TLegend; -class TPaveLabel; - -class MHArray; - -class MHOnSubtraction : public MH -{ -private: - TString fHistogramType; - TH1D* fChiSquareHisto; - TH1D* fSignificanceHisto; - TH1D* fSummedAlphaPlots; - MHArray *fThetaHistoArray; - TLegend *fThetaLegend; - Double_t fMaxSignif; - Double_t fMaxRedChiSq; - Double_t fSlope; // slope for exponential fit - Int_t fThetaBin; - Int_t fSigniPlotIndex; - Int_t fSigniPlotColor; - - Double_t fSignificance; - - Bool_t CalcAET(TH3D *histon, MParList *parlist, const Bool_t Draw); - - Bool_t FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa, - Double_t &lowerBin, Double_t &upperBin, - Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE, - TString funcName = ""); - - Bool_t ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa, - Double_t &lowerBin, Double_t &upperBin, - Double_t &gammaSignal, Double_t &errorGammaSignal, - Double_t &off, Double_t &errorOff, - Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE, - TString funcName = "", TPad *drawPad = 0, Int_t drawBase = 0); - - - -public: - MHOnSubtraction(const char *name=NULL, const char *title=NULL); - ~MHOnSubtraction(); - - Double_t CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta); - Double_t GetSignificance() { return fSignificance; }; - - void SetExpoSlope(Double_t slope) { fSlope = slope; } - - Bool_t Calc(MParList *parlist, const Bool_t Draw); - - Bool_t Calc(TH3 *histon, MParList *parlist, const Bool_t Draw); - Bool_t TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t Draw, - TPad *drawPad = 0, Int_t drawBase = 0); - Bool_t Calc(TH1 *histon, MParList *parlist, const Bool_t Draw, - Float_t signalRegion = 0); - - void Draw(Option_t *option=""); - TObject *DrawClone(Option_t *option="") const; - - ClassDef(MHOnSubtraction, 1) //Extracts gamma signals from pure ON-data -}; - -#endif -