Changeset 4889 for trunk/MagicSoft/Mars/mhist

Timestamp:

09/08/04 18:49:00 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHEffectiveOnTime.cc (modified) (8 diffs)
• MHEffectiveOnTime.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc

@@ -42 +42 @@
 #include <TCanvas.h>
 #include <TMinuit.h>
+#include <TRandom.h>
 
 #include "MTime.h"
@@ -62 +63 @@
 //
 MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
-    : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60)
+    : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60)
 {
     //
@@ -141 +142 @@
 }
 
-Double_t testval = 0;
-Double_t testerr = 0;
+// FIXME: Just for a preliminary check
+static Double_t testval = 0;
+static Double_t testerr = 0;
 
 // --------------------------------------------------------------------------
@@ -185 +187 @@
 Bool_t MHEffectiveOnTime::Finalize()
 {
-    Fit();
+    FitThetaBins();
+    Calc();
+
     fIsFinalized = kTRUE;
+
     return kTRUE;
 }
 
-void MHEffectiveOnTime::Fit()
-{
+void MHEffectiveOnTime::FitThetaBins()
+{
+    const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
+
     // nbins = number of Theta bins
     const Int_t nbins = fHTimeDiff.GetNbinsY();
 
+    TH1D *h=0;
     for (int i=1; i<=nbins; i++)
     {
         //        TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
-        TH1D *h = fHTimeDiff.ProjectionX("CalcTheta", i, i, "E");
-
-        const Double_t Nm = h->Integral();
-
-        if (Nm<=0)
+        h = fHTimeDiff.ProjectionX(name, i, i, "E");
+
+        Double_t res[7];
+        if (!FitH(h, res))
             continue;
 
-        // 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);
-
-        // Nmdel = Nm * binwidth,  with Nm = number of observed events
-        const Double_t Nmdel = h->Integral("width");
-
-        //
-        // Setup Poisson function for the fit:
-        // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
-        //
-        // parameter 0 = lambda
-        // parameter 1 = N0*del      
-        //
-        TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
-        func.SetParNames("lambda", "N0", "del");
-
-        func.SetParameter(0, 100);       // Hz
-        func.SetParameter(1, Nm);
-        func.FixParameter(2, Nmdel/Nm);
-
-        // 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(&func, "0IRQ");
-
-        const Double_t chi2   = func.GetChisquare();
-        const Int_t    NDF    = func.GetNDF();
-
-        // was fit successful ?
-        if (NDF>0 && chi2<2.5*NDF)
-        {
-            const Double_t lambda = func.GetParameter(0);
-            const Double_t N0     = func.GetParameter(1);
-            const Double_t prob   = func.GetProb();
-
-            /*
-             *fLog << all << "Nm/lambda=" << Nm/lambda << "  chi2/NDF=";
-             *fLog << (NDF ? chi2/NDF : 0.0) << "  lambda=";
-             *fLog << lambda << "  N0=" << N0 << endl;
-             */
-
-            Double_t emat[2][2];
-            gMinuit->mnemat((Double_t*)emat, 2);
-
-            const Double_t dldl   = emat[0][0];
-            //const Double_t dN0dN0 = emat[1][1];
-
-            const Double_t teff   = Nm/lambda;
-            const Double_t dteff  = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
-
-            const Double_t dl     = TMath::Sqrt(dldl);
-
-            //const Double_t kappa  = Nm/N0;
-            //const Double_t Rdead  = 1.0 - kappa;
-            //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
-
-            // the effective on time is Nm/lambda
-            fHEffOn.SetBinContent(i, teff);
-            fHEffOn.SetBinError  (i, dteff);
-
-            // plot chi2-probability of fit
-            fHProb.SetBinContent(i, prob*100);
-
-            // plot chi2/NDF of fit
-            fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0);
-
-            // lambda of fit
-            fHLambda.SetBinContent(i, lambda);
-            fHLambda.SetBinError  (i,     dl);
-
-            // N0 of fit
-            fHN0.SetBinContent(i, N0);
-
-            // Rdead (from fit) is the fraction from real time lost by the dead time
-            //fHRdead.SetBinContent(i, Rdead);
-            //fHRdead.SetBinError  (i,dRdead);
-        }
-
+        // the effective on time is Nm/lambda
+        fHEffOn.SetBinContent(i, res[0]);
+        fHEffOn.SetBinError  (i, res[1]);
+
+        // plot chi2-probability of fit
+        fHProb.SetBinContent(i, res[2]);
+
+        // plot chi2/NDF of fit
+        fHChi2.SetBinContent(i, res[3]);
+
+        // lambda of fit
+        fHLambda.SetBinContent(i, res[4]);
+        fHLambda.SetBinError  (i, res[5]);
+
+        // N0 of fit
+        fHN0.SetBinContent(i, res[6]);
+
+        // Rdead (from fit) is the fraction from real time lost by the dead time
+        //fHRdead.SetBinContent(i, Rdead);
+        //fHRdead.SetBinError  (i,dRdead);
+    }
+
+    // Histogram is reused via gROOT->FindObject()
+    // Need to be deleted only once
+    if (h)
         delete h;
-    }
-}
-
+}
+
+Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const
+{
+    const Double_t Nm = h->Integral();
+
+    // FIXME: Do fit only if contents of bin has changed
+    if (Nm<=0)
+        return kFALSE;
+
+    // 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.05, 0.95 };
+    Double_t yq[2];
+    h->GetQuantiles(2, yq, xq);
+
+    // Nmdel = Nm * binwidth,  with Nm = number of observed events
+    const Double_t Nmdel = h->Integral("width");
+
+    //
+    // Setup Poisson function for the fit:
+    // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
+    //
+    // parameter 0 = lambda
+    // parameter 1 = N0*del
+    //
+    TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
+    func.SetParNames("lambda", "N0", "del");
+
+    func.SetParameter(0, 100);       // Hz
+    func.SetParameter(1, Nm);
+    func.FixParameter(2, Nmdel/Nm);
+
+    // 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(&func, "0IRQ");
+
+    const Double_t chi2   = func.GetChisquare();
+    const Int_t    NDF    = func.GetNDF();
+
+    // was fit successful ?
+    if (NDF<=0 || chi2>=2.5*NDF)
+        return kFALSE;
+
+    const Double_t lambda = func.GetParameter(0);
+    const Double_t N0     = func.GetParameter(1);
+    const Double_t prob   = func.GetProb();
+
+    /*
+     *fLog << all << "Nm/lambda=" << Nm/lambda << "  chi2/NDF=";
+     *fLog << (NDF ? chi2/NDF : 0.0) << "  lambda=";
+     *fLog << lambda << "  N0=" << N0 << endl;
+     */
+
+    Double_t emat[2][2];
+    gMinuit->mnemat((Double_t*)emat, 2);
+
+    const Double_t dldl   = emat[0][0];
+    //const Double_t dN0dN0 = emat[1][1];
+
+    const Double_t teff   = Nm/lambda;
+    const Double_t dteff  = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
+
+    const Double_t dl     = TMath::Sqrt(dldl);
+
+    //const Double_t kappa  = Nm/N0;
+    //const Double_t Rdead  = 1.0 - kappa;
+    //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
+
+    // the effective on time is Nm/lambda
+    res[0] = teff;
+    res[1] = dteff;
+
+    // plot chi2-probability of fit
+    res[2] = prob*100;
+
+    // plot chi2/NDF of fit
+    res[3] = NDF ? chi2/NDF : 0.0;
+
+    // lambda of fit
+    res[4] = lambda;
+    res[5] = dl;
+
+    // N0 of fit
+    res[6] = N0;
+
+    // Rdead (from fit) is the fraction from real time lost by the dead time
+    //fHRdead.SetBinContent(i, Rdead);
+    //fHRdead.SetBinError  (i,dRdead);
+
+    return kTRUE;
+}
 
 void MHEffectiveOnTime::Paint(Option_t *opt)
