Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2163) +++ trunk/MagicSoft/Mars/Changelog (revision 2164) @@ -1,3 +1,16 @@ -*-*- END OF LINE -*-*- + + 2003/06/06: Robert Wagner + + * mhist/MHOnSubtraction.[h,cc] + - Class for extracting a gamma signal from on data only. Works + on fully differential data in Alpha, Energy and Theta as well + as on single Alpha plots. Experimental version, expect + functionality but code still optimized for debugging purposes + + * mhist/MHAlphaEnergyTheta.cc + - Fill signed alpha value instead of absolute value + + 2003/06/06: Wolfgang Wittek @@ -19,4 +32,6 @@ * mhist/Makefile, HistLinkDef.h : added new class. + + 2003/06/05: Thomas Bretz Index: trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 2163) +++ trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 2164) @@ -114,7 +114,7 @@ MHillasSrc &hil = *(MHillasSrc*)par; - fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), + fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), fMcEvt->GetTelescopeTheta()*kRad2Deg, w); - + return kTRUE; } Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2164) +++ trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2164) @@ -0,0 +1,1430 @@ +/* ======================================================================== *\ +! +! * +! * 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 a MHArray of MH3 histograms in // +// energy, one per theta bin. // +// // +////////////////////////////////////////////////////////////////////////////// + +// This part of MARS is code still evolving. Please do not remove commentary +// parts and clearly state where and how you changed the code, if you did so. + +#include "MHOnSubtraction.h" + +#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); + + +// -------------------------------------------------------------------------- +// +// 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. +// +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))))); +} + +// ----------------------------------------------------------------------- +// +// 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), " << flush; + *fLog << dbg << + 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); + 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); + 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)alphaHisto.GetBinContent(centerBin-i) < outerEvents/outerBins) { + lowerSignalEdge = centerBin-i+1; + break; + } + } + if (centerBinupperSignalEdge) upperSignalEdge=centerBin; + + double nOnInt = 0; + for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++) { + nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins; + } + + double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) * (outerEvents/outerBins); + + sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); + + *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 + + double gausMean = gp->GetParameter(1); + double 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 scalingFactor = + (alphaHisto.GetXaxis()->GetBinLowEdge(alphaHisto.GetNbinsX()+1) + - alphaHisto.GetXaxis()->GetBinLowEdge(1)) / alphaHisto.GetNbinsX(); + + double nOnInt = (gp->Integral(-signalRegionFactor*gausSigma+gausMean, + signalRegionFactor*gausSigma+gausMean) / scalingFactor); + double 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(lowerBin); + upperBin = alphaHisto.GetBinLowEdge(upperBin)+alphaHisto.GetBinWidth(upperBin); + + sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); + + *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, the fit from MHOnSubtraction::FitHistogram is 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) +{ + 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; +// if (draw) { + 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; + + if (gausPol) { + double gausMean = gausPol->GetParameter(1); + double 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 nOn = 0; + double 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) { + // 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]; + + for (Int_t i=0; i < maxBin; i++) { + double nOn = 0; double eNOn = 0; + double nOff = 0; double 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); + } + + TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600); + c3->cd(); + TGraph* gr = new TGraph(totalBins-2,alpha,signi); + gr->Draw("ACP"); + gr->GetXaxis()->SetTitle("|#alpha_{max}|"); + gr->GetYaxis()->SetTitle("Significance"); + gr->SetMarkerStyle(2); + gr->SetMarkerSize(1); + gr->SetMarkerColor(4); + gr->SetLineColor(4); + gr->SetLineWidth(2); + gr->SetTitle("Significance vs ALPHA integration range"); + gr->Draw("LP"); + } //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", 40.5, -0.5, 41); + 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 expoSlope = e->GetParameter(1); + + *fLog << inf << "Determined slope for energy distribution is " + << expoSlope << endl; + + double 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"); + binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1), + aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1)); + parlist->AddToList(binsmh3x); + + MBinning *binsmh3y = new MBinning("BinningMH3Y"); + 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); + + 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]"); + signalTH2DHist.Sumw2(); + + offTH2DHist.Reset(); + offTH2DHist.SetTitle("Off Signal"); + offTH2DHist.SetXTitle("E [GeV]"); + offTH2DHist.SetYTitle("theta [deg]"); + offTH2DHist.Sumw2(); + + signifTH2DHist.Reset(); + signifTH2DHist.SetTitle("Significance"); + signifTH2DHist.SetXTitle("E [GeV]"); + signifTH2DHist.SetYTitle("theta [deg]"); + 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 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) +{ + // 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); + + offHisto = new MH3(1); // Off(E) + offHisto->SetName("MHOnSubtractionOff"); + offHisto->SetTitle("Off Signal"); + parlist->AddToList(offHisto); + + signifHisto = new MH3(1); // Significance(E) + signifHisto->SetName("MHOnSubtractionSignif"); + signifHisto->SetTitle("Significance"); + parlist->AddToList(signifHisto); + } + 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", 40.5, -0.5, 41); + 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]; + + 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); + + 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 2164) +++ trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h (revision 2164) @@ -0,0 +1,85 @@ +#ifndef MARS_MHOnSubtraction +#define MARS_MHOnSubtraction + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef ROOT_TH1 +#include "TH1.h" +#endif +#ifndef ROOT_TH3 +#include "TH3.h" +#endif +#ifndef ROOT_TGraph +#include "TGraph.h" +#endif +#ifndef ROOT_TPaveLabel +#include "TPaveLabel.h" +#endif + + + +class TH2D; +class TPad; +class TLegend; +class TString; + +class MParList; +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; + +public: + MHOnSubtraction(const char *name=NULL, const char *title=NULL); + ~MHOnSubtraction(); + + virtual Bool_t SetupFill(const MParList *pliston); + virtual Bool_t Fill(const MParContainer *pcont, const Stat_t w); + +/* const TH3D *GetHist() { return fH; } */ +/* const TH3D *GetHist() const { return fH; } */ + + Double_t CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta); + + void SetExpoSlope(Double_t slope) { fSlope = slope; } + + 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 = ""); + + Bool_t Calc(MParList *parlist, const Bool_t Draw); + Bool_t CalcAET(TH3D *histon, 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 + Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 2163) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 2164) @@ -59,4 +59,5 @@ MHSigmabarTheta.cc \ MHSigmaTheta.cc \ + MHOnSubtraction.cc \ MHTrigLvl0.cc