Index: /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5005) +++ /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5006) @@ -46,4 +46,6 @@ #include "MMath.h" +#include "MLogManip.h" + ClassImp(MAlphaFitter); @@ -86,10 +88,17 @@ Double_t bgmin=fBgMin; Double_t bgmax=fBgMax; - Byte_t polynom=fPolynom; + Byte_t polynom=fPolynomOrder; // 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()); + + if (fCoefficients.GetSize() != polynom+3) + { + fCoefficients.Set(polynom+3); + fCoefficients.Reset(); + } + else + func.SetParameters(fCoefficients.GetArray()); func.FixParameter(1, 0); @@ -154,6 +163,6 @@ h.Fit(&func, "NQI", "", 0, sigmax); - fChiSqSignal = func.GetChisquare()/func.GetNDF(); - fSigmaGaus = func.GetParameter(2); + fChiSqSignal = func.GetChisquare()/func.GetNDF(); + fCoefficients.Set(func.GetNpar(), func.GetParameters()); //const Bool_t ok = NDF>0 && chi2<2.5*NDF; @@ -167,16 +176,20 @@ } // ------------------------------------ - const Double_t s = func.Integral(0, sigint)/alphaw; + //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); + //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; + + fIntegralMax = h.GetBinLowEdge(bin+1); + fEventsBackground = func.Integral(0, fIntegralMax)/alphaw; + fEventsSignal = h.Integral(0, bin); + fEventsExcess = fEventsSignal-fEventsBackground; + fSignificance = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground); return kTRUE; @@ -186,6 +199,6 @@ { 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, + fSignificance, fSigInt, GetGausSigma(), + (int)fEventsExcess, fIntegralMax, fChiSqBg, fChiSqSignal)); Index: /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5005) +++ /trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5006) @@ -4,4 +4,8 @@ #ifndef MARS_MParContainer #include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayD +#include #endif @@ -15,15 +19,19 @@ Float_t fBgMin; Float_t fBgMax; - Int_t fPolynom; + Int_t fPolynomOrder; Double_t fSignificance; - Double_t fExcessEvents; + Double_t fEventsExcess; + Double_t fEventsSignal; + Double_t fEventsBackground; + Double_t fChiSqSignal; Double_t fChiSqBg; - Double_t fSigmaGaus; Double_t fIntegralMax; + TArrayD fCoefficients; + public: - MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynom(1) + MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1) { } @@ -38,22 +46,30 @@ MAlphaFitter &f = static_cast(o); - f.fSigInt = fSigInt; - f.fSigMax = fSigMax; - f.fBgMin = fBgMin; - f.fBgMax = fBgMax; - f.fPolynom = fPolynom; + f.fSigInt = fSigInt; + f.fSigMax = fSigMax; + f.fBgMin = fBgMin; + f.fBgMax = fBgMax; + f.fPolynomOrder = fPolynomOrder; } - void SetSignalIntegralMax(Float_t s) { fSigInt = s; } - void SetSignalFitMax(Float_t s) { fSigMax = s; } - void SetBackgroundFitMin(Float_t s) { fBgMin = s; } - void SetBackgroundFitMax(Float_t s) { fBgMax = s; } - void SetNumPolynom(Int_t s) { fPolynom = s; } + void SetSignalIntegralMax(Float_t s) { fSigInt = s; } + void SetSignalFitMax(Float_t s) { fSigMax = s; } + void SetBackgroundFitMin(Float_t s) { fBgMin = s; } + void SetBackgroundFitMax(Float_t s) { fBgMax = s; } + void SetPolynomOrder(Int_t s) { fPolynomOrder = s; } - Double_t GetExcessEvents() const { return fExcessEvents; } - Double_t GetSignificance() const { return fSignificance; } - Double_t GetChiSqSignal() const { return fChiSqSignal; } - Double_t GetChiSqBg() const { return fChiSqBg; } - Double_t GetSigmaGaus() const { return fSigmaGaus; } + Double_t GetEventsExcess() const { return fEventsExcess; } + Double_t GetEventsSignal() const { return fEventsSignal; } + Double_t GetEventsBackground() const { return fEventsBackground; } + + Double_t GetSignificance() const { return fSignificance; } + Double_t GetChiSqSignal() const { return fChiSqSignal; } + Double_t GetChiSqBg() const { return fChiSqBg; } + + Double_t GetGausSigma() const { return fCoefficients[2]; } + Double_t GetGausMu() const { return fCoefficients[1]; } + Double_t GetGausA() const { return fCoefficients[0]; } + Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; } + const TArrayD &GetCoefficients() const { return fCoefficients; } void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const; Index: /trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5005) +++ /trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5006) @@ -175,6 +175,6 @@ if (fit.Fit(*h)) { - fHEnergy.SetBinContent(i, fit.GetExcessEvents()); - fHEnergy.SetBinError(i, fit.GetExcessEvents()*0.2); + fHEnergy.SetBinContent(i, fit.GetEventsExcess()); + fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2); } delete h; @@ -194,6 +194,6 @@ if (fit.Fit(*h)) { - fHTheta.SetBinContent(i, fit.GetExcessEvents()); - fHTheta.SetBinError(i, fit.GetExcessEvents()*0.2); + fHTheta.SetBinContent(i, fit.GetEventsExcess()); + fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2); } delete h; @@ -334,8 +334,8 @@ // Fill histogram // - fHTime.SetBinContent(n+1, fit.GetExcessEvents()); - fHTime.SetBinError(n+1, fit.GetExcessEvents()*0.1); - - *fLog << all << *fTimeEffOn << ": " << fit.GetExcessEvents() << endl; + fHTime.SetBinContent(n+1, fit.GetEventsExcess()); + fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1); + + *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl; rebin = steps; @@ -353,6 +353,6 @@ { alpha = (*fMatrix)[fMap[0]]; - energy = 1000; //(*fMatrix)[fMap[1]]; - theta = 0; //(*fMatrix)[fMap[2]]; + energy = 1000; + theta = 0; } else @@ -375,42 +375,4 @@ fHAlphaTime.Fill(TMath::Abs(alpha), w); - /* - // N_gamma vs Energy and Theta - const Double_t Ngam = fHUnfold.GetBinContent(m,n); - // A_eff vs Energy and Theta - const Double_t Aeff = fHAeff.GetBinContent(m,n); - // T_eff vs Theta - const Double_t Effon = teff.GetBinContent(n); - - const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam; - const Double_t c2 = teff.GetBinError(n) /Effon; - const Double_t c3 = fHAeff.GetBinError(m,n) /Aeff; - - const Double_t cont = Ngam / (DeltaE * Effon * Aeff); - const Double_t dcont = sqrt(c1*c1 + c2*c2 + c3*c3); - - // - // Out of Range - // - const Bool_t oor = Ngam<=0 || DeltaE<=0 || Effon<=0 || Aeff<=0; - - if (oor) - *fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") "; - - if (Ngam<=0) - *fLog << " Ngam=0"; - if (DeltaE<=0) - *fLog << " DeltaE=0"; - if (Effon<=0) - *fLog << " Effon=0"; - if (Aeff<=0) - *fLog << " Aeff=0"; - - if (oor) - *fLog << endl; - - fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont); - fHFlux.SetBinError(m,n, oor ? 1e-20 : dcont*cont); - */ return kTRUE; } @@ -487,10 +449,12 @@ if (o==(TString)"energy") { - if (fHEnergy.GetEntries()>0) + /* + if (fHEnergy.GetEntries()>10) { gPad->SetLogx(); gPad->SetLogy(); } - FitEnergySpec(kTRUE); + FitEnergySpec(kTRUE);*/ + } Index: /trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5005) +++ /trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5006) @@ -131,5 +131,7 @@ #include "MAstroSky2Local.h" #include "MStatusDisplay.h" + #include "MMath.h" +#include "MAlphaFitter.h" #include "MBinning.h" @@ -847,5 +849,12 @@ delete h0; - TH1 *h=0; + MAlphaFitter fit; + fit.SetSignalIntegralMax(sigint); + fit.SetSignalFitMax(sigmax); + fit.SetBackgroundFitMin(bgmin); + fit.SetBackgroundFitMax(bgmax); + fit.SetPolynomOrder(polynom); + + TH1D *h=0; for (int ix=1; ix<=nx; ix++) for (int iy=1; iy<=ny; iy++) @@ -855,85 +864,29 @@ h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy); + if (h->GetBinContent(1)==0) + continue; + + + h->Scale(entries/h->GetEntries()); + + if (!fit.Fit(*h)) + continue; + const Double_t alpha0 = h->GetBinContent(1); - - // Check for the regios which is not filled... - if (alpha0==0) - continue; - - h->Scale(entries/h->GetEntries()); - if (alpha0>maxalpha0) maxalpha0=alpha0; - - // First fit a polynom in the off region - func.FixParameter(0, 0); - func.FixParameter(2, 1); - func.ReleaseParameter(3); - - for (int i=5; iFit(&func, "N0Q", "", bgmin, bgmax); - - h4a.Fill(func.GetChisquare()); - //h5a.Fill(func.GetProb()); - - //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1)); - //func.SetParLimits(2, 5, 90); - - func.ReleaseParameter(0); - //func.ReleaseParameter(1); - func.ReleaseParameter(2); - func.FixParameter(3, func.GetParameter(3)); - for (int i=5; iFit(&func, "N0Q", "", 0, sigmax); - - TArrayD p(func.GetNpar(), func.GetParameters()); - + // Fill results into some histograms - h0a.Fill(p[0]); - h0b.Fill(p[3]); - h1.Fill(p[1]); - h2.Fill(p[2]); + h0a.Fill(fit.GetGausA()); // gaus-A + h0b.Fill(fit.GetCoefficient(3)); // 3 + h1.Fill(fit.GetGausMu()); // mu + h2.Fill(fit.GetGausSigma()); // sigma-gaus if (polynom>1) - h3.Fill(p[5]); - h4b.Fill(func.GetChisquare()); - //h5b.Fill(func.GetProb()); - - // Implementing the integral as analytical function - // gives the same result in the order of 10e-5 - // and it is not faster at all... - //const Bool_t ok = /*p[0]>=0 && /*p[0]1.75;// && p[2]<88.5; - /* - if (p[0]<-fabs(alpha0-p[3])*7 && p[2]GetBin(ix, iy); @@ -945,14 +898,12 @@ h6.Fill(sig); - if (sig>maxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)maxs) { maxs = sig; maxx = ix; maxy = iy; - maxpar = p; + maxpar = fit.GetCoefficients(); } } - - *fLog << inf << "done." << endl; h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);