Changeset 5692 for trunk/MagicSoft/Mars/mhflux

Timestamp:

01/03/05 12:02:16 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MAlphaFitter.cc (modified) (7 diffs)
• MAlphaFitter.h (modified) (6 diffs)
• MHAlpha.cc (modified) (22 diffs)
• MHAlpha.h (modified) (6 diffs)

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

@@ -41 +41 @@
 #include <TF1.h>
 #include <TH1.h>
+#include <TH3.h>
+
+#include <TRandom.h>
 
 #include <TLatex.h>
@@ -51 +54 @@
 
 using namespace std;
+
+void MAlphaFitter::Clear(Option_t *o)
+{
+    fSignificance=0;
+    fEventsExcess=0;
+    fEventsSignal=0;
+    fEventsBackground=0;
+
+    fChiSqSignal=0;
+    fChiSqBg=0;
+    fIntegralMax=0;
+
+    fCoefficients.Reset();
+}
 
 // --------------------------------------------------------------------------
@@ -84 +101 @@
 Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
 {
+    Clear();
+    if (h.GetEntries()==0)
+        return kFALSE;
+
     Double_t sigmax=fSigMax;
     Double_t bgmin =fBgMin;
@@ -175 +196 @@
      //fSignificance = MMath::SignificanceLiMaSigned(s, b);
 
-
     const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);
 
@@ -190 +210 @@
 }
 
+Bool_t MAlphaFitter::Fit(TH1D &hon, TH1D &hof, Bool_t paint)
+{
+    /*
+     Clear();
+     if (hon.GetEntries()==0)
+     return kFALSE;
+     */
+
+    TH1D h(hon);
+    h.Add(&hof, -1);
+
+    MAlphaFitter fit(*this);
+    fit.SetPolynomOrder(1);
+
+    if (!fit.Fit(h, paint))
+        return kFALSE;
+
+    fChiSqSignal = fit.GetChiSqSignal();
+    fChiSqBg     = fit.GetChiSqBg();
+    fCoefficients = fit.GetCoefficients();
+
+    const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);
+
+    fIntegralMax      = hon.GetBinLowEdge(bin+1);
+    fEventsBackground = hof.Integral(0, bin);
+    fEventsSignal     = hon.Integral(0, bin);
+    fEventsExcess     = fEventsSignal-fEventsBackground;
+    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);
+
+    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
+    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;
+}
+
 void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
 {
-    TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f  \\omega=%.1f\\circ  E=%d  (\\alpha<%.1f\\circ)  (\\chi_{b}^{2}=%.1f  \\chi_{s}^{2}=%.1f)",
+    TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f  \\omega=%.1f\\circ  E=%d  (\\alpha<%.1f\\circ)  (\\chi_{b}^{2}/ndf=%.1f  \\chi_{s}^{2}/ndf=%.1f  c_{0}=%.1f)",
                            fSignificance, GetGausSigma(),
                            (int)fEventsExcess, fIntegralMax,
-                           fChiSqBg, fChiSqSignal));
+                           fChiSqBg, fChiSqSignal, fCoefficients[3]));
 
     text.SetBit(TLatex::kTextNDC);
@@ -210 +310 @@
     f.fBgMin        = fBgMin;
     f.fBgMax        = fBgMax;
+    f.fScaleMin     = fScaleMin;
+    f.fScaleMax     = fScaleMax;
     f.fPolynomOrder = fPolynomOrder;
+    f.fScaleMode    = fScaleMode;
+    f.fScaleUser    = fScaleUser;
     f.fCoefficients.Set(fCoefficients.GetSize());
     f.fCoefficients.Reset();
@@ -227 +331 @@
     *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.";   break;
