Ignore:
Timestamp:
09/14/04 16:54:56 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r5002 r5006  
    4646#include "MMath.h"
    4747
     48#include "MLogManip.h"
     49
    4850ClassImp(MAlphaFitter);
    4951
    … …  
    8688    Double_t bgmin=fBgMin;
    8789    Double_t bgmax=fBgMax;
    88     Byte_t polynom=fPolynom;
     90    Byte_t polynom=fPolynomOrder;
    8991
    9092    // Implementing the function yourself is only about 5% faster
    9193    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    9294    //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
    93     TArrayD maxpar(func.GetNpar());
     95
     96    if (fCoefficients.GetSize() != polynom+3)
     97    {
     98        fCoefficients.Set(polynom+3);
     99        fCoefficients.Reset();
     100    }
     101    else
     102        func.SetParameters(fCoefficients.GetArray());
    94103
    95104    func.FixParameter(1, 0);
    … …  
    154163    h.Fit(&func, "NQI", "", 0, sigmax);
    155164
    156     fChiSqSignal  = func.GetChisquare()/func.GetNDF();
    157     fSigmaGaus    = func.GetParameter(2);
     165    fChiSqSignal = func.GetChisquare()/func.GetNDF();
     166    fCoefficients.Set(func.GetNpar(), func.GetParameters());
    158167
    159168    //const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    … …  
    167176    }
    168177    // ------------------------------------
    169     const Double_t s  = func.Integral(0, sigint)/alphaw;
     178    //const Double_t s = func.Integral(0, sigint)/alphaw;
    170179    func.SetParameter(0, 0);
    171180    func.SetParameter(2, 1);
    172     const Double_t b  = func.Integral(0, sigint)/alphaw;
    173 
    174     fSignificance = MMath::SignificanceLiMaSigned(s, b);
     181    //const Double_t b = func.Integral(0, sigint)/alphaw;
     182
     183    //fSignificance = MMath::SignificanceLiMaSigned(s, b);
    175184    //exc = s-b;
    176185
    177186    const Double_t uin = 1.25*sigint;
    178187    const Int_t    bin = h.GetXaxis()->FindFixBin(uin);
    179     fIntegralMax  = h.GetBinLowEdge(bin+1);
    180     fExcessEvents = h.Integral(0, bin)-func.Integral(0, fIntegralMax)/alphaw;
     188
     189    fIntegralMax      = h.GetBinLowEdge(bin+1);
     190    fEventsBackground = func.Integral(0, fIntegralMax)/alphaw;
     191    fEventsSignal     = h.Integral(0, bin);
     192    fEventsExcess     = fEventsSignal-fEventsBackground;
     193    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
    181194
    182195    return kTRUE;
    … …  
    186199{
    187200    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)",
    188                            fSignificance, fSigInt, fSigmaGaus,
    189                            (int)fExcessEvents, fIntegralMax,
     201                           fSignificance, fSigInt, GetGausSigma(),
     202                           (int)fEventsExcess, fIntegralMax,
    190203                           fChiSqBg, fChiSqSignal));
    191204

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

    r5002 r5006  
    44#ifndef MARS_MParContainer
    55#include "MParContainer.h"
     6#endif
     7
     8#ifndef ROOT_TArrayD
     9#include <TArrayD.h>
    610#endif
    711
    … …  
    1519    Float_t fBgMin;
    1620    Float_t fBgMax;
    17     Int_t   fPolynom;
     21    Int_t   fPolynomOrder;
    1822
    1923    Double_t fSignificance;
    20     Double_t fExcessEvents;
     24    Double_t fEventsExcess;
     25    Double_t fEventsSignal;
     26    Double_t fEventsBackground;
     27
    2128    Double_t fChiSqSignal;
    2229    Double_t fChiSqBg;
    23     Double_t fSigmaGaus;
    2430    Double_t fIntegralMax;
    2531
     32    TArrayD fCoefficients;
     33
    2634public:
    27     MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynom(1)
     35    MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynomOrder(1)
    2836    {
    2937    }
    … …  
    3846        MAlphaFitter &f = static_cast<MAlphaFitter&>(o);
    3947
    40         f.fSigInt  = fSigInt;
    41         f.fSigMax  = fSigMax;
    42         f.fBgMin   = fBgMin;
    43         f.fBgMax   = fBgMax;
    44         f.fPolynom = fPolynom;
     48        f.fSigInt       = fSigInt;
     49        f.fSigMax       = fSigMax;
     50        f.fBgMin        = fBgMin;
     51        f.fBgMax        = fBgMax;
     52        f.fPolynomOrder = fPolynomOrder;
    4553    }
    4654
    47     void SetSignalIntegralMax(Float_t s) { fSigInt = s; }
    48     void SetSignalFitMax(Float_t s)      { fSigMax = s; }
    49     void SetBackgroundFitMin(Float_t s)  { fBgMin = s; }
    50     void SetBackgroundFitMax(Float_t s)  { fBgMax = s; }
    51     void SetNumPolynom(Int_t s)          { fPolynom = s; }
     55    void SetSignalIntegralMax(Float_t s)   { fSigInt      = s; }
     56    void SetSignalFitMax(Float_t s)        { fSigMax      = s; }
     57    void SetBackgroundFitMin(Float_t s)    { fBgMin        = s; }
     58    void SetBackgroundFitMax(Float_t s)    { fBgMax        = s; }
     59    void SetPolynomOrder(Int_t s)          { fPolynomOrder = s; }
    5260
    53     Double_t GetExcessEvents() const { return fExcessEvents; }
    54     Double_t GetSignificance() const { return fSignificance; }
    55     Double_t GetChiSqSignal() const  { return fChiSqSignal; }
    56     Double_t GetChiSqBg() const      { return fChiSqBg; }
    57     Double_t GetSigmaGaus() const    { return fSigmaGaus; }
     61    Double_t GetEventsExcess() const       { return fEventsExcess; }
     62    Double_t GetEventsSignal() const       { return fEventsSignal; }
     63    Double_t GetEventsBackground() const   { return fEventsBackground; }
     64
     65    Double_t GetSignificance() const       { return fSignificance; }
     66    Double_t GetChiSqSignal() const        { return fChiSqSignal; }
     67    Double_t GetChiSqBg() const            { return fChiSqBg; }
     68
     69    Double_t GetGausSigma() const          { return fCoefficients[2]; }
     70    Double_t GetGausMu() const             { return fCoefficients[1]; }
     71    Double_t GetGausA() const              { return fCoefficients[0]; }
     72    Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
     73    const TArrayD &GetCoefficients() const { return fCoefficients; }
    5874
    5975    void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const;

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

    r5002 r5006  
    175175        if (fit.Fit(*h))
    176176        {
    177             fHEnergy.SetBinContent(i, fit.GetExcessEvents());
    178             fHEnergy.SetBinError(i, fit.GetExcessEvents()*0.2);
     177            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
     178            fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
    179179        }
    180180        delete h;
    … …  
    194194        if (fit.Fit(*h))
    195195        {
    196             fHTheta.SetBinContent(i, fit.GetExcessEvents());
    197             fHTheta.SetBinError(i, fit.GetExcessEvents()*0.2);
     196            fHTheta.SetBinContent(i, fit.GetEventsExcess());
     197            fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
    198198        }
    199199        delete h;
    … …  
    334334    // Fill histogram
    335335    //
    336     fHTime.SetBinContent(n+1, fit.GetExcessEvents());
    337     fHTime.SetBinError(n+1, fit.GetExcessEvents()*0.1);
    338 
    339     *fLog << all << *fTimeEffOn << ": " << fit.GetExcessEvents() << endl;
     336    fHTime.SetBinContent(n+1, fit.GetEventsExcess());
     337    fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);
     338
     339    *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;
    340340
    341341    rebin = steps;
    … …  
    353353    {
    354354        alpha  = (*fMatrix)[fMap[0]];
    355         energy = 1000; //(*fMatrix)[fMap[1]];
    356         theta  =    0; //(*fMatrix)[fMap[2]];
     355        energy = 1000;
     356        theta  =    0;
    357357    }
    358358    else
    … …  
    375375    fHAlphaTime.Fill(TMath::Abs(alpha), w);
    376376
    377     /*
    378      // N_gamma vs Energy and Theta
    379      const Double_t Ngam  = fHUnfold.GetBinContent(m,n);
    380      // A_eff   vs Energy and Theta
    381      const Double_t Aeff  = fHAeff.GetBinContent(m,n);
    382      // T_eff   vs Theta
    383      const Double_t Effon = teff.GetBinContent(n);
    384 
    385      const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam;
    386      const Double_t c2 = teff.GetBinError(n)      /Effon;
    387      const Double_t c3 = fHAeff.GetBinError(m,n)  /Aeff;
    388 
    389      const Double_t cont  = Ngam / (DeltaE * Effon * Aeff);
    390      const Double_t dcont = sqrt(c1*c1 + c2*c2 + c3*c3);
    391 
    392      //
    393      // Out of Range
    394      //
    395      const Bool_t oor = Ngam<=0 || DeltaE<=0 || Effon<=0 || Aeff<=0;
    396 
    397      if (oor)
    398      *fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") ";
    399 
    400      if (Ngam<=0)
    401      *fLog << " Ngam=0";
    402      if (DeltaE<=0)
    403      *fLog << " DeltaE=0";
    404      if (Effon<=0)
    405      *fLog << " Effon=0";
    406      if (Aeff<=0)
    407      *fLog << " Aeff=0";
    408 
    409      if (oor)
    410      *fLog << endl;
    411 
    412      fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont);
    413      fHFlux.SetBinError(m,n,   oor ? 1e-20 : dcont*cont);
    414      */
    415377    return kTRUE;
    416378}
    … …  
    487449    if (o==(TString)"energy")
    488450    {
    489         if (fHEnergy.GetEntries()>0)
     451        /*
     452        if (fHEnergy.GetEntries()>10)
    490453        {
    491454            gPad->SetLogx();
    492455            gPad->SetLogy();
    493456        }
    494         FitEnergySpec(kTRUE);
     457        FitEnergySpec(kTRUE);*/
     458
    495459    }
    496460

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

    r4999 r5006  
    131131#include "MAstroSky2Local.h"
    132132#include "MStatusDisplay.h"
     133
    133134#include "MMath.h"
     135#include "MAlphaFitter.h"
    134136
    135137#include "MBinning.h"
    … …  
    847849    delete h0;
    848850
    849     TH1 *h=0;
     851    MAlphaFitter fit;
     852    fit.SetSignalIntegralMax(sigint);
     853    fit.SetSignalFitMax(sigmax);
     854    fit.SetBackgroundFitMin(bgmin);
     855    fit.SetBackgroundFitMax(bgmax);
     856    fit.SetPolynomOrder(polynom);
     857
     858    TH1D *h=0;
    850859    for (int ix=1; ix<=nx; ix++)
    851860        for (int iy=1; iy<=ny; iy++)
    … …  
    855864            h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
    856865
     866            if (h->GetBinContent(1)==0)
     867                continue;
     868
     869
     870            h->Scale(entries/h->GetEntries());
     871
     872            if (!fit.Fit(*h))
     873                continue;
     874
    857875            const Double_t alpha0 = h->GetBinContent(1);
    858 
    859             // Check for the regios which is not filled...
    860             if (alpha0==0)
    861                 continue;
    862 
    863             h->Scale(entries/h->GetEntries());
    864 
    865876            if (alpha0>maxalpha0)
    866877                maxalpha0=alpha0;
    867 
    868             // First fit a polynom in the off region
    869             func.FixParameter(0, 0);
    870             func.FixParameter(2, 1);
    871             func.ReleaseParameter(3);
    872 
    873             for (int i=5; i<func.GetNpar(); i++)
    874                 func.ReleaseParameter(i);
    875 
    876             h->Fit(&func, "N0Q", "", bgmin, bgmax);
    877 
    878             h4a.Fill(func.GetChisquare());
    879             //h5a.Fill(func.GetProb());
    880 
    881             //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    882             //func.SetParLimits(2, 5, 90);
    883 
    884             func.ReleaseParameter(0);
    885             //func.ReleaseParameter(1);
    886             func.ReleaseParameter(2);
    887             func.FixParameter(3, func.GetParameter(3));
    888             for (int i=5; i<func.GetNpar(); i++)
    889                 func.FixParameter(i, func.GetParameter(i));
    890 
    891             // Do not allow signals smaller than the background
    892             const Double_t A  = alpha0-func.GetParameter(3);
    893             const Double_t dA = TMath::Abs(A);
    894             func.SetParLimits(0, -dA*4, dA*4);
    895             func.SetParLimits(2, 0, 90);
    896 
    897             // Now fit a gaus in the on region on top of the polynom
    898             func.SetParameter(0, A);
    899             func.SetParameter(2, sigmax*0.75);
    900 
    901             h->Fit(&func, "N0Q", "", 0, sigmax);
    902 
    903             TArrayD p(func.GetNpar(), func.GetParameters());
    904 
     878 
    905879            // Fill results into some histograms
    906             h0a.Fill(p[0]);
    907             h0b.Fill(p[3]);
    908             h1.Fill(p[1]);
    909             h2.Fill(p[2]);
     880            h0a.Fill(fit.GetGausA()); // gaus-A
     881            h0b.Fill(fit.GetCoefficient(3)); // 3
     882            h1.Fill(fit.GetGausMu());  // mu
     883            h2.Fill(fit.GetGausSigma());  // sigma-gaus
    910884            if (polynom>1)
    911                 h3.Fill(p[5]);
    912             h4b.Fill(func.GetChisquare());
    913             //h5b.Fill(func.GetProb());
    914 
    915             // Implementing the integral as analytical function
    916             // gives the same result in the order of 10e-5
    917             // and it is not faster at all...
    918             //const Bool_t ok = /*p[0]>=0 && /*p[0]<alpha0*2 &&*/ p[2]>1.75;// && p[2]<88.5;
    919             /*
    920             if (p[0]<-fabs(alpha0-p[3])*7 && p[2]<alphaw/3)
    921             {
    922                 func.SetParameter(0, alpha0-p[3]);
    923                 cout << p[2] << " " << p[0] << " " << alpha0-p[3] << endl;
    924             }
    925             */
    926 
    927             // The fitted function returned units of
    928             // counts bin binwidth. To get the correct number
    929             // of events we must adapt the functions by dividing
    930             // the result of the integration by the bin-width
    931             const Double_t s = func.Integral(0, sigint)/w;
    932 
    933             func.SetParameter(0, 0);
    934             func.SetParameter(2, 1);
    935 
    936             const Double_t b   = func.Integral(0, sigint)/w;
    937             const Double_t sig = SignificanceLiMa(s, b);
     885                h3.Fill(fit.GetCoefficient(5));
     886            h4b.Fill(fit.GetChiSqSignal());
     887
     888            const Double_t sig = fit.GetSignificance();
     889            const Double_t b   = fit.GetEventsBackground();
     890            const Double_t s   = fit.GetEventsSignal();
    938891
    939892            const Int_t n = hist->GetBin(ix, iy);
    … …  
    945898                h6.Fill(sig);
    946899
    947             if (sig>maxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
     900            if (sig>maxs)
    948901            {
    949902                maxs = sig;
    950903                maxx = ix;
    951904                maxy = iy;
    952                 maxpar = p;
     905                maxpar = fit.GetCoefficients();
    953906            }
    954907        }
    955 
    956     *fLog << inf << "done." << endl;
    957908
    958909    h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
Note: See TracChangeset for help on using the changeset viewer.