Changeset 1293 for trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc

Timestamp:

04/24/02 13:47:29 (23 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc (modified) (6 diffs)

trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc

@@ -26 +26 @@
 //////////////////////////////////////////////////////////////////////////////
 //                                                                          //
-//  MHEffOnTimeTheta                                                       //
+//  MHEffOnTimeTheta                                                        //
 //                                                                          //
+//  calculates the effective on time for each bin in Theta                  //
 //                                                                          //
 //////////////////////////////////////////////////////////////////////////////
 
 #include "MHEffOnTimeTheta.h"
+
+#include <TStyle.h>
 
 #include <TF1.h>
@@ -47 +50 @@
 ClassImp(MHEffOnTimeTheta);
 
-
 // --------------------------------------------------------------------------
 //
-// Default Constructor. It sets name and title only. Typically you won't
-// need to change this.
+// Default Constructor. It sets name and title of the histograms. 
 //
 MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title)
-    : fHist()
+    : fHEffOn()
 {
     //
@@ -62 +63 @@
     fTitle = title ? title : "1-D histogram of Eff On Time";
 
-    fHist.SetName("EffOn");
-    fHist.SetTitle("Effective On Time Vs. Theta");
-
-    fHist.SetDirectory(NULL);
-
-    fHist.SetXTitle("t_{eff} [s]");
-    fHist.SetYTitle("\\Theta [\\circ]");
-}
-
+    // effective on time versus theta
+    fHEffOn.SetName("EffOn");
+    fHEffOn.SetTitle("Effective On Time vs. Theta");
+
+    fHEffOn.SetDirectory(NULL);
+
+    fHEffOn.SetXTitle("\\Theta [\\circ]");
+    fHEffOn.SetYTitle("t-eff [s]");
+
+    // chi2/NDF versus theta
+    fHChi2.SetName("Chi2/NDF");
+    fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Theta");
+
+    fHChi2.SetDirectory(NULL);
+
+    fHChi2.SetXTitle("\\Theta [\\circ]");
+    fHChi2.SetYTitle("chi2/NDF");
+
+    // lambda versus theta
+    fHLambda.SetName("lambda");
+    fHLambda.SetTitle("lambda of OnTimeFit vs. Theta");
+
+    fHLambda.SetDirectory(NULL);
+
+    fHLambda.SetXTitle("\\Theta [\\circ]");
+    fHLambda.SetYTitle("\\lambda [Hz]");
+
+    // N0del versus theta
+    fHN0del.SetName("N0del");
+    fHN0del.SetTitle("N0del of OnTimeFit vs. Theta");
+
+    fHN0del.SetDirectory(NULL);
+
+    fHN0del.SetXTitle("\\Theta [\\circ]");
+    fHN0del.SetYTitle("N0del");
+}
+
+// -----------------------------------------------------------------------
+//
+// Calculate the effective on time by fitting the distribution of
+// time differences
+//
 void MHEffOnTimeTheta::Calc(TH2D *hist)
 {
+    // nbins = number of Theta bins
     const Int_t nbins = hist->GetNbinsY();
 
-    for (int i=1; i<nbins; i++)
+    for (int i=1; i<=nbins; i++)
     {
-        //char txt[100];
-        //sprintf(txt, "Name%d", 100*i);
-        //new TCanvas(txt, "Title");
-
-        TH1D &h = *hist->ProjectionX("dTime-Distribution for fixed Time", i, i);
-
-        //hist->Draw();
-        //gPad->SetLogy();
-
-        Double_t Nmdel = h.Integral("width");
-        Double_t mean  = h.GetMean();
-
-        TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)",
-                 mean, hist->GetXaxis()->GetXmax());
-
-        func.SetParameter(0, 100); // [Hz]
-        func.SetParameter(1, Nmdel);
-
-        func.SetParLimits(0, 0, 1000);    // [Hz]
-        func.SetParLimits(1, 0, 2*Nmdel);
-
-        func.SetParName(0, "lambda");
-        func.SetParName(1, "Nmdel");
-
-        h.Fit("Poisson", "RN");
-
-        //func.SetRange(0, 0.1); // Range of Drawing
-        //func.DrawCopy("same");
-
-        //cout << func.GetParameter(0) << " " << func.GetParameter(1) << endl;
-
-        Double_t lambda = func.GetParameter(0);
-        //Double_t a      = func.GetParameter(1);
-
-        //cout << "t_eff = " << h->Integral()/lambda << "  T(last)=" << time.GetTimeLo()*0.0001 << endl;
-
-        fHist.SetBinContent(i, h.Integral()/lambda);
+        //        TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
+        TH1D &h = *hist->ProjectionX("Calc-theta", i, i, "E");
+
+        // Nmdel = Nm * binwidth,  with Nm = number of observed events
+        const Double_t Nmdel = h.Integral("width");
+        const Double_t Nm    = h.Integral();
+        //        Double_t mean  = h->GetMean();
+
+        //...................................................
+        // determine range (yq[0], yq[1]) of time differences 
+        // where fit should be performed;
+        // require a fraction >=xq[0] of all entries to lie below yq[0]
+        //     and a fraction <=xq[1] of all entries to lie below yq[1];  
+        // within the range (yq[0], yq[1]) there must be no empty bin;
+        // choose pedestrian approach as long as GetQuantiles is not available
+
+        Double_t xq[2] = { 0.15, 0.95 };
+        Double_t yq[2];
+
+        // GetQuantiles doesn't seem to be available in root 3.01/06
+        // h->GetQuantiles(2,yq,xq);
+
+        const Double_t sumtot = h.Integral();
+        const Int_t    jbins  = h.GetNbinsX();
+
+        if (sumtot > 0.0)
+        {
+            Double_t sum1 = 0.0;
+            yq[0] = h.GetBinLowEdge(jbins+1);
+            for (int j=1; j<=jbins; j++)
+            {
+                if (sum1 >= xq[0]*sumtot)
+                {
+                    yq[0] = h.GetBinLowEdge(j);
+                    break;
+                }
+                sum1 += h.GetBinContent(j);
+            }
+
+            Double_t sum2 = 0.0;
+            yq[1] = h.GetBinLowEdge(jbins+1);
+            for (int j=1; j<=jbins; j++)
+            {
+                Double_t content = h.GetBinContent(j);
+                sum2 += content;
+                if (sum2 >= xq[1]*sumtot || content == 0.0)
+                {
+                    yq[1] = h.GetBinLowEdge(j);
+                    break;
+                }
+            }
+
+            //...................................................
+
+            // parameter 0 = lambda
+            // parameter 1 = N0*del        with N0 = ideal number of events
+            //                             and del = bin width of time difference
+            TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", yq[0], yq[1]);
+
+            func.SetParameter(0, 100); // [Hz]
+            func.SetParameter(1, Nmdel);
+
+            func.SetParLimits(0, 0, 1000);    // [Hz]
+            func.SetParLimits(1, 0, 10*Nmdel);
+
+            func.SetParName(0, "lambda");
+            func.SetParName(1, "Nmdel");
+
+            // options : 0  do not plot the function
+            //           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
+            h.Fit("Poisson", "0IRQ");
+
+            // gPad->SetLogy();
+            // gStyle->SetOptFit(1011);
+            // h->GetXaxis()->SetTitle("time difference [s]");
+            // h->GetYaxis()->SetTitle("Counts");
+            // h->DrawCopy();
+
+            // func.SetRange(yq[0], yq[1]); // Range of Drawing
+            // func.DrawCopy("same");
+
+            const Double_t lambda = func.GetParameter(0);
+            const Double_t N0del  = func.GetParameter(1);
+            const Double_t chi2   = func.GetChisquare();
+            const Int_t    NDF    = func.GetNDF();
+
+            // was fit successful ?
+            if (NDF>0 && chi2<2.5*NDF)
+            {
+                // the effective on time is Nm/lambda
+                fHEffOn.SetBinContent(i, Nm/lambda);
+
+                // plot chi2/NDF of fit
+                fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0);
+
+                // lambda of fit
+                fHLambda.SetBinContent(i, lambda);
+
+                // N0del of fit
+                fHN0del.SetBinContent(i, N0del);
+
+                delete &h;
+                continue;
+            }
+        }
+
+        fHEffOn.SetBinContent (i, 1.e-20);
+        fHChi2.SetBinContent  (i, 1.e-20);
+        fHLambda.SetBinContent(i, 1.e-20);
+        fHN0del.SetBinContent (i, 1.e-20);
 
         delete &h;