+    case kLeastSquare: *fLog << "least square."; break;
+    case kUserScale:   *fLog << "user def (" << fScaleUser << ")"; break;
+    }
+    *fLog << endl;
+
+    if (TString(o).Contains("result"))
+    {
+        *fLog << "Result:" << endl;
+        *fLog << " - Significance           " << fSignificance << endl;
+        *fLog << " - Excess Events          " << fEventsExcess << endl;
+        *fLog << " - Signal Events          " << fEventsSignal << endl;
+        *fLog << " - Background Events      " << fEventsBackground << endl;
+        *fLog << " - Chi^2/ndf (Signal)     " << fChiSqSignal << endl;
+        *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl;
+        *fLog << " - Signal integrated      " << fIntegralMax << "°" << endl;
+    }
+}
+
+Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
+{
+    const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
+    TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E");
+    h->SetDirectory(0);
+
+    const Bool_t rc = Fit(*h, paint);
+    delete h;
+    return rc;
+}
+
+Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
+{
+    const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000)));
+    TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E");
+    h->SetDirectory(0);
+
+    const Bool_t rc = Fit(*h, paint);
+    delete h;
+    return rc;
+}
+
+Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
+{
+    const TString name(Form("TempAlpha%06d", gRandom->Integer(1000000)));
+    TH1D *h = hon.ProjectionZ(name, -1, 9999, -1, 9999, "E");
+    h->SetDirectory(0);
+
+    const Bool_t rc = Fit(*h, paint);
+    delete h;
+    return rc;
+}
+
+Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
+{
+    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
+    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
+    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
+    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
+    h1->SetDirectory(0);
+    h0->SetDirectory(0);
+
+    Scale(*h0, *h1);
+
+    const Bool_t rc = Fit(*h1, *h0, paint);
+
+    delete h0;
+    delete h1;
+
+    return rc;
+}
+
+Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
+{
+    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
+    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
+
+    TH1D *h1 = hon.ProjectionZ(name1, bin, bin, -1, 9999, "E");
+    TH1D *h0 = hof.ProjectionZ(name0, bin, bin, -1, 9999, "E");
+    h1->SetDirectory(0);
+    h0->SetDirectory(0);
+
+    Scale(*h0, *h1);
+
+    const Bool_t rc = Fit(*h1, *h0, paint);
+
+    delete h0;
+    delete h1;
+
+    return rc;
+}
+
+Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
+{
+    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
+    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
+
+    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E");
+    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E");
+    h1->SetDirectory(0);
+    h0->SetDirectory(0);
+
+    Scale(*h0, *h1);
+
+    const Bool_t rc = Fit(*h1, *h0, paint);
+
+    delete h0;
+    delete h1;
+
+    return rc;
+}
+
+void MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
+{
+    Float_t scaleon = 1;
+    Float_t scaleof = 1;
+    switch (fScaleMode)
+    {
+    case kNone:
+        return;
+
+    case kEntries:
+        scaleon = on.GetEntries();
+        scaleof = of.GetEntries();
+        break;
+
+    case kIntegral:
+        scaleon = on.Integral();
+        scaleof = of.Integral();
+        break;
+
+    case kOffRegion:
+        {
+            const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
+            const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
+            scaleon = on.Integral(min, max);
+            scaleof = of.Integral(min, max);
+        }
+        break;
+
+    case kUserScale:
+        scaleon = fScaleUser;
+        break;
+
+    default:
+        return;
+    }
+
+    if (scaleof!=0)
+        of.Scale(scaleon/scaleof);
+}

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

@@ -15 +15 @@
 
 class TH1D;
+class TH3D;
 
 class MAlphaFitter : public MParContainer
 {
+public:
+    enum ScaleMode_t {
+        kNone,
+        kEntries,
+        kIntegral,
+        kOffRegion,
+        kLeastSquare,
+        kUserScale
+    };
 private:
     Float_t fSigInt;
@@ -23 +33 @@
     Float_t fBgMin;
     Float_t fBgMax;
+    Float_t fScaleMin;
+    Float_t fScaleMax;
     Int_t   fPolynomOrder;
 
@@ -38 +50 @@
     TF1 *fFunc;
 
+    ScaleMode_t fScaleMode;
+    Double_t fScaleUser;
+
 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), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90))
+    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(kEntries), fScaleUser(1)
     {
         fName  = name  ? name  : "MAlphaFitter";
@@ -58 +73 @@
     }
 
-    void Print(Option_t *o=0) const;
+    void Clear(Option_t *o="");
+    void Print(Option_t *o="") const;
     void Copy(TObject &o) const;
     /*
@@ -67 +83 @@
     */
 
+    void SetScaleUser(Float_t scale)       { fScaleUser = scale; fScaleMode=kUserScale; }
+    void SetScaleMode(ScaleMode_t mode)    { fScaleMode    = mode; }
     void SetSignalIntegralMax(Float_t s)   { fSigInt       = s; }
     void SetSignalFitMax(Float_t s)        { fSigMax       = s; }
     void SetBackgroundFitMin(Float_t s)    { fBgMin        = s; }
     void SetBackgroundFitMax(Float_t s)    { fBgMax        = s; }
