Changeset 7602 for trunk/MagicSoft/Mars/mhflux

Timestamp:

03/14/06 16:07:39 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc

@@ -474 +474 @@
 
     // FIXME: Do fit only if contents of bin has changed
-    if (Nm<=0)
+    if (Nm<=0 || h->GetBinContent(1)<=0)
         return kFALSE;
 
@@ -483 +483 @@
     // 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.6, 0.99 };
+    Double_t xq[2] = { 0.6, 0.95 }; // previously 0.99
     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)");
-    //func.SetParNames("lambda", "N0", "del");
+    //
+    // Setup exponential function for the fit:
+    //
+    // parameter 0 = rate [Hz]
+    // parameter 1 = normalization
+    //
+    TF1 func("Exp", " exp([1]-[0]*x)");
 
     func.SetParameter(0, 200);       // Hz
-    func.SetParameter(1, Nm);
-    func.FixParameter(2, Nmdel/Nm);
+    func.SetParameter(1, log(h->GetBinContent(1)));       // Hz
 
     // options : N  do not store the function, do not draw
@@ -508 +502 @@
     //           R  use the range specified in the function range
     //           Q  quiet mode
-    h->Fit(&func, "NIQE", "", yq[0], yq[1]);
+    //           L  Use log-likelihood (better for low statistics)
+    h->Fit(&func, "NIQEL", "", h->GetBinLowEdge(3)/*yq[0]*/, yq[1]);
 
     const Double_t chi2 = func.GetChisquare();
+    const Double_t prob = func.GetProb();
     const Int_t    NDF  = func.GetNDF();
 
     // was fit successful ?
-    const Bool_t ok = NDF>0 && chi2<3*NDF;
+    const Bool_t ok = prob>0.001; //NDF>0 && chi2<3*NDF;
 
     if (paint)
@@ -523 +519 @@
     }
 
+    // The effective on time is the "real rate" (slope of the exponential)
+    // divided by the total number of events (histogram integral including 
+    // under- and overflows)
     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);
+    const Double_t dldl   = func.GetParError(0)*func.GetParError(0);
+    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);
 
     // the effective on time is Nm/lambda
@@ -563 +544 @@
     // Chi2
     res[6] = chi2;
-
-    // Rdead (from fit) is the fraction from real time lost by the dead time
-    //fHRdead.SetBinContent(i, Rdead);
-    //fHRdead.SetBinError  (i,dRdead);
 
     return ok;
@@ -583 +560 @@
     fHThetaNDF.Reset();
 
+    // Use a random name to make sure the object is unique
     const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
 
@@ -594 +572 @@
         h = fH2DeltaT.ProjectionX(name, i, i, "E");
 
-        Double_t res[7];
-        if (!FitH(h, res))
+        Double_t res[7] = {0, 0, 0, 0, 0, 0, 0};
+        //if (!FitH(h, res))
+        //    continue;
+        FitH(h, res);
+
+        if (res[0]==0)
             continue;