Changeset 1292 for trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc

Timestamp:

04/24/02 12:29:54 (22 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.cc

@@ -28 +28 @@
 //  MHAlphaEnergyTime                                                       //
 //                                                                          //
+//  3D-histogram in alpha, E-est and time                                   //
 //                                                                          //
 //////////////////////////////////////////////////////////////////////////////
@@ -34 +35 @@
 
 #include <TCanvas.h>
+
+#include <math.h>
 
 #include "MHillasSrc.h"
@@ -50 +53 @@
 // --------------------------------------------------------------------------
 //
-// Default Constructor. It sets name and title only. Typically you won't
-// need to change this.
-//
-#include <iostream.h>
+// Default Constructor. It sets name and title of the histogram. 
+//
 MHAlphaEnergyTime::MHAlphaEnergyTime(const char *name, const char *title) 
   : fHist()
@@ -65 +66 @@
     fHist.SetDirectory(NULL);
 
-    fHist.GetXaxis()->SetTitle("\\alpha [\\circ]");
-    fHist.GetYaxis()->SetTitle("E_{est} [GeV]");
-    fHist.GetZaxis()->SetTitle("t [s]");
-}
-
+    fHist.SetTitle("3D-plot of alpha, E-est, time");
+    fHist.SetXTitle("\\alpha [\\circ]");
+    fHist.SetYTitle("E-est [GeV]            ");
+    fHist.SetZTitle("time [s]");
+}
+
+// --------------------------------------------------------------------------
+//
+// Set binnings and prepare filling of the histogram
+// 
 Bool_t MHAlphaEnergyTime::SetupFill(const MParList *plist)
 {
@@ -102 +108 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Fill the histogram
+// 
 Bool_t MHAlphaEnergyTime::Fill(const MParContainer *par)
 {
@@ -110 +120 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Draw the histogram
+// 
 void MHAlphaEnergyTime::Draw(Option_t *opt)
 {
     if (!gPad)
-        MakeDefCanvas("AlphaEnergyTime", "Distrib of alpha, E, t");
+        MakeDefCanvas("AlphaEnergyTime", fTitle);
 
     gPad->Divide(2,2);
@@ -120 +134 @@
 
     gPad->cd(1);
-    h = fHist.Project3D("x");
+    h = fHist.Project3D("ex");
+
+    h->SetTitle("Distribution of \\alpha [\\circ]");
+    h->SetXTitle("\\alpha [\\circ]");
+    h->SetYTitle("Counts");
+
     h->Draw(opt);
     h->SetBit(kCanDelete);
 
     gPad->cd(2);
-    h = fHist.Project3D("y");
-    h->Draw(opt);
-    h->SetBit(kCanDelete);
+    h = fHist.Project3D("ey");
+
+    h->SetTitle("Distribution of E-est [GeV]");
+    h->SetXTitle("E-est [GeV]            ");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+    gPad->SetLogx();
 
     gPad->cd(3);
-    h = fHist.Project3D("z");
+    h = fHist.Project3D("ez");
+
+    h->SetTitle("Distribution of time [s]");
+    h->SetXTitle("time [s]");
+    h->SetYTitle("Counts");
+
     h->Draw(opt);
     h->SetBit(kCanDelete);
@@ -139 +169 @@
     gPad->Modified();
     gPad->Update();
-}
-
+
+}
+
+// --------------------------------------------------------------------------
+//
+// Draw copies of the histogram
+// 
 TObject *MHAlphaEnergyTime::DrawClone(Option_t *opt) const
 {
-    TCanvas *c = MakeDefCanvas("AlphaEnergyTime", "Distrib of alpha, E, t");
-    c->Divide(2, 2);
+    TCanvas &c = *MakeDefCanvas("AlphaEnergyTime", fTitle);
+
+    c.Divide(2, 2);
 
     gROOT->SetSelectedPad(NULL);
 
-    //
-    // FIXME: ProjectionX,Y is not const within root
-    //
     TH1 *h;
 
-    c->cd(1);
-    h = ((TH3D*)(&fHist))->Project3D("x");
-    h->Draw(opt);
-    h->SetBit(kCanDelete);
-
-    c->cd(2);
-    h = ((TH3D*)(&fHist))->Project3D("y");
-    h->Draw(opt);
-    h->SetBit(kCanDelete);
-
-    c->cd(3);
-    h = ((TH3D*)(&fHist))->Project3D("z");
-    h->Draw(opt);
-    h->SetBit(kCanDelete);
-
-    c->cd(4);
+    c.cd(1);
+    h = ((TH3D*)(&fHist))->Project3D("ex");
+
+    h->SetTitle("Distribution of \\alpha [\\circ]");
+    h->SetXTitle("\\alpha [\\circ]");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+
+    c.cd(2);
+    h = ((TH3D*)(&fHist))->Project3D("ey");
+
+    h->SetTitle("Distribution of E-est [GeV]");
+    h->SetXTitle("E-est [GeV]            ");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+    gPad->SetLogx();
+
+    c.cd(3);
+    h = ((TH3D*)(&fHist))->Project3D("ez");
+
+    h->SetTitle("Distribution of time [s]");
+    h->SetXTitle("time [s]");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+
+    c.cd(4);
     ((TH3D&)fHist).DrawCopy(opt);
 
-    c->Modified();
-    c->Update();
-
-    return c;
-}
-
-void MHAlphaEnergyTime::Substract(const TH3D *h1, const TH3D *h2)
-{
-    MH::SetBinning(&fHist, (TH1*)h1);
-
-    fHist.Sumw2();
+    c.Modified();
+    c.Update();
+
+    return &c;
+}
+
+// --------------------------------------------------------------------------
+//
+// Calculate the histogram as the difference of two histograms :
+//          fHist(gamma) = h1(source) - h2(antisource)
+// 
+void MHAlphaEnergyTime::Subtract(const TH3D *h1, const TH3D *h2)
+{
+  //    MH::SetBinning(&fHist, (TH1*)h1);
+
+  //    fHist.Sumw2();
     fHist.Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME!
 }
 
-void MHAlphaEnergyTime::SetAlphaRange(Axis_t lo, Axis_t up)
-{
+// --------------------------------------------------------------------------
+//
+// Integrate fHist(gamma) in the alpha range (lo, up)
+// 
+TH2D *MHAlphaEnergyTime::GetAlphaProjection(Axis_t lo, Axis_t up)
+{
+    if (up < lo)
+    {
+        *fLog << err << fName << ": Alpha projection not possible: lo=" << lo << " up=" << up << endl;
+        return NULL;
+    }
+
     TAxis &axe = *fHist.GetXaxis();
 
-    //
-    // FIXME: ROOT Binning??? of projection
-    // root 3.02: SetRangeUser
-    axe.SetRange(axe.FindFixBin(lo), axe.FindFixBin(up));
-}
-
-TH2D *MHAlphaEnergyTime::GetAlphaProjection(Axis_t lo, Axis_t up)
-{
-    SetAlphaRange(lo, up);
-    return (TH2D*)fHist.Project3D("yz");
-}
+    Int_t ilo = axe.FindFixBin(lo);
+    Int_t iup = axe.FindFixBin(up);
+
+    const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo);
+    const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo;
+
+    const Double_t epsup1 = up-axe.GetBinLowEdge(iup);
+    const Double_t epsup2 = axe.GetBinUpEdge(iup)-up;
+
+    const Double_t epslo = epslo1<epslo2 ? epslo1 : epslo2;
+    const Double_t epsup = epsup1<epsup2 ? epsup1 : epsup2;
+
+    if (epslo1>epslo2)
+        ilo++;
+
+    if (epsup1<epsup2)
+        iup--;
+
+    if (epslo>0.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup))
+    {
+        *fLog << err << fName << ": binning is not adequate for the requested projection:" << endl;
+        *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl;
+        *fLog << " epslo = " << epslo << endl;
+        *fLog << " epsup = " << epsup << endl;
+        *fLog << " dwl   = " << axe.GetBinWidth(ilo) << endl;
+        *fLog << " dwu   = " << axe.GetBinWidth(iup) << endl;
+        return NULL;
+    }
+
+    axe.SetRange(ilo, iup);
+
+    TH2D &h2D = *(TH2D *)fHist.Project3D("ezy");
+
+    h2D.SetTitle("2D-plot  of time vs. E-est");
+    h2D.SetXTitle("E-est [GeV]            ");
+    h2D.SetYTitle("time [s]");
+
+    return &h2D;
+}
+
+//---------------------------------------------------------
+//
+// Draw the projected histogram
+//
+TH2D *MHAlphaEnergyTime::DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt)
+{
+    TH2D *h2D = GetAlphaProjection(lo, up);
+
+    if (!h2D)
+        return NULL;
+
+    char txt[100];
+    sprintf(txt, "No.of Gammas vs. E-est and Time (%.1f < alpha < %.1f deg)", lo, up);
+
+    //    TCanvas *c = MakeDefCanvas("AlphaEnergyTime", "2D histogram of gamma signal in energy and time");
+    TCanvas &c = *MakeDefCanvas("AlphaEnergyTime", txt);
+
+    c.Divide(2, 2);
+
+    gROOT->SetSelectedPad(NULL);
+
+    TH1 *h;
+
+    c.cd(1);
+    h = h2D->ProjectionX("Eest", -1, 9999, "E");
+    h->SetTitle("Distribution of E-est [GeV]");
+    h->SetXTitle("E-est [GeV]            ");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+    gPad->SetLogx();
+
+    c.cd(2);
+    h = h2D->ProjectionY("time", -1, 9999, "E");
+    h->SetTitle("Distribution of time [s]");
+    h->SetXTitle("time [s]");
+    h->SetYTitle("Counts");
+
+    h->Draw(opt);
+    h->SetBit(kCanDelete);
+
+    c.cd(3);
+
+    h2D->DrawCopy(opt);
+    gPad->SetLogx();
+
+    c.Modified();
+    c.Update();
+
+    return h2D;
+}
+