+    void SetScaleMin(Float_t s)            { fScaleMin     = s; }
+    void SetScaleMax(Float_t s)            { fScaleMax     = s; }
     void SetPolynomOrder(Int_t s)          { fPolynomOrder = s; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", s));
         gROOT->GetListOfFunctions()->Remove(fFunc);
@@ -93 +113 @@
 
     Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
+    Bool_t Fit(TH1D &on, TH1D &off, Bool_t paint=kFALSE);
+    Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
+    {
+        return off ? Fit(on, *off, paint) : Fit(on, paint);
+    }
+
+    Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE);
+    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 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 FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE)
+    {
+        return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint);
+    }
+    Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
+    {
+        return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint);
+    }
+    Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
+    {
+        return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint);
+    }
+
+    void Scale(TH1D &off, const TH1D &on) const;
 
     ClassDef(MAlphaFitter, 1)

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -75 +75 @@
 //
 MHAlpha::MHAlpha(const char *name, const char *title)
-    : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0),
-    fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0)
+    : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
+    fPointPos(0), fTimeEffOn(0), fTime(0),
+    fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE),
+    //fEnergyMin(-1), fEnergyMax(-1), fSizeMin(-1), fSizeMax(-1),
+    fMatrix(0)
 {
     //
@@ -100 +103 @@
 
 
-    fHEnergy.SetName("Energy");
-    fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
-    fHEnergy.SetXTitle("E_{est} [GeV]");
+    //fHEnergy.SetName("Energy");
+    //fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
+    //fHEnergy.SetXTitle("E_{est} [GeV]");
     fHEnergy.SetYTitle("N_{exc}");
     fHEnergy.SetDirectory(NULL);
@@ -184 +187 @@
     for (int i=1; i<=n; i++)
     {
-        TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":"");
-        if (fit.Fit(*h))
+        if (fit.FitEnergy(fHAlpha, fOffData, i))
         {
             fHEnergy.SetBinContent(i, fit.GetEventsExcess());
             fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
         }
-        delete h;
     }
 }
@@ -203 +204 @@
     for (int i=1; i<=n; i++)
     {
-        TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":"");
-        if (fit.Fit(*h))
+        if (fit.FitTheta(fHAlpha, fOffData, i))
         {
             fHTheta.SetBinContent(i, fit.GetEventsExcess());
             fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
         }
-        delete h;
     }
 }
@@ -219 +218 @@
     fHTheta.Reset();
 
+    if (fName!=(TString)"MHAlphaOff")
+    {
+        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff");
+        if (!hoff)
+            *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
+        else
+        {
+            *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
+            SetOffData(*hoff);
+        }
+    }
+
+    fHillas = 0;
+    /*
+    if (fSizeMin>=0 || fSizeMax>=0)
+    {
+        fHillas = (MHillas*)pl->FindObject("MHillas");
+        if (!fHillas)
+        {
+            *fLog << warn << "Size cut set, but MHillas not found... abort." << endl;
+            return kFALSE;
+        }
+    }
+    */
     fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
     if (!fEnergy)
-        *fLog << warn << "MEnergyEst not found... ignored." << endl;
+    { /*
+        if (fEnergyMin>=0 || fEnergyMax>=0)
+        {
+            *fLog << warn << "Energy cut set, but MEnergyEst not found... abort." << endl;
+            return kFALSE;
+        } */
+
+        *fLog << warn << "MEnergyEst not found... " << flush;
+
+        if (!fHillas)
+            fHillas = (MHillas*)pl->FindObject("MHillas");
+        if (fHillas)
+            *fLog << "using SIZE instead!" << endl;
+        else
+            *fLog << "ignored." << endl;
+
+        fHEnergy.SetName("Size");
+        fHEnergy.SetTitle(" N_{exc} vs. Size ");
+        fHEnergy.SetXTitle("Size [\\gamma]");
+    }
+    else
+    {
+        fHEnergy.SetName("Energy");
+        fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
+        fHEnergy.SetXTitle("E_{est} [GeV]");
+    }
 
     fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