@@ -320 +359 @@
 
         if (!fIsFinalized)
-            Fit();
+            FitThetaBins();
     }
     if (o==(TString)"paint")
@@ -386 +425 @@
 void MHEffectiveOnTime::Calc()
 {
-    const Double_t val = fHEffOn.Integral();
-
-    Double_t error = 0;
-    for (int i=0; i<fHEffOn.GetXaxis()->GetNbins(); i++)
-        error += fHEffOn.GetBinError(i);
+    TH1D *h = fHTimeDiff.ProjectionX("", -1, 99999, "E");
+    h->SetDirectory(0);
+
+    Double_t res[7];
+    Bool_t rc = FitH(h, res);
+    delete h;
+
+    if (!rc)
+        return;
+
+    const Double_t val   = res[0];
+    const Double_t error = res[1];
 
     fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0);
@@ -403 +449 @@
     MTime now(*fTime);
     now.AddMilliSeconds(fInterval*1000);
+
     *fLog <<all << now << " - ";// << fLastTime-fTime;
     *fLog << Form("T_{eff} = %.1fs \\pm %.1fs",
@@ -435 +482 @@
 
         *fTime = *time;
-        // fTime->MinusNull()
+
+        // Make this just a ns before the first event
+        fTime->Minus1ns();
     }
 
     if (*fTime+fInterval<*time)
     {
-        Fit();
+        FitThetaBins();
         Calc();
     }

trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h

@@ -24 +24 @@
 {
 private:
-    MPointingPos *fPointPos;   //!
-    MTime         fLastTime;   //!
+    MPointingPos   *fPointPos;   //!
+    MTime           fLastTime;   //!
 
-    Double_t      fEffOnTime0; //!
-    Double_t      fEffOnErr0;  //!
+    Double_t        fEffOnTime0; //!
+    Double_t        fEffOnErr0;  //!
 
-    MTime          *fTime;
-    MParameterDerr *fParam;
+    MTime          *fTime;       //!
+    MParameterDerr *fParam;      //!
 
     TH2D fHTimeDiff;
@@ -44 +44 @@
     Int_t fInterval;
 
-    void Fit();
+    Bool_t FitH(TH1D *h, Double_t *res) const;
+    void FitThetaBins();
     void Calc();