Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.cc (revision 6634) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.cc (revision 6634) @@ -0,0 +1,846 @@ +/* ======================================================================== *\ +! +! * +! * 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 express +! * or implied warranty. +! * +! +! +! Author(s): Marcols Lopez, August 2004 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHMyFindSignificance +// +// +// Non: # events in the signal region (0,fAlphasig) of the ON-Data hist. +// Nbg: # events in the signal region (0,fAlphasig) of the OFF-Data hist. +// Nex = Non - Nbg +// +///////////////////////////////////////////////////////////////////////////// +#include "MHMyFindSignificanceONOFF.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + + +ClassImp(MHMyFindSignificanceONOFF); + +using namespace std; + + +// -------------------------------------------------------------------------- +// +// Constructor +// +MHMyFindSignificanceONOFF::MHMyFindSignificanceONOFF(const char *name, const char *title) + : fNon(0.0), fdNon(0.0), fNbg(0.0), fdNbg(0.0), fNex(0.0), fdNex(0.0), + fSigLiMa(0.0), + fNbgFit(0.0), fdNbgFit(0.0), fNexFit(0.0), fdNexFit(0.0), fGamma(0.0), + fNoff(0.0), fSigLiMaFit(0.0), + fBestSigma(0.0), fBestAlphaCut(0.0), fBestExcess(0.0) +{ + fName = name ? name : "MHMyFindSignificanceONOFF"; + fTitle = title ? title : "Find Significance in alpha plot"; + + fHistON = NULL; + fHistOrigON = NULL; + fHistOFF = NULL; + fHistOrigOFF = NULL; + fHistOFFNormalized = NULL; + + // allow rebinning of the alpha plot + fRebin = kTRUE; + + // allow reducing the degree of the polynomial + fReduceDegree = kTRUE; +} + + + + +// --------------------------------------------------------------------------- +// +// Calculate +// +// - fNon: number of events in the signal region from ON data +// +// and 3 different estimation of the background: +// +// - fNbgOFF: number of events in the signal region from OFF data +// - fNbgFitOFF: estimation of # events in the signal region from a fit of +// the OFF data in the Bkg region. +// - fNbgFitON: estimation of # events in the signal region from a fit of +// the ON data in the Bkg region. +// - fNex +// +Bool_t MHMyFindSignificanceONOFF::FindSigmaONOFF(TH1 *histON, TH1 *histOFF, + Double_t alphasig, Double_t alphamin, Double_t alphamax, Int_t degree) +{ + + fHistOrigON = histON; + fHistOrigOFF = histOFF; + + fAlphasig = alphasig; + fAlphaminON = alphamin; + fAlphamaxON = alphamax; + fAlphaminOFF = 0; // <-------------------- + fAlphamaxOFF = alphamax; + fDegree = degree; + + + + + // ----------------------------------------------------------------- + // + // Calculate the normalization factor in range (alphanin,alphamax) + // + // FIXME: This can be done directly with: + // fFindsigON.GetNbtot()/fFindsigOFF.GetNbtot() + // but for that the polynomial for OFF-Data has to be fitted + // (alphamin, alphamax) instead of (0,alphamax) + // + CalcNormalizationFactor(fHistOrigON, fHistOrigOFF, fAlphaminON, fAlphamaxON); + + + // ----------------------------------------------------------------------- + // + // Use MHFindSignificance to find the # evts in the signal and Bkg region + // and to fit the histograms. + // The ON hist is fitted to a polynomial in reg. (fAlphaminON,fAlphamaxON) + // and to a gaussian+polynomial. + // The OFF hist is only fitted to a polynomial in(fAlphaminOFF,fAlphamaxOFF + // + fFindsigON.SetReduceDegree(fReduceDegree); + fFindsigOFF.SetReduceDegree(fReduceDegree); + fFindsigON.SetRebin(fRebin); + fFindsigOFF.SetRebin(fRebin); + + const Bool_t rc1 = fFindsigON.FindSigma(fHistOrigON, fAlphaminON, fAlphamaxON, fDegree, fAlphasig, 0,1,0); + + const Bool_t rc2 = fFindsigOFF.FindSigma(fHistOrigOFF, fAlphaminOFF,fAlphamaxOFF, fDegree, fAlphasig, 0,0,0); + + if (rc1==kFALSE || rc2==kFALSE) + { + *fLog << err << dbginf << "MHFindSignificance retuns kFALSE." << endl; + return kFALSE; + } + + + // Check consistency + if (fabs(fFindsigON.GetAlphasi() - fFindsigOFF.GetAlphasi()) > fEps) + { + *fLog << "fAlphasiOFF (" << fFindsigOFF.GetAlphasi()<< ") is not equal to fAlphasi (" + << fFindsigON.GetAlphasi() << "), this is something that should not happen" + << endl; + //return kFALSE; It might happen in pathological cases (very few OFF) + } + + + // + // Get results of the fits + // + fHistON = fFindsigON.GetHist(); // don't use the original hist to draw + fHistOFF = fFindsigOFF.GetHist(); + // From ON-Data + fNon = fFindsigON.GetNon(); + fdNon = fFindsigON.GetdNon(); + // From OFF-Data + fNbg = fFindsigOFF.GetNon(); // counts in signal reg. OFF hist NO + fdNbg = fFindsigOFF.GetdNon(); //normalized !!! + // From OFF-Data polynomial fit + fNbgFit = fFindsigOFF.GetNbg(); + fdNbgFit = fFindsigOFF.GetdNbg(); + fGamma = fFindsigOFF.GetGamma() * fNormFactor; //<----------- + fNoff = fFindsigOFF.GetNoff(); + + + // + // Compare the counter number of bg events, NbgMea, and the fitted, NbgFit + // + if (fabs(fNbg - fNbgFit) > 0.1 * fNbgFit) + { + *fLog << err << dbginf << "number of OFF events and Fitted number of OFF events in signal region do not agree (within 10 %)" << endl; + + *fLog << "fNbg = " << fNbg << " ; fNbgFit = " << fNbgFit << endl; + } + + + + // --------------------------------------------------- + // + // calculate the number of excess events in the signal region + // + fNex = fNon - fNbg * fNormFactor; + fNexFit = fNon - fNbgFit * fNormFactor; + + + // + // Calculate Significance + // + // Significance computed using counted number of OFF-data events in signal + // region (fNbg) and normalization factor (fNormFactor). + // This is the strictly speaking the significance in Li&Ma paper... + // + if ( !SigmaLiMa(fNon, fNbg, fNormFactor, &fSigLiMa) ) + { + *fLog << err << dbginf << "SigmaLiMa failed" << endl; + return kFALSE; + } + // + // Significance computed using effective number of OFF events in signal + // region (fNoff) and gamma factor (fGama). + // This is Wolfgang approach to the calulation of significance + // using Li&Ma formula and estimated OFF events from polynomial fit. + // + if ( !SigmaLiMa(fNon, fNoff, fGamma , &fSigLiMaFit) ) + { + *fLog << err << dbginf << "SigmaLiMaFit failed" << endl; + return kFALSE; + } + + + // + // calculate the error of the number of excess events + // using fitted quantities and counted quantities + // + fdNex = fNex / fSigLiMa; + fdNexFit = fNexFit / fSigLiMaFit; + + + + // + // Create histogram of OFF data scaled by the normalization factor + // + NormalizedHistOFF(); + + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Create a new histogram to store the hist of OFF data scale with by the +// normalization factor +// +void MHMyFindSignificanceONOFF::NormalizedHistOFF() +{ + + // + // We need the following factor in case that one of the histogram + // (ON or OFF) is rebin when calculatin the sifnigicance + // + const Double_t BinWidthRatioONOFF = fHistON->GetBinWidth(1)/fHistOFF->GetBinWidth(1) ; + + + // + // Clone OFF histogram and scale it + // + fHistOFFNormalized = (TH1*)fHistOFF->Clone(); + fHistOFFNormalized->SetName(Form("%s (Normalized)",fHistOFF->GetName())); + + fHistOFFNormalized->Scale(fNormFactor*BinWidthRatioONOFF); + + + // + // Scale the fit function of OFF data + // + TF1* fPolyOFFNormalized = (TF1*)fFindsigOFF.GetPoly()->Clone(); + + fPolyOFFNormalized->SetLineColor(2); + fPolyOFFNormalized->SetLineWidth(2); + + + const Int_t npar = fPolyOFFNormalized->GetNpar(); + Double_t Parameter = 0.0; + Double_t ParameterError = 0.0; + for(Int_t i=0; iGetParameter(i); + ParameterError = fNormFactor*BinWidthRatioONOFF *fPolyOFFNormalized->GetParError(i); + + fPolyOFFNormalized -> SetParameter(i, Parameter); + fPolyOFFNormalized -> SetParError(i,ParameterError); + } + + TList *funclistOFF = fHistOFFNormalized->GetListOfFunctions(); + funclistOFF->Clear(); // delete the functions cloned from fHistOFF + funclistOFF->Add(fPolyOFFNormalized); + + +} + + + +// -------------------------------------------------------------------------- +// +// Draw hte Alpha plot for the ON data, and on top of it, the normalized Alpha +// plot for the OFF data. +// +void MHMyFindSignificanceONOFF::Draw(Option_t* option) +{ + + if(!fHistON || !fHistOFF) + return; + + + + // + // If a canvas exists use it, otherwise Creates a new canvas + // + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + // pad->Divide(2,1); + + AppendPad(""); + + pad->cd(1); + + // + // Draw histograms and fit functions + // + fHistON->SetTitle(""); + fHistON->SetMarkerStyle(21); + fHistON->Draw("p"); + + fHistOFFNormalized->SetLineColor(2); + fHistOFFNormalized->SetMarkerStyle(25); + fHistOFFNormalized->SetMarkerColor(2); + fHistOFFNormalized->Draw("Psame"); + + TList *funclistON = fHistON->GetListOfFunctions(); + TF1* polyON = (TF1*)funclistON->FindObject("Poly"); + TF1* gaussON = (TF1*)funclistON->FindObject("PolyGauss"); + TF1* backgON = (TF1*)funclistON->FindObject("Backg"); + + if(gaussON) + { + gaussON->SetLineColor(4); + gaussON->SetLineWidth(4); + } + if(backgON) + funclistON->Remove(backgON); + if(polyON) + { + polyON->SetLineStyle(2); + polyON->SetLineColor(2); + polyON->SetLineWidth(2); + } + + + // + // Draw results onto the figure + // + DrawResults(option); + + + gPad->Modified(); + gPad->Update(); + + + + + // // +// // Draw the difference of ON-OFF +// // +// pad->cd(2); +// TH1D* h = (TH1D*)fHistON->Clone(); +// h->GetListOfFunctions()->Clear(); +// h->SetTitle("Alpha Plot ON-OFF"); +// h->Add(fHistOFFNormalized,-1.); +// h->SetFillColor(17); +// h->Draw("hist"); + +// TLine* line = new TLine(0, 0, 90, 0); +// line->Draw(); + +// DrawResults("short"); + +} + + + +// -------------------------------------------------------------------------- +// +// Draw the results on top the histograms pad +// +void MHMyFindSignificanceONOFF::DrawResults(Option_t* option) +{ + + TString opt(option); + + TPaveText *pt; + char tx[100]; + + + if(opt.Contains("short",TString::kIgnoreCase)) + { + pt= new TPaveText(0.24, 0.59, 0.86, 0.87, "NDC"); + + + // sprintf(tx, " Nbg measured (Normalized) = %8.1f #pm %8.1f", +// fNbg * fNormFactor, fdNbg * fNormFactor); +// pt->AddText(tx); + +// sprintf(tx, " Nex (ON - OFF measured) = %8.1f #pm %8.1f",fNex,fdNex); +// pt->AddText(tx); + +// Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0; +// sprintf(tx, " Significance = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMa, ratio); +// pt->AddText(tx); + // sprintf(tx, " Measured OFF background:"); +// TText *t7 = pt->AddText(tx); +// t7->SetTextSize(0.02); +// t7->SetTextColor(8); + +// sprintf(tx, "Results for |alpha|< %6.2f [\\circ] using fit to OFF:", +// fFindsigON.GetAlphasi()); +// TText *t6 = pt->AddText(tx); +// t6->SetTextSize(0.03); +// t6->SetTextColor(8); + + // sprintf(tx, "Noff = %6.1f", fNoff); +// pt->AddText(tx); + +// sprintf(tx, "Nbg = %8.1f", +// fNbgFit*fNormFactor); +// pt->AddText(tx); + +// sprintf(tx, "Nex = %8.1f", fNexFit); +// pt->AddText(tx); + + + sprintf(tx, "T_{obs} = 64 min."); + pt->AddText(tx); + + sprintf(tx, "Zenith Angle = 7-17 #circ"); + pt->AddText(tx); + + sprintf(tx, "Significance = %6.2f #sigma", fSigLiMaFit); + TText *t8 =pt->AddText(tx); + //t8->SetTextSize(0.1); + t8->SetTextColor(8); + } + + else if(opt.Contains("all",TString::kIgnoreCase)) + { + pt = new TPaveText(0.20, 0.33, 0.88, 0.90, "NDC"); + + // + // Fit to OFF data + // + sprintf(tx, "Results of polynomial fit to OFF (order %2d) :", fDegree); + TText *t1 = pt->AddText(tx); + t1->SetTextSize(0.03); + t1->SetTextColor(2); + + sprintf(tx, " (%6.2f< |alpha| <%6.2f [\\circ])", fAlphaminOFF, + fAlphamaxOFF); + pt->AddText(tx); + + sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f", + fFindsigOFF.GetChisq(), fFindsigOFF.GetNdf(), fFindsigOFF.GetProb()); + pt->AddText(tx); + + sprintf(tx, " Nbgtot(fit)= %8.1f #pm %8.1f", + fFindsigOFF.GetNbgtotFitted(), fFindsigOFF.GetdNbgtotFitted()); + pt->AddText(tx); + + sprintf(tx, " Nbgtot(meas) = %8.1f", fFindsigOFF.GetNbgtot()); + pt->AddText(tx); + + + // + // Results + // + sprintf(tx, "Results for |alpha|< %6.2f [\\circ] :", + fFindsigON.GetAlphasi()); + TText *t6 = pt->AddText(tx); + t6->SetTextSize(0.03); + t6->SetTextColor(8); + + sprintf(tx, " Non = %8.1f #pm %8.1f", fNon, fdNon); + pt->AddText(tx); + + + // PRINT INFORMATION ABOUT MEASURED QUANTITIES + sprintf(tx, " Measured OFF background:"); + TText *t7 = pt->AddText(tx); + t7->SetTextSize(0.02); + t7->SetTextColor(8); + + sprintf(tx, " Nbg(meas)*NormFactor = %8.1f #pm %8.1f", + fNbg*fNormFactor, fdNbg*fNormFactor); + pt->AddText(tx); + + sprintf(tx, " Nex (meas) = %8.1f #pm %8.1f", fNex, fdNex); + pt->AddText(tx); + + sprintf(tx," Normalization Factor (= Non/Noff) = %4.4f",fNormFactor); + pt->AddText(tx); + + Double_t ratio = fNbg>0.0 ? fNex/(fNbg*fNormFactor) : 0.0; + sprintf(tx, " Significance(mea) = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMa, ratio); + pt->AddText(tx); + + + // PRINT INFORMATION ABOUT FITTED QUANTITIES + sprintf(tx, " Fitted OFF background:"); + TText *t8 = pt->AddText(tx); + t8->SetTextSize(0.02); + t8->SetTextColor(8); + + sprintf(tx, " Nbg(fit)*NormFactor = %8.1f #pm %8.1f", + fNbgFit*fNormFactor, fdNbgFit*fNormFactor); + pt->AddText(tx); + + sprintf(tx, " Nex(fit) = %8.1f #pm %8.1f", fNexFit, fdNexFit); + pt->AddText(tx); + + sprintf(tx, " Effective Noff (i.e. fNoff) = %6.1f, Gamma = %4.4f", + fNoff, fGamma); + pt->AddText(tx); + + ratio = fNbgFit>0.0 ? fNexFit/(fNbgFit*fNormFactor):0.0; + sprintf(tx, " Significance(fit) = %6.2f, Nex/(Nbg*NormFactor) = %6.2f", fSigLiMaFit, ratio); + pt->AddText(tx); + + + // + // Fit to ON data (polynomial+Gaus) + // + sprintf(tx, "Results of (polynomial+Gauss) fit to ON:"); + TText *t9 = pt->AddText(tx); + t9->SetTextSize(0.03); + t9->SetTextColor(4); + + //sprintf(tx, " chi2 = %8.2f, Ndof = %4d, Prob = %6.2f", + //fGChisq, fGNdf, fGProb); + //pt->AddText(tx); + + sprintf(tx, " Sigma of Gauss = %8.1f #pm %8.1f [\\circ]", + fFindsigON.GetSigmaGauss(), fFindsigON.GetdSigmaGauss()); + pt->AddText(tx); + + sprintf(tx, " total no.of excess events = %8.1f #pm %8.1f", + fFindsigON.GetNexGauss(), fFindsigON.GetdNexGauss()); + pt->AddText(tx); + } + else + return; + + + pt->SetFillStyle(0); + pt->SetBorderSize(0); + pt->SetTextAlign(12); + + pt->Draw(); +} + + +// -------------------------------------------------------------------------- +// +// Set flag fRebin +// +void MHMyFindSignificanceONOFF::SetRebin(Bool_t b) +{ + fRebin = b; + + *fLog << "MHMyFindSignificanceONOFF::SetRebin; flag fRebin set to " + << (b? "kTRUE" : "kFALSE") << endl; +} + +// -------------------------------------------------------------------------- +// +// Set flag fReduceDegree +// +void MHMyFindSignificanceONOFF::SetReduceDegree(Bool_t b) +{ + fReduceDegree = b; + + *fLog << "MHMyFindSignificanceONOFF::SetReduceDegree; flag fReduceDegree set to " + << (b? "kTRUE" : "kFALSE") << endl; +} + +// --------------------------------------------------------------------------- +// +// Given an histogram, count the number of events within a given range. +// Is equivaltent to TH1::Integral(binlo,binhi) but making sure that the bins +// are fully contain in the interval (within a given tolerance). +// +Int_t MHMyFindSignificanceONOFF::CountEventsInRange(TH1* hist, Double_t alphalo, Double_t alphaup, Double_t* numevents, Double_t* numeventserror) +{ + + const Int_t nbins = hist->GetNbinsX(); + const Double_t binwidth = hist->GetBinWidth(1); + + Double_t nevt = 0; + Double_t nevtErr = 0; + + // + // Loop over all the bins + // + for(Int_t i=1; i<=nbins; i++) + { + Double_t xlo = hist->GetBinLowEdge(i); + Double_t xup = hist->GetBinLowEdge(i+1); + + // bin must be completely contained in the signal region + if( xlo >= (alphalo-fEps) && xup <= (alphaup+fEps)) + { + // check all bins all of same width + const Double_t width = fabs(xup-xlo); + if (fabs(width-binwidth) > fEps) + { + *fLog << err << GetName() << " alpha plot has variable binning, which is not allowed...exit." << endl; + return kFALSE; + } + + + nevt += hist->GetBinContent(i); + nevtErr += hist->GetBinError(i) * hist->GetBinError(i); + } + } + + nevtErr = TMath::Sqrt(nevtErr); + + + // Return results + *numevents = nevt; + *numeventserror = nevtErr; +} + + +// --------------------------------------------------------------------------- +// +// Counts the number of events in the Bkg region in the ON_Data and OFF-Data +// histograms. Then compute the normalization factor as: +// NormFactor = NoffON / NoffNOFF +// +Int_t MHMyFindSignificanceONOFF::CalcNormalizationFactor(TH1* hon, TH1* hoff, Double_t alphalo, Double_t alphaup) +{ + Double_t nevton = 0; + Double_t nevtoff = 0; + Double_t dnevton = 0; + Double_t dnevtoff = 0; + + // counts events in the region (alphalo,alphaup) for the ON-Data hist + CountEventsInRange(hon, alphalo, alphaup,&nevton,&dnevton); + + // counts events in the region (alphalo,alphaup) for the OFF-Data hist + CountEventsInRange(hoff, alphalo, alphaup,&nevtoff,&dnevtoff); + + + // Normalization factor and its error + fNormFactor = nevton/nevtoff; + + fNormFactorError = 1./nevton + 1./nevtoff; + fNormFactorError = TMath::Sqrt(fNormFactorError); + fNormFactorError *= fNormFactor; // <-------------- + + + *fLog << "non = " << nevton << " noff = " << nevtoff << endl; + *fLog << "Normalization factor = " << fNormFactor << " +- " + << fNormFactorError << endl; + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// SigmaLiMa +// +// calculates the significance according to Li & Ma +// ApJ 272 (1983) 317 +// +Bool_t MHMyFindSignificanceONOFF::SigmaLiMa(Double_t non, Double_t noff, + Double_t gamma, Double_t *siglima) +{ + if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0) + { + *siglima = 0.0; + return kFALSE; + } + + Double_t help1 = non * log( (1.0+gamma)*non / (gamma*(non+noff)) ); + Double_t help2 = noff * log( (1.0+gamma)*noff / ( non+noff ) ); + *siglima = sqrt( 2.0 * (help1+help2) ); + + Double_t nex = non - gamma*noff; + if (nex < 0.0) + *siglima = - *siglima; + + //*fLog << "MHMyFindSignificanceONOFF::SigmaLiMa; non, noff, gamma, *siglima = " + // << non << ", " << noff << ", " << gamma << ", " << *siglima << endl; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Sigma Vs Alpha +// +// +Bool_t MHMyFindSignificanceONOFF::SigmaVsAlpha(TH1 *histON, TH1 *histOFF, Double_t alphamin, Double_t alphamax, Int_t degree, Bool_t draw) +{ + + + //-------------------------------------------- + // Create histogram + + //Int_t nsteps = 15; + Int_t nsteps = histON->FindBin(40); + + TH1D* fSigVsAlpha = new TH1D("SigVsAlpha","Sigma vs Alpha", nsteps, 0.0, 30); + MH::SetBinning(fSigVsAlpha, histON->GetXaxis()); + fSigVsAlpha->SetXTitle("upper edge of signal region in |alpha| [\\circ]"); + fSigVsAlpha->SetYTitle("Significance of gamma signal"); + + + TGraph* grSigVsAlpha = new TGraph; + TGraph* grNexVsAlpha = new TGraph; + + + + + //-------------------------------------------- + // loop over different signal regions + for (Int_t i=1; i<=nsteps; i++) + { + //Double_t alphasig = fSigVsAlpha->GetBinCenter(i); + //Double_t alphasig = fSigVsAlpha->GetBinLowEdge(i+1); + Double_t alphasig = histON->GetBinLowEdge(i+1); + + + gLog.SetNullOutput(1); + //FindSigmaONOFF(histON, histOFF, alphasig, alphamin, alphamax, degree); + MHMyFindSignificanceONOFF findsig; + findsig.SetRebin(fRebin); + findsig.SetReduceDegree(fReduceDegree); + + findsig.FindSigmaONOFF( histON, histOFF, alphasig, alphamin, alphamax, degree); + const Double_t sig = findsig.GetSigLiMa(); + const Double_t Nex = findsig.GetNex(); + + fSigVsAlpha->SetBinContent(i, sig); + + if( sig>fBestSigma ) + { + fBestSigma = sig; + fBestAlphaCut = alphasig; + fBestExcess = Nex; + } + + grSigVsAlpha->SetPoint(i-1,alphasig,sig); + grNexVsAlpha->SetPoint(i-1,alphasig,Nex); + + cout << "****************** " << i<< " " << alphasig << " " << sig<< " " << Nex << "*************"<< endl; + gLog.SetNullOutput(0); + } + + +// fSigVsAlpha->DrawCopy(); + + + if (!draw) + return kTRUE; + + // ----------------------------------------------------------------- + // + // Draw + // + TCanvas *c = new TCanvas("SigVsAlpha", "Sigma vs Alpha", 600, 600); + c->SetGridx(); + c->SetGridy(); + + Double_t margin = 0.12; // give more space for drawing the axis titles + c->SetLeftMargin(margin); + c->SetRightMargin(margin); + + + // + // Draw Sig Vs Alpha + // + grSigVsAlpha->SetMarkerStyle(8); + //grSigVsAlpha->SetMarkerColor(2); + //grSigVsAlpha->SetLineColor(2); + grSigVsAlpha->Draw("apl"); + grSigVsAlpha->GetHistogram()->SetXTitle("Alpha Cut [\\circ]"); + grSigVsAlpha->GetHistogram()->SetYTitle("Significance"); + grSigVsAlpha->GetHistogram()->GetYaxis()->SetTitleOffset(1.3); + + + // + // Draw second graph on a transparent pad + // + TPad *overlay = new TPad("overlay","",0,0,1,1); + overlay->SetFillStyle(4000); + overlay->SetFillColor(0); + overlay->SetFrameFillStyle(4000); + overlay->SetLeftMargin(margin); + overlay->SetRightMargin(margin); + overlay->Draw(); + overlay->cd(); + + grNexVsAlpha->SetMarkerStyle(8); + grNexVsAlpha->SetMarkerColor(kRed); + grNexVsAlpha->SetLineColor(kRed); + grNexVsAlpha->Draw("apl"); + + // Trick to not draw the graph Y-Axis + grNexVsAlpha->GetHistogram()->GetYaxis()->SetNdivisions(0); + + + // + // Draw an axis on the right side + // + c->Update(); //<-- Imprescindible for updating the gPad coordinates + + TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(), + gPad->GetUxmax(), gPad->GetUymax(), + gPad->GetUymin(),gPad->GetUymax() ,510,"+L"); + axis->SetTitle("Excess events"); + axis->SetTitleOffset(1.3); + //axis->SetLineColor(kRed); + axis->SetLabelColor(kRed); + axis->SetTextColor(kRed); + axis->Draw(); + + + c->Modified(); + c->Update(); + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.h (revision 6634) +++ /trunk/MagicSoft/Mars/mtemp/mucm/classes/MHMyFindSignificanceONOFF.h (revision 6634) @@ -0,0 +1,138 @@ +#ifndef MARS_MHMyFindSignificanceONOFF +#define MARS_MHMyFindSignificanceONOFF + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef MARS_MHFindSignificance +#include "MHFindSignificance.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class TF1; +class TH1; +class TCanvas; + +class MHMyFindSignificanceONOFF : public MH +{ + private: + + TH1 *fHistOrigON; // original plot of |alpha| (0.0 to 90.0 degrees) + TH1 *fHistON; // copy of fHistOrig or rebinned histogram + TH1 *fHistOrigOFF; // original plot of |alpha| (0.0 to 90.0 degrees) + TH1 *fHistOFF; // copy of fHistOrig or rebinned histogram + TH1 *fHistOFFNormalized; // fHistOFF normalized (contents and errors) with + + MHFindSignificance fFindsigON; // Will contain the fit of the on data + MHFindSignificance fFindsigOFF; + + + Double_t fAlphasig; + Double_t fAlphaminON; + Double_t fAlphamaxON; + Double_t fAlphaminOFF; + Double_t fAlphamaxOFF; + Int_t fDegree; + + Double_t fNormFactor; // Normalization factor between ON and OFF hists. + Double_t fNormFactorError; + + // Couted quatities + Double_t fNon; + Double_t fdNon; + Double_t fNbg; // Counts in the signal region of OFF hist, NO normalized !!! + Double_t fdNbg; + Double_t fNex; + Double_t fdNex; + Double_t fSigLiMa; // significance of gamma signal according to Li & Ma + + // from the polynomial fit of the OFF data + Double_t fNbgFit; + Double_t fdNbgFit; + Double_t fNexFit; + Double_t fdNexFit; + Double_t fGamma; + Double_t fNoff; + Double_t fSigLiMaFit; // significance of gamma signal according to Li & Ma + + Bool_t fRebin; // if true : allow rebinning of the alpha plot + Bool_t fReduceDegree; // if true : allow reducing of the order of the polynomial + + const static Double_t fEps = 1.e-4; // tolerance for floating point + // comparisons + + + Double_t fBestSigma; + Double_t fBestAlphaCut; + Double_t fBestExcess; + + void NormalizedHistOFF(); + void DrawResults(Option_t* option="all"); + +public: + + MHMyFindSignificanceONOFF(const char *name=NULL, const char *title=NULL); + + Bool_t FindSigmaONOFF(TH1 *histON, TH1 *histOFF, Double_t alphasig, Double_t alphamin, Double_t alphamax, Int_t degree); + + Int_t CountEventsInRange(TH1* hist, Double_t alphalo, Double_t alphaup, Double_t* numevents, Double_t* numeventserror); + + Int_t CalcNormalizationFactor(TH1* hon, TH1* hoff, Double_t alphalo, Double_t alphaup); + + Bool_t SigmaLiMa(Double_t non, Double_t noff, Double_t gamma, + Double_t *siglima); + + + Double_t GetNormFactor() const { return fNormFactor; } + Double_t GetNon() const { return fNon; } + Double_t GetNbg() const { return fNbg*fNormFactor; } // <------!!! + Double_t GetNex() const { return fNex; } + Double_t GetdNex() const { return fdNex; } + + Double_t GetSigLiMa() const { return fSigLiMa; } + Double_t GetNbgFit() const { return fNbgFit*fNormFactor; } //<------!!! + Double_t GetNexFit() const { return fNexFit; } + Double_t GetSigLiMaFit() const { return fSigLiMaFit; } + + + void SetRebin(Bool_t b=kTRUE); + void SetReduceDegree(Bool_t b=kTRUE); + + + void Draw(Option_t* option="all"); + + Bool_t SigmaVsAlpha(TH1 *histON, TH1 *histOFF, Double_t alphamin, Double_t alphamax, Int_t degree, Bool_t draw=kTRUE); + + + + Double_t GetBestSigma() const {return fBestSigma;} + Double_t GetBestAlphaCut() const {return fBestAlphaCut;} + Double_t GetBestExcess() const {return fBestExcess;} + + ClassDef(MHMyFindSignificanceONOFF, 1) // Determine significance from ON and OFF data +}; + +#endif + + + + + + + + + + + + + + + + + + + +