Changeset 6988 for trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

Timestamp:

05/02/05 10:17:41 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (modified) (39 diffs)

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -81 +81 @@
 //
 MHAlpha::MHAlpha(const char *name, const char *title)
-    : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
+    : fNameParameter("MHillasSrc"), fParameter(0),
+    fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
     fPointPos(0), fTimeEffOn(0), fTime(0),
     fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE),
@@ -93 +94 @@
     fTitle = title ? title : "Alpha plot";
 
-    fHAlpha.SetName("Alpha");
-    fHAlpha.SetTitle("Alpha");
-    fHAlpha.SetXTitle("\\Theta [deg]");
-    //fHAlpha.SetYTitle("E_{est} [GeV]");
-    fHAlpha.SetZTitle("|\\alpha| [\\circ]");
-    fHAlpha.SetDirectory(NULL);
-    fHAlpha.UseCurrentStyle();
+    fHist.SetName("Alpha");
+    fHist.SetTitle("Alpha");
+    fHist.SetXTitle("\\Theta [deg]");
+    //fHist.SetYTitle("E_{est} [GeV]");
+    fHist.SetZTitle("|\\alpha| [\\circ]");
+    fHist.SetDirectory(NULL);
+    fHist.UseCurrentStyle();
 
     // Main histogram
-    fHAlphaTime.SetName("Alpha");
-    fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
-    fHAlphaTime.SetYTitle("Counts");
-    fHAlphaTime.UseCurrentStyle();
-    fHAlphaTime.SetDirectory(NULL);
+    fHistTime.SetName("Alpha");
+    fHistTime.SetXTitle("|\\alpha| [\\circ]");
+    fHistTime.SetYTitle("Counts");
+    fHistTime.UseCurrentStyle();
+    fHistTime.SetDirectory(NULL);
 
 
@@ -142 +143 @@
     binse.Apply(fHEnergy);
     binst.Apply(fHTheta);
-    binsa.Apply(fHAlphaTime);
-
-    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
-}
-
-void MHAlpha::FitEnergySpec(Bool_t paint)
-{
-    /*
-    TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
-    f.SetParameter(0,  -2.1);
-    f.SetParameter(1,  1400); // 50
-    f.SetParameter(2,  fHEnergy.Integral()); // 50
-    f.SetLineWidth(2);
-
-    fHEnergy.Fit(&f,  "0QI", "", 100, 2400); // Integral? 90, 2800
-
-    const Double_t chi2 = f.GetChisquare();
-    const Int_t    NDF  = f.GetNDF();
-
-    // was fit successful ?
-    const Bool_t ok = NDF>0 && chi2<2.5*NDF;
-    f.SetLineColor(ok ? kGreen : kRed);
-
-    if (paint)
-    {
-        f.Paint("same");
-
-        if (ok)
-        {
-        TString s;
-        s += Form("\\alpha=%.1f ",   f.GetParameter(0));
-        s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
-        s += Form("N=%.1e",          f.GetParameter(2));
-
-        TLatex text(0.25, 0.94, s.Data());
-        text.SetBit(TLatex::kTextNDC);
-        text.SetTextSize(0.04);
-        text.Paint();
-        }
-    }*/
+    binsa.Apply(fHistTime);
+
+    MH::SetBinning(&fHist, &binst, &binse, &binsa);
+}
+
+Double_t MHAlpha::GetVal() const
+{
+    return static_cast<const MHillasSrc*>(fParameter)->GetAlpha();
 }
 
@@ -189 +158 @@
     MAlphaFitter fit(fFit);
 
-    const Int_t n = fHAlpha.GetNbinsY();
+    const Int_t n = fHist.GetNbinsY();
 
     fHEnergy.SetEntries(0);
