#ifndef MARS_MAlphaFitter #define MARS_MAlphaFitter #ifndef MARS_MParContainer #include "MParContainer.h" #endif #ifndef ROOT_TArrayD #include #endif #ifndef ROOT_TF1 #include #endif class TH1D; class TH3D; class MAlphaFitter : public MParContainer { public: enum ScaleMode_t { kNone, // No scaling kEntries, // scale by the number of entries in on and off kIntegral, // scale by the integral in on and off kOffRegion, // scale by the integral between fScaleMin, fScaleMax in on and off kBackground, // scale by the integral between fBgMin, fBgMax in on and off kLeastSquare, // not yet implemented kUserScale // scale by fixed factor set by SetScaleUser }; enum Strategy_t { kSignificance, kSignificanceChi2, kSignificanceLogExcess, kSignificanceExcess, kExcess, kGaussSigma, kWeakSource }; enum SignalFunc_t { kGauss, kThetaSq }; private: // Fitting Setup Float_t fSigInt; // minimum of range to fit the signal Float_t fSigMax; // maximum of range to fit the signal Float_t fBgMin; // minimum of range to fit the background Float_t fBgMax; // minimum of range to fit the background Float_t fScaleMin; // minimum of range to determin the scale factor of the background Float_t fScaleMax; // maximum of range to determin the scale factor of the background Int_t fPolynomOrder; // order of polyom to be fitted to the background Bool_t fFitBackground; // Backround fit: yes/no SignalFunc_t fSignalFunc; // Type of signal function // Result Double_t fSignificance; // significance of an unknown signal (Li/Ma 17) Double_t fSignificanceExc; // significance of a known excess (Li/Ma 5) Double_t fEventsExcess; // calculated number of excess events (signal-bg) Double_t fEventsSignal; // calculated number of signal events Double_t fEventsBackground; // calculated number of bg events (fScaleFactor already applied) Double_t fChiSqSignal; // Reduced (chi^2/NDF) chisq of signal fit Double_t fChiSqBg; // Reduced (chi^2/NDF) chisq of bg fit Double_t fIntegralMax; // Calculated bin border to which it was integrated Double_t fScaleFactor; // Scale factor determined for off-data TArrayD fCoefficients; // Fit result TArrayD fErrors; // Fit errors // Function TF1 *fFunc; // fit function (gauss + polynom) // Scaling setup ScaleMode_t fScaleMode; // scaling mode Double_t fScaleUser; // user scale factor // 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: // Implementing the function yourself is only about 5% faster MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fFitBackground(kTRUE), fSignalFunc(kGauss), 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) { fName = name ? name : "MAlphaFitter"; fTitle = title ? title : "Fit alpha"; fFunc->SetName("Dummy"); gROOT->GetListOfFunctions()->Remove(fFunc); Clear(); } MAlphaFitter(const MAlphaFitter &f) : fFunc(0) { f.Copy(*this); } ~MAlphaFitter() { delete fFunc; } // TObject void Clear(Option_t *o=""); void Print(Option_t *o="") const; //*MENU* void Copy(TObject &o) const; // Setter void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; } void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; } void SetMinimizationStrategy(Strategy_t mode) { fStrategy = mode; } 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 SetScaleMin(Float_t s) { fScaleMin = s; } void SetScaleMax(Float_t s) { fScaleMax = s; } void SetPolynomOrder(Int_t s) { if (s==fPolynomOrder) return; fPolynomOrder = s; SetSignalFunction(fSignalFunc); } void SetSignalFunction(SignalFunc_t func) { delete fFunc; switch (func) { case kGauss: fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", fPolynomOrder)); break; case kThetaSq: // if (fPolynomOrder==0) // fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + pol0(3)"); // else // { if (fPolynomOrder>0) fPolynomOrder = 1; fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + expo(3)"); // } break; } fSignalFunc=func; 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; } // Getter Double_t GetSignalIntegralMax() const { return fSigInt; } 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 GetSignificanceExc() const { return fSignificanceExc; } Double_t GetChiSqSignal() const { return fChiSqSignal; } Double_t GetChiSqBg() const { return fChiSqBg; } Double_t GetScaleFactor() const { return fScaleFactor; } Double_t GetMinimizationValue() const; ScaleMode_t GetScaleMode() const { return fScaleMode; } 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; } const TArrayD &GetErrors() const { return fErrors; } Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; } Double_t CalcUpperLimit() const; // Interface to fit Bool_t Fit(TH1D &h, Bool_t paint=kFALSE); Bool_t Fit(const TH1D &on, const TH1D &off, Double_t alpha, Bool_t paint=kFALSE); Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE) { return off ? Fit(on, *off, alpha, paint) : Fit(on, paint); } Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE) { return off ? Fit(on, *off, 1, paint) : Fit(on, paint); } Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE) { const Double_t alpha = off ? Scale(*off, on) : 1; return off ? Fit(on, *off, alpha, paint) : Fit(on, paint); } Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE); Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE); Bool_t FitTheta(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE); //Bool_t FitTime(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE); Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE); Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE); Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE); //Bool_t FitTime(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE); Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE) { return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint); } Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE) { return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint); } Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE) { return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint); }/* Bool_t FitTime(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE) { return off ? FitTime(on, *off, bin, paint) : FitTime(on, bin, paint); }*/ Double_t Scale(TH1D &off, const TH1D &on) const; // Interface to result void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035, Bool_t draw=kFALSE) const; void DrawResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const { PaintResult(x, y, size, kTRUE); } // MTask Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); ClassDef(MAlphaFitter, 3) }; #endif