@@ -230 +278 @@
     if (!fTimeEffOn)
         *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;
+    else
+        *fTimeEffOn = MTime(); // FIXME: How to do it different?
 
     fTime = (MTime*)pl->FindObject("MTime");
@@ -246 +296 @@
 
     MBinning binst, binse, binsa;
-    binst.SetEdges(fHAlpha, 'x');
-    binse.SetEdges(fHAlpha, 'y');
-    binsa.SetEdges(fHAlpha, 'z');
-
-    MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
-    if (fPointPos && bins)
-        binst.SetEdges(*bins);
-    if (!fPointPos)
-        binst.SetEdges(1, 0, 90);
-
-    bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
-    if (fEnergy && bins)
-        binse.SetEdges(*bins);
-    if (!fEnergy)
-        binse.SetEdges(1, 10, 100000);
-
-    bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
-    if (bins)
-        binsa.SetEdges(*bins);
+    binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
+    binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
+    binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');
+
+    if (!fOffData)
+    {
+        MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
+        if (fPointPos && bins)
+            binst.SetEdges(*bins);
+        if (!fPointPos)
+            binst.SetEdges(1, 0, 60);
+
+        bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
+        if ((fEnergy||fHillas) && bins)
+            binse.SetEdges(*bins);
+        if (!fEnergy && !fHillas)
+            binse.SetEdges(1, 10, 100000);
+
+        bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
+        if (bins)
+            binsa.SetEdges(*bins);
+    }
 
     binse.Apply(fHEnergy);
@@ -321 +374 @@
     }
 
+    // Check new 'last time'
+    MTime *time = final ? fTime : fTimeEffOn;
+
+    if (time->GetAxisTime()<=fHTime.GetXaxis()->GetXmax())
+    {
+        *fLog << warn << "WARNING - New time-stamp " << *time << " lower" << endl;
+        *fLog << "than upper edge of histogram... skipped." << endl;
+        *fLog << "This should not happen. Maybe you started you eventloop with" << endl;
+        *fLog << "an already initialized time stamp MTimeEffectiveOnTime?" << endl;
+        rebin++;
+        return;
+    }
+
+    // Fit time histogram
     MAlphaFitter fit(fFit);
-    if (!fit.Fit(fHAlphaTime))
+
+    TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
+    const Bool_t rc = fit.Fit(fHAlphaTime, h);
+    if (h)
+        delete h;
+
+    if (!rc)
         return;
 
@@ -334 +407 @@
     // Prepare Histogram
     //
-
-    MTime *time = final ? fTime : fTimeEffOn;
     if (final)
         time->Plus1ns();
@@ -363 +434 @@
 {
     Double_t alpha, energy, theta;
+    Double_t size=-1;
 
     if (fMatrix)
     {
         alpha  = (*fMatrix)[fMap[0]];
-        energy = 1000; 
-        theta  =    0; 
+        energy = fMap[1]<0 ? -1 : (*fMatrix)[fMap[1]];
+        size   = fMap[2]<0 ? -1 : (*fMatrix)[fMap[2]];
+        //<0 ? 1000 : (*fMatrix)[fMap[1]];
+        theta  = 0;
+
+        if (energy<0)
+            energy=size;
+        if (size<0)
+            size=energy;
+
+        if (energy<0 && size<0)
+            energy = size = 1000;
     }
     else
@@ -380 +462 @@
 
         alpha  = hil->GetAlpha();
-        energy = fEnergy   ? fEnergy->GetEnergy() : 1000;
+        if (fHillas)
+            size = fHillas->GetSize();
+        energy = fEnergy   ? fEnergy->GetEnergy() : (fHillas?fHillas->GetSize():1000);
         theta  = fPointPos ? fPointPos->GetZd()   : 0;
     }
 
+    // enhance histogram if necessary
     UpdateAlphaTime();
 
+    // Fill histograms
     fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
-    fHAlphaTime.Fill(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)
+        fHAlphaTime.Fill(TMath::Abs(alpha), w);
 
     return kTRUE;
@@ -441 +538 @@
 
         padsave->cd(1);
