Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4887) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4889) @@ -42,4 +42,5 @@ #include #include +#include #include "MTime.h" @@ -62,5 +63,5 @@ // MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title) - : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60) + : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60) { // @@ -141,6 +142,7 @@ } -Double_t testval = 0; -Double_t testerr = 0; +// FIXME: Just for a preliminary check +static Double_t testval = 0; +static Double_t testerr = 0; // -------------------------------------------------------------------------- @@ -185,116 +187,153 @@ Bool_t MHEffectiveOnTime::Finalize() { - Fit(); + FitThetaBins(); + Calc(); + fIsFinalized = kTRUE; + return kTRUE; } -void MHEffectiveOnTime::Fit() -{ +void MHEffectiveOnTime::FitThetaBins() +{ + const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999)); + // nbins = number of Theta bins const Int_t nbins = fHTimeDiff.GetNbinsY(); + TH1D *h=0; for (int i=1; i<=nbins; i++) { // TH1D &h = *hist->ProjectionX("Calc-theta", i, i); - TH1D *h = fHTimeDiff.ProjectionX("CalcTheta", i, i, "E"); - - const Double_t Nm = h->Integral(); - - if (Nm<=0) + h = fHTimeDiff.ProjectionX(name, i, i, "E"); + + Double_t res[7]; + if (!FitH(h, res)) continue; - // 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); - - // Nmdel = Nm * binwidth, with Nm = number of observed events - const Double_t Nmdel = h->Integral("width"); - - // - // Setup Poisson function for the fit: - // lambda [Hz], N0 = ideal no of evts, del = bin width of dt - // - // parameter 0 = lambda - // parameter 1 = N0*del - // - TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]); - func.SetParNames("lambda", "N0", "del"); - - func.SetParameter(0, 100); // Hz - func.SetParameter(1, Nm); - func.FixParameter(2, Nmdel/Nm); - - // 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(&func, "0IRQ"); - - const Double_t chi2 = func.GetChisquare(); - const Int_t NDF = func.GetNDF(); - - // was fit successful ? - if (NDF>0 && chi2<2.5*NDF) - { - const Double_t lambda = func.GetParameter(0); - const Double_t N0 = func.GetParameter(1); - const Double_t prob = func.GetProb(); - - /* - *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF="; - *fLog << (NDF ? chi2/NDF : 0.0) << " lambda="; - *fLog << lambda << " N0=" << N0 << endl; - */ - - Double_t emat[2][2]; - gMinuit->mnemat((Double_t*)emat, 2); - - const Double_t dldl = emat[0][0]; - //const Double_t dN0dN0 = emat[1][1]; - - const Double_t teff = Nm/lambda; - const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm); - - const Double_t dl = TMath::Sqrt(dldl); - - //const Double_t kappa = Nm/N0; - //const Double_t Rdead = 1.0 - kappa; - //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm); - - // the effective on time is Nm/lambda - fHEffOn.SetBinContent(i, teff); - fHEffOn.SetBinError (i, dteff); - - // plot chi2-probability of fit - fHProb.SetBinContent(i, prob*100); - - // plot chi2/NDF of fit - fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0); - - // lambda of fit - fHLambda.SetBinContent(i, lambda); - fHLambda.SetBinError (i, dl); - - // N0 of fit - fHN0.SetBinContent(i, N0); - - // Rdead (from fit) is the fraction from real time lost by the dead time - //fHRdead.SetBinContent(i, Rdead); - //fHRdead.SetBinError (i,dRdead); - } - + // the effective on time is Nm/lambda + fHEffOn.SetBinContent(i, res[0]); + fHEffOn.SetBinError (i, res[1]); + + // plot chi2-probability of fit + fHProb.SetBinContent(i, res[2]); + + // plot chi2/NDF of fit + fHChi2.SetBinContent(i, res[3]); + + // lambda of fit + fHLambda.SetBinContent(i, res[4]); + fHLambda.SetBinError (i, res[5]); + + // N0 of fit + fHN0.SetBinContent(i, res[6]); + + // Rdead (from fit) is the fraction from real time lost by the dead time + //fHRdead.SetBinContent(i, Rdead); + //fHRdead.SetBinError (i,dRdead); + } + + // Histogram is reused via gROOT->FindObject() + // Need to be deleted only once + if (h) delete h; - } -} - +} + +Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const +{ + const Double_t Nm = h->Integral(); + + // FIXME: Do fit only if contents of bin has changed + if (Nm<=0) + return kFALSE; + + // 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.05, 0.95 }; + Double_t yq[2]; + h->GetQuantiles(2, yq, xq); + + // Nmdel = Nm * binwidth, with Nm = number of observed events + const Double_t Nmdel = h->Integral("width"); + + // + // Setup Poisson function for the fit: + // lambda [Hz], N0 = ideal no of evts, del = bin width of dt + // + // parameter 0 = lambda + // parameter 1 = N0*del + // + TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]); + func.SetParNames("lambda", "N0", "del"); + + func.SetParameter(0, 100); // Hz + func.SetParameter(1, Nm); + func.FixParameter(2, Nmdel/Nm); + + // 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(&func, "0IRQ"); + + const Double_t chi2 = func.GetChisquare(); + const Int_t NDF = func.GetNDF(); + + // was fit successful ? + if (NDF<=0 || chi2>=2.5*NDF) + return kFALSE; + + const Double_t lambda = func.GetParameter(0); + const Double_t N0 = func.GetParameter(1); + const Double_t prob = func.GetProb(); + + /* + *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF="; + *fLog << (NDF ? chi2/NDF : 0.0) << " lambda="; + *fLog << lambda << " N0=" << N0 << endl; + */ + + Double_t emat[2][2]; + gMinuit->mnemat((Double_t*)emat, 2); + + const Double_t dldl = emat[0][0]; + //const Double_t dN0dN0 = emat[1][1]; + + const Double_t teff = Nm/lambda; + const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm); + + const Double_t dl = TMath::Sqrt(dldl); + + //const Double_t kappa = Nm/N0; + //const Double_t Rdead = 1.0 - kappa; + //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm); + + // the effective on time is Nm/lambda + res[0] = teff; + res[1] = dteff; + + // plot chi2-probability of fit + res[2] = prob*100; + + // plot chi2/NDF of fit + res[3] = NDF ? chi2/NDF : 0.0; + + // lambda of fit + res[4] = lambda; + res[5] = dl; + + // N0 of fit + res[6] = N0; + + // Rdead (from fit) is the fraction from real time lost by the dead time + //fHRdead.SetBinContent(i, Rdead); + //fHRdead.SetBinError (i,dRdead); + + return kTRUE; +} void MHEffectiveOnTime::Paint(Option_t *opt) @@ -320,5 +359,5 @@ if (!fIsFinalized) - Fit(); + FitThetaBins(); } if (o==(TString)"paint") @@ -386,9 +425,16 @@ void MHEffectiveOnTime::Calc() { - const Double_t val = fHEffOn.Integral(); - - Double_t error = 0; - for (int i=0; iGetNbins(); i++) - error += fHEffOn.GetBinError(i); + TH1D *h = fHTimeDiff.ProjectionX("", -1, 99999, "E"); + h->SetDirectory(0); + + Double_t res[7]; + Bool_t rc = FitH(h, res); + delete h; + + if (!rc) + return; + + const Double_t val = res[0]; + const Double_t error = res[1]; fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0); @@ -403,4 +449,5 @@ MTime now(*fTime); now.AddMilliSeconds(fInterval*1000); + *fLog <MinusNull() + + // Make this just a ns before the first event + fTime->Minus1ns(); } if (*fTime+fInterval<*time) { - Fit(); + FitThetaBins(); Calc(); } Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4887) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4889) @@ -24,12 +24,12 @@ { private: - MPointingPos *fPointPos; //! - MTime fLastTime; //! + MPointingPos *fPointPos; //! + MTime fLastTime; //! - Double_t fEffOnTime0; //! - Double_t fEffOnErr0; //! + Double_t fEffOnTime0; //! + Double_t fEffOnErr0; //! - MTime *fTime; - MParameterDerr *fParam; + MTime *fTime; //! + MParameterDerr *fParam; //! TH2D fHTimeDiff; @@ -44,5 +44,6 @@ Int_t fInterval; - void Fit(); + Bool_t FitH(TH1D *h, Double_t *res) const; + void FitThetaBins(); void Calc();