@@ -119 +229 @@
 }
 
+// -------------------------------------------------------------------------
+//
+// Set the binnings and prepare the filling of the histograms
+//
 Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist)
 {
@@ -128 +242 @@
    }
 
-   SetBinning(&fHist, bins);
+   SetBinning(&fHEffOn,  bins);
+   SetBinning(&fHChi2,   bins);
+   SetBinning(&fHLambda, bins);
+   SetBinning(&fHN0del,  bins);
+
+   fHEffOn.Sumw2();
+   fHChi2.Sumw2();
+   fHLambda.Sumw2();
+   fHN0del.Sumw2();
 
    return kTRUE;
 }
 
+// -------------------------------------------------------------------------
+//
+// Dummy Fill
+// without it get the error message :
+// Error: MHEffOnTimeTime() no default constructor FILE:macros/wowflux.C LINE:359
+// *** Interpreter error recovered ***
+//
+Bool_t MHEffOnTimeTheta::Fill(const MParContainer *par)
+{
+  return kTRUE; 
+}     
+// -------------------------------------------------------------------------
+//
+// Draw a copy of the histogram
+//
 TObject *MHEffOnTimeTheta::DrawClone(Option_t *opt) const
 {
-    TCanvas *c = MakeDefCanvas("EffOnTimeTheta", "t_eff vs. Theta");
+    TCanvas &c = *MakeDefCanvas("EffOnTimeTheta", "Results from on time fit vs. Theta");
+
+    c.Divide(2, 2);
 
     gROOT->SetSelectedPad(NULL);
 
-    ((TH2&)fHist).DrawCopy(opt);
-
-    c->Modified();
-    c->Update();
-
-    return c;
-}
-
+    c.cd(1);
+    ((TH2*)&fHEffOn)->DrawCopy(opt);
+
+    c.cd(2);
+    ((TH2*)&fHChi2)->DrawCopy(opt);
+
+    c.cd(3);
+    ((TH2*)&fHLambda)->DrawCopy(opt);
+
+    c.cd(4);
+    ((TH2*)&fHN0del)->DrawCopy(opt);
+
+    c.Modified();
+    c.Update();
+
+    return &c;
+}
+
+// -------------------------------------------------------------------------
+//
+// Draw the histogram
+//
 void MHEffOnTimeTheta::Draw(Option_t *opt)
 {
     if (!gPad)
-        MakeDefCanvas("EffOnTimeTheta", "t_eff vs. Theta");
-
-    fHist.Draw(opt);
+        MakeDefCanvas("EffOnTimeTheta", "Results from on time fit vs. Theta");
+
+    gPad->Divide(2,2);
+
+    gPad->cd(1);
+    fHEffOn.Draw(opt);
+
+    gPad->cd(2);
+    fHChi2.Draw(opt);
+
+    gPad->cd(3);
+    fHLambda.Draw(opt);
+
+    gPad->cd(4);
+    fHN0del.Draw(opt);
 
     gPad->Modified();
@@ -159 +324 @@
 
 
-Bool_t MHEffOnTimeTheta::Fill(const MParContainer *par)
-{
-    return kTRUE;
-}
-
+
+
+
+
+
+
+
+
+
+
+