-        fHAlpha.ProjectionZ(fNameProjAlpha);
+        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
+        if (fOffData)
+        {
+            TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
+            fFit.Scale(*hoff, *hon);
+
+            hon->SetMaximum();
+            hon->SetMaximum(TMath::Max(hon->GetMaximum(), hoff->GetMaximum())*1.05);
+
+            if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
+            {
+                h0->Reset();
+                h0->Add(hoff, hon, -1);
+                const Float_t min = h0->GetMinimum()*1.05;
+                hon->SetMinimum(min<0 ? min : 0);
+            }
+        }
         FitEnergyBins();
         FitThetaBins();
@@ -447 +560 @@
 
     if (o==(TString)"alpha")
-        if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha)))
+        if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
         {
             // Do not store the 'final' result in fFit
             MAlphaFitter fit(fFit);
-            fit.Fit(*h0, kTRUE);
+            TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff");
+            fit.Fit(*h0, hoff, kTRUE);
             fit.PaintResult();
         }
@@ -459 +573 @@
 
     if (o==(TString)"theta")
+    {
+        TH1 *h = (TH1*)gPad->FindObject("Alpha_x");
+        if (h)
+        {
+            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x");
+            h2->SetDirectory(0);
+            h2->Scale(fHTheta.Integral()/h2->Integral());
+            h->Reset();
+            h->Add(h2);
+            delete h2;
+        }
         PaintText(fHTheta.Integral(), 0);
+    }
 
     if (o==(TString)"energy")
     {
+        TH1 *h = (TH1*)gPad->FindObject("Alpha_y");
+        if (h)
+        {
+            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y");
+            h2->SetDirectory(0);
+            h2->Scale(fHEnergy.Integral()/h2->Integral());
+            h->Reset();
+            h->Add(h2);
+            delete h2;
+        }
+
         if (fHEnergy.GetMaximum()>1)
         {
@@ -469 +606 @@
         }
         FitEnergySpec(kTRUE);
-
     }
 
@@ -494 +630 @@
     pad->cd(1);
     gPad->SetBorderMode(0);
-    h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E");
+
+    h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
     h->SetBit(TH1::kNoTitle);
+    h->SetStats(kTRUE);
     h->SetXTitle("\\alpha [\\circ]");
     h->SetYTitle("Counts");
@@ -501 +639 @@
     h->SetMarkerStyle(kFullDotMedium);
     h->SetBit(kCanDelete);
-    h->Draw();
+    h->Draw("");
+
+    if (fOffData)
+    {
+        h->SetMarkerColor(kGreen);
+
+        h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
+        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
+        h->SetXTitle("\\alpha [\\circ]");
+        h->SetYTitle("Counts");
+        h->SetDirectory(NULL);
+        h->SetMarkerStyle(kFullDotMedium);
+        h->SetBit(kCanDelete);
+        h->SetMarkerColor(kRed);
+        h->Draw("same");
+
+        h = (TH1D*)h->Clone("ProjAlphaOnOff");
+        h->SetBit(TH1::kNoTitle);
+        h->SetXTitle("\\alpha [\\circ]");
+        h->SetYTitle("Counts");
+        h->SetDirectory(NULL);
+        h->SetMarkerStyle(kFullDotMedium);
+        h->SetBit(kCanDelete);
+        h->SetMarkerColor(kBlue);
+        h->Draw("same");
+
+        TLine lin;
+        lin.SetLineStyle(kDashed);
+        lin.DrawLine(0, 0, 90, 0);
+    }
+
     // After the Alpha-projection has been drawn. Fit the histogram
     // and paint the result into this pad
@@ -511 +679 @@
         gPad->SetBorderMode(0);
         fHEnergy.Draw();
+
         AppendPad("energy");
+
+        h = (TH1D*)fHAlpha.Project3D("y");
+        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
+        h->SetXTitle("E [GeV]");
+        h->SetYTitle("Counts");
+        h->SetDirectory(NULL);
+        h->SetMarkerStyle(kFullDotMedium);
+        h->SetBit(kCanDelete);
+        h->SetMarkerColor(kCyan);
+        h->SetLineColor(kCyan);
+        h->Draw("Psame");
     }
     else
@@ -531 +711 @@
         gPad->SetBorderMode(0);
         fHTheta.Draw();
+
         AppendPad("theta");
