Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5012) @@ -86,24 +86,15 @@ Double_t sigint=fSigInt; Double_t sigmax=fSigMax; - Double_t bgmin=fBgMin; - Double_t bgmax=fBgMax; - Byte_t polynom=fPolynomOrder; + Double_t bgmin =fBgMin; + Double_t bgmax =fBgMax; - // Implementing the function yourself is only about 5% faster - TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90); - //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90); + //TF1 fFunc("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90); - if (fCoefficients.GetSize() != polynom+3) - { - fCoefficients.Set(polynom+3); - fCoefficients.Reset(); - } - else - func.SetParameters(fCoefficients.GetArray()); + //fFunc->SetParameters(fCoefficients.GetArray()); - func.FixParameter(1, 0); - func.FixParameter(4, 0); - func.SetParLimits(2, 0, 90); - func.SetParLimits(3, -1, 1); + fFunc->FixParameter(1, 0); + fFunc->FixParameter(4, 0); + fFunc->SetParLimits(2, 0, 90); + fFunc->SetParLimits(3, -1, 1); const Double_t alpha0 = h.GetBinContent(1); @@ -115,10 +106,10 @@ // First fit a polynom in the off region - func.FixParameter(0, 0); - func.FixParameter(2, 1); - func.ReleaseParameter(3); + fFunc->FixParameter(0, 0); + fFunc->FixParameter(2, 1); + fFunc->ReleaseParameter(3); - for (int i=5; iGetNpar(); i++) + fFunc->ReleaseParameter(i); // options : N do not store the function, do not draw @@ -126,34 +117,34 @@ // R use the range specified in the function range // Q quiet mode - h.Fit(&func, "NQI", "", bgmin, bgmax); + h.Fit(fFunc, "NQI", "", bgmin, bgmax); - fChiSqBg = func.GetChisquare()/func.GetNDF(); + fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF(); // ------------------------------------ if (paint) { - func.SetRange(0, 90); - func.SetLineColor(kRed); - func.SetLineWidth(2); - func.Paint("same"); + fFunc->SetRange(0, 90); + fFunc->SetLineColor(kRed); + fFunc->SetLineWidth(2); + fFunc->Paint("same"); } // ------------------------------------ - func.ReleaseParameter(0); + fFunc->ReleaseParameter(0); //func.ReleaseParameter(1); - func.ReleaseParameter(2); - func.FixParameter(3, func.GetParameter(3)); - for (int i=5; iReleaseParameter(2); + fFunc->FixParameter(3, fFunc->GetParameter(3)); + for (int i=5; iGetNpar(); i++) + fFunc->FixParameter(i, fFunc->GetParameter(i)); // Do not allow signals smaller than the background - const Double_t A = alpha0-func.GetParameter(3); + const Double_t A = alpha0-fFunc->GetParameter(3); const Double_t dA = TMath::Abs(A); - func.SetParLimits(0, -dA*4, dA*4); - func.SetParLimits(2, 0, 90); + fFunc->SetParLimits(0, -dA*4, dA*4); + fFunc->SetParLimits(2, 0, 90); // Now fit a gaus in the on region on top of the polynom - func.SetParameter(0, A); - func.SetParameter(2, sigmax*0.75); + fFunc->SetParameter(0, A); + fFunc->SetParameter(2, sigmax*0.75); // options : N do not store the function, do not draw @@ -161,8 +152,8 @@ // R use the range specified in the function range // Q quiet mode - h.Fit(&func, "NQI", "", 0, sigmax); + h.Fit(fFunc, "NQI", "", 0, sigmax); - fChiSqSignal = func.GetChisquare()/func.GetNDF(); - fCoefficients.Set(func.GetNpar(), func.GetParameters()); + fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF(); + fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters()); //const Bool_t ok = NDF>0 && chi2<2.5*NDF; @@ -171,16 +162,14 @@ if (paint) { - func.SetLineColor(kGreen); - func.SetLineWidth(2); - func.Paint("same"); + fFunc->SetLineColor(kGreen); + fFunc->SetLineWidth(2); + fFunc->Paint("same"); } // ------------------------------------ - //const Double_t s = func.Integral(0, sigint)/alphaw; - func.SetParameter(0, 0); - func.SetParameter(2, 1); - //const Double_t b = func.Integral(0, sigint)/alphaw; - - //fSignificance = MMath::SignificanceLiMaSigned(s, b); - //exc = s-b; + const Double_t s = fFunc->Integral(0, sigint)/alphaw; + fFunc->SetParameter(0, 0); + fFunc->SetParameter(2, 1); + const Double_t b = fFunc->Integral(0, sigint)/alphaw; + fSignificance = MMath::SignificanceLiMaSigned(s, b); const Double_t uin = 1.25*sigint; @@ -188,8 +177,8 @@ fIntegralMax = h.GetBinLowEdge(bin+1); - fEventsBackground = func.Integral(0, fIntegralMax)/alphaw; + fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw; fEventsSignal = h.Integral(0, bin); fEventsExcess = fEventsSignal-fEventsBackground; - fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground); + //fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground); return kTRUE; Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5012) @@ -10,4 +10,8 @@ #endif +#ifndef ROOT_TF1 +#include +#endif + class TH1D; @@ -15,4 +19,6 @@ { private: + TF1 *fFunc; + Float_t fSigInt; Float_t fSigMax; @@ -33,11 +39,16 @@ public: - MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1) + // Implementing the function yourself is only about 5% faster + MAlphaFitter() : fFunc(new TF1("", "gaus(0) + pol1(3)", 0, 90)), fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1), fCoefficients(3+fPolynomOrder+1) { } - MAlphaFitter(const MAlphaFitter &f) + MAlphaFitter(const MAlphaFitter &f) : fFunc(0) { f.Copy(*this); + } + ~MAlphaFitter() + { + delete fFunc; } @@ -51,4 +62,10 @@ f.fBgMax = fBgMax; f.fPolynomOrder = fPolynomOrder; + f.fCoefficients.Set(fCoefficients.GetSize()); + f.fCoefficients.Reset(); + + TF1 *fcn = f.fFunc; + f.fFunc = new TF1(*fFunc); + delete fcn; } @@ -57,5 +74,5 @@ void SetBackgroundFitMin(Float_t s) { fBgMin = s; } void SetBackgroundFitMax(Float_t s) { fBgMax = s; } - void SetPolynomOrder(Int_t s) { fPolynomOrder = s; } + void SetPolynomOrder(Int_t s) { fPolynomOrder = s; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", s)); fCoefficients.Set(3+s+1); fCoefficients.Reset(); } Double_t GetEventsExcess() const { return fEventsExcess; } @@ -76,4 +93,5 @@ Bool_t Fit(TH1D &h, Bool_t paint=kFALSE); + ClassDef(MAlphaFitter, 1) }; Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5012) @@ -140,12 +140,19 @@ void MHAlpha::FitEnergySpec(Bool_t paint) { - TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]"); - f.SetParameter(0, -2.0); - f.SetParameter(1, 500); // 50 - f.SetParameter(2, 1); // 50 + /* + TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");// + f.SetParameter(0, -2.1); + f.SetParameter(1, 1400); // 50 + f.SetParameter(2, fHEnergy.Integral()); // 50 f.SetLineWidth(2); - f.SetLineColor(kGreen); - - fHEnergy.Fit(&f, "0QI", "", 90, 1900); // Integral? + + fHEnergy.Fit(&f, "0QI", "", 100, 2400); // Integral? 90, 2800 + + const Double_t chi2 = f.GetChisquare(); + const Int_t NDF = f.GetNDF(); + + // was fit successful ? + const Bool_t ok = NDF>0 && chi2<2.5*NDF; + f.SetLineColor(ok ? kGreen : kRed); if (paint) @@ -153,12 +160,17 @@ f.Paint("same"); - TLatex text(0.2, 0.94, Form("\\alpha=%.1f E_{0}=%.1fGeV N=%.1f", - f.GetParameter(0), - f.GetParameter(1), - 0/*f.GetParameter(2)*/)); + if (ok) + { + TString s; + s += Form("\\alpha=%.1f ", f.GetParameter(0)); + s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000); + s += Form("N=%.1e", f.GetParameter(2)); + + TLatex text(0.25, 0.94, s.Data()); text.SetBit(TLatex::kTextNDC); text.SetTextSize(0.04); text.Paint(); - } + } + }*/ } @@ -449,5 +461,4 @@ if (o==(TString)"energy") { - /* if (fHEnergy.GetEntries()>10) { @@ -455,5 +466,5 @@ gPad->SetLogy(); } - FitEnergySpec(kTRUE);*/ + FitEnergySpec(kTRUE); } @@ -525,150 +536,4 @@ } -// -------------------------------------------------------------------------- -// -// This is a preliminary implementation of a alpha-fit procedure for -// all possible source positions. It will be moved into its own -// more powerfull class soon. -// -// The fit function is "gaus(0)+pol2(3)" which is equivalent to: -// [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2 -// or -// A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2 -// -// Parameter [1] is fixed to 0 while the alpha peak should be -// symmetric around alpha=0. -// -// Parameter [4] is fixed to 0 because the first derivative at -// alpha=0 should be 0, too. -// -// In a first step the background is fitted between bgmin and bgmax, -// while the parameters [0]=0 and [2]=1 are fixed. -// -// In a second step the signal region (alphaGetBinWidth(1); - - // Check for the regios which is not filled... - if (alpha0==0) - return kFALSE; //*fLog << warn << "Histogram empty." << endl; - - // First fit a polynom in the off region - func.FixParameter(0, 0); - func.FixParameter(2, 1); - func.ReleaseParameter(3); - - for (int i=5; i0 && chi2<2.5*NDF; - - // ------------------------------------ - if (paint) - { - func.SetLineColor(kGreen); - func.SetLineWidth(2); - func.Paint("same"); - } - // ------------------------------------ - const Double_t s = func.Integral(0, sigint)/alphaw; - func.SetParameter(0, 0); - func.SetParameter(2, 1); - const Double_t b = func.Integral(0, sigint)/alphaw; - - fSignificance = MMath::SignificanceLiMaSigned(s, b); - //exc = s-b; - - const Double_t uin = 1.25*sigint; - const Int_t bin = h.GetXaxis()->FindFixBin(uin); - fIntegralMax = h.GetBinLowEdge(bin+1); - fExcessEvents = h.Integral(0, bin)-func.Integral(0, fIntegralMax)/alphaw; - - return kTRUE; -} - -void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const -{ - TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ) \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}=%.1f \\chi_{s}^{2}=%.1f)", - fSignificance, fSigInt, fSigmaGaus, - (int)fExcessEvents, fIntegralMax, - fChiSqBg, fChiSqSignal)); - - text.SetBit(TLatex::kTextNDC); - text.SetTextSize(size); - text.Paint(); -} -*/ Bool_t MHAlpha::Finalize() { Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 5012) @@ -84,5 +84,5 @@ { private: - MAlphaFitter fFit; + MAlphaFitter fFit; //! SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT) TH3D fHAlpha; // Alpha vs. theta and energy Index: trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 5012) @@ -188,5 +188,5 @@ fHThetaLambda.SetTitle("Slope (Rate) vs Theta"); fHThetaLambda.SetXTitle("\\Theta [\\circ]"); - fHThetaLambda.SetYTitle("\\labda"); + fHThetaLambda.SetYTitle("\\lambda [s^{-1}]"); fHThetaLambda.UseCurrentStyle(); fHThetaLambda.SetDirectory(NULL); @@ -197,5 +197,5 @@ fHTimeLambda.SetTitle("Slope (Rate) vs Time"); fHTimeLambda.SetXTitle("\\Time [\\circ]"); - fHTimeLambda.SetYTitle("\\lambda"); + fHTimeLambda.SetYTitle("\\lambda [s^{-1}]"); fHTimeLambda.UseCurrentStyle(); fHTimeLambda.SetDirectory(NULL); @@ -318,5 +318,5 @@ // parameter 1 = N0*del // - static TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)"); + TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)"); //func.SetParNames("lambda", "N0", "del"); @@ -791,5 +791,5 @@ AppendPad("time"); fHTimeLambda.Draw("same"); - DrawRightAxis("\\lambda"); + DrawRightAxis("\\lambda [s^{-1}]"); pad->cd(2); @@ -821,4 +821,4 @@ AppendPad("theta"); fHThetaLambda.Draw("same"); - DrawRightAxis("\\lambda"); -} + DrawRightAxis("\\lambda [s^{-1}]"); +} Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5012) @@ -808,10 +808,5 @@ const Double_t w = fHist.GetZaxis()->GetBinWidth(1); - // xmin, xmax, npar - //TF1 func("MyFunc", fcn, 0, 90, 6); - // Implementing the function yourself is only about 5% faster - TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90); - //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90); - TArrayD maxpar(func.GetNpar()); + TArrayD maxpar; /* func.SetParName(0, "A"); @@ -819,9 +814,4 @@ * func.SetParName(2, "sigma"); */ - - func.FixParameter(1, 0); - func.FixParameter(4, 0); - func.SetParLimits(2, 0, 90); - func.SetParLimits(3, -1, 1); const Int_t nx = hist->GetXaxis()->GetNbins(); @@ -864,5 +854,5 @@ h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy); - if (h->GetBinContent(1)==0) + if (h->GetEntries()==0) continue; @@ -988,6 +978,6 @@ result->Draw(); - TF1 f1("", func.GetExpFormula(), 0, 90); - TF1 f2("", func.GetExpFormula(), 0, 90); + TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90); + TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90); f1.SetParameters(maxpar.GetArray()); f2.SetParameters(maxpar.GetArray()); Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 5011) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 5012) @@ -61,5 +61,5 @@ TH1 *GetHistByName(const TString name) { return &fHist; } - void FitSignificance(Float_t sigint=15, Float_t sigmax=70, Float_t bgmin=40, Float_t bgmax=70, Byte_t polynom=1); //*MENU* + void FitSignificance(Float_t sigint=10, Float_t sigmax=75, Float_t bgmin=45, Float_t bgmax=85, Byte_t polynom=2); //*MENU* void FitSignificanceStd() { FitSignificance(); } //*MENU*