Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7719) +++ trunk/MagicSoft/Mars/Changelog (revision 7720) @@ -43,4 +43,12 @@ * mreport/MReportFileReadCC.cc: - always output the file format version + + * mhflux/MAlphaFitter.[h,cc]: + - implemented fitting of the off-data for background determination + - the result values are not yet accessible + + * mhflux/MHThetaSqN.[h,cc]: + - some improvements regarding the treatment of the signal-regions + - allow to set a different cut-level for off-cuts Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 7719) +++ trunk/MagicSoft/Mars/NEWS (revision 7720) @@ -4,6 +4,15 @@ - merpp: better handling of problems with the TH and TD part of the - CC-REPORT for files older than 200507190 and 200412210 - respectively + CC-REPORT for files older than 200507190 and 200412210 + respectively + + - ganymed: implemented a new class (MHThetaSqN) which allows to use more + than one off-source region in wobble-mode. To use it add the following to + your ganymed_wobble.rc: + + MJCut.NameHist: MHThetaSqN (switch the new feature on) + + MHThetaSqN.NumOffSourcePos: 3 (define the number of off-regions) + + MHThetaSqN.DoOffCut: Yes,No (switch on/off the off-cut) + + Cut1.ThetaCut: None + + MHThetaSqN.SignificanceCutLevel: 2.0 (increase off-cut by 2.0/1.7) Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 7719) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 7720) @@ -38,4 +38,8 @@ // + Double_t fSignificanceExc; // significance of a known excess // +// Version 3: +// ---------- +// + TArrayD fErrors; // errors of coefficients +// // ////////////////////////////////////////////////////////////////////////////// @@ -75,52 +79,15 @@ fCoefficients.Reset(); + fErrors.Reset(); } // -------------------------------------------------------------------------- // -// 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 @@ -138,4 +105,7 @@ fFunc->SetParameter(i, 0); + if (!fFitBackground) + return kTRUE; + // options : N do not store the function, do not draw // I use integral of function in bin rather than value at bin center @@ -143,64 +113,40 @@ // Q quiet mode // E Perform better Errors estimation using Minos technique - if (fFitBackground) - { - h.Fit(fFunc, "NQI", "", bgmin, bgmax); - fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF(); - } - + h.Fit(fFunc, "NQIE", "", fBgMin, fBgMax); + fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF(); + + fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters()); + fErrors.Set(fFunc->GetNpar()); + for (int i=3; iGetNpar(); i++) + fErrors[i] = fFunc->GetParError(i); // ------------------------------------ - if (paint && fFitBackground) - { - fFunc->SetRange(0, 90); - fFunc->SetLineColor(kRed); - fFunc->SetLineWidth(2); + + if (paint) + { + if (paint==2) + { + fFunc->SetLineColor(kBlack); + fFunc->SetLineWidth(1); + } + else + { + fFunc->SetRange(0, 90); + fFunc->SetLineColor(kRed); + fFunc->SetLineWidth(2); + } fFunc->Paint("same"); } - // ------------------------------------ - - fFunc->ReleaseParameter(0); // It is also released by SetParLimits later on - //func.ReleaseParameter(1); // It is also released by SetParLimits later on - fFunc->ReleaseParameter(2); - for (int i=3; iGetNpar(); i++) - fFunc->FixParameter(i, fFunc->GetParameter(i)); - - // Do not allow signals smaller than the background - const Double_t s = fSignalFunc==kGauss ? fFunc->GetParameter(3) : TMath::Exp(fFunc->GetParameter(3)); - const Double_t A = alpha0-s; - const Double_t dA = TMath::Abs(A); - fFunc->SetParLimits(0, -dA*4, dA*4); // SetParLimits also releases the parameter - fFunc->SetParLimits(2, 0, 90); // SetParLimits also releases the parameter - - // Now fit a gaus in the on region on top of the polynom - fFunc->SetParameter(0, A); - fFunc->SetParameter(2, sigmax*0.75); - - // options : N do not store the function, do not draw - // I use integral of function in bin rather than value at bin center - // R use the range specified in the function range - // Q quiet mode - // E Perform better Errors estimation using Minos technique - h.Fit(fFunc, "NQI", "", 0, sigmax); - - fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF(); - fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters()); - - //const Bool_t ok = NDF>0 && chi2<2.5*NDF; - - // ------------------------------------ - if (paint) - { - fFunc->SetLineColor(kGreen); - fFunc->SetLineWidth(2); - fFunc->Paint("same"); - } - // ------------------------------------ - - //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw; - fFunc->SetParameter(0, 0); - fFunc->SetParameter(2, 1); - //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw; - //fSignificance = MMath::SignificanceLiMaSigned(s, b); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the result of the fit and set the corresponding data members +// +void MAlphaFitter::FitResult(const TH1D &h) +{ + const Double_t alphaw = h.GetXaxis()->GetBinWidth(1); const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999); @@ -218,4 +164,166 @@ if (fEventsExcess<0) fEventsExcess=0; +} + +// -------------------------------------------------------------------------- +// +// 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 (alphaReleaseParameter(0); // It is also released by SetParLimits later on + //func.ReleaseParameter(1); // It is also released by SetParLimits later on + fFunc->ReleaseParameter(2); + for (int i=3; iGetNpar(); i++) + fFunc->FixParameter(i, fFunc->GetParameter(i)); + + + // Do not allow signals smaller than the background + const Double_t alpha0 = h.GetBinContent(1); + const Double_t s = fSignalFunc==kGauss ? fFunc->GetParameter(3) : TMath::Exp(fFunc->GetParameter(3)); + const Double_t A = alpha0-s; + //const Double_t dA = TMath::Abs(A); + //fFunc->SetParLimits(0, -dA*4, dA*4); // SetParLimits also releases the parameter + fFunc->SetParLimits(2, 0, 90); // SetParLimits also releases the parameter + + // Now fit a gaus in the on region on top of the polynom + fFunc->SetParameter(0, A); + fFunc->SetParameter(2, fSigMax*0.75); + + // options : N do not store the function, do not draw + // I use integral of function in bin rather than value at bin center + // R use the range specified in the function range + // Q quiet mode + // E Perform better Errors estimation using Minos technique + h.Fit(fFunc, "NQI", "", 0, fSigMax); + + fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF(); + fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters()); + for (int i=0; i<3; i++) + fErrors[i] = fFunc->GetParError(i); + //const Bool_t ok = NDF>0 && chi2<2.5*NDF; + + // ------------------------------------ + if (paint) + { + fFunc->SetLineColor(kGreen); + fFunc->SetLineWidth(2); + fFunc->Paint("same"); + } + // ------------------------------------ + + //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw; + fFunc->SetParameter(0, 0); + fFunc->SetParameter(2, 1); + //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw; + //fSignificance = MMath::SignificanceLiMaSigned(s, b); + + // Calculate the fit result and set the corresponding data members + FitResult(h); + + return kTRUE; +} + +Double_t MAlphaFitter::DoOffFit(const TH1D &hon, const TH1D &hof, Bool_t paint) +{ + if (fSignalFunc!=kThetaSq) + return 0; + + // ---------------------------------------------------------------------------- + + const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999); + + + MAlphaFitter fit(*this); + fit.EnableBackgroundFit(); + fit.SetBackgroundFitMin(0); + + // produce a histogram containing the off-samples from on-source and + // off-source in the off-source region and the on-data in the source-region + TH1D h(hof); + h.Add(&hon); + h.Scale(0.5); + for (int i=1; i<=bin+3; i++) + { + h.SetBinContent(i, hof.GetBinContent(i)); + h.SetBinError( i, hof.GetBinError(i)); + } + + // Now fit the off-data + if (!fit.FitOff(h, paint?2:0)) // FIXME: Show fit! + return -1; + + // Calculate fit-result + fit.FitResult(h); + + // Do a gaussian error propagation to calculated the error of + // the background estimated from the fit + const Double_t ea = fit.GetErrors()[3]; + const Double_t eb = fit.GetErrors()[4]; + const Double_t a = fit.GetCoefficients()[3]; + const Double_t b = fit.GetCoefficients()[4]; + + const Double_t t = fIntegralMax; + + const Double_t ex = TMath::Exp(t*b); + const Double_t eab = TMath::Exp(a)/b; + + const Double_t eA = ex-1; + const Double_t eB = t*ex - eA/b; + + const Double_t w = h.GetXaxis()->GetBinWidth(1); + + // Error of estimated background + const Double_t er = TMath::Abs(eab)*TMath::Hypot(eA*ea, eB*eb)/w; + + // Calculate arbitrary scale factor from propagated error from the + // condistion: sqrt(alpha*background) = est.background/est.error + // const Double_t bg = hof.Integral(1, bin); + // const Double_t sc = bg * er*er / (fit2.GetEventsBackground()*fit2.GetEventsBackground()); + // Assuming that bg and fit2.GetEventsBackground() are rather identical: + const Double_t sc = er*er / fit.GetEventsBackground(); + + return sc; + /* + cout << MMath::SignificanceLiMaSigned(hon.Integral(1, bin), fit.GetEventsBackground()/sc, sc) << " "; + cout << sc << " "; + cout << fit.fChiSqBg << endl; + */ return kTRUE; @@ -230,5 +338,4 @@ MAlphaFitter fit(*this); fit.SetPolynomOrder(0); - if (alpha<=0 || !fit.Fit(h, paint)) return kFALSE; @@ -237,4 +344,14 @@ fChiSqBg = fit.GetChiSqBg(); fCoefficients = fit.GetCoefficients(); + fErrors = fit.GetErrors(); + + + // ---------------------------------------------------------------------------- + + const Double_t scale = DoOffFit(hon, hof, paint); + if (scale<0) + return kFALSE; + + // ---------------------------------------------------------------------------- const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999); @@ -321,4 +438,6 @@ f.fCoefficients.Set(fCoefficients.GetSize()); f.fCoefficients.Reset(); + f.fErrors.Set(fCoefficients.GetSize()); + f.fErrors.Reset(); // Result Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 7719) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 7720) @@ -66,4 +66,5 @@ TArrayD fCoefficients; // Fit result + TArrayD fErrors; // Fit errors // Function @@ -76,4 +77,8 @@ // Minimization strategy Strategy_t fStrategy; // How to calc minimization value + + Double_t DoOffFit(const TH1D &hon, const TH1D &hof, Bool_t paint); + Bool_t FitOff(TH1D &h, Int_t paint); + void FitResult(const TH1D &h); public: @@ -82,5 +87,5 @@ fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fFitBackground(kTRUE), fSignalFunc(kGauss), - fCoefficients(3+fPolynomOrder+1), + fCoefficients(3+fPolynomOrder+1), fErrors(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance) @@ -150,6 +155,10 @@ fFunc->SetName("Dummy"); gROOT->GetListOfFunctions()->Remove(fFunc); + fCoefficients.Set(3+fPolynomOrder+1); fCoefficients.Reset(); + + fErrors.Set(3+fPolynomOrder+1); + fErrors.Reset(); } void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; } @@ -176,4 +185,5 @@ Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; } const TArrayD &GetCoefficients() const { return fCoefficients; } + const TArrayD &GetErrors() const { return fErrors; } Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; } @@ -232,5 +242,5 @@ Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); - ClassDef(MAlphaFitter, 2) + ClassDef(MAlphaFitter, 3) }; Index: trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc (revision 7719) +++ trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc (revision 7720) @@ -67,5 +67,5 @@ // MHThetaSqN::MHThetaSqN(const char *name, const char *title) - : MHAlpha(name, title), fDisp(0), fSrcPosCam(0), + : MHAlpha(name, title), fDisp(0), fSrcPosCam(0), fSignificanceCutLevel(1.7), fNumBinsSignal(3), fNumBinsTotal(75), fNumOffSourcePos(3), fDoOffCut(kTRUE) { @@ -175,4 +175,6 @@ fFit.SetScaleUser(1./fNumOffSourcePos); + fThetaSqCut = fSignificanceCutLevel*fFit.GetSignalIntegralMax()/1.7; + return kTRUE; } @@ -215,12 +217,147 @@ const Float_t rad = TMath::TwoPi()/n; - for (UInt_t i=0; i(src0).Rotate(i*rad); - const Double_t d = (src-org).Mod2(); - +/* + for (UInt_t i=0; i(src0).Rotate(i*rad); + SetVal((src-org).Mod2()); + + if (!MHAlpha::Fill(NULL, weight)) + return kFALSE; + } + */ + + // Calculate distance (theta-sq) to all (off-)source regions + TArrayD dist(n); + for (UInt_t i=0; i(src0).Rotate(i*rad); + dist[i] = (src-org).Mod2(); + } + + // Processing off-data + // Check if event's origin is in the on-regions + if (!fOffData && fDoOffCut && dist[0](src0).Rotate(i*rad); + dist[i] = (src-org).Mod2(); + } + + for (UInt_t i=0; i(src0).Rotate(i*rad+TMath::Pi()); + const Double_t d1 = (src1-org).Mod2(); + + if (d1(src0).Rotate(i*rad); + const Double_t d0 = (src-org).Mod2(); + + SetVal(d0); + + if (!MHAlpha::Fill(NULL, weight)) + return kFALSE; + } + */ + +/* + for (UInt_t i=0; i(src0).Rotate(i*rad); + const Double_t d = (src-org).Mod2(); + + // off: is in src region on: is in off regions + /// if (!fOffData && i==0) || (fOffData && i!=0) + if ((bool)fOffData ^ (i==0) ) { if (d