+
+        h = (TH1D*)fHAlpha.Project3D("x");
+        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
+        h->SetXTitle("\\theta [\\circ]");
+        h->SetYTitle("Counts");
+        h->SetDirectory(NULL);
+        h->SetMarkerStyle(kFullDotMedium);
+        h->SetBit(kCanDelete);
+        h->SetMarkerColor(kCyan);
+        h->SetLineColor(kCyan);
+        h->Draw("Psame");
     }
     else
         delete pad->GetPad(4);
-
-}
+}
+
+void MHAlpha::DrawAll()
+{
+    // FIXME: Do in Paint if special option given!
+    TCanvas *c = new TCanvas;
+    Int_t n = fHAlpha.GetNbinsY();
+    Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
+    c->Divide(nc, nc, 0, 0);
+
+    // Do not store the 'final' result in fFit
+    MAlphaFitter fit(fFit);
+
+    for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
+    {
+        c->cd(i);
+
+        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
+        hon->SetBit(TH1::kNoTitle);
+        hon->SetStats(kTRUE);
+        hon->SetXTitle("\\alpha [\\circ]");
+        hon->SetYTitle("Counts");
+        hon->SetDirectory(NULL);
+        hon->SetMarkerStyle(kFullDotMedium);
+        hon->SetBit(kCanDelete);
+        hon->Draw("");
+
+        TH1D *hof = 0;
+
+        if (fOffData)
+        {
+            hon->SetMarkerColor(kGreen);
+
+            hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
+            hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
+            hof->SetXTitle("\\alpha [\\circ]");
+            hof->SetYTitle("Counts");
+            hof->SetDirectory(NULL);
+            hof->SetMarkerStyle(kFullDotMedium);
+            hof->SetBit(kCanDelete);
+            hof->SetMarkerColor(kRed);
+            hof->Draw("same");
+
+            fit.Scale(*hof, *hon);
+
+            hon->SetMaximum();
+            hon->SetMaximum(TMath::Max(hon->GetMaximum(), hof->GetMaximum())*1.05);
+
+            TH1D *diff = new TH1D(*hon);
+            diff->Add(hof, -1);
+            diff->SetBit(TH1::kNoTitle);
+            diff->SetXTitle("\\alpha [\\circ]");
+            diff->SetYTitle("Counts");
+            diff->SetDirectory(NULL);
+            diff->SetMarkerStyle(kFullDotMedium);
+            diff->SetBit(kCanDelete);
+            diff->SetMarkerColor(kBlue);
+            diff->Draw("same");
+
+            TLine lin;
+            lin.SetLineStyle(kDashed);
+            lin.DrawLine(0, 0, 90, 0);
+
+            const Float_t min = diff->GetMinimum()*1.05;
+            hon->SetMinimum(min<0 ? min : 0);
+        }
+
+        if (hof ? fit.Fit(*hon, *hof) : fit.Fit(*hon))
+                *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
+        /*
+        if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
+        {
+            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
+            fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
+        }*/
+    }
+
+}
+
 
 Bool_t MHAlpha::Finalize()
 {
-    // Store the final result in fFit
-    TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
-    h->SetDirectory(0);
-    fFit.Print();
-    Bool_t rc = fFit.Fit(*h);
-    delete h;
-    if (!rc)
+    //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
+    //h->SetDirectory(0);
+    //Bool_t rc = fFit.Fit(*h);
+    //delete h;
+    if (!fFit.FitAlpha(fHAlpha, fOffData))
     {
         *fLog << warn << "Histogram empty." << endl;
@@ -552 +819 @@
     }
 
+    // Store the final result in fFit
+    fFit.Print("result");
+
     if (fResult)
         fResult->SetVal(fFit.GetSignificance());
 
-    FitEnergyBins();
-    FitThetaBins();
-    UpdateAlphaTime(kTRUE);
-    MH::RemoveFirstBin(fHTime);
+    if (!fSkipHistEnergy)
+    {
+        *fLog << inf << "Processing energy bins..." << endl;
+        FitEnergyBins();
+    }
+    if (!fSkipHistTheta)
+    {
+        *fLog << inf << "Processing theta bins..." << endl;
+        FitThetaBins();
+    }
+    if (!fSkipHistTime)
+    {
+        *fLog << inf << "Processing time bins..." << endl;
+        UpdateAlphaTime(kTRUE);
+        MH::RemoveFirstBin(fHTime);
+    }
 
     return kTRUE;
@@ -572 +854 @@
 // will take the values from the matrix instead of the containers.
 //
-void MHAlpha::InitMapping(MHMatrix *mat)
+// It takes fSkipHist* into account!
+//
+void MHAlpha::InitMapping(MHMatrix *mat, Int_t type)
 {
     if (fMatrix)
@@ -580 +864 @@
 
     fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
-    //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
+    fMap[1] = -1;
+    fMap[2] = -1;
+    fMap[3] = -1;
+    fMap[4] = -1;
+
+    if (!fSkipHistEnergy)
+        if (type==0)
+        {
+            fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
+            fMap[2] = -1;
+        }
+        else
+        {
+            fMap[1] = -1;
+            fMap[2] = fMatrix->AddColumn("MHillas.fSize");
+        }
+
+    if (!fSkipHistTheta)
+        fMap[3] = fMatrix->AddColumn("MPointingPos.fZd");
+
+   // if (!fSkipHistTime)
+   //     fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime");
 }
 

trunk/MagicSoft/Mars/mhflux/MHAlpha.h

@@ -21 +21 @@
 class MParameterD;
 class MEnergyEst;
+class MHillas;
 class MHMatrix;
 class MPointingPos;
@@ -84 +85 @@
 {
 private:
+    const TH3D *fOffData;
+
     MAlphaFitter fFit;          // SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT)
 
@@ -94 +97 @@
     MParameterD  *fResult;      //!
     MEnergyEst   *fEnergy;      //!
+    MHillas      *fHillas;      //!
     MPointingPos *fPointPos;    //!
 
@@ -100 +104 @@
     MTime    fLastTime;         //! Last fTimeEffOn
 
-    const TString fNameProjAlpha;  //! This should make sure, that gROOT doen't confuse the projection with something else
+    //Float_t fEnergyMin;
+    //Float_t fEnergyMax;
+    //Float_t fSizeMin;
+    //Float_t fSizeMax;
+
+    Bool_t fSkipHistTime;
+    Bool_t fSkipHistTheta;
+    Bool_t fSkipHistEnergy;
+
+    //const TString fNameProjAlpha;  //! This should make sure, that gROOT doen't confuse the projection with something else
 
     MHMatrix *fMatrix;          //!
-    Int_t fMap[2];              //!
+    Int_t fMap[5];              //!
 
     void FitEnergySpec(Bool_t paint=kFALSE);
@@ -115 +128 @@
     void PaintText(Double_t val, Double_t error) const;
 
+    Int_t DistancetoPrimitive(Int_t px, Int_t py) { return 0; }
+
 public:
     MHAlpha(const char *name=NULL, const char *title=NULL);
@@ -125 +140 @@
     const MAlphaFitter &GetAlphaFitter() const { return fFit; }
 
+    void SetOffData(const MHAlpha &h)
+    {
+        fOffData = &h.fHAlpha;
+    }
+/*
+    void SetSizeCuts(Float_t min, Float_t max)   { fSizeMin=min; fSizeMax=max; }
+    void SetSizeMin(Float_t min)                 { fSizeMin=min; }
+    void SetSizeMax(Float_t max)                 { fSizeMax=max; }
+    void SetEnergyCuts(Float_t min, Float_t max) { fEnergyMin=min; fEnergyMax=max; }
+    void SetEnergyMin(Float_t min)               { fEnergyMin=min; }
+    void SetEnergyMax(Float_t max)               { fEnergyMax=max; }
+
+    void SetCuts(const MHAlpha &h) {
+        fSizeMin = h.fSizeMin; fEnergyMin = h.fEnergyMin;
+        fSizeMax = h.fSizeMax; fEnergyMax = h.fEnergyMax;
+    }
+    */
+
+    void SkipHistTime(Bool_t b=kTRUE)   { fSkipHistTime=b; }
+    void SkipHistTheta(Bool_t b=kTRUE)  { fSkipHistTheta=b; }
+    void SkipHistEnergy(Bool_t b=kTRUE) { fSkipHistEnergy=b; }
+
     void Paint(Option_t *opt="");
     void Draw(Option_t *option="");
 
-    void InitMapping(MHMatrix *mat);
+    void DrawAll(); //*MENU*
+
+    void InitMapping(MHMatrix *mat, Int_t type=0);
     void StopMapping();