Changeset 6983 for trunk/MagicSoft/Mars/mhflux

Timestamp:

04/28/05 11:10:12 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (modified) (5 diffs)

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -45 +45 @@
 #include "MMcEvt.hxx"
 
-#include "MEnergyEst.h"
 #include "MBinning.h"
 #include "MParList.h"
@@ -156 +155 @@
     return kTRUE;
 }
-
+/*
+#include <TSpline.h>
+Double_t x[] = {33, 75, 153, 343, 745, 1617, 3509, 7175};
+Double_t y[] = {2,  24, 302, 333, 132,   53,   11,    1};
+Int_t n = 8;
+TSpline3 spline("", x, y, n);
+*/
 // --------------------------------------------------------------------------
 //
@@ -180 +185 @@
     // lg(N0) = 4.3*2.2-6
 
-    const Double_t n  = 2; //pow(10, -4.3*(log10(etru)-2.2));
-
-    fChisq += log10(etru)<2.2? res*res*n/2:res*res;
-    fBias  += log10(etru)<2.2? res*sqrt(n/2):res;
+    const Double_t n  = 2.;///spline.Eval(etru); //pow(10, -4.3*(log10(etru)-2.2));
+    if (n<=0)
+        return kCONTINUE;
+
+    fChisq += log10(etru)>0 ? res*res*n/2   : res*res;
+    fBias  += log10(etru)>0 ? res*sqrt(n/2) : res;
 
     return kTRUE;
 }
-
+/*
+void MHEnergyEst::CalcChisq(Double_t &chisq, Double_t &prob) const
+{
+    TH1D *h1 = (TH1D*)fHEnergy.Project3D("dum2_ex");
+    TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum2_ey");
+
+    Int_t df = 0;
+    chisq    = 0;
+    prob     = 0;
+
+    for (Int_t i=1; i<=h1->GetNbinsX(); i++)
+    {
+        if (h1->GetBinContent(i)>0 && h2->GetBinContent(i)>0)
+        {
+            const Double_t bin1 = log10(h1->GetBinContent(i));
+            const Double_t bin2 = log10(h2->GetBinContent(i));
+
+            const Double_t temp = bin1-bin2;
+
+            chisq += temp*temp/(bin1+bin2);
+            df ++;
+        }
+    }
+  
+    prob = TMath::Prob(0.5*chisq, Int_t(0.5*df));
+
+    chisq /= df;
+  
+    delete h1;
+    delete h2;
+}
+*/
 Bool_t MHEnergyEst::Finalize()
 {
+    //Double_t chisq, prob;
+    //CalcChisq(chisq, prob);
+  
     fChisq /= GetNumExecutions();
     fBias  /= GetNumExecutions();
 
+    fResult->SetVal(TMath::Sqrt(fChisq));
+
+/*
     Double_t sigma = TMath::Sqrt(fChisq);//+TMath::Abs(fBias);
     fResult->SetVal(sigma);
-
+*/
     Print();
 
@@ -225 +269 @@
     TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey");
 
-    h1->Copy(hist);
-    hist.Divide(h2);
+    h2->Copy(hist);
+    hist.Divide(h1);
 
     delete h1;
@@ -233 +277 @@
     hist.SetNameTitle("EnergyRatio", "Ratio between estimated and monte carlo energy");
     hist.SetXTitle("E [GeV]");
-    hist.SetYTitle("N_{est}/N_{mc} [1]");
+    hist.SetYTitle("N_{mc}/N_{est} [1]");
     hist.SetDirectory(0);
 }