Changeset 5012 for trunk


Ignore:
Timestamp:
09/14/04 19:04:01 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r5006 r5012  
    8686    Double_t sigint=fSigInt;
    8787    Double_t sigmax=fSigMax;
    88     Double_t bgmin=fBgMin;
    89     Double_t bgmax=fBgMax;
    90     Byte_t polynom=fPolynomOrder;
     88    Double_t bgmin =fBgMin;
     89    Double_t bgmax =fBgMax;
    9190
    92     // Implementing the function yourself is only about 5% faster
    93     TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    94     //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
     91    //TF1 fFunc("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90);
    9592
    96     if (fCoefficients.GetSize() != polynom+3)
    97     {
    98         fCoefficients.Set(polynom+3);
    99         fCoefficients.Reset();
    100     }
    101     else
    102         func.SetParameters(fCoefficients.GetArray());
     93    //fFunc->SetParameters(fCoefficients.GetArray());
    10394
    104     func.FixParameter(1, 0);
    105     func.FixParameter(4, 0);
    106     func.SetParLimits(2, 0, 90);
    107     func.SetParLimits(3, -1, 1);
     95    fFunc->FixParameter(1, 0);
     96    fFunc->FixParameter(4, 0);
     97    fFunc->SetParLimits(2, 0, 90);
     98    fFunc->SetParLimits(3, -1, 1);
    10899
    109100    const Double_t alpha0 = h.GetBinContent(1);
    … …  
    115106
    116107    // First fit a polynom in the off region
    117     func.FixParameter(0, 0);
    118     func.FixParameter(2, 1);
    119     func.ReleaseParameter(3);
     108    fFunc->FixParameter(0, 0);
     109    fFunc->FixParameter(2, 1);
     110    fFunc->ReleaseParameter(3);
    120111
    121     for (int i=5; i<func.GetNpar(); i++)
    122         func.ReleaseParameter(i);
     112    for (int i=5; i<fFunc->GetNpar(); i++)
     113        fFunc->ReleaseParameter(i);
    123114
    124115    // options : N  do not store the function, do not draw
    … …  
    126117    //           R  use the range specified in the function range
    127118    //           Q  quiet mode
    128     h.Fit(&func, "NQI", "", bgmin, bgmax);
     119    h.Fit(fFunc, "NQI", "", bgmin, bgmax);
    129120
    130     fChiSqBg = func.GetChisquare()/func.GetNDF();
     121    fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
    131122
    132123    // ------------------------------------
    133124    if (paint)
    134125    {
    135         func.SetRange(0, 90);
    136         func.SetLineColor(kRed);
    137         func.SetLineWidth(2);
    138         func.Paint("same");
     126        fFunc->SetRange(0, 90);
     127        fFunc->SetLineColor(kRed);
     128        fFunc->SetLineWidth(2);
     129        fFunc->Paint("same");
    139130    }
    140131    // ------------------------------------
    141132
    142     func.ReleaseParameter(0);
     133    fFunc->ReleaseParameter(0);
    143134    //func.ReleaseParameter(1);
    144     func.ReleaseParameter(2);
    145     func.FixParameter(3, func.GetParameter(3));
    146     for (int i=5; i<func.GetNpar(); i++)
    147         func.FixParameter(i, func.GetParameter(i));
     135    fFunc->ReleaseParameter(2);
     136    fFunc->FixParameter(3, fFunc->GetParameter(3));
     137    for (int i=5; i<fFunc->GetNpar(); i++)
     138        fFunc->FixParameter(i, fFunc->GetParameter(i));
    148139
    149140    // Do not allow signals smaller than the background
    150     const Double_t A  = alpha0-func.GetParameter(3);
     141    const Double_t A  = alpha0-fFunc->GetParameter(3);
    151142    const Double_t dA = TMath::Abs(A);
    152     func.SetParLimits(0, -dA*4, dA*4);
    153     func.SetParLimits(2, 0, 90);
     143    fFunc->SetParLimits(0, -dA*4, dA*4);
     144    fFunc->SetParLimits(2, 0, 90);
    154145
    155146    // Now fit a gaus in the on region on top of the polynom
    156     func.SetParameter(0, A);
    157     func.SetParameter(2, sigmax*0.75);
     147    fFunc->SetParameter(0, A);
     148    fFunc->SetParameter(2, sigmax*0.75);
    158149
    159150    // options : N  do not store the function, do not draw
    … …  
    161152    //           R  use the range specified in the function range
    162153    //           Q  quiet mode
    163     h.Fit(&func, "NQI", "", 0, sigmax);
     154    h.Fit(fFunc, "NQI", "", 0, sigmax);
    164155
    165     fChiSqSignal = func.GetChisquare()/func.GetNDF();
    166     fCoefficients.Set(func.GetNpar(), func.GetParameters());
     156    fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
     157    fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());
    167158
    168159    //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    … …  
    171162    if (paint)
    172163    {
    173         func.SetLineColor(kGreen);
    174         func.SetLineWidth(2);
    175         func.Paint("same");
     164        fFunc->SetLineColor(kGreen);
     165        fFunc->SetLineWidth(2);
     166        fFunc->Paint("same");
    176167    }
    177168    // ------------------------------------
    178     //const Double_t s = func.Integral(0, sigint)/alphaw;
    179     func.SetParameter(0, 0);
    180     func.SetParameter(2, 1);
    181     //const Double_t b = func.Integral(0, sigint)/alphaw;
    182 
    183     //fSignificance = MMath::SignificanceLiMaSigned(s, b);
    184     //exc = s-b;
     169    const Double_t s = fFunc->Integral(0, sigint)/alphaw;
     170    fFunc->SetParameter(0, 0);
     171    fFunc->SetParameter(2, 1);
     172    const Double_t b = fFunc->Integral(0, sigint)/alphaw;
     173    fSignificance = MMath::SignificanceLiMaSigned(s, b);
    185174
    186175    const Double_t uin = 1.25*sigint;
    … …  
    188177
    189178    fIntegralMax      = h.GetBinLowEdge(bin+1);
    190     fEventsBackground = func.Integral(0, fIntegralMax)/alphaw;
     179    fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
    191180    fEventsSignal     = h.Integral(0, bin);
    192181    fEventsExcess     = fEventsSignal-fEventsBackground;
    193     fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
     182    //fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
    194183
    195184    return kTRUE;

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

    r5006 r5012  
    1010#endif
    1111
     12#ifndef ROOT_TF1
     13#include <TF1.h>
     14#endif
     15
    1216class TH1D;
    1317
    … …  
    1519{
    1620private:
     21    TF1 *fFunc;
     22
    1723    Float_t fSigInt;
    1824    Float_t fSigMax;
    … …  
    3339
    3440public:
    35     MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1)
     41    // Implementing the function yourself is only about 5% faster
     42    MAlphaFitter() : fFunc(new TF1("", "gaus(0) + pol1(3)", 0, 90)), fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1), fCoefficients(3+fPolynomOrder+1)
    3643    {
    3744    }
    3845
    39     MAlphaFitter(const MAlphaFitter &f)
     46    MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
    4047    {
    4148        f.Copy(*this);
     49    }
     50    ~MAlphaFitter()
     51    {
     52        delete fFunc;
    4253    }
    4354
    … …  
    5162        f.fBgMax        = fBgMax;
    5263        f.fPolynomOrder = fPolynomOrder;
     64        f.fCoefficients.Set(fCoefficients.GetSize());
     65        f.fCoefficients.Reset();
     66
     67        TF1 *fcn = f.fFunc;
     68        f.fFunc = new TF1(*fFunc);
     69        delete fcn;
    5370    }
    5471
    … …  
    5774    void SetBackgroundFitMin(Float_t s)    { fBgMin        = s; }
    5875    void SetBackgroundFitMax(Float_t s)    { fBgMax        = s; }
    59     void SetPolynomOrder(Int_t s)          { fPolynomOrder = s; }
     76    void SetPolynomOrder(Int_t s)          { fPolynomOrder = s; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", s)); fCoefficients.Set(3+s+1); fCoefficients.Reset(); }
    6077
    6178    Double_t GetEventsExcess() const       { return fEventsExcess; }
    … …  
    7693
    7794    Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
     95
    7896    ClassDef(MAlphaFitter, 1)
    7997};

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

    r5006 r5012  
    140140void MHAlpha::FitEnergySpec(Bool_t paint)
    141141{
    142     TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]");
    143     f.SetParameter(0, -2.0);
    144     f.SetParameter(1,  500); // 50
    145     f.SetParameter(2,    1); // 50
     142    /*
     143    TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
     144    f.SetParameter(0,  -2.1);
     145    f.SetParameter(1,  1400); // 50
     146    f.SetParameter(2,  fHEnergy.Integral()); // 50
    146147    f.SetLineWidth(2);
    147     f.SetLineColor(kGreen);
    148 
    149     fHEnergy.Fit(&f, "0QI", "", 90, 1900); // Integral?
     148
     149    fHEnergy.Fit(&f,  "0QI", "", 100, 2400); // Integral? 90, 2800
     150
     151    const Double_t chi2 = f.GetChisquare();
     152    const Int_t    NDF  = f.GetNDF();
     153
     154    // was fit successful ?
     155    const Bool_t ok = NDF>0 && chi2<2.5*NDF;
     156    f.SetLineColor(ok ? kGreen : kRed);
    150157
    151158    if (paint)
    … …  
    153160        f.Paint("same");
    154161
    155         TLatex text(0.2, 0.94, Form("\\alpha=%.1f E_{0}=%.1fGeV N=%.1f",
    156                                     f.GetParameter(0),
    157                                     f.GetParameter(1),
    158                                     0/*f.GetParameter(2)*/));
     162        if (ok)
     163        {
     164        TString s;
     165        s += Form("\\alpha=%.1f ",   f.GetParameter(0));
     166        s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
     167        s += Form("N=%.1e",          f.GetParameter(2));
     168
     169        TLatex text(0.25, 0.94, s.Data());
    159170        text.SetBit(TLatex::kTextNDC);
    160171        text.SetTextSize(0.04);
    161172        text.Paint();
    162     }
     173        }
     174    }*/
    163175}
    164176
    … …  
    449461    if (o==(TString)"energy")
    450462    {
    451         /*
    452463        if (fHEnergy.GetEntries()>10)
    453464        {
    … …  
    455466            gPad->SetLogy();
    456467        }
    457         FitEnergySpec(kTRUE);*/
     468        FitEnergySpec(kTRUE);
    458469
    459470    }
    … …  
    525536}
    526537
    527 // --------------------------------------------------------------------------
    528 //
    529 // This is a preliminary implementation of a alpha-fit procedure for
    530 // all possible source positions. It will be moved into its own
    531 // more powerfull class soon.
    532 //
    533 // The fit function is "gaus(0)+pol2(3)" which is equivalent to:
    534 //   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
    535 // or
    536 //   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
    537 //
    538 // Parameter [1] is fixed to 0 while the alpha peak should be
    539 // symmetric around alpha=0.
    540 //
    541 // Parameter [4] is fixed to 0 because the first derivative at
    542 // alpha=0 should be 0, too.
    543 //
    544 // In a first step the background is fitted between bgmin and bgmax,
    545 // while the parameters [0]=0 and [2]=1 are fixed.
    546 //
    547 // In a second step the signal region (alpha<sigmax) is fittet using
    548 // the whole function with parameters [1], [3], [4] and [5] fixed.
    549 //
    550 // The number of excess and background events are calculated as
    551 //   s = int(hist,    0, 1.25*sigint)
    552 //   b = int(pol2(3), 0, 1.25*sigint)
    553 //
    554 // The Significance is calculated using the Significance() member
    555 // function.
    556 //
    557 /*
    558 Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
    559 {
    560     Double_t sigint=fSigInt;
    561     Double_t sigmax=fSigMax;
    562     Double_t bgmin=fBgMin;
    563     Double_t bgmax=fBgMax;
    564     Byte_t polynom=fPolynom;
    565 
    566     // Implementing the function yourself is only about 5% faster
    567     TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    568     //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
    569     TArrayD maxpar(func.GetNpar());
    570 
    571     func.FixParameter(1, 0);
    572     func.FixParameter(4, 0);
    573     func.SetParLimits(2, 0, 90);
    574     func.SetParLimits(3, -1, 1);
    575 
    576     const Double_t alpha0 = h.GetBinContent(1);
    577     const Double_t alphaw = h.GetXaxis()->GetBinWidth(1);
    578 
    579     // Check for the regios which is not filled...
    580     if (alpha0==0)
    581         return kFALSE; //*fLog << warn << "Histogram empty." << endl;
    582 
    583     // First fit a polynom in the off region
    584     func.FixParameter(0, 0);
    585     func.FixParameter(2, 1);
    586     func.ReleaseParameter(3);
    587 
    588     for (int i=5; i<func.GetNpar(); i++)
    589         func.ReleaseParameter(i);
    590 
    591     // options : N  do not store the function, do not draw
    592     //           I  use integral of function in bin rather than value at bin center
    593     //           R  use the range specified in the function range
    594     //           Q  quiet mode
    595     h.Fit(&func, "NQI", "", bgmin, bgmax);
    596 
    597     fChiSqBg = func.GetChisquare()/func.GetNDF();
    598 
    599     // ------------------------------------
    600     if (paint)
    601     {
    602         func.SetRange(0, 90);
    603         func.SetLineColor(kRed);
    604         func.SetLineWidth(2);
    605         func.Paint("same");
    606     }
    607     // ------------------------------------
    608 
    609     func.ReleaseParameter(0);
    610     //func.ReleaseParameter(1);
    611     func.ReleaseParameter(2);
    612     func.FixParameter(3, func.GetParameter(3));
    613     for (int i=5; i<func.GetNpar(); i++)
    614         func.FixParameter(i, func.GetParameter(i));
    615 
    616     // Do not allow signals smaller than the background
    617     const Double_t A  = alpha0-func.GetParameter(3);
    618     const Double_t dA = TMath::Abs(A);
    619     func.SetParLimits(0, -dA*4, dA*4);
    620     func.SetParLimits(2, 0, 90);
    621 
    622     // Now fit a gaus in the on region on top of the polynom
    623     func.SetParameter(0, A);
    624     func.SetParameter(2, sigmax*0.75);
    625 
    626     // options : N  do not store the function, do not draw
    627     //           I  use integral of function in bin rather than value at bin center
    628     //           R  use the range specified in the function range
    629     //           Q  quiet mode
    630     h.Fit(&func, "NQI", "", 0, sigmax);
    631 
    632     fChiSqSignal  = func.GetChisquare()/func.GetNDF();
    633     fSigmaGaus    = func.GetParameter(2);
    634 
    635     //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    636 
    637     // ------------------------------------
    638     if (paint)
    639     {
    640         func.SetLineColor(kGreen);
    641         func.SetLineWidth(2);
    642         func.Paint("same");
    643     }
    644     // ------------------------------------
    645     const Double_t s   = func.Integral(0, sigint)/alphaw;
    646     func.SetParameter(0, 0);
    647     func.SetParameter(2, 1);
    648     const Double_t b   = func.Integral(0, sigint)/alphaw;
    649 
    650     fSignificance = MMath::SignificanceLiMaSigned(s, b);
    651     //exc = s-b;
    652 
    653     const Double_t uin = 1.25*sigint;
    654     const Int_t    bin = h.GetXaxis()->FindFixBin(uin);
    655     fIntegralMax  = h.GetBinLowEdge(bin+1);
    656     fExcessEvents = h.Integral(0, bin)-func.Integral(0, fIntegralMax)/alphaw;
    657 
    658     return kTRUE;
    659 }
    660 
    661 void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
    662 {
    663     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)",
    664                            fSignificance, fSigInt, fSigmaGaus,
    665                            (int)fExcessEvents, fIntegralMax,
    666                            fChiSqBg, fChiSqSignal));
    667 
    668     text.SetBit(TLatex::kTextNDC);
    669     text.SetTextSize(size);
    670     text.Paint();
    671 }
    672 */
    673538Bool_t MHAlpha::Finalize()
    674539{

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

    r5002 r5012  
    8484{
    8585private:
    86     MAlphaFitter fFit;
     86    MAlphaFitter fFit;          //! SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT)
    8787
    8888    TH3D fHAlpha;               // Alpha vs. theta and energy

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

    r4999 r5012  
    188188    fHThetaLambda.SetTitle("Slope (Rate) vs Theta");
    189189    fHThetaLambda.SetXTitle("\\Theta [\\circ]");
    190     fHThetaLambda.SetYTitle("\\labda");
     190    fHThetaLambda.SetYTitle("\\lambda [s^{-1}]");
    191191    fHThetaLambda.UseCurrentStyle();
    192192    fHThetaLambda.SetDirectory(NULL);
    … …  
    197197    fHTimeLambda.SetTitle("Slope (Rate) vs Time");
    198198    fHTimeLambda.SetXTitle("\\Time [\\circ]");
    199     fHTimeLambda.SetYTitle("\\lambda");
     199    fHTimeLambda.SetYTitle("\\lambda [s^{-1}]");
    200200    fHTimeLambda.UseCurrentStyle();
    201201    fHTimeLambda.SetDirectory(NULL);
    … …  
    318318    // parameter 1 = N0*del
    319319    //
    320     static TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)");
     320    TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)");
    321321    //func.SetParNames("lambda", "N0", "del");
    322322
    … …  
    791791    AppendPad("time");
    792792    fHTimeLambda.Draw("same");
    793     DrawRightAxis("\\lambda");
     793    DrawRightAxis("\\lambda [s^{-1}]");
    794794
    795795    pad->cd(2);
    … …  
    821821    AppendPad("theta");
    822822    fHThetaLambda.Draw("same");
    823     DrawRightAxis("\\lambda");
    824 }
     823    DrawRightAxis("\\lambda [s^{-1}]");
     824}

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

    r5006 r5012  
    808808    const Double_t w = fHist.GetZaxis()->GetBinWidth(1);
    809809
    810     //                      xmin, xmax, npar
    811     //TF1 func("MyFunc", fcn, 0, 90, 6);
    812     // Implementing the function yourself is only about 5% faster
    813     TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    814     //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
    815     TArrayD maxpar(func.GetNpar());
     810    TArrayD maxpar;
    816811
    817812    /*  func.SetParName(0, "A");
    … …  
    819814     *  func.SetParName(2, "sigma");
    820815    */
    821 
    822     func.FixParameter(1, 0);
    823     func.FixParameter(4, 0);
    824     func.SetParLimits(2, 0, 90);
    825     func.SetParLimits(3, -1, 1);
    826816
    827817    const Int_t nx = hist->GetXaxis()->GetNbins();
    … …  
    864854            h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
    865855
    866             if (h->GetBinContent(1)==0)
     856            if (h->GetEntries()==0)
    867857                continue;
    868858
    … …  
    988978        result->Draw();
    989979
    990         TF1 f1("", func.GetExpFormula(), 0, 90);
    991         TF1 f2("", func.GetExpFormula(), 0, 90);
     980        TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
     981        TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    992982        f1.SetParameters(maxpar.GetArray());
    993983        f2.SetParameters(maxpar.GetArray());

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

    r4999 r5012  
    6161    TH1 *GetHistByName(const TString name) { return &fHist; }
    6262
    63     void FitSignificance(Float_t sigint=15, Float_t sigmax=70, Float_t bgmin=40, Float_t bgmax=70, Byte_t polynom=1); //*MENU*
     63    void FitSignificance(Float_t sigint=10, Float_t sigmax=75, Float_t bgmin=45, Float_t bgmax=85, Byte_t polynom=2); //*MENU*
    6464    void FitSignificanceStd() { FitSignificance(); } //*MENU*
    6565
Note: See TracChangeset for help on using the changeset viewer.