Changeset 6924 for trunk/MagicSoft/Mars/mhflux

Timestamp:

04/11/05 16:15:08 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

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

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -78 +78 @@
     fHResolution.SetDirectory(NULL);
     fHResolution.SetName("EnergyRes");
-    fHResolution.SetTitle("Histogram in \\frac{\\Delta E}{E_{mc}} vs E_{est} and E_{mc}");
+    fHResolution.SetTitle("Histogram in \\frac{\\Delta lg(E)}{lg(E_{mc})} vs E_{est} and E_{mc}");
     fHResolution.SetXTitle("E_{est} [GeV]");
     fHResolution.SetYTitle("E_{mc} [GeV]");
-    fHResolution.SetZTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
+    fHResolution.SetZTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
 
     fHImpact.SetDirectory(NULL);
     fHImpact.SetName("Impact");
-    fHImpact.SetTitle("\\frac{\\Delta E}{E} vs Impact parameter");
+    fHImpact.SetTitle("\\frac{\\Delta lg(E)}{lg(E)} vs Impact parameter");
     fHImpact.SetXTitle("Impact parameter [m]");
-    fHImpact.SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
+    fHImpact.SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
 
     fHEnergy.Sumw2();
@@ -97 +97 @@
     binst.SetEdgesCos(50, 0, 60);
     binsi.SetEdges(10, 0, 400);
-    binsr.SetEdges(10, 0, 1);
+    binsr.SetEdges(50, -1.3, 1.3);
 
     SetBinning(&fHEnergy,     &binse, &binse, &binst);
@@ -164 +164 @@
     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   = (eest-etru)/etru;
+    const Double_t res   = (log10(eest)-log10(etru));///log10(etru);
+    const Double_t resE   = (eest-etru)/etru;
 
     fHEnergy.Fill(eest, etru, theta, w);
-    fHResolution.Fill(eest, etru, TMath::Abs(res), w);
-    fHImpact.Fill(imp, TMath::Abs(res), w);
-
-    fChisq += TMath::Abs(res);//*res;
+    fHResolution.Fill(eest, etru, resE, w);
+    fHImpact.Fill(imp, resE, w);
+
+    fChisq += res*res;
     fBias  += res;
 
@@ -181 +182 @@
     fBias  /= GetNumExecutions();
 
-    Double_t res = fChisq; //TMath::Sqrt(fChisq - fBias*fBias);
-
-    fResult->SetVal(TMath::IsNaN(res)?0:res);/// GetNumExecutions());
-
-    *fLog << all << "Mean Energy Resoltuion: " << Form("%.1f%%", fResult->GetVal()*100) << endl;
-    *fLog << all << "Energy Bias at:         " << Form("%.1f%%", fBias*100) << endl;
+    Double_t sigma = TMath::Sqrt(fChisq);//+TMath::Abs(fBias);
+    fResult->SetVal(sigma);
+
+    *fLog << all << "Mean log10(Energy) Resoltion: " << Form("%.1f%%", TMath::Sqrt(fChisq-fBias*fBias)*100) << endl;
+    *fLog << all << "Mean log10(Energy) Bias:      " << Form("%.1f%%", fBias*100) << endl;
 
     return kTRUE;
@@ -239 +239 @@
             }
 
-            //pad->GetPad(1)->GetPad(2)->cd(2);
+            pad->GetPad(1)->GetPad(2)->cd(2);
+            /*=*/fHResolution.ProjectionZ("Resolution");
             ///*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e");
         }
@@ -294 +295 @@
     gPad->SetBorderMode(0);
     gPad->SetLogx();
+    gPad->SetGridx();
+    gPad->SetGridy();
 
     //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
@@ -301 +304 @@
     h2->SetBit(kCanDelete);
     h2->SetFillColor(kBlue);
+    h2->SetLineColor(kRed);
 
     TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, 9999, "s");
@@ -310 +314 @@
     h1->SetStats(kFALSE);
 
- 
-    //h1->Draw("E3");
-    h2->Draw();
+    h2->Draw("");
+    h1->Draw("E3same");
     h1->Draw("Chistsame");
 
@@ -346 +349 @@
     gPad->SetBorderMode(0);
 
+    gPad->SetGridx();
+    gPad->SetGridy();
     gPad->SetLogx();
     h = (TH1D*)fHEnergy.Project3D("ey");
@@ -376 +381 @@
     pad3->Divide(2, 1, 1e-10, 1e-10);
     pad3->cd(1);
-    gPad->SetBorderMode(0);/*
+    gPad->SetBorderMode(0);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    h = fHEnergy.Project3D("ez");
+    h->SetTitle("Zenith Angle Distribution");
+    h->SetBit(TH1::kNoStats);
+    h->SetDirectory(NULL);
+    h->SetXTitle("\\Theta [\\circ]");
+    h->SetBit(kCanDelete);
+    h->Draw();
+
+    pad3->cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    /*
     h = fHImpact.ProjectionX("Impact", -1, 9999, "e");
     h->SetBit(TH1::kNoStats);
@@ -384 +404 @@
     h->SetBit(kCanDelete);
     h->Draw();*/
-    h = fHEnergy.Project3D("ez");
-    h->SetTitle("Distribution of Theta");
-    h->SetBit(TH1::kNoStats);
+    // ----------------------------------------
+    h = fHResolution.ProjectionZ("Resolution");
+    //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetYTitle("Counts");
+    h->SetTitle("Distribution of \\Delta lg(E) / lg(E_{mc})");
     h->SetDirectory(NULL);
-    h->SetXTitle("\\Theta [\\circ]");
     h->SetBit(kCanDelete);
-    h->Draw();
-
-    pad3->cd(2);
-    gPad->SetBorderMode(0);
+    h->GetXaxis()->SetRangeUser(-1.3, 1.3);
+    h->Draw("");
+    // ----------------------------------------
 
     pad->cd(2);
@@ -413 +433 @@
     h = MakePlot(fHResolution, "zy");
     h->SetXTitle("E_{mc} [GeV]");
-    h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
-    h->SetTitle("Energy resolution \\Delta E / E vs Monte Carlo energy E_{mc}");
-    h->SetMinimum(0);
-    h->SetMaximum(1);
+    h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetTitle("Energy resolution \\Delta lg(E) / lg(E) vs Monte Carlo energy E_{mc}");
+    h->SetMinimum(-1.3);
+    h->SetMaximum(1.3);
 
     pad2->cd(3);
     h = MakePlot(fHResolution, "zx");
     h->SetXTitle("E_{est} [GeV]");
-    h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
-    h->SetTitle("Energy Resolution \\Delta E / E vs estimated Energy E_{est}");
-    h->SetMinimum(0);
-    h->SetMaximum(1);
+    h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetTitle("Energy Resolution \\Delta lg(E) / lg(E) vs estimated Energy E_{est}");
+    h->SetMinimum(-1.3);
+    h->SetMaximum(1.3);
 }