Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc (revision 1292) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc (revision 1293) @@ -26,10 +26,13 @@ ////////////////////////////////////////////////////////////////////////////// // // -// MHEffOnTimeTheta // +// MHEffOnTimeTheta // // // +// calculates the effective on time for each bin in Theta // // // ////////////////////////////////////////////////////////////////////////////// #include "MHEffOnTimeTheta.h" + +#include #include @@ -47,12 +50,10 @@ ClassImp(MHEffOnTimeTheta); - // -------------------------------------------------------------------------- // -// Default Constructor. It sets name and title only. Typically you won't -// need to change this. +// Default Constructor. It sets name and title of the histograms. // MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title) - : fHist() + : fHEffOn() { // @@ -62,56 +63,165 @@ fTitle = title ? title : "1-D histogram of Eff On Time"; - fHist.SetName("EffOn"); - fHist.SetTitle("Effective On Time Vs. Theta"); - - fHist.SetDirectory(NULL); - - fHist.SetXTitle("t_{eff} [s]"); - fHist.SetYTitle("\\Theta [\\circ]"); -} - + // effective on time versus theta + fHEffOn.SetName("EffOn"); + fHEffOn.SetTitle("Effective On Time vs. Theta"); + + fHEffOn.SetDirectory(NULL); + + fHEffOn.SetXTitle("\\Theta [\\circ]"); + fHEffOn.SetYTitle("t-eff [s]"); + + // chi2/NDF versus theta + fHChi2.SetName("Chi2/NDF"); + fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Theta"); + + fHChi2.SetDirectory(NULL); + + fHChi2.SetXTitle("\\Theta [\\circ]"); + fHChi2.SetYTitle("chi2/NDF"); + + // lambda versus theta + fHLambda.SetName("lambda"); + fHLambda.SetTitle("lambda of OnTimeFit vs. Theta"); + + fHLambda.SetDirectory(NULL); + + fHLambda.SetXTitle("\\Theta [\\circ]"); + fHLambda.SetYTitle("\\lambda [Hz]"); + + // N0del versus theta + fHN0del.SetName("N0del"); + fHN0del.SetTitle("N0del of OnTimeFit vs. Theta"); + + fHN0del.SetDirectory(NULL); + + fHN0del.SetXTitle("\\Theta [\\circ]"); + fHN0del.SetYTitle("N0del"); +} + +// ----------------------------------------------------------------------- +// +// Calculate the effective on time by fitting the distribution of +// time differences +// void MHEffOnTimeTheta::Calc(TH2D *hist) { + // nbins = number of Theta bins const Int_t nbins = hist->GetNbinsY(); - for (int i=1; iProjectionX("dTime-Distribution for fixed Time", i, i); - - //hist->Draw(); - //gPad->SetLogy(); - - Double_t Nmdel = h.Integral("width"); - Double_t mean = h.GetMean(); - - TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", - mean, hist->GetXaxis()->GetXmax()); - - func.SetParameter(0, 100); // [Hz] - func.SetParameter(1, Nmdel); - - func.SetParLimits(0, 0, 1000); // [Hz] - func.SetParLimits(1, 0, 2*Nmdel); - - func.SetParName(0, "lambda"); - func.SetParName(1, "Nmdel"); - - h.Fit("Poisson", "RN"); - - //func.SetRange(0, 0.1); // Range of Drawing - //func.DrawCopy("same"); - - //cout << func.GetParameter(0) << " " << func.GetParameter(1) << endl; - - Double_t lambda = func.GetParameter(0); - //Double_t a = func.GetParameter(1); - - //cout << "t_eff = " << h->Integral()/lambda << " T(last)=" << time.GetTimeLo()*0.0001 << endl; - - fHist.SetBinContent(i, h.Integral()/lambda); + // TH1D &h = *hist->ProjectionX("Calc-theta", i, i); + TH1D &h = *hist->ProjectionX("Calc-theta", i, i, "E"); + + // Nmdel = Nm * binwidth, with Nm = number of observed events + const Double_t Nmdel = h.Integral("width"); + const Double_t Nm = h.Integral(); + // Double_t mean = h->GetMean(); + + //................................................... + // determine range (yq[0], yq[1]) of time differences + // where fit should be performed; + // require a fraction >=xq[0] of all entries to lie below yq[0] + // and a fraction <=xq[1] of all entries to lie below yq[1]; + // within the range (yq[0], yq[1]) there must be no empty bin; + // choose pedestrian approach as long as GetQuantiles is not available + + Double_t xq[2] = { 0.15, 0.95 }; + Double_t yq[2]; + + // GetQuantiles doesn't seem to be available in root 3.01/06 + // h->GetQuantiles(2,yq,xq); + + const Double_t sumtot = h.Integral(); + const Int_t jbins = h.GetNbinsX(); + + if (sumtot > 0.0) + { + Double_t sum1 = 0.0; + yq[0] = h.GetBinLowEdge(jbins+1); + for (int j=1; j<=jbins; j++) + { + if (sum1 >= xq[0]*sumtot) + { + yq[0] = h.GetBinLowEdge(j); + break; + } + sum1 += h.GetBinContent(j); + } + + Double_t sum2 = 0.0; + yq[1] = h.GetBinLowEdge(jbins+1); + for (int j=1; j<=jbins; j++) + { + Double_t content = h.GetBinContent(j); + sum2 += content; + if (sum2 >= xq[1]*sumtot || content == 0.0) + { + yq[1] = h.GetBinLowEdge(j); + break; + } + } + + //................................................... + + // parameter 0 = lambda + // parameter 1 = N0*del with N0 = ideal number of events + // and del = bin width of time difference + TF1 func("Poisson", "[1] * [0] * exp(-[0] *x)", yq[0], yq[1]); + + func.SetParameter(0, 100); // [Hz] + func.SetParameter(1, Nmdel); + + func.SetParLimits(0, 0, 1000); // [Hz] + func.SetParLimits(1, 0, 10*Nmdel); + + func.SetParName(0, "lambda"); + func.SetParName(1, "Nmdel"); + + // options : 0 do not plot the function + // 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("Poisson", "0IRQ"); + + // gPad->SetLogy(); + // gStyle->SetOptFit(1011); + // h->GetXaxis()->SetTitle("time difference [s]"); + // h->GetYaxis()->SetTitle("Counts"); + // h->DrawCopy(); + + // func.SetRange(yq[0], yq[1]); // Range of Drawing + // func.DrawCopy("same"); + + const Double_t lambda = func.GetParameter(0); + const Double_t N0del = func.GetParameter(1); + const Double_t chi2 = func.GetChisquare(); + const Int_t NDF = func.GetNDF(); + + // was fit successful ? + if (NDF>0 && chi2<2.5*NDF) + { + // the effective on time is Nm/lambda + fHEffOn.SetBinContent(i, Nm/lambda); + + // plot chi2/NDF of fit + fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0); + + // lambda of fit + fHLambda.SetBinContent(i, lambda); + + // N0del of fit + fHN0del.SetBinContent(i, N0del); + + delete &h; + continue; + } + } + + fHEffOn.SetBinContent (i, 1.e-20); + fHChi2.SetBinContent (i, 1.e-20); + fHLambda.SetBinContent(i, 1.e-20); + fHN0del.SetBinContent (i, 1.e-20); delete &h; @@ -119,4 +229,8 @@ } +// ------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histograms +// Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist) { @@ -128,29 +242,80 @@ } - SetBinning(&fHist, bins); + SetBinning(&fHEffOn, bins); + SetBinning(&fHChi2, bins); + SetBinning(&fHLambda, bins); + SetBinning(&fHN0del, bins); + + fHEffOn.Sumw2(); + fHChi2.Sumw2(); + fHLambda.Sumw2(); + fHN0del.Sumw2(); return kTRUE; } +// ------------------------------------------------------------------------- +// +// Dummy Fill +// without it get the error message : +// Error: MHEffOnTimeTime() no default constructor FILE:macros/wowflux.C LINE:359 +// *** Interpreter error recovered *** +// +Bool_t MHEffOnTimeTheta::Fill(const MParContainer *par) +{ + return kTRUE; +} +// ------------------------------------------------------------------------- +// +// Draw a copy of the histogram +// TObject *MHEffOnTimeTheta::DrawClone(Option_t *opt) const { - TCanvas *c = MakeDefCanvas("EffOnTimeTheta", "t_eff vs. Theta"); + TCanvas &c = *MakeDefCanvas("EffOnTimeTheta", "Results from on time fit vs. Theta"); + + c.Divide(2, 2); gROOT->SetSelectedPad(NULL); - ((TH2&)fHist).DrawCopy(opt); - - c->Modified(); - c->Update(); - - return c; -} - + c.cd(1); + ((TH2*)&fHEffOn)->DrawCopy(opt); + + c.cd(2); + ((TH2*)&fHChi2)->DrawCopy(opt); + + c.cd(3); + ((TH2*)&fHLambda)->DrawCopy(opt); + + c.cd(4); + ((TH2*)&fHN0del)->DrawCopy(opt); + + c.Modified(); + c.Update(); + + return &c; +} + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// void MHEffOnTimeTheta::Draw(Option_t *opt) { if (!gPad) - MakeDefCanvas("EffOnTimeTheta", "t_eff vs. Theta"); - - fHist.Draw(opt); + MakeDefCanvas("EffOnTimeTheta", "Results from on time fit vs. Theta"); + + gPad->Divide(2,2); + + gPad->cd(1); + fHEffOn.Draw(opt); + + gPad->cd(2); + fHChi2.Draw(opt); + + gPad->cd(3); + fHLambda.Draw(opt); + + gPad->cd(4); + fHN0del.Draw(opt); gPad->Modified(); @@ -159,7 +324,13 @@ -Bool_t MHEffOnTimeTheta::Fill(const MParContainer *par) -{ - return kTRUE; -} - + + + + + + + + + + +