Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2169) +++ trunk/MagicSoft/Mars/Changelog (revision 2170) @@ -1,3 +1,9 @@ -*-*- END OF LINE -*-*- + + 2003/06/11: Robert Wagner + + * mhist/MHOnSubtraction.[h,cc] + - Some bugfixes, e.g. concerning binning of result histograms + - Improvements in output Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2169) +++ trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2170) @@ -26,15 +26,19 @@ ////////////////////////////////////////////////////////////////////////////// // // -// MHOnSubtraction // +// 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. // +// 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 remove commentary -// parts and clearly state where and how you changed the code, if you did so. +// This part of MARS is code still evolving. Please do not change the code +// without prior feedback by the author. #include "MHOnSubtraction.h" @@ -64,7 +68,5 @@ // Default Constructor. // -MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0), - fMaxRedChiSq(0), - fSlope(20.0) +MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0), fMaxRedChiSq(0), fSlope(20.0) { fName = name ? name : "MHOnSubtraction"; @@ -96,20 +98,18 @@ // 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) { + if (nOn<=0) *fLog << warn << "Got " << nOn << " total events, " << flush; - } - - if (nOff<=0) { - *fLog << warn << "Got " << nOff << " background 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))))); @@ -117,4 +117,7 @@ // ----------------------------------------------------------------------- +// +// 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 @@ -177,23 +180,19 @@ *fLog << dbg << "Plot contains " << - outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " << flush; - *fLog << dbg << + 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); @@ -203,7 +202,8 @@ alphaHisto.GetListOfFunctions()->Add(po); alphaHisto.SetLineColor(97); - alphaHisto.DrawCopy(); //rwagner - po->Draw("SAME"); - + if (draw) { + alphaHisto.DrawCopy(); //rwagner + po->Draw("SAME"); + } sigLiMa = 0; lowerBin = 1; @@ -213,4 +213,5 @@ return kFALSE; //No gaus fit applied } + if (outerEvents/outerBins < 2.5) { *fLog << warn << "Only " << /*setprecision(2) <<*/ outerEvents/outerBins @@ -224,6 +225,8 @@ po->SetParameter(0,outerEvents/outerBins); alphaHisto.GetListOfFunctions()->Add(po); - alphaHisto.DrawCopy(); //rwagner - po->Draw("SAME"); + if (draw) { + alphaHisto.DrawCopy(); //rwagner + po->Draw("SAME"); + } Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0); @@ -231,8 +234,8 @@ alphaHisto.GetNbinsX()- centerBin : centerBin; - Int_t lowerSignalEdge = centerBin; for (Int_t i=3; i < maxBin; i++) { - if ((double)alphaHisto.GetBinContent(centerBin-i) < outerEvents/outerBins) { + if ((double)alphaHisto.GetBinContent(centerBin-i) + < outerEvents/outerBins) { lowerSignalEdge = centerBin-i+1; break; @@ -243,5 +246,6 @@ Int_t upperSignalEdge = centerBin; for (Int_t i=3; i < maxBin; i++) { - if ((double)alphaHisto.GetBinContent(centerBin+i) <= outerEvents/outerBins) { + if ((double)alphaHisto.GetBinContent(centerBin+i) + <= outerEvents/outerBins) { upperSignalEdge=centerBin+i-1; break; @@ -251,18 +255,23 @@ double nOnInt = 0; - for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++) { + for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++) nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins; - } - double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) * (outerEvents/outerBins); + double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) + * (outerEvents/outerBins); sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); - - *fLog << inf << "Estimated significance is " << sigLiMa << " sigma " << endl; + + 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; + << lowerBin << " < ALPHA < " << upperBin + << " (Most likely you wanna ignore this)" << endl; return kFALSE; //No gaus fit applied - } @@ -325,6 +334,10 @@ sigLiMa = CalcSignificance(nOnInt, nOffInt, 1); - - *fLog << inf << "Fit estimates significance to be " << sigLiMa << " sigma " << endl; +// if (sigLiMa < 3) +// alphaHisto.SetLineColor(2); + + + *fLog << inf << "Fit estimates significance to be " + << sigLiMa << " sigma " << endl; *fLog << inf << "Fit yields signal region to be " @@ -336,8 +349,13 @@ } + + + + // ----------------------------------------------------------------------- // // Does the actual extraction of the gamma signal. For performance -// reasons, the fit from MHOnSubtraction::FitHistogram is used and not redone. +// 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 < @@ -347,8 +365,9 @@ // 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) + 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); @@ -358,12 +377,9 @@ *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); - -// } + 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; } @@ -414,5 +430,4 @@ if (pol) { - po = pol; @@ -449,6 +464,5 @@ } //for bin eNOn = sqrt(eNOn); - - + // Evaluate background @@ -463,8 +477,7 @@ eNOff = sqrt(fabs(nOff)); - if (nOn==0) { // there should not be a negative number of signal events + 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; @@ -472,5 +485,4 @@ } - *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", "; @@ -498,8 +510,8 @@ pt->SetBorderSize(0); pt->SetTextAlign(12); - pt->Draw(); + pt->Draw(""); } - if (draw) { + if (draw||drawPad) { // Plot significance vs alpha_cut @@ -507,10 +519,10 @@ 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 nOn = 0; double eNOn = 0; @@ -529,19 +541,36 @@ alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2; signi[i] = CalcSignificance(nOn, nOff, 1); + maxSigni = maxSigni > signi[i] ? maxSigni : signi[i]; } - - TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600); - c3->cd(); + + 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); - gr->Draw("ACP"); - gr->GetXaxis()->SetTitle("|#alpha_{max}|"); - gr->GetYaxis()->SetTitle("Significance"); + 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); - gr->SetLineColor(4); - gr->SetLineWidth(2); + gr->SetMarkerColor(4+fSigniPlotIndex); + gr->SetLineColor(4+fSigniPlotIndex); + gr->SetLineWidth(1); gr->SetTitle("Significance vs ALPHA integration range"); gr->Draw("LP"); + + fSigniPlotIndex++; + } //draw @@ -724,5 +753,5 @@ // So now we can indeed extract the signal. - cout << "STARTIGN REAL CALC "<< endl; +// cout << "STARTIGN REAL CALC "<< endl; sumSigLiMa = 0; @@ -757,5 +786,5 @@ //fit exponential function to data - TF1 *e = new TF1("expF","expo",0,5); + TF1* e = new TF1("expF","expo",0,5); e->SetLineWidth(3); e->SetLineColor(thetaBin); @@ -847,9 +876,11 @@ // ------------------------------------- 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)); + 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)); @@ -894,4 +925,5 @@ signalTH2DHist.SetXTitle("E [GeV]"); signalTH2DHist.SetYTitle("theta [deg]"); + signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); signalTH2DHist.Sumw2(); @@ -900,4 +932,5 @@ offTH2DHist.SetXTitle("E [GeV]"); offTH2DHist.SetYTitle("theta [deg]"); + offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); offTH2DHist.Sumw2(); @@ -906,4 +939,5 @@ signifTH2DHist.SetXTitle("E [GeV]"); signifTH2DHist.SetYTitle("theta [deg]"); + signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis()); signifTH2DHist.Sumw2(); @@ -1005,20 +1039,13 @@ c4->cd((energyBins+3)*(thetaBin-1)+2); - sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2); TPad* tp = (TPad*)(gROOT->FindObject(hname)); - tp->SetLogx(); - + tp->SetLogx(); signalTH1DHist.SetLineColor(2); - signalTH1DHist.DrawCopy(); - + signalTH1DHist.DrawCopy(); offTH1DHist.SetLineColor(4); - offTH1DHist.DrawCopy("SAME"); - - + offTH1DHist.DrawCopy("SAME"); c4->Update(); - - - + signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries()); offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries()); @@ -1032,6 +1059,7 @@ } + c4->Update(); - + return kTRUE; } @@ -1084,4 +1112,7 @@ Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase) { + + fSigniPlotColor = 0; + // Analyze aeHisto // ------------------------------------- @@ -1113,9 +1144,13 @@ 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) @@ -1123,4 +1158,7 @@ signifHisto->SetTitle("Significance"); parlist->AddToList(signifHisto); + signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()), + aeHisto->GetYaxis()); + } if (fChiSquareHisto==0x0) @@ -1142,4 +1180,6 @@ TH1F* alphaHisto[energyBins]; + + fSigniPlotIndex = 0; // cf. ExtractSignal for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) { @@ -1238,5 +1278,5 @@ sigLiMa, minLowerBin, maxUpperBin, gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit, - funcName); + funcName, drawPad, drawBase); sumSigLiMa += sigLiMa; @@ -1304,7 +1344,8 @@ // // 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) + const Bool_t draw, Float_t signalRegion) { // Find signal region and estimate significance @@ -1323,5 +1364,6 @@ ExtractSignal(*aHisto, sigLiMa, lowerBin, upperBin, - gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3.0, kTRUE); + gammaSignal, errorGammaSignal, off, errorOff, + (Float_t)3.0, draw); *fLog << inf << "Signal is " Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h (revision 2169) +++ trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h (revision 2170) @@ -17,4 +17,8 @@ #include "TPaveLabel.h" #endif +#ifndef ROOT_TFile +#include "TFile.h" +#endif + @@ -41,4 +45,22 @@ Double_t fSlope; // slope for exponential fit Int_t fThetaBin; + Int_t fSigniPlotIndex; + Int_t fSigniPlotColor; + + 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: @@ -55,22 +77,10 @@ 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 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);