Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4974) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4975) @@ -48,5 +48,5 @@ // // For each "time-bin" the histogram is fitted and the resulting effective -// on-time is stored in the fHEffOnTime histogram. Each entry in this +// on-time is stored in the fHTimeEffOn histogram. Each entry in this // histogram is the effective observation time between the upper and // lower edges of the bins. @@ -83,4 +83,5 @@ #include +#include #include #include @@ -131,79 +132,101 @@ // effective on time versus theta - fHEffOnTheta.SetName("EffOnTheta"); - fHEffOnTheta.SetTitle("Effective On Time T_{eff}"); - fHEffOnTheta.SetXTitle("\\Theta [\\circ]"); - fHEffOnTheta.SetYTitle("T_{eff} [s]"); - fHEffOnTheta.UseCurrentStyle(); - fHEffOnTheta.SetDirectory(NULL); - fHEffOnTheta.GetYaxis()->SetTitleOffset(1.2); + fHThetaEffOn.SetName("EffOnTheta"); + fHThetaEffOn.SetTitle("Effective On Time T_{eff}"); + fHThetaEffOn.SetXTitle("\\Theta [\\circ]"); + fHThetaEffOn.SetYTitle("T_{eff} [s]"); + fHThetaEffOn.UseCurrentStyle(); + fHThetaEffOn.SetDirectory(NULL); + fHThetaEffOn.GetYaxis()->SetTitleOffset(1.2); + fHThetaEffOn.GetYaxis()->SetTitleColor(kBlue); + fHThetaEffOn.SetLineColor(kBlue); //fHEffOn.Sumw2(); // effective on time versus time - fHEffOnTime.SetName("EffOnTime"); - fHEffOnTime.SetTitle("Effective On Time T_{eff}"); - fHEffOnTime.SetXTitle("Time"); - fHEffOnTime.SetYTitle("T_{eff} [s]"); - fHEffOnTime.UseCurrentStyle(); - fHEffOnTime.SetDirectory(NULL); - fHEffOnTime.GetYaxis()->SetTitleOffset(1.2); - fHEffOnTime.GetXaxis()->SetLabelSize(0.033); - fHEffOnTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); - fHEffOnTime.GetXaxis()->SetTimeDisplay(1); - fHEffOnTime.Sumw2(); + fHTimeEffOn.SetName("EffOnTime"); + fHTimeEffOn.SetTitle("Effective On Time T_{eff}"); + fHTimeEffOn.SetXTitle("Time"); + fHTimeEffOn.SetYTitle("T_{eff} [s]"); + fHTimeEffOn.UseCurrentStyle(); + fHTimeEffOn.SetDirectory(NULL); + fHTimeEffOn.GetYaxis()->SetTitleOffset(1.2); + fHTimeEffOn.GetXaxis()->SetLabelSize(0.033); + fHTimeEffOn.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fHTimeEffOn.GetXaxis()->SetTimeDisplay(1); + fHTimeEffOn.GetYaxis()->SetTitleColor(kBlue); + fHTimeEffOn.SetLineColor(kBlue); // chi2 probability - fHProbTheta.SetName("ProbTheta"); - fHProbTheta.SetTitle("\\chi^{2} Probability of Fit"); - fHProbTheta.SetXTitle("\\Theta [\\circ]"); - fHProbTheta.SetYTitle("p [%]"); - fHProbTheta.UseCurrentStyle(); - fHProbTheta.SetDirectory(NULL); - fHProbTheta.GetYaxis()->SetTitleOffset(1.2); - fHProbTheta.SetMaximum(101); + fHThetaProb.SetName("ProbTheta"); + fHThetaProb.SetTitle("\\chi^{2} Probability of Fit"); + fHThetaProb.SetXTitle("\\Theta [\\circ]"); + fHThetaProb.SetYTitle("p [%]"); + fHThetaProb.UseCurrentStyle(); + fHThetaProb.SetDirectory(NULL); + fHThetaProb.GetYaxis()->SetTitleOffset(1.2); + fHThetaProb.SetMaximum(101); + fHThetaProb.GetYaxis()->SetTitleColor(kBlue); + fHThetaProb.SetLineColor(kBlue); // chi2 probability - fHProbTime.SetName("ProbTime"); - fHProbTime.SetTitle("\\chi^{2} Probability of Fit"); - fHProbTime.SetXTitle("Time"); - fHProbTime.SetYTitle("p [%]"); - fHProbTime.UseCurrentStyle(); - fHProbTime.SetDirectory(NULL); - fHProbTime.GetYaxis()->SetTitleOffset(1.2); - fHProbTime.GetXaxis()->SetLabelSize(0.033); - fHProbTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); - fHProbTime.GetXaxis()->SetTimeDisplay(1); - fHProbTime.SetMaximum(101); + fHTimeProb.SetName("ProbTime"); + fHTimeProb.SetTitle("\\chi^{2} Probability of Fit"); + fHTimeProb.SetXTitle("Time"); + fHTimeProb.SetYTitle("p [%]"); + fHTimeProb.UseCurrentStyle(); + fHTimeProb.SetDirectory(NULL); + fHTimeProb.GetYaxis()->SetTitleOffset(1.2); + fHTimeProb.GetXaxis()->SetLabelSize(0.033); + fHTimeProb.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fHTimeProb.GetXaxis()->SetTimeDisplay(1); + fHTimeProb.SetMaximum(101); + fHTimeProb.GetYaxis()->SetTitleColor(kBlue); + fHTimeProb.SetLineColor(kBlue); // lambda versus theta - fHLambda.SetName("lambda"); - fHLambda.SetTitle("lambda of Effective On Time Fit"); - fHLambda.SetXTitle("\\Theta [\\circ]"); - fHLambda.SetYTitle("\\lambda [Hz]"); - fHLambda.UseCurrentStyle(); - fHLambda.SetDirectory(NULL); - // fHLambda.Sumw2(); - - // N0 versus theta - fHN0.SetName("N0"); - fHN0.SetTitle("Ideal number of events N_{0}"); - fHN0.SetXTitle("\\Theta [\\circ]"); - fHN0.SetYTitle("N_{0}"); - fHN0.UseCurrentStyle(); - fHN0.SetDirectory(NULL); - //fHN0del.Sumw2(); - - // chi2/NDF versus theta + fHThetaLambda.SetName("LambdaTheta"); + fHThetaLambda.SetTitle("Slope (Rate) vs Theta"); + fHThetaLambda.SetXTitle("\\Theta [\\circ]"); + fHThetaLambda.SetYTitle("S"); + fHThetaLambda.UseCurrentStyle(); + fHThetaLambda.SetDirectory(NULL); + fHThetaLambda.SetLineColor(kGreen); + + // lambda versus time + fHTimeLambda.SetName("LambdaTime"); + fHTimeLambda.SetTitle("Slope (Rate) vs Time"); + fHTimeLambda.SetXTitle("\\Time [\\circ]"); + fHTimeLambda.SetYTitle("S"); + fHTimeLambda.UseCurrentStyle(); + fHTimeLambda.SetDirectory(NULL); + fHTimeLambda.GetYaxis()->SetTitleOffset(1.2); + fHTimeLambda.GetXaxis()->SetLabelSize(0.033); + fHTimeLambda.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fHTimeLambda.GetXaxis()->SetTimeDisplay(1); + fHTimeLambda.SetLineColor(kGreen); + + // NDoF versus theta + fHThetaNDF.SetName("NDofTheta"); + fHThetaNDF.SetTitle("Number of Degrees of freedom vs Theta"); + fHThetaNDF.SetXTitle("\\Theta [\\circ]"); + fHThetaNDF.SetYTitle("NDoF [#]"); + fHThetaNDF.UseCurrentStyle(); + fHThetaNDF.SetDirectory(NULL); + fHThetaNDF.SetLineColor(kGreen); + + // NDoF versus time /* - fHChi2.SetName("Chi2/NDF"); - fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit"); - fHChi2.SetXTitle("\\Theta [\\circ]"); - fHChi2.SetYTitle("\\chi^{2}/NDF"); - fHChi2.UseCurrentStyle(); - fHChi2.SetDirectory(NULL); - fHChi2.GetYaxis()->SetTitleOffset(1.2); - //fHChi2.Sumw2(); - */ - + fHTimeNDF.SetName("NDofTime"); + fHTimeNDF.SetTitle("Number of Degrees of freedom vs Time"); + fHTimeNDF.SetXTitle("Time"); + fHTimeNDF.SetYTitle("NDoF [#]"); + fHTimeNDF.UseCurrentStyle(); + fHTimeNDF.SetDirectory(NULL); + fHTimeNDF.GetYaxis()->SetTitleOffset(1.2); + fHTimeNDF.GetXaxis()->SetLabelSize(0.033); + fHTimeNDF.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fHTimeNDF.GetXaxis()->SetTimeDisplay(1); + fHTimeNDF.SetLineColor(kBlue); + */ // setup binning MBinning btheta("BinningTheta"); @@ -217,8 +240,8 @@ btime.Apply(fH1DeltaT); - btheta.Apply(fHEffOnTheta); - btheta.Apply(fHLambda); - btheta.Apply(fHN0); - btheta.Apply(fHProbTheta); + btheta.Apply(fHThetaEffOn); + btheta.Apply(fHThetaLambda); + btheta.Apply(fHThetaNDF); + btheta.Apply(fHThetaProb); //btheta.Apply(fHChi2); } @@ -251,8 +274,8 @@ if (binstheta) { - binstheta->Apply(fHEffOnTheta); - binstheta->Apply(fHLambda); - binstheta->Apply(fHN0); - binstheta->Apply(fHProbTheta); + binstheta->Apply(fHThetaEffOn); + binstheta->Apply(fHThetaLambda); + binstheta->Apply(fHThetaNDF); + binstheta->Apply(fHThetaProb); //binstheta->Apply(fHChi2); } @@ -295,6 +318,6 @@ // parameter 1 = N0*del // - TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]); - func.SetParNames("lambda", "N0", "del"); + static TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)"); + //func.SetParNames("lambda", "N0", "del"); func.SetParameter(0, 100); // Hz @@ -302,9 +325,9 @@ func.FixParameter(2, Nmdel/Nm); - // options : 0 do not plot the function + // options : N do not store the function, do not draw // I use integral of function in bin rather than value at bin center // R use the range specified in the function range // Q quiet mode - h->Fit(&func, "0IRQ"); + h->Fit(&func, "NIQ", "", yq[0], yq[1]); const Double_t chi2 = func.GetChisquare(); @@ -325,5 +348,5 @@ const Double_t lambda = func.GetParameter(0); - const Double_t N0 = func.GetParameter(1); + //const Double_t N0 = func.GetParameter(1); const Double_t prob = func.GetProb(); @@ -362,6 +385,6 @@ res[4] = dl; - // N0 of fit - res[5] = N0; + // NDoF of fit + res[5] = NDF; // Rdead (from fit) is the fraction from real time lost by the dead time @@ -379,8 +402,8 @@ void MHEffectiveOnTime::FitThetaBins() { - fHEffOnTheta.Reset(); - fHProbTheta.Reset(); - fHLambda.Reset(); - fHN0.Reset(); + fHThetaEffOn.Reset(); + fHThetaProb.Reset(); + fHThetaLambda.Reset(); + fHThetaNDF.Reset(); const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999)); @@ -400,9 +423,9 @@ // the effective on time is Nm/lambda - fHEffOnTheta.SetBinContent(i, res[0]); - fHEffOnTheta.SetBinError (i, res[1]); + fHThetaEffOn.SetBinContent(i, res[0]); + fHThetaEffOn.SetBinError (i, res[1]); // plot chi2-probability of fit - fHProbTheta.SetBinContent(i, res[2]); + fHThetaProb.SetBinContent(i, res[2]); // plot chi2/NDF of fit @@ -410,9 +433,9 @@ // lambda of fit - fHLambda.SetBinContent(i, res[3]); - fHLambda.SetBinError (i, res[4]); - - // N0 of fit - fHN0.SetBinContent(i, res[5]); + fHThetaLambda.SetBinContent(i, res[3]); + fHThetaLambda.SetBinError (i, res[4]); + + // NDoF of fit + fHThetaNDF.SetBinContent(i, res[5]); // Rdead (from fit) is the fraction from real time lost by the dead time @@ -449,20 +472,27 @@ // Get number of bins - const Int_t n = fHEffOnTime.GetNbinsX(); + const Int_t n = fHTimeEffOn.GetNbinsX(); // Enhance binning MBinning bins; - bins.SetEdges(fHEffOnTime, 'x'); + bins.SetEdges(fHTimeEffOn, 'x'); bins.AddEdge(fLastTime.GetAxisTime()); - bins.Apply(fHEffOnTime); - bins.Apply(fHProbTime); + bins.Apply(fHTimeEffOn); + bins.Apply(fHTimeProb); + bins.Apply(fHTimeLambda); + //bins.Apply(fHTimeNDF); // // Fill histogram // - fHEffOnTime.SetBinContent(n, res[0]); - fHEffOnTime.SetBinError(n, res[1]); - - fHProbTime.SetBinContent(n, res[2]); + fHTimeEffOn.SetBinContent(n, res[0]); + fHTimeEffOn.SetBinError(n, res[1]); + + fHTimeProb.SetBinContent(n, res[2]); + + fHTimeLambda.SetBinContent(n, res[3]); + fHTimeLambda.SetBinError(n, res[4]); + + //fHTimeNDF.SetBinContent(n, res[5]); // @@ -500,6 +530,8 @@ MBinning bins; bins.SetEdges(1, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime()); - bins.Apply(fHEffOnTime); - bins.Apply(fHProbTime); + bins.Apply(fHTimeEffOn); + bins.Apply(fHTimeProb); + bins.Apply(fHTimeLambda); + //bins.Apply(fHTimeNDF); fParam->SetVal(0, 0); @@ -540,6 +572,8 @@ FitThetaBins(); FitTimeBin(); - MH::RemoveFirstBin(fHEffOnTime); - MH::RemoveFirstBin(fHProbTime); + MH::RemoveFirstBin(fHTimeEffOn); + MH::RemoveFirstBin(fHTimeProb); + MH::RemoveFirstBin(fHTimeLambda); + //MH::RemoveFirstBin(fHTimeNDF); fIsFinalized = kTRUE; @@ -558,4 +592,49 @@ text.SetTextSize(0.04); text.Paint(); +} + +void MHEffectiveOnTime::PaintProb(TH1 &h) const +{ + Double_t sum = 0; + Int_t n = 0; + for (int i=0; i0) + { + sum += h.GetBinContent(i+1); + n++; + } + + if (n==0) + return; + + TLatex text(0.45, 0.94, Form("\\bar{p} = %.1f%% (n=%d)", sum/n, n)); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); +} + +void MHEffectiveOnTime::UpdateRightAxis(TH1 &h) +{ + cout << "Update... " << flush; + const Double_t max = h.GetMaximum()*1.1; + if (max==0) + return; + + h.SetNormFactor(h.Integral()*gPad->GetUymax()/max); + + TGaxis *axis = (TGaxis*)gPad->FindObject("RightAxis"); + if (!axis) + return; + + cout << axis << " "; + cout << "X: " << gPad->GetUxmax() << " "; + cout << "Y: " << gPad->GetUymin() << "/" << gPad->GetUymax() << " "; + cout << "max=" << max << endl; + + axis->SetX1(gPad->GetUxmax()); + axis->SetX2(gPad->GetUxmax()); + axis->SetY1(gPad->GetUymin()); + axis->SetY2(gPad->GetUymax()); + axis->SetWmax(max); } @@ -621,9 +700,29 @@ } + if (o==(TString)"timendf") + { + // UpdateRightAxis(fHTimeNDF); + // FIXME: first bin? + PaintProb(fHTimeProb); + } + + if (o==(TString)"thetandf") + { + UpdateRightAxis(fHThetaNDF); + // FIXME: first bin? + PaintProb(fHThetaProb); + } + h=0; if (o==(TString)"theta") - h = &fHEffOnTheta; + { + h = &fHThetaEffOn; + UpdateRightAxis(fHThetaLambda); + } if (o==(TString)"time") - h = &fHEffOnTime; + { + h = &fHTimeEffOn; + UpdateRightAxis(fHTimeLambda); + } if (!h) @@ -635,4 +734,18 @@ PaintText(h->Integral(), error); +} + +void MHEffectiveOnTime::DrawRightAxis(const char *title) +{ + TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(), + gPad->GetUxmax(), gPad->GetUymax(), + 0, 1, 510, "+L"); + axis->SetName("RightAxis"); + axis->SetTitle(title); + axis->SetTitleOffset(0.9); + axis->SetTextColor(kGreen); + axis->CenterTitle(); + axis->SetBit(kCanDelete); + axis->Draw(); } @@ -669,11 +782,15 @@ pad->GetPad(1)->cd(2); gPad->SetBorderMode(0); - fHProbTime.Draw(); - AppendPad("paint"); + fHTimeProb.Draw(); + AppendPad("timendf"); + //fHTimeNDF.Draw("same"); + //DrawRightAxis("NDF"); pad->GetPad(1)->cd(3); gPad->SetBorderMode(0); - fHEffOnTime.Draw(); + fHTimeEffOn.Draw(); AppendPad("time"); + fHTimeLambda.Draw("same"); + DrawRightAxis("Slope"); pad->cd(2); @@ -695,9 +812,14 @@ pad->GetPad(2)->cd(2); gPad->SetBorderMode(0); - fHProbTheta.Draw(); + fHThetaProb.Draw(); + AppendPad("thetandf"); + fHThetaNDF.Draw("same"); + DrawRightAxis("NDF"); pad->GetPad(2)->cd(3); gPad->SetBorderMode(0); - fHEffOnTheta.Draw(); + fHThetaEffOn.Draw(); AppendPad("theta"); -} + fHThetaLambda.Draw("same"); + DrawRightAxis("Slope"); +} Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4974) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4975) @@ -13,4 +13,7 @@ #ifndef ROOT_TH2 #include +#endif +#ifndef ROOT_TF1 +#include #endif @@ -33,14 +36,15 @@ TH1D fH1DeltaT; //! Counts vs Delta T (for a time interval) - TH1D fHEffOnTheta; // Effective On time versus Theta - TH1D fHEffOnTime; // Effective On time versus Time + TH1D fHThetaEffOn; // Effective On time versus Theta + TH1D fHThetaProb; // Chisq prob fit of Effective On time versus Theta + TH1D fHThetaNDF; // NDF vs Theta + TH1D fHThetaLambda; // Slope (rate) vs Theta - TH1D fHProbTheta; // Chisq prob fit of Effective On time versus Theta - TH1D fHProbTime; // Chisq prob fit of Effective On time versus Time + TH1D fHTimeEffOn; // Effective On time versus Time + TH1D fHTimeProb; // Chisq prob fit of Effective On time versus Time + //TH1D fHTimeNDF; // NDF vs Time + TH1D fHTimeLambda; // Slope (rate) vs Time - TH1D fHN0; - TH1D fHLambda; - - Bool_t fIsFinalized; // Flag telling you whether fHEffOnTheta is the final result + Bool_t fIsFinalized; // Flag telling you whether fHThetaEffOn is the final result Int_t fNumEvents; // Number of events to be used for a bin in time @@ -52,5 +56,8 @@ void FitThetaBins(); void FitTimeBin(); + void PaintProb(TH1 &h) const; void PaintText(Double_t val, Double_t error) const; + void DrawRightAxis(const char *title); + void UpdateRightAxis(TH1 &h); public: @@ -63,6 +70,6 @@ Bool_t Finalize(); - const TH1D &GetHEffOnTheta() const { return fHEffOnTheta; } - const TH1D &GetHEffOnTime() const { return fHEffOnTime; } + const TH1D &GetHEffOnTheta() const { return fHThetaEffOn; } + const TH1D &GetHEffOnTime() const { return fHTimeEffOn; } void Draw(Option_t *option="");