Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4887) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4887) @@ -0,0 +1,451 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 8/2002 +! Author(s): Wolfgang Wittek, 1/2002 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHEffectiveOnTime +// +// Fills a 2-D histogram Time-difference vs. Theta +// +// From this histogram the effective on-time is determined by a fit. +// The result of the fit (see Fit()) and the fit-parameters (like chi^2) +// are stored in corresponding histograms +// +////////////////////////////////////////////////////////////////////////////// + +#include "MHEffectiveOnTime.h" + +#include +#include +#include +#include + +#include "MTime.h" +#include "MParameters.h" +#include "MPointingPos.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHEffectiveOnTime); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title) + : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHEffectiveOnTime"; + fTitle = title ? title : "2-D histogram in Theta and time difference"; + + // Main histogram + fHTimeDiff.SetXTitle("\\Delta t [s]"); + fHTimeDiff.SetYTitle("\\Theta [\\circ]"); + fHTimeDiff.UseCurrentStyle(); + fHTimeDiff.SetDirectory(NULL); + + // effective on time versus theta + fHEffOn.SetName("EffOnTime"); + fHEffOn.SetTitle("Effective On Time T_{eff}"); + fHEffOn.SetXTitle("\\Theta [\\circ]"); + fHEffOn.SetYTitle("T_{eff} [s]"); + fHEffOn.UseCurrentStyle(); + fHEffOn.SetDirectory(NULL); + fHEffOn.GetYaxis()->SetTitleOffset(1.2); + //fHEffOn.Sumw2(); + + // chi2/NDF versus theta + 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(); + + // chi2 probability + fHProb.SetName("Prob"); + fHProb.SetTitle("\\chi^{2} Probability of Fit"); + fHProb.SetXTitle("\\Theta [\\circ]"); + fHProb.SetYTitle("p [%]"); + fHProb.UseCurrentStyle(); + fHProb.SetDirectory(NULL); + fHProb.GetYaxis()->SetTitleOffset(1.2); + fHProb.SetMaximum(101); + //fHChi2.Sumw2(); + + // 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(); + + // setup binning + MBinning btheta("BinningTheta"); + btheta.SetEdgesCos(51, 0, 60); + + MBinning btime("BinningTimeDiff"); + btime.SetEdges(50, 0, 0.1); + + MH::SetBinning(&fHTimeDiff, &btime, &btheta); + + btheta.Apply(fHEffOn); + btheta.Apply(fHChi2); + btheta.Apply(fHLambda); + btheta.Apply(fHN0); + btheta.Apply(fHProb); +} + +Double_t testval = 0; +Double_t testerr = 0; + +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histogram +// +Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist) +{ + fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); + if (!fPointPos) + { + *fLog << err << dbginf << "MPointingPos not found... aborting." << endl; + return kFALSE; + } + + // FIXME: Remove const-qualifier from base-class! + fTime = (MTime*)const_cast(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime"); + if (!fTime) + return kFALSE; + fParam = (MParameterDerr*)const_cast(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime"); + if (!fParam) + return kFALSE; + + const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningTimeDiff"); + const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binstheta || !binsdtime) + *fLog << warn << dbginf << "At least one MBinning not found... ignored." << endl; + else + { + SetBinning(&fHTimeDiff, binsdtime, binstheta); + + binstheta->Apply(fHEffOn); + binstheta->Apply(fHChi2); + binstheta->Apply(fHLambda); + binstheta->Apply(fHN0); + binstheta->Apply(fHProb); + } + + return kTRUE; +} + +Bool_t MHEffectiveOnTime::Finalize() +{ + Fit(); + fIsFinalized = kTRUE; + return kTRUE; +} + +void MHEffectiveOnTime::Fit() +{ + // nbins = number of Theta bins + const Int_t nbins = fHTimeDiff.GetNbinsY(); + + 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) + 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); + } + + delete h; + } +} + + +void MHEffectiveOnTime::Paint(Option_t *opt) +{ + TVirtualPad *padsave = gPad; + + TString o(opt); + if (o==(TString)"fit") + { + TH1D *h0=0; + + padsave->cd(1); + if ((h0 = (TH1D*)gPad->FindObject("TimeDiff"))) + { + h0 = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E"); + if (h0->GetEntries()>0) + gPad->SetLogy(); + } + + padsave->cd(2); + if ((h0 = (TH1D*)gPad->FindObject("Theta"))) + fHTimeDiff.ProjectionY("Theta", -1, 9999, "E"); + + if (!fIsFinalized) + Fit(); + } + if (o==(TString)"paint") + { + Double_t error = 0; + for (int i=0; iGetNbins(); i++) + error += fHEffOn.GetBinError(i); + + TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", fHEffOn.Integral(), error)); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); + } + gPad = padsave; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHEffectiveOnTime::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad("fit"); + + pad->Divide(2,2); + + TH1 *h; + + pad->cd(1); + gPad->SetBorderMode(0); + h = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E"); + h->SetTitle("Distribution of \\Delta t [s]"); + h->SetXTitle("\\Delta t [s]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->Draw(); + + pad->cd(2); + gPad->SetBorderMode(0); + h = fHTimeDiff.ProjectionY("Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta [\\circ]"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->GetYaxis()->SetTitleOffset(1.1); + h->Draw(); + + pad->cd(3); + gPad->SetBorderMode(0); + fHEffOn.Draw(); + AppendPad("paint"); + + pad->cd(4); + gPad->SetBorderMode(0); + fHProb.Draw(); +} + +void MHEffectiveOnTime::Calc() +{ + const Double_t val = fHEffOn.Integral(); + + Double_t error = 0; + for (int i=0; iGetNbins(); i++) + error += fHEffOn.GetBinError(i); + + fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0); + fParam->SetReadyToSave(); + + testval += fParam->GetVal(); + testerr += fParam->GetErr(); + + fEffOnTime0 = val; + fEffOnErr0 = error; + + MTime now(*fTime); + now.AddMilliSeconds(fInterval*1000); + *fLog <GetVal(), fParam->GetErr()); + *fLog << Form(" %.1f %.1f %.1f %.1f", + val, testval, error, testerr) << endl; + + fTime->AddMilliSeconds(fInterval*1000); +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w) +{ + const MTime *time = dynamic_cast(par); + if (!time) + { + *fLog << err << "ERROR - MHEffectiveOnTimeTime::Fill without argument or container doesn't inherit from MTime... abort." << endl; + return kFALSE; + } + + if (*fTime==MTime()) + { + *fLog << all << *time << " - " << *fTime << " " << TMath::Floor(*time) << endl; + fEffOnTime0 = 0; + fEffOnErr0 = 0; + + fParam->SetVal(0, 0); + fParam->SetReadyToSave(); + + *fTime = *time; + // fTime->MinusNull() + } + + if (*fTime+fInterval<*time) + { + Fit(); + Calc(); + } + + fHTimeDiff.Fill(*time-fLastTime, fPointPos->GetZd(), w); + fLastTime = *time; + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4887) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4887) @@ -0,0 +1,70 @@ +#ifndef MARS_MHEffectiveOnTime +#define MARS_MHEffectiveOnTime + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef MARS_MTime +#include "MTime.h" +#endif +#ifndef ROOT_TH1 +#include +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MTime; +class MPointingPos; +class MParameterDerr; + +class MParList; + +class MHEffectiveOnTime : public MH +{ +private: + MPointingPos *fPointPos; //! + MTime fLastTime; //! + + Double_t fEffOnTime0; //! + Double_t fEffOnErr0; //! + + MTime *fTime; + MParameterDerr *fParam; + + TH2D fHTimeDiff; + TH1D fHEffOn; + TH1D fHChi2; + TH1D fHProb; + TH1D fHN0; + TH1D fHLambda; + + Bool_t fIsFinalized; + + Int_t fInterval; + + void Fit(); + void Calc(); + +public: + MHEffectiveOnTime(const char *name=NULL, const char *title=NULL); + + void SetInterval(Int_t i) { fInterval=i; } + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + const TH2D *GetHist() { return &fHTimeDiff; } + const TH2D *GetHist() const { return &fHTimeDiff; } + + TH1 *GetHistByName(const TString name) { return &fHTimeDiff; } + + void Draw(Option_t *option=""); + void Paint(Option_t *opt=""); + + ClassDef(MHEffectiveOnTime, 1) // 2D-histogram to determin effective on-time vs. theta +}; + +#endif +