Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 6949) @@ -347,12 +347,12 @@ { *fLog << "Result:" << endl; - *fLog << " - Significance " << fSignificance << endl; - *fLog << " - Excess Events " << fEventsExcess << endl; - *fLog << " - Signal Events " << fEventsSignal << endl; - *fLog << " - Background Events " << fEventsBackground << endl; - *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl; - *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl; - *fLog << " - Signal integrated " << fIntegralMax << "°" << endl; - *fLog << " - Scale Factor (Off) " << fScaleFactor << endl; + *fLog << " - Significance (Li/Ma) " << fSignificance << endl; + *fLog << " - Excess Events " << fEventsExcess << endl; + *fLog << " - Signal Events " << fEventsSignal << endl; + *fLog << " - Background Events " << fEventsBackground << endl; + *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl; + *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl; + *fLog << " - Signal integrated up to " << fIntegralMax << "°" << endl; + *fLog << " - Scale Factor (Off) " << fScaleFactor << endl; } } Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 6949) @@ -86,5 +86,5 @@ void Clear(Option_t *o=""); - void Print(Option_t *o="") const; + void Print(Option_t *o="") const; //*MENU* void Copy(TObject &o) const; /* Index: trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 6949) @@ -660,5 +660,5 @@ { TVirtualPad *p=gPad->GetPad(x+1)->GetPad(y+1); - if (!(st = (TPaveStats*)p->GetPrimitive("stats"))) + if (!(st = dynamic_cast(p->GetPrimitive("stats")))) continue; Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 6938) +++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 6949) @@ -36,6 +36,8 @@ #include -#include +#include +#include #include +#include #include "MString.h" @@ -78,14 +80,14 @@ fHResolution.SetDirectory(NULL); fHResolution.SetName("EnergyRes"); - fHResolution.SetTitle("Histogram in \\frac{\\Delta lg(E)}{lg(E_{mc})} vs E_{est} and E_{mc}"); + fHResolution.SetTitle("Histogram in \\Delta lg(E) vs E_{est} and E_{mc}"); fHResolution.SetXTitle("E_{est} [GeV]"); fHResolution.SetYTitle("E_{mc} [GeV]"); - fHResolution.SetZTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); + fHResolution.SetZTitle("lg(E_{est}) - lg(E_{mc})"); fHImpact.SetDirectory(NULL); fHImpact.SetName("Impact"); - fHImpact.SetTitle("\\frac{\\Delta lg(E)}{lg(E)} vs Impact parameter"); + fHImpact.SetTitle("\\Delta lg(E) vs Impact parameter"); fHImpact.SetXTitle("Impact parameter [m]"); - fHImpact.SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); + fHImpact.SetYTitle("lg(E_{est}) - lg(E_{mc})"); fHEnergy.Sumw2(); @@ -147,4 +149,5 @@ fChisq = 0; + fBias = 0; fHEnergy.Reset(); fHImpact.Reset(); @@ -164,6 +167,6 @@ 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; + const Double_t res = log10(eest)-log10(etru);///log10(etru); + const Double_t resE = (eest-etru)/etru; fHEnergy.Fill(eest, etru, theta, w); @@ -171,6 +174,14 @@ fHImpact.Fill(imp, resE, w); - fChisq += res*res; - fBias += res; + //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; //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; return kTRUE; @@ -192,7 +203,18 @@ 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); + + if (TMath::IsNaN(res)) + { + *fLog << all << "MHEnergyEst::Print: ERROR - Resolution is NaN (not a number)." << endl; + return; + } + *fLog << all; - *fLog << "Mean log10(Energy) Resoltion: " << Form("%5.2f%%", TMath::Sqrt(fChisq-fBias*fBias)) << endl; - *fLog << "Mean log10(Energy) Bias: " << Form("%5.2f%%", fBias) << endl; + *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; } @@ -246,6 +268,16 @@ pad->GetPad(1)->GetPad(2)->cd(2); - /*=*/fHResolution.ProjectionZ("Resolution"); - ///*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e"); + hx = (TH1D*)fHResolution.ProjectionZ("Resolution"); + TPaveStats *stats = dynamic_cast(hx->FindObject("stats")); + if (stats) + { + stats->SetBit(BIT(17)); // TStyle.cxx: kTakeStyle=BIT(17) + stats->SetX1NDC(0.63); + stats->SetY1NDC(0.68); + } + + hx->Fit("gaus", "Q"); + gPad=NULL; + gStyle->SetOptFit(101); } } @@ -325,5 +357,5 @@ // h1->Draw("Chistsame"); - return h1; + return h2; } @@ -415,9 +447,10 @@ //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->SetTitle("Distribution of \\Delta lg(E)"); h->SetDirectory(NULL); h->SetBit(kCanDelete); h->GetXaxis()->SetRangeUser(-1.3, 1.3); h->Draw(""); + //h->Fit("gaus"); // ---------------------------------------- @@ -440,6 +473,6 @@ h = MakePlot(fHResolution, "zy"); h->SetXTitle("E_{mc} [GeV]"); - 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->SetYTitle("lg(E_{est}) - lg(E_{mc})"); + h->SetTitle("Energy resolution \\Delta lg(E) vs Monte Carlo energy E_{mc}"); h->SetMinimum(-1.3); h->SetMaximum(1.3); @@ -448,6 +481,6 @@ h = MakePlot(fHResolution, "zx"); h->SetXTitle("E_{est} [GeV]"); - 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->SetYTitle("lg(E_{est}) - lg(E_{mc})"); + h->SetTitle("Energy Resolution \\Delta lg(E) vs estimated Energy E_{est}"); h->SetMinimum(-1.3); h->SetMaximum(1.3);