Ignore:
Timestamp:
02/07/09 20:40:28 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

    r9281 r9302  
    7171using namespace std;
    7272
     73// --------------------------------------------------------------------------
     74//
     75// Default constructor.
     76//
     77MAlphaFitter::MAlphaFitter(const char *name, const char *title) : fSigInt(15),
     78    fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80),
     79    fPolynomOrder(2), fFitBackground(kTRUE), fFunc(0),
     80    fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
     81{
     82    fName  = name  ? name  : "MAlphaFitter";
     83    fTitle = title ? title : "Fit alpha";
     84
     85    SetSignalFunction(kGauss);
     86
     87    Clear();
     88}
     89
     90// --------------------------------------------------------------------------
     91//
     92// Destructor
     93//
     94MAlphaFitter::~MAlphaFitter()
     95{
     96    delete fFunc;
     97}
     98
     99// --------------------------------------------------------------------------
     100//
     101// Re-initializes fFunc either according to SignalFunc_t
     102//
     103void MAlphaFitter::SetSignalFunction(SignalFunc_t func)
     104{
     105    if (gROOT->GetListOfFunctions()->FindObject(""))
     106    {
     107        gLog << err << "MAlphaFitter::SetSignalFunction -- '' found!" << endl;
     108        return;
     109    }
     110
     111    delete fFunc;
     112    fFunc = 0;
     113
     114    switch (func)
     115    {
     116    case kGauss:
     117        fFunc=new TF1("", MString::Format("gaus(0) + pol%d(3)", fPolynomOrder));
     118        break;
     119    case kThetaSq:
     120        if (fPolynomOrder>0)
     121            fPolynomOrder = 1;
     122        fFunc=new TF1("", "[0]*exp(-0.5*((sqrt(x)-[1])/[2])^2) + expo(3)");
     123        break;
     124    }
     125
     126    fSignalFunc=func;
     127
     128    fFunc->SetName("Dummy");
     129    gROOT->GetListOfFunctions()->Remove(fFunc);
     130
     131    fCoefficients.Set(3+fPolynomOrder+1);
     132    fCoefficients.Reset();
     133
     134    fErrors.Set(3+fPolynomOrder+1);
     135    fErrors.Reset();
     136}
     137
     138// --------------------------------------------------------------------------
     139//
     140// Reset variables which belong to results. Reset the arrays.
     141//
    73142void MAlphaFitter::Clear(Option_t *o)
    74143{
    … …  
    86155    fCoefficients.Reset();
    87156    fErrors.Reset();
     157}
     158
     159// --------------------------------------------------------------------------
     160//
     161// Returns fFunc->Eval(d) or 0 if fFunc==NULL
     162//
     163Double_t MAlphaFitter::Eval(Double_t d) const
     164{
     165    return fFunc ? fFunc->Eval(d) : 0;
    88166}
    89167
    … …  
    294372    h.SetDirectory(0);
    295373    h.Add(&hon);
     374
    296375    h.Scale(0.5);
    297376    for (int i=1; i<=bin+3; i++)
    … …  
    310389    // Do a gaussian error propagation to calculated the error of
    311390    // the background estimated from the fit
    312     const Double_t ea = fit.GetErrors()[3];
    313     const Double_t eb = fit.GetErrors()[4];
    314     const Double_t a  = fit.GetCoefficients()[3];
    315     const Double_t b  = fit.GetCoefficients()[4];
     391    const Double_t ea = fit.fErrors[3];
     392    const Double_t eb = fit.fErrors[4];
     393    const Double_t a  = fit.fCoefficients[3];
     394    const Double_t b  = fit.fCoefficients[4];
    316395
    317396    const Double_t t  = fIntegralMax;
    … …  
    333412    // const Double_t sc = bg * er*er / (fit2.GetEventsBackground()*fit2.GetEventsBackground());
    334413    // Assuming that bg and fit2.GetEventsBackground() are rather identical:
    335     const Double_t sc = er*er / fit.GetEventsBackground();
     414    const Double_t sc = er*er / fit.fEventsBackground;
     415
     416
    336417    /*
    337418     cout << MMath::SignificanceLiMaSigned(hon.Integral(1, bin), fit.GetEventsBackground()/sc, sc) << " ";
    … …  
    354435        return kFALSE;
    355436
    356     fChiSqSignal  = fit.GetChiSqSignal();
    357     fChiSqBg      = fit.GetChiSqBg();
    358     fCoefficients = fit.GetCoefficients();
    359     fErrors       = fit.GetErrors();
     437    fChiSqSignal  = fit.fChiSqSignal;
     438    fChiSqBg      = fit.fChiSqBg;
     439    fCoefficients = fit.fCoefficients;
     440    fErrors       = fit.fErrors;
    360441
    361442    // ----------------------------------------------------------------------------
    … …  
    472553
    473554    // Function
    474     TF1 *fcn = f.fFunc;
     555    delete f.fFunc;
     556
    475557    f.fFunc = new TF1(*fFunc);
    476558    f.fFunc->SetName("Dummy");
    477559    gROOT->GetListOfFunctions()->Remove(f.fFunc);
    478     delete fcn;
    479560}
    480561
    … …  
    527608{
    528609    const TString name(MString::Format("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
     610
    529611    TH1D *h = hon.ProjectionZ(name, 0, hon.GetNbinsX()+1, bin, bin, "E");
     612    h->SetDirectory(0);
     613
     614    const Bool_t rc = Fit(*h, paint);
     615
     616    delete h;
     617
     618    return rc;
     619}
     620
     621Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
     622{
     623    const TString name(MString::Format("TempAlphaTheta%06d", gRandom->Integer(1000000)));
     624
     625    TH1D *h = hon.ProjectionZ(name, bin, bin, 0, hon.GetNbinsY()+1, "E");
     626    h->SetDirectory(0);
     627
     628    const Bool_t rc = Fit(*h, paint);
     629
     630    delete h;
     631
     632    return rc;
     633}
     634/*
     635Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint)
     636{
     637    const TString name(Form("TempAlphaTime%06d", gRandom->Integer(1000000)));
     638
     639    hon.GetZaxis()->SetRange(bin,bin);
     640    TH1D *h = (TH1D*)hon.Project3D("ye");
     641    hon.GetZaxis()->SetRange(-1,-1);
     642
    530643    h->SetDirectory(0);
    531644
    … …  
    534647    return rc;
    535648}
    536 
    537 Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
    538 {
    539     const TString name(MString::Format("TempAlphaTheta%06d", gRandom->Integer(1000000)));
    540     TH1D *h = hon.ProjectionZ(name, bin, bin, 0, hon.GetNbinsY()+1, "E");
    541     h->SetDirectory(0);
    542 
    543     const Bool_t rc = Fit(*h, paint);
    544     delete h;
    545     return rc;
    546 }
    547 /*
    548 Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint)
    549 {
    550     const TString name(Form("TempAlphaTime%06d", gRandom->Integer(1000000)));
    551 
    552     hon.GetZaxis()->SetRange(bin,bin);
    553     TH1D *h = (TH1D*)hon.Project3D("ye");
    554     hon.GetZaxis()->SetRange(-1,-1);
    555 
    556     h->SetDirectory(0);
    557 
    558     const Bool_t rc = Fit(*h, paint);
    559     delete h;
    560     return rc;
    561 }
    562649*/
    563650Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
    564651{
    565652    const TString name(MString::Format("TempAlpha%06d", gRandom->Integer(1000000)));
     653
    566654    TH1D *h = hon.ProjectionZ(name, 0, hon.GetNbinsX()+1, 0, hon.GetNbinsY()+1, "E");
    567655    h->SetDirectory(0);
    568656
    569657    const Bool_t rc = Fit(*h, paint);
     658
    570659    delete h;
     660
    571661    return rc;
    572662}
    … …  
    578668
    579669    TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, bin, bin, "E");
     670    h1->SetDirectory(0);
     671
    580672    TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, bin, bin, "E");
    581     h1->SetDirectory(0);
    582673    h0->SetDirectory(0);
    583674
    … …  
    596687
    597688    TH1D *h1 = hon.ProjectionZ(name1, bin, bin, 0, hon.GetNbinsY()+1, "E");
     689    h1->SetDirectory(0);
     690
    598691    TH1D *h0 = hof.ProjectionZ(name0, bin, bin, 0, hof.GetNbinsY()+1, "E");
    599     h1->SetDirectory(0);
    600692    h0->SetDirectory(0);
    601693
    … …  
    638730
    639731    TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, 0, hon.GetNbinsY()+1, "E");
     732    h1->SetDirectory(0);
     733
    640734    TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, 0, hof.GetNbinsY()+1, "E");
    641     h1->SetDirectory(0);
    642735    h0->SetDirectory(0);
    643736
    … …  
    656749
    657750    TH1D *h1 = hon.ProjectionZ(name1, 0, hon.GetNbinsX()+1, bin, bin, "E");
     751    h1->SetDirectory(0);
     752
    658753    TH1D *h0 = hof.ProjectionZ(name0, 0, hof.GetNbinsX()+1, bin, bin, "E");
    659     h1->SetDirectory(0);
    660754    h0->SetDirectory(0);
    661755

  • trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

    r8989 r9302  
    1010#endif
    1111
    12 #ifndef ROOT_TF1
    13 #include <TF1.h>
    14 #endif
    15 
     12class TF1;
    1613class TH1D;
    1714class TH3D;
    … …  
    8582public:
    8683    // 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 
     84    MAlphaFitter(const char *name=0, const char *title=0);
    10385    MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
    10486    {
    10587        f.Copy(*this);
    10688    }
    107     ~MAlphaFitter()
    108     {
    109         delete fFunc;
    110     }
     89    ~MAlphaFitter();
    11190
    11291    // TObject
    … …  
    134113        SetSignalFunction(fSignalFunc);
    135114    }
    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     }
     115    void SetSignalFunction(SignalFunc_t func);
     116
    165117    void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; }
    166118
    … …  
    193145    const TArrayD &GetCoefficients() const { return fCoefficients; }
    194146    const TArrayD &GetErrors() const       { return fErrors; }
    195     Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; }
     147    Double_t Eval(Double_t d) const;
    196148
    197149    Double_t CalcUpperLimit() const;
Note: See TracChangeset for help on using the changeset viewer.