@@ -197 +166 @@
     for (int i=1; i<=n; i++)
     {
-        if (fit.FitEnergy(fHAlpha, fOffData, i))
+        if (fit.FitEnergy(fHist, fOffData, i))
         {
             fHEnergy.SetBinContent(i, fit.GetEventsExcess());
@@ -220 +189 @@
     MAlphaFitter fit(fFit);
 
-    const Int_t n = fHAlpha.GetNbinsX();
+    const Int_t n = fHist.GetNbinsX();
 
     fHTheta.SetEntries(0);
@@ -226 +195 @@
     for (int i=1; i<=n; i++)
     {
-        if (fit.FitTheta(fHAlpha, fOffData, i))
+        if (fit.FitTheta(fHist, fOffData, i))
         {
             fHTheta.SetBinContent(i, fit.GetEventsExcess());
@@ -237 +206 @@
 }
 
+Bool_t MHAlpha::GetParameter(const MParList &pl)
+{
+    fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc");
+    if (fParameter)
+        return kTRUE;
+
+    *fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl;
+    return kFALSE;
+}
+
 Bool_t MHAlpha::SetupFill(const MParList *pl)
 {
-    fHAlpha.Reset();
+    fHist.Reset();
     fHEnergy.Reset();
     fHTheta.Reset();
     fHTime.Reset();
 
-    if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
-    {
-        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha");
+    const TString off(Form("%sOff", ClassName()));
+    if (fName!=off && fOffData==NULL)
+    {
+        const TString desc(Form("%s [%s] found... using ", off.Data(), ClassName()));
+        MHAlpha *hoff = (MHAlpha*)pl->FindObject(off, ClassName());
         if (!hoff)
-            *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
+            *fLog << inf << "No " << desc << "current data only!" << endl;
         else
         {
-            *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
+            *fLog << inf << desc << "on-off mode!" << endl;
             SetOffData(*hoff);
         }
     }
+
+    if (!GetParameter(*pl))
+        return kFALSE;
 
     fHillas = 0;
@@ -271 +255 @@
         fHEnergy.SetTitle(" N_{exc} vs. Size ");
         fHEnergy.SetXTitle("Size [phe]");
-        fHAlpha.SetYTitle("Size [phe]");
+        fHist.SetYTitle("Size [phe]");
     }
     else
@@ -277 +261 @@
         fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
         fHEnergy.SetXTitle("E_{est} [GeV]");
-        fHAlpha.SetYTitle("E_{est} [GeV]");
+        fHist.SetYTitle("E_{est} [GeV]");
     }
 
@@ -306 +290 @@
 
     MBinning binst, binse, binsa;
-    binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
-    binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
-    binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');
+    binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
+    binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
+    binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
     if (!fOffData)
     {
@@ -324 +308 @@
                 binse.SetEdges(*pl, "BinningSize");
 
-        binsa.SetEdges(*pl, "BinningAlpha");
+        binsa.SetEdges(*pl, Form("Binning%s", ClassName()+2));
     }
     else
@@ -330 +314 @@
         fHEnergy.SetTitle(fOffData->GetTitle());
         fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
-        fHAlpha.SetYTitle(fOffData->GetYaxis()->GetTitle());
+        fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
     }
 
     binse.Apply(fHEnergy);
     binst.Apply(fHTheta);
-    binsa.Apply(fHAlphaTime);
-    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
+    binsa.Apply(fHistTime);
+    MH::SetBinning(&fHist, &binst, &binse, &binsa);
 
     MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
@@ -402 +386 @@
 
     TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
-    const Bool_t rc = fit.ScaleAndFit(fHAlphaTime, h);
+    const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
     if (h)
         delete h;
@@ -410 +394 @@
 
     // Reset Histogram
-    fHAlphaTime.Reset();
-    fHAlphaTime.SetEntries(0);
+    fHistTime.Reset();
+    fHistTime.SetEntries(0);
 
     //
@@ -469 +453 @@
     else
     {
-        const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
-        if (!par)
-        {
-            *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
-            return kFALSE;
-        }
-
-        alpha  = hil->GetAlpha();
+        alpha  = GetVal();
 
         if (fHillas)
@@ -491 +468 @@
 
     // Fill histograms
-    fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
+    fHist.Fill(theta, energy, TMath::Abs(alpha), w);
 
     // Check cuts
@@ -503 +480 @@
 
     if (!fSkipHistTime)
-        fHAlphaTime.Fill(TMath::Abs(alpha), w);
+        fHistTime.Fill(TMath::Abs(alpha), w);
 
     return kTRUE;
@@ -557 +534 @@
 
         padsave->cd(1);
-        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
+        TH1D *hon = fHist.ProjectionZ("Proj");
         if (fOffData)
         {
-            TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
+            TH1D *hoff = fOffData->ProjectionZ("ProjOff");
             const Double_t alpha = fFit.Scale(*hoff, *hon);
 
@@ -569 +546 @@
             hoff->SetMaximum(alpha);
 
-            if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
+            if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
             {
                 h0->Reset();
@@ -583 +560 @@
     }
 
-    if (o==(TString)"alpha")
-        if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
+    if (o==(TString)"variable")
+        if ((h0 = (TH1D*)gPad->FindObject("Proj")))
         {
             // Check whether we are dealing with on-off analysis
-            TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff");
+            TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
 
             // BE CARFEULL: This is a really weird workaround!
@@ -603 +580 @@
     if (o==(TString)"theta")
     {
-        TH1 *h = (TH1*)gPad->FindObject("Alpha_x");
+        TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
         if (h)
         {
-            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x");
+            TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
             h2->SetDirectory(0);
             h2->Scale(fHTheta.Integral()/h2->Integral());
@@ -618 +595 @@
     if (o==(TString)"energy")
     {
-        TH1 *h = (TH1*)gPad->FindObject("Alpha_y");
+        TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
         if (h)
         {
-            TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y");
+            TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
             h2->SetDirectory(0);
             h2->Scale(fHEnergy.Integral()/h2->Integral());
@@ -634 +611 @@
             gPad->SetLogy();
         }
-        FitEnergySpec(kTRUE);
     }
 
@@ -660 +636 @@
     gPad->SetBorderMode(0);
 
-    h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
+    h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E");
     h->SetBit(TH1::kNoTitle);
     h->SetStats(kTRUE);
-    h->SetXTitle("\\alpha [\\circ]");
+    h->SetXTitle(fHist.GetZaxis()->GetTitle());
     h->SetYTitle("Counts");
     h->SetDirectory(NULL);
@@ -674 +650 @@
         h->SetMarkerColor(kGreen);
 
-        h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
+        h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E");
         h->SetBit(TH1::kNoTitle);
-        h->SetXTitle("\\alpha [\\circ]");
+        h->SetXTitle(fHist.GetZaxis()->GetTitle());
         h->SetYTitle("Counts");
         h->SetDirectory(NULL);
@@ -684 +660 @@
         h->Draw("same");
 
-        h = (TH1D*)h->Clone("ProjAlphaOnOff");
+        h = (TH1D*)h->Clone("ProjOnOff");
         h->SetBit(TH1::kNoTitle);
-        h->SetXTitle("\\alpha [\\circ]");
+        h->SetXTitle(fHist.GetZaxis()->GetTitle());
         h->SetYTitle("Counts");
         h->SetDirectory(NULL);
@@ -696 +672 @@
         TLine lin;
         lin.SetLineStyle(kDashed);
-        lin.DrawLine(0, 0, 90, 0);
+        lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
     }
 
     // After the Alpha-projection has been drawn. Fit the histogram
     // and paint the result into this pad
-    AppendPad("alpha");
+    AppendPad("variable");
 
     if (fHEnergy.GetNbinsX()>1)
@@ -713 +689 @@
         AppendPad("energy");
 
-        h = (TH1D*)fHAlpha.Project3D("y");
+        h = (TH1D*)fHist.Project3D("y");
         h->SetBit(TH1::kNoTitle|TH1::kNoStats);
         h->SetXTitle("E [GeV]");
@@ -749 +725 @@
         AppendPad("theta");
 
-        h = (TH1D*)fHAlpha.Project3D("x");
+        h = (TH1D*)fHist.Project3D("x");
         h->SetBit(TH1::kNoTitle|TH1::kNoStats);
         h->SetXTitle("\\theta [\\circ]");
@@ -768 +744 @@
     // FIXME: Do in Paint if special option given!
     TCanvas *c = new TCanvas;
-    Int_t n = fHAlpha.GetNbinsY();
+    Int_t n = fHist.GetNbinsY();
     Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
     c->Divide(nc, nc, 1e-10, 1e-10);
@@ -775 +751 @@
     MAlphaFitter fit(fFit);
 
-    for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
+    for (int i=1; i<=fHist.GetNbinsY(); i++)
     {
         c->cd(i);
 
-        TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
+        TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E");
         hon->SetBit(TH1::kNoTitle);
         hon->SetStats(kTRUE);
-        hon->SetXTitle("\\alpha [\\circ]");
+        hon->SetXTitle(fHist.GetZaxis()->GetTitle());
         hon->SetYTitle("Counts");
         hon->SetDirectory(NULL);
@@ -796 +772 @@
             hon->SetMarkerColor(kGreen);
 
-            hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
+            hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E");
             hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
-            hof->SetXTitle("\\alpha [\\circ]");
+            hof->SetXTitle(fHist.GetZaxis()->GetTitle());
             hof->SetYTitle("Counts");
             hof->SetDirectory(NULL);
@@ -814 +790 @@
             diff->Add(hof, -1);
             diff->SetBit(TH1::kNoTitle);
-            diff->SetXTitle("\\alpha [\\circ]");
+            diff->SetXTitle(fHist.GetZaxis()->GetTitle());
             diff->SetYTitle("Counts");
             diff->SetDirectory(NULL);
@@ -824 +800 @@
             TLine lin;
             lin.SetLineStyle(kDashed);
-            lin.DrawLine(0, 0, 90, 0);
+            lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
 
             const Float_t min = diff->GetMinimum()*1.05;
@@ -833 +809 @@
             *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
         /*
-        if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
+        if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
         {
             fHEnergy.SetBinContent(i, fit.GetEventsExcess());
@@ -845 +821 @@
 Bool_t MHAlpha::Finalize()
 {
-    //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
+    //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
     //h->SetDirectory(0);
     //Bool_t rc = fFit.Fit(*h);
     //delete h;
 
-    if (!fFit.FitAlpha(fHAlpha, fOffData))
+    if (!fFit.FitAlpha(fHist, fOffData))
     {
         *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
@@ -900 +876 @@
     fMatrix = mat;
 
-    fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
+    fMap[0] = fMatrix->AddColumn(GetParameterRule());
     fMap[1] = -1;
     fMap[2] = -1;