Changeset 7151 for trunk/MagicSoft/Mars/mhflux

Timestamp:

06/13/05 10:20:58 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MAlphaFitter.cc (modified) (11 diffs)
• MAlphaFitter.h (modified) (6 diffs)
• MHAlpha.cc (modified) (4 diffs)
• MHAlpha.h (modified) (1 diff)
• MHThetaSq.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

@@ -135 +135 @@
 
     for (int i=5; i<fFunc->GetNpar(); i++)
-        fFunc->ReleaseParameter(i);
+        if (fFitBackground)
+            fFunc->ReleaseParameter(i);
+        else
+            fFunc->SetParameter(i, 0);
 
     // options : N  do not store the function, do not draw
@@ -142 +145 @@
     //           Q  quiet mode
     //           E  Perform better Errors estimation using Minos technique
-    h.Fit(fFunc, "NQI", "", bgmin, bgmax);
-
-    fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
+    if (fFitBackground)
+    {
+        h.Fit(fFunc, "NQI", "", bgmin, bgmax);
+        fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
+    }
+
 
     // ------------------------------------
-    if (paint)
+    if (paint && fFitBackground)
     {
         fFunc->SetRange(0, 90);
@@ -209 +215 @@
     fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
 
+    if (TMath::IsNaN(fSignificance))
+        fSignificance=0;
+
     if (fEventsExcess<0)
         fEventsExcess=0;
@@ -240 +249 @@
     fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);
 
+    if (TMath::IsNaN(fSignificance))
+        fSignificance=0;
     if (fEventsExcess<0)
         fEventsExcess=0;
-/*
-    TF1 func("", "gaus(0)+pol0(3)", 0., 90.);
-
-    const Double_t A  = fEventsSignal/bin;
-    const Double_t dA = TMath::Abs(A);
-    func.SetParLimits(0, -dA*4, dA*4);
-    func.SetParLimits(2, 0, 90);
-    func.SetParLimits(3, -dA, dA);
-
-    func.SetParameter(0, A);
-    func.FixParameter(1, 0);
-    func.SetParameter(2, fSigMax*0.75);
-    func.SetParameter(3, 0);
-
-    // options : N  do not store the function, do not draw
-    //           I  use integral of function in bin rather than value at bin center
-    //           R  use the range specified in the function range
-    //           Q  quiet mode
-    //           E  Perform better Errors estimation using Minos technique
-    TH1D h(hon);
-    h.Add(&hof, -1);
-    h.Fit(&func, "NQI", "", 0, 90);
-
-    fChiSqSignal = func.GetChisquare()/func.GetNDF();
-
-    const Int_t bin1 = h.GetXaxis()->FindFixBin(func.GetParameter(2)*2);
-
-    fChiSqBg = 0;
-    for (int i=bin1; i<=h.GetNbinsX(); i++)
-    {
-        const Float_t val = h.GetBinContent(i);
-        fChiSqBg = val*val;
-    }
-    if (fChiSqBg>0)
-        fChiSqBg /= h.GetNbinsX()+1-bin1;
-
-    fCoefficients.Set(func.GetNpar(), func.GetParameters());
-
-    // ------------------------------------
-    if (paint)
-    {
-        func.SetLineColor(kBlue);
-        func.SetLineWidth(2);
-        func.Paint("same");
-    }
-    // ------------------------------------
-  */
+
     return kTRUE;
 }
@@ -323 +288 @@
     f.fScaleMax     = fScaleMax;
     f.fPolynomOrder = fPolynomOrder;
+    f.fFitBackground= fFitBackground;
+    f.fSignalFunc   = fSignalFunc;
     f.fScaleMode    = fScaleMode;
     f.fScaleUser    = fScaleUser;
