/* ======================================================================== *\ ! ! * ! * 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(); }