Changeset 7817 for trunk/MagicSoft/Mars

Timestamp:

07/29/06 11:20:14 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhflux/MHEffectiveOnTime.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2006/07/29 Thomas Bretz
+
+   * mhflux/MHEffectiveOnTime.cc:
+     - changed the fit such that initial values are calculated
+       automatically now instead of using build in values. This
+       makes the fit more indepedant of the underlaying rates.
+       With a test of ~350 sequences in the test database it could
+       be shown that the new fit gives the same result +/-1sek.
+       The highest deviation was +5s the lowest -10s.
+     - the number of the first bin used in the fit became a variable
+     - A limit of 15kHz was set for the rate
+     - sanity checkes for lambda==0 added (possible division by zero)
+     - increased class version number by one
+     - made functions derived from MH private
+     - to fit the resulting "gammas" use fFirstBin=1 and fNumEvents=120
+
+
+
  2006/07/28 Daniela Dorner
 

trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc

@@ -76 +76 @@
 //
 //
+//  Class version 2:
+//  ----------------
+//   +  UInt_t fFirstBin;
+//   +  UInt_t fNumEvents;
+//   -  Int_t fNumEvents;
+//
+//
 // ==========================================================================
 // Dear Colleagues,
@@ -286 +293 @@
 MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
     : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), 
-    fNumEvents(200*60)/*, fNameProjDeltaT(Form("DeltaT_%p", this)),
-    fNameProjTheta(Form("Theta_%p", this))*/
+    fNumEvents(200*60), fFirstBin(3)
+    //fNumEvents(2*60), fFirstBin(1)
 {
     //
@@ -488 +495 @@
 
     //
+    // Determine a good starting value for the slope
+    //
+    const TAxis  &axe = *h->GetXaxis();
+    const UInt_t ibin = axe.FindFixBin(yq[1]);
+    const Double_t x1 = axe.GetBinCenter(ibin<=fFirstBin?fFirstBin+1:ibin);
+    const Double_t x0 = axe.GetBinCenter(fFirstBin);
+    const Double_t y1 = h->GetBinContent(ibin)>1 ? TMath::Log(h->GetBinContent(ibin)) : 0;
+    const Double_t y0 = TMath::Log(h->GetBinContent(fFirstBin));
+
+    // Estimated slope
+    const Float_t m = -(y1-y0)/(x1-x0);
+
+    //
     // Setup exponential function for the fit:
     //
@@ -495 +515 @@
     TF1 func("Exp", " exp([1]-[0]*x)");
 
-    func.SetParameter(0, 200);       // Hz
+    func.SetParameter(0, m);       // Hz
     func.SetParameter(1, log(h->GetBinContent(1)));       // Hz
+
+    // We set a limit to make sure that almost empty histograms which
+    // are fitted dont't produce hang ups or crashes
+    func.SetParLimits(0, 0, 15000); // Hz
 
     // options : N  do not store the function, do not draw
@@ -503 +527 @@
     //           Q  quiet mode
     //           L  Use log-likelihood (better for low statistics)
-    h->Fit(&func, "NIQEL", "", h->GetBinLowEdge(3)/*yq[0]*/, yq[1]);
+    h->Fit(&func, "NIQEL", "", h->GetBinLowEdge(fFirstBin)/*yq[0]*/, yq[1]);
 
     const Double_t chi2 = func.GetChisquare();
@@ -524 +548 @@
     const Double_t lambda = func.GetParameter(0);
     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 teff   = lambda==0 ? 0 : Nm / lambda;
+    const Double_t dteff  = lambda==0 ? 0 : teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
     const Double_t dl     = TMath::Sqrt(dldl);
 
@@ -715 +739 @@
     // histogram - if it fails try again when 1.6% more events available
     //
-    const Int_t n = (Int_t)fH1DeltaT.GetEntries();
+    const UInt_t n = (UInt_t)fH1DeltaT.GetEntries();
     if (n>=fNumEvents && n%(fNumEvents/60)==0)
         FitTimeBin();
@@ -993 +1017 @@
     DrawRightAxis("\\lambda [s^{-1}]");
 }
+
+// --------------------------------------------------------------------------
+//
+// The following resources are available:
+//
+//    MHEffectiveOnTime.FistBin:   3
+//    MHEffectiveOnTime.NuMEvents: 12000
+//
+Int_t MHEffectiveOnTime::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "FirstBin", print))
+    {
+        rc = kTRUE;
+        SetFirstBin(GetEnvValue(env, prefix, "FirstBin", (Int_t)fFirstBin));
+    }
+    if (IsEnvDefined(env, prefix, "NumEvents", print))
+    {
+        rc = kTRUE;
+        SetNumEvents(GetEnvValue(env, prefix, "NumEvents", (Int_t)fNumEvents));
+    }
+    return rc;
+}