Changeset 7409 for trunk/MagicSoft/Mars/mhflux

Timestamp:

11/18/05 17:15:30 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

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

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -35 +35 @@
 #include "MHEnergyEst.h"
 
+#include <TF1.h>
 #include <TLine.h>
 #include <TCanvas.h>
@@ -79 +80 @@
     fHResolution.SetDirectory(NULL);
     fHResolution.SetName("EnergyRes");
-    fHResolution.SetTitle("Histogram in \\Delta lg(E) vs E_{est} and E_{mc}");
+    fHResolution.SetTitle("Histogram in \\Delta E/E vs E_{est} and E_{mc}");
     fHResolution.SetXTitle("E_{est} [GeV]");
     fHResolution.SetYTitle("E_{mc} [GeV]");
-    fHResolution.SetZTitle("lg(E_{est}) - lg(E_{mc})");
+    fHResolution.SetZTitle("E_{est}/E_{mc}-1");
 
     fHImpact.SetDirectory(NULL);
     fHImpact.SetName("Impact");
-    fHImpact.SetTitle("\\Delta lg(E) vs Impact parameter");
+    fHImpact.SetTitle("\\Delta E/E vs Impact parameter");
     fHImpact.SetXTitle("Impact parameter [m]");
-    fHImpact.SetYTitle("lg(E_{est}) - lg(E_{mc})");
+    fHImpact.SetYTitle("E_{est}/E_{mc}-1");
 
     fHEnergy.Sumw2();
@@ -95 +96 @@
 
     MBinning binsi, binse, binst, binsr;
-    binse.SetEdgesLog(15, 10, 1000000);
+    binse.SetEdgesLog(25, 10, 1000000);
     binst.SetEdgesASin(51, -0.005, 0.505);
     binsi.SetEdges(10, 0, 400);
-    binsr.SetEdges(50, -1.3, 1.3);
+    binsr.SetEdges(75, -1.75, 1.75);
 
     SetBinning(&fHEnergy,     &binse, &binse, &binst);
@@ -149 +150 @@
     fChisq = 0;
     fBias  = 0;
+
     fHEnergy.Reset();
     fHImpact.Reset();
@@ -155 +157 @@
     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);
-*/
+
 // --------------------------------------------------------------------------
 //
@@ -172 +168 @@
     const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
     const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
-    const Double_t res   = log10(eest)-log10(etru);///log10(etru);
     const Double_t resE  = (eest-etru)/etru;
 
@@ -179 +174 @@
     fHImpact.Fill(imp, resE, w);
 
-    //if (etru>125 && etru<750)
-    //    return kCONTINUE;
-
-    // lg(N)  = 4.3*lg(E)-6
-    // lg(N0) = 4.3*2.2-6
-
-    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;
+    // For the fit we use a different quantity
+    //const Double_t res = TMath::Log10(eest/etru);
+    const Double_t res = eest-etru;
+
+    fChisq += res*res;
+    fBias  += 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);
-*/
+    fResult->SetVal(fChisq);
+
     Print();
 
@@ -247 +198 @@
 void MHEnergyEst::Print(Option_t *o) const
 {
-    const Double_t res  =   TMath::Sqrt(fChisq-fBias*fBias);
-    const Double_t resp =   TMath::Power(10,  res)-1;
-    const Double_t resm = 1-TMath::Power(10, -res);
-
+    //    const Double_t resp =   TMath::Power(10,  res)-1;
+    //    const Double_t resm = 1-TMath::Power(10, -res);
+
+    const Double_t res = TMath::Sqrt(fChisq-fBias*fBias);
     if (TMath::IsNaN(res))
     {
@@ -257 +208 @@
     }
 
-    *fLog << all;
-    *fLog << "Mean log10(Energy) Resoltion: +/- " << Form("%4.2f", res) << endl;
-    *fLog << "Mean log10(Energy) Bias:         "  << Form("%+4.2f", fBias) << endl;
-    *fLog << "Asymmetric Energy  Resoltion:   +"  << Form("%4.1f%%", resp*100) << " -" << Form("%4.1f%%", resm*100) << endl;
+    TH1D *h = (TH1D*)fHResolution.ProjectionZ("Resolution");
+    h->Fit("gaus", "Q0");
+
+    TF1 *f = h->GetFunction("gaus");
+
+    *fLog << all << "F=" << fChisq << endl;
+    *fLog << "Results from Histogram:" << endl;
+    *fLog << " Mean  of Delta E/E: " << Form("%+4.2f%%", 100*h->GetMean()) << endl;
+    *fLog << " RMS   of Delta E/E: " << Form("%4.2f%%",  100*h->GetRMS()) << endl;
+    *fLog << "Results from Histogram-Fit:" << endl;
+    *fLog << " Bias  of Delta E/E: " << Form("%+4.2f%%", 100*f->GetParameter(1)) << endl;
+    *fLog << " Sigma of Delta E/E: " << Form("%4.2f%%",  100*f->GetParameter(2)) << endl;
+    *fLog << " Res   of Delta E/E: " << Form("%4.2f%%",  100*TMath::Hypot(f->GetParameter(1), f->GetParameter(2))) << endl;
+
+    delete h;
 }
 
@@ -411 +373 @@
     h1->SetBit(kCanDelete);
     h1->SetLineWidth(2);
-    h1->SetLineColor(kRed);
+    h1->SetLineColor(kBlue);
     h1->SetFillStyle(4000);
     h1->SetStats(kFALSE);
 
     h2->Draw("");
-    h1->Draw("E3 hist C same");
+    h1->Draw("E0 hist C same");
 //    h1->Draw("Chistsame");
 
@@ -509 +471 @@
     //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
     h->SetYTitle("Counts");
-    h->SetTitle("Distribution of \\Delta lg(E)");
+    h->SetTitle("Distribution of \\Delta E/E");
     h->SetDirectory(NULL);
     h->SetBit(kCanDelete);
@@ -535 +497 @@
     h = MakePlot(fHResolution, "zy");
     h->SetXTitle("E_{mc} [GeV]");
-    h->SetYTitle("lg(E_{est}) - lg(E_{mc})");
-    h->SetTitle("Energy resolution \\Delta lg(E) vs Monte Carlo energy E_{mc}");
+    h->SetYTitle("(E_{est}/E_{mc}-1");
+    h->SetTitle("Energy resolution \\Delta E/E vs Monte Carlo Energy E_{mc}");
     h->SetMinimum(-1.3);
     h->SetMaximum(1.3);
@@ -543 +505 @@
     h = MakePlot(fHResolution, "zx");
     h->SetXTitle("E_{est} [GeV]");
-    h->SetYTitle("lg(E_{est}) - lg(E_{mc})");
-    h->SetTitle("Energy Resolution \\Delta lg(E) vs estimated Energy E_{est}");
+    h->SetYTitle("(E_{est}/E_{mc}-1");
+    h->SetTitle("Energy resolution \\Delta E/E vs estimated Energy E_{est}");
     h->SetMinimum(-1.3);
     h->SetMaximum(1.3);