@@ -350 +317 @@
 {
     *fLog << GetDescriptor() << ": Fitting..." << endl;
-    *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
     *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
-    *fLog << " ...polynom order " << fPolynomOrder << endl;
-    *fLog << " ...scale mode: ";
-    switch (fScaleMode)
-    {
-    case kNone:        *fLog << "none.";         break;
-    case kEntries:     *fLog << "entries.";      break;
-    case kIntegral:    *fLog << "integral.";     break;
-    case kOffRegion:   *fLog << "off region (integral between " << fScaleMin << " and " << fScaleMax << ")"; break;
-    case kBackground:  *fLog << "background (integral between " << fBgMin    << " and " << fBgMax    << ")"; break;
-    case kLeastSquare: *fLog << "least square (N/A)"; break;
-    case kUserScale:   *fLog << "user def (" << fScaleUser << ")"; break;
+    *fLog << " ...signal function: ";
+    switch (fSignalFunc)
+    {
+    case kGauss:       *fLog << "gauss(x)";        break;
+    case kThetaSq:     *fLog << "gauss(sqrt(x))";  break;
     }
     *fLog << endl;
+    if (!fFitBackground)
+        *fLog << " ...no background." << endl;
+    else
+    {
+        *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
+        *fLog << " ...polynom order " << fPolynomOrder << endl;
+        *fLog << " ...scale mode: ";
+        switch (fScaleMode)
+        {
+        case kNone:        *fLog << "none.";         break;
+        case kEntries:     *fLog << "entries.";      break;
+        case kIntegral:    *fLog << "integral.";     break;
+        case kOffRegion:   *fLog << "off region (integral between " << fScaleMin << " and " << fScaleMax << ")"; break;
+        case kBackground:  *fLog << "background (integral between " << fBgMin    << " and " << fBgMax    << ")"; break;
+        case kLeastSquare: *fLog << "least square (N/A)"; break;
+        case kUserScale:   *fLog << "user def (" << fScaleUser << ")"; break;
+        }
+        *fLog << endl;
+    }
 
     if (TString(o).Contains("result"))
@@ -402 +381 @@
 }
 
+Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint)
+{
+    const TString name(Form("TempAlphaTime%06d", gRandom->Integer(1000000)));
+
+    hon.GetZaxis()->SetRange(bin,bin);
+    TH1D *h = (TH1D*)hon.Project3D("ye");
+    hon.GetZaxis()->SetRange(-1,9999);
+
+    h->SetDirectory(0);
+
+    const Bool_t rc = Fit(*h, paint);
+    delete h;
+    return rc;
+}
+
 Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
 {
@@ -448 +442 @@
     return rc;
 }
-
+/*
+Bool_t MAlphaFitter::FitTime(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
+{
+    const TString name1(Form("TempAlphaTime%06d_on",  gRandom->Integer(1000000)));
+    const TString name0(Form("TempAlphaTime%06d_off", gRandom->Integer(1000000)));
+
+    hon.GetZaxis()->SetRange(bin,bin);
+    TH1D *h1 = (TH1D*)hon.Project3D("ye");
+    hon.GetZaxis()->SetRange(-1,9999);
+    h1->SetDirectory(0);
+
+    hof.GetZaxis()->SetRange(bin,bin);
+    TH1D *h0 = (TH1D*)hof.Project3D("ye");
+    hof.GetZaxis()->SetRange(-1,9999);
+    h0->SetDirectory(0);
+
+    const Bool_t rc = ScaleAndFit(*h1, h0, paint);
+
+    delete h0;
+    delete h1;
+
+    return rc;
+}
+*/
 Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
 {
@@ -541 +558 @@
     case kExcess:
         return -GetEventsExcess();
+    case kGaussSigma:
+        return GetGausSigma();
     }
     return 0;
@@ -601 +620 @@
         if (txt==(TString)"excess")
             fStrategy = kExcess;
+        if (txt==(TString)"gausssigma" || txt==(TString)"gaussigma")
+            fStrategy = kGaussSigma;
         rc = kTRUE;
     }
@@ -625 +646 @@
         rc = kTRUE;
     }
+    if (IsEnvDefined(env, prefix, "Signalfunction", print))
+    {
+        TString txt = GetEnvValue(env, prefix, "SignalFunction", "");
+        txt = txt.Strip(TString::kBoth);
+        txt.ToLower();
+        if (txt==(TString)"gauss" || txt==(TString)"gaus")
+            SetSignalFunction(kGauss);
+        if (txt==(TString)"thetasq")
+            SetSignalFunction(kThetaSq);
+        rc = kTRUE;
+    }
     if (IsEnvDefined(env, prefix, "Scale", print))
     {

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

@@ -34 +34 @@
         kSignificanceLogExcess,
         kSignificanceExcess,
-        kExcess
+        kExcess,
+        kGaussSigma
     };
+    enum SignalFunc_t {
+        kGauss, kThetaSq
+    };
+
 private:
-    Float_t fSigInt;
-    Float_t fSigMax;
-    Float_t fBgMin;
-    Float_t fBgMax;
-    Float_t fScaleMin;
-    Float_t fScaleMax;
-    Int_t   fPolynomOrder;
-
-    Double_t fSignificance;
-    Double_t fEventsExcess;
-    Double_t fEventsSignal;
-    Double_t fEventsBackground; // fScaleFactor already applied
-
-    Double_t fChiSqSignal;
-    Double_t fChiSqBg;
-    Double_t fIntegralMax;
-    Double_t fScaleFactor;  // Scale factor for off-data
-
-    TArrayD fCoefficients;
-
-    TF1 *fFunc;
-
-    ScaleMode_t fScaleMode;
-    Double_t fScaleUser;
-
-    Strategy_t fStrategy;
+    // 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 signal
+    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 for off-data
+
+    TArrayD fCoefficients;      // Fit result
+
+    // 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
 
 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), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
+    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),
+        fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)),
+        fScaleMode(kOffRegion), fScaleUser(1), fStrategy(kSignificance)
     {
         fName  = name  ? name  : "MAlphaFitter";
@@ -86 +102 @@
     }
 
+    // TObject
     void Clear(Option_t *o="");
     void Print(Option_t *o="") const; //*MENU*
     void Copy(TObject &o) const;
-    /*
-    TObject *Clone(const char *name) const
-    {
-        return new MAlphaFitter(*this);
-    }
-    */
-
+
+    // Setter
     void SetScaleUser(Float_t scale)       { fScaleUser = scale; fScaleMode=kUserScale; }
     void SetScaleMode(ScaleMode_t mode)    { fScaleMode    = mode; }
@@ -105 +117 @@
     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; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", 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:
+            fFunc=new TF1 ("", Form("[0]*exp(-((sqrt(x)-[1])/[2])^2) + pol%d(3)", fPolynomOrder));
+            break;
+        }
+        fSignalFunc=func;
         gROOT->GetListOfFunctions()->Remove(fFunc);
         fFunc->SetName("Dummy");
-        fCoefficients.Set(3+s+1); fCoefficients.Reset(); }
-
+        fCoefficients.Set(3+fPolynomOrder+1);
+        fCoefficients.Reset();
+    }
+    void EnableBackgroundFit(Bool_t b=kTRUE) { fFitBackground=b; }
+
+    // Getter
     Double_t GetSignalIntegralMax() const  { return fSigInt; }
 
@@ -127 +164 @@
     Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
     const TArrayD &GetCoefficients() const { return fCoefficients; }
-
-    void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const;
-
+    Double_t Eval(Double_t d) const { return fFunc ? fFunc->Eval(d) : 0; }
+
+    // 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);
@@ -149 +186 @@
     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)
@@ -165 +204 @@
     {
         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) const;
+
+    // MTask
     Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
 

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -292 +292 @@
     if (!fResult)
         return kFALSE;
+    fSigma = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "GaussSigma");
+    if (!fSigma)
+        return kFALSE;
 
     //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
@@ -469 +472 @@
             size = fHillas->GetSize();
         energy = fEnergy   ? fEnergy->GetVal() : (fHillas?fHillas->GetSize():1000);
-        theta  = fPointPos ? fPointPos->GetZd()   : 0;
-    }
-
-    //if (size>0)
-    //    alpha /= (2.4 + 1.13*(log10((energy-14)/0.37)-5)*(log10((energy-14)/0.37)-5))/15;
+        theta  = fPointPos ? fPointPos->GetZd() : 0;
+    }
 
     // enhance histogram if necessary
@@ -480 +480 @@
     // Fill histograms
     fHist.Fill(theta, energy, TMath::Abs(alpha), w);
-
-    // Check cuts
-    /*
-    if ( (fEnergyMin>=0 && energy<fEnergyMin) ||
-         (fEnergyMax>=0 && energy>fEnergyMax) ||
-         (fSizeMin>=0   && size  <fSizeMin)   ||
-         (fSizeMax>=0   && size  >fSizeMin) )
-         return kTRUE;
-         */
 
     if (!fSkipHistTime)
@@ -851 +842 @@
 
     fResult->SetVal(fFit.GetMinimizationValue());
+    fSigma->SetVal(fFit.GetGausSigma());
 
     if (!fSkipHistEnergy)

trunk/MagicSoft/Mars/mhflux/MHAlpha.h

@@ -44 +44 @@
 
     MParameterD   *fResult;     //!
+    MParameterD   *fSigma;      //!
     MParameterD   *fEnergy;     //!
     MHillas       *fHillas;     //!

trunk/MagicSoft/Mars/mhflux/MHThetaSq.h

@@ -24 +24 @@
 
     Bool_t SetupFill(const MParList *pl);
-    void InitMapping(MHMatrix *mat, Int_t type=0);
 
     // MParContainer
@@ -31 +30 @@
 public:
     MHThetaSq(const char *name=NULL, const char *title=NULL);
+
+    void InitMapping(MHMatrix *mat, Int_t type=0);
 
     void SetNumBinsSignal(UInt_t n) { fNumBinsSignal=TMath::Max(n, 1U); }