| 1 | #ifndef MARS_MAlphaFitter
|
|---|
| 2 | #define MARS_MAlphaFitter
|
|---|
| 3 |
|
|---|
| 4 | #ifndef MARS_MParContainer
|
|---|
| 5 | #include "MParContainer.h"
|
|---|
| 6 | #endif
|
|---|
| 7 |
|
|---|
| 8 | #ifndef ROOT_TArrayD
|
|---|
| 9 | #include <TArrayD.h>
|
|---|
| 10 | #endif
|
|---|
| 11 |
|
|---|
| 12 | #ifndef ROOT_TF1
|
|---|
| 13 | #include <TF1.h>
|
|---|
| 14 | #endif
|
|---|
| 15 |
|
|---|
| 16 | class TH1D;
|
|---|
| 17 | class TH3D;
|
|---|
| 18 |
|
|---|
| 19 | class MAlphaFitter : public MParContainer
|
|---|
| 20 | {
|
|---|
| 21 | public:
|
|---|
| 22 | enum ScaleMode_t {
|
|---|
| 23 | kNone, // No scaling
|
|---|
| 24 | kEntries, // scale by the number of entries in on and off
|
|---|
| 25 | kIntegral, // scale by the integral in on and off
|
|---|
| 26 | kOffRegion, // scale by the integral between fScaleMin, fScaleMax in on and off
|
|---|
| 27 | kBackground, // scale by the integral between fBgMin, fBgMax in on and off
|
|---|
| 28 | kLeastSquare, // not yet implemented
|
|---|
| 29 | kUserScale // scale by fixed factor set by SetScaleUser
|
|---|
| 30 | };
|
|---|
| 31 | enum Strategy_t {
|
|---|
| 32 | kSignificance,
|
|---|
| 33 | kSignificanceChi2,
|
|---|
| 34 | kSignificanceLogExcess,
|
|---|
| 35 | kSignificanceExcess,
|
|---|
| 36 | kExcess,
|
|---|
| 37 | kGaussSigma,
|
|---|
| 38 | kWeakSource,
|
|---|
| 39 | kWeakSourceLogExcess
|
|---|
| 40 | };
|
|---|
| 41 | enum SignalFunc_t {
|
|---|
| 42 | kGauss, kThetaSq
|
|---|
| 43 | };
|
|---|
| 44 |
|
|---|
| 45 | private:
|
|---|
| 46 | // Fitting Setup
|
|---|
| 47 | Float_t fSigInt; // minimum of range to fit the signal
|
|---|
| 48 | Float_t fSigMax; // maximum of range to fit the signal
|
|---|
| 49 | Float_t fBgMin; // minimum of range to fit the background
|
|---|
| 50 | Float_t fBgMax; // minimum of range to fit the background
|
|---|
| 51 | Float_t fScaleMin; // minimum of range to determin the scale factor of the background
|
|---|
| 52 | Float_t fScaleMax; // maximum of range to determin the scale factor of the background
|
|---|
| 53 | Int_t fPolynomOrder; // order of polyom to be fitted to the background
|
|---|
| 54 | Bool_t fFitBackground; // Backround fit: yes/no
|
|---|
| 55 | SignalFunc_t fSignalFunc; // Type of signal function
|
|---|
| 56 | // Result
|
|---|
| 57 | Double_t fSignificance; // significance of an unknown signal (Li/Ma 17)
|
|---|
| 58 | Double_t fErrorExcess; // Simple error propagation
|
|---|
| 59 | Double_t fEventsExcess; // calculated number of excess events (signal-bg)
|
|---|
| 60 | Double_t fEventsSignal; // calculated number of signal events
|
|---|
| 61 | Double_t fEventsBackground; // calculated number of bg events (fScaleFactor already applied)
|
|---|
| 62 |
|
|---|
| 63 | Double_t fChiSqSignal; // Reduced (chi^2/NDF) chisq of signal fit
|
|---|
| 64 | Double_t fChiSqBg; // Reduced (chi^2/NDF) chisq of bg fit
|
|---|
| 65 | Double_t fIntegralMax; // Calculated bin border to which it was integrated
|
|---|
| 66 | Double_t fScaleFactor; // Scale factor determined for off-data
|
|---|
| 67 |
|
|---|
| 68 | TArrayD fCoefficients; // Fit result
|
|---|
| 69 | TArrayD fErrors; // Fit errors
|
|---|
| 70 |
|
|---|
| 71 | // Function
|
|---|
| 72 | TF1 *fFunc; // fit function (gauss + polynom)
|
|---|
| 73 |
|
|---|
| 74 | // Scaling setup
|
|---|
| 75 | ScaleMode_t fScaleMode; // scaling mode
|
|---|
| 76 | Double_t fScaleUser; // user scale factor
|
|---|
| 77 |
|
|---|
| 78 | // Minimization strategy
|
|---|
| 79 | Strategy_t fStrategy; // How to calc minimization value
|
|---|
| 80 |
|
|---|
| 81 | Double_t DoOffFit(const TH1D &hon, const TH1D &hof, Bool_t paint);
|
|---|
| 82 | Bool_t FitOff(TH1D &h, Int_t paint);
|
|---|
| 83 | void FitResult(const TH1D &h);
|
|---|
| 84 |
|
|---|
| 85 | public:
|
|---|
| 86 | // Implementing the function yourself is only about 5% faster
|
|---|
| 87 | MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15),
|
|---|
| 88 | fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80),
|
|---|
| 89 | fPolynomOrder(2), fFitBackground(kTRUE), fSignalFunc(kGauss),
|
|---|
| 90 | fCoefficients(3+fPolynomOrder+1), fErrors(3+fPolynomOrder+1),
|
|---|
| 91 | fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)),
|
|---|
| 92 | fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
|
|---|
| 93 | {
|
|---|
| 94 | fName = name ? name : "MAlphaFitter";
|
|---|
| 95 | fTitle = title ? title : "Fit alpha";
|
|---|
| 96 |
|
|---|
| 97 | fFunc->SetName("Dummy");
|
|---|
| 98 | gROOT->GetListOfFunctions()->Remove(fFunc);
|
|---|
| 99 |
|
|---|
| 100 | Clear();
|
|---|
| 101 | }
|
|---|
| 102 |
|
|---|
| 103 | MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
|
|---|
| 104 | {
|
|---|
| 105 | f.Copy(*this);
|
|---|
| 106 | }
|
|---|
| 107 | ~MAlphaFitter()
|
|---|
| 108 | {
|
|---|
| 109 | delete fFunc;
|
|---|
| 110 | }
|
|---|
| 111 |
|
|---|
| 112 | // TObject
|
|---|
| 113 | void Clear(Option_t *o="");
|
|---|
| 114 | void Print(Option_t *o="") const; //*MENU*
|
|---|
| 115 | void Copy(TObject &o) const;
|
|---|
| 116 |
|
|---|
| 117 | // Setter
|
|---|
| 118 | void SetScaleUser(Float_t scale) { fScaleUser = scale; fScaleMode=kUserScale; }
|
|---|
| 119 | void SetScaleMode(ScaleMode_t mode) { fScaleMode = mode; }
|
|---|
| 120 | void SetMinimizationStrategy(Strategy_t mode) { fStrategy = mode; }
|
|---|
| 121 | void SetSignalIntegralMax(Float_t s) { fSigInt = s; }
|
|---|
| 122 | void SetSignalFitMax(Float_t s) { fSigMax = s; }
|
|---|
| 123 | void SetBackgroundFitMin(Float_t s) { fBgMin = s; }
|
|---|
| 124 | void SetBackgroundFitMax(Float_t s) { fBgMax = s; }
|
|---|
| 125 | void SetScaleMin(Float_t s) { fScaleMin = s; }
|
|---|
| 126 | void SetScaleMax(Float_t s) { fScaleMax = s; }
|
|---|
| 127 | void SetPolynomOrder(Int_t s)
|
|---|
| 128 | {
|
|---|
| 129 | if (s==fPolynomOrder)
|
|---|
| 130 | return;
|
|---|
| 131 |
|
|---|
| 132 | fPolynomOrder = s;
|
|---|
| 133 |
|
|---|
| 134 | SetSignalFunction(fSignalFunc);
|
|---|
| 135 | }
|
|---|
| 136 | void SetSignalFunction(SignalFunc_t func)
|
|---|
| 137 | {
|
|---|
| 138 | delete fFunc;
|
|---|
| 139 | switch (func)
|
|---|
| 140 | {
|
|---|
| 141 | case kGauss:
|
|---|
| 142 | fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", fPolynomOrder));
|
|---|
| 143 | break;
|
|---|
| 144 | case kThetaSq:
|
|---|
| 145 | // if (fPolynomOrder==0)
|
|---|
| 146 | // fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + pol0(3)");
|
|---|
| 147 | // else
|
|---|
| 148 | // {
|
|---|
| 149 | if (fPolynomOrder>0)
|
|---|
| 150 | fPolynomOrder = 1;
|
|---|
| 151 | fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + expo(3)");
|
|---|
| 152 | // }
|
|---|
| 153 | break;
|
|---|
| 154 | }
|
|---|
| 155 | fSignalFunc=func;
|
|---|
| 156 | fFunc->SetName("Dummy");
|
|---|
| 157 | gROOT->GetListOfFunctions()->Remove(fFunc);
|
|---|
| 158 |
|
|---|
| 159 | fCoefficients.Set(3+fPolynomOrder+1);
|
|---|
| 160 | fCoefficients.Reset();
|
|---|
| 161 |
|
|---|
| 162 | fErrors.Set(3+fPolynomOrder+1);
|
|---|
| 163 | fErrors.Reset();
|
|---|
| 164 | }
|
|---|
| 165 | void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; }
|
|---|
| 166 |
|
|---|
| 167 | void FixScale()
|
|---|
| 168 | {
|
|---|
| 169 | fScaleMode = kUserScale;
|
|---|
| 170 | fScaleUser = fScaleFactor;
|
|---|
| 171 | }
|
|---|
| 172 |
|
|---|
| 173 | // Getter
|
|---|
| 174 | Double_t GetSignalIntegralMax() const { return fSigInt; }
|
|---|
| 175 |
|
|---|
| 176 | Double_t GetEventsExcess() const { return fEventsExcess; }
|
|---|
| 177 | Double_t GetEventsSignal() const { return fEventsSignal; }
|
|---|
| 178 | Double_t GetEventsBackground() const { return fEventsBackground; }
|
|---|
| 179 |
|
|---|
| 180 | Double_t GetSignificance() const { return fSignificance; }
|
|---|
| 181 | Double_t GetErrorExcess() const { return fErrorExcess; }
|
|---|
| 182 | Double_t GetChiSqSignal() const { return fChiSqSignal; }
|
|---|
| 183 | Double_t GetChiSqBg() const { return fChiSqBg; }
|
|---|
| 184 | Double_t GetScaleFactor() const { return fScaleFactor; }
|
|---|
| 185 | Double_t GetMinimizationValue() const;
|
|---|
| 186 |
|
|---|
| 187 | ScaleMode_t GetScaleMode() const { return fScaleMode; }
|
|---|
| 188 |
|
|---|
| 189 | Double_t GetGausSigma() const { return fCoefficients[2]; }
|
|---|
| 190 | Double_t GetGausMu() const { return fCoefficients[1]; }
|
|---|
| 191 | Double_t GetGausA() const { return fCoefficients[0]; }
|
|---|
| 192 | Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
|
|---|
| 193 | const TArrayD &GetCoefficients() const { return fCoefficients; }
|
|---|
| 194 | const TArrayD &GetErrors() const { return fErrors; }
|
|---|
| 195 | Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; }
|
|---|
| 196 |
|
|---|
| 197 | Double_t CalcUpperLimit() const;
|
|---|
| 198 |
|
|---|
| 199 | // Interface to fit
|
|---|
| 200 | Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
|
|---|
| 201 | Bool_t Fit(const TH1D &on, const TH1D &off, Double_t alpha, Bool_t paint=kFALSE);
|
|---|
| 202 | Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE)
|
|---|
| 203 | {
|
|---|
| 204 | return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
|
|---|
| 205 | }
|
|---|
| 206 | Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
|
|---|
| 207 | {
|
|---|
| 208 | return off ? Fit(on, *off, 1, paint) : Fit(on, paint);
|
|---|
| 209 | }
|
|---|
| 210 | Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
|
|---|
| 211 | {
|
|---|
| 212 | const Double_t alpha = off ? Scale(*off, on) : 1;
|
|---|
| 213 | return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
|
|---|
| 214 | }
|
|---|
| 215 |
|
|---|
| 216 | Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE);
|
|---|
| 217 | Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 218 | Bool_t FitTheta(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 219 | //Bool_t FitTime(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 220 |
|
|---|
| 221 | Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE);
|
|---|
| 222 | Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 223 | Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 224 | //Bool_t FitTime(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
|
|---|
| 225 |
|
|---|
| 226 | Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE)
|
|---|
| 227 | {
|
|---|
| 228 | return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint);
|
|---|
| 229 | }
|
|---|
| 230 | Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
|---|
| 231 | {
|
|---|
| 232 | return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint);
|
|---|
| 233 | }
|
|---|
| 234 | Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
|---|
| 235 | {
|
|---|
| 236 | return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint);
|
|---|
| 237 | }/*
|
|---|
| 238 | Bool_t FitTime(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
|
|---|
| 239 | {
|
|---|
| 240 | return off ? FitTime(on, *off, bin, paint) : FitTime(on, bin, paint);
|
|---|
| 241 | }*/
|
|---|
| 242 |
|
|---|
| 243 | Double_t Scale(TH1D &off, const TH1D &on) const;
|
|---|
| 244 |
|
|---|
| 245 | Bool_t ApplyScaling(const TH3D &hon, TH3D &hof, UInt_t bin) const;
|
|---|
| 246 | Bool_t ApplyScaling(const TH3D &hon, TH3D &hof) const;
|
|---|
| 247 |
|
|---|
| 248 | // Interface to result
|
|---|
| 249 | void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035, Bool_t draw=kFALSE) const;
|
|---|
| 250 | void DrawResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const { PaintResult(x, y, size, kTRUE); }
|
|---|
| 251 |
|
|---|
| 252 | // MTask
|
|---|
| 253 | Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
|
|---|
| 254 |
|
|---|
| 255 | ClassDef(MAlphaFitter, 4)
|
|---|
| 256 | };
|
|---|
| 257 |
|
|---|
| 258 | #endif
|
|---|
