Index: trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 5065) +++ trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 5066) @@ -20,7 +20,4 @@ #pragma link C++ class MHFlux; #pragma link C++ class MHSigmaTheta; -#pragma link C++ class MHEffOnTime+; -#pragma link C++ class MHTimeDiffTime+; -#pragma link C++ class MHTimeDiffTheta+; #pragma link C++ class MHThetabarTime+; #pragma link C++ class MHThetabarTheta+; Index: trunk/MagicSoft/Mars/mhist/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlpha.cc (revision 5065) +++ (revision ) @@ -1,334 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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, 3/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHAlpha -// -// Create a single Alpha-Plot. The alpha-plot is fitted online. You can -// check the result when it is filles in the MStatusDisplay -// For more information see MHFalseSource::FitSignificance -// -// For convinience (fit) the output significance is stored in a -// container in the parlisrt -// -// PRELIMINARY! -// -////////////////////////////////////////////////////////////////////////////// -#include "MHAlpha.h" - -#include -#include -#include -#include -#include -#include -#include - -#include "MGeomCam.h" -#include "MSrcPosCam.h" -#include "MHillasSrc.h" -#include "MTime.h" -#include "MObservatory.h" -#include "MPointingPos.h" -#include "MAstroSky2Local.h" -#include "MStatusDisplay.h" -#include "MParameters.h" -#include "MHMatrix.h" - -#include "MMath.h" -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHAlpha); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor -// -MHAlpha::MHAlpha(const char *name, const char *title) - : fAlphaCut(12.5), fBgMean(55), fResult(0), fMatrix(0) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHAlpha"; - fTitle = title ? title : "Alpha plot"; - - fHist.SetDirectory(NULL); - - fHist.SetName("Alpha"); - fHist.SetTitle("Alpha"); - fHist.SetXTitle("|\\alpha| [\\circ]"); - fHist.SetYTitle("Counts"); - fHist.UseCurrentStyle(); - - MBinning binsa; - binsa.SetEdges(18, 0, 90); - binsa.Apply(fHist); - - fSigInt=15; - fSigMax=75; - fBgMin=45; - fBgMax=90; - fPolynom=1; -} - -Bool_t MHAlpha::SetupFill(const MParList *pl) -{ - fHist.Reset(); - - fResult = (MParameterD*)pl->FindObject("Significance", "MParameterD"); - if (!fResult) - *fLog << warn << "Significance [MParameterD] not found... ignored." << endl; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the histogram. For details see the code or the class description -// -Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w) -{ - Double_t alpha; - - if (fMatrix) - alpha = (*fMatrix)[fMap]; - else - { - const MHillasSrc *hil = dynamic_cast(par); - if (!par) - { - *fLog << err << dbginf << "MHillasSrc not found... abort." << endl; - return kFALSE; - } - - alpha = hil->GetAlpha(); - } - - fHist.Fill(TMath::Abs(alpha), w); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Update the projections -// -void MHAlpha::Paint(Option_t *opt) -{ - DoFit(); -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHAlpha::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - fHist.Draw(); - - AppendPad(""); -} - -// -------------------------------------------------------------------------- -// -// This is a preliminary implementation of a alpha-fit procedure for -// all possible source positions. It will be moved into its own -// more powerfull class soon. -// -// The fit function is "gaus(0)+pol2(3)" which is equivalent to: -// [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2 -// or -// A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2 -// -// Parameter [1] is fixed to 0 while the alpha peak should be -// symmetric around alpha=0. -// -// Parameter [4] is fixed to 0 because the first derivative at -// alpha=0 should be 0, too. -// -// In a first step the background is fitted between bgmin and bgmax, -// while the parameters [0]=0 and [2]=1 are fixed. -// -// In a second step the signal region (alphaGetBinContent(1); - const Double_t alphaw = h->GetXaxis()->GetBinWidth(1); - - // Check for the regios which is not filled... - if (alpha0==0) - return kFALSE; //*fLog << warn << "Histogram empty." << endl; - - // First fit a polynom in the off region - func.FixParameter(0, 0); - func.FixParameter(2, 1); - func.ReleaseParameter(3); - - for (int i=5; iFit(&func, "N0Q", "", bgmin, bgmax); - - const Double_t chisq1 = func.GetChisquare()/func.GetNDF(); - - // ------------------------------------ - if (paint) - { - func.SetRange(0, 90); - func.SetLineColor(kRed); - func.Paint("same"); - } - // ------------------------------------ - - func.ReleaseParameter(0); - //func.ReleaseParameter(1); - func.ReleaseParameter(2); - func.FixParameter(3, func.GetParameter(3)); - for (int i=5; iFit(&func, "N0Q", "", 0, sigmax); - - const Double_t chisq2 = func.GetChisquare()/func.GetNDF(); - const Double_t sigmagaus = func.GetParameter(2); - - // ------------------------------------ - if (paint) - { - func.SetLineColor(kGreen); - func.Paint("same"); - } - // ------------------------------------ - - const Double_t s = func.Integral(0, sigint)/alphaw; - func.SetParameter(0, 0); - func.SetParameter(2, 1); - const Double_t b = func.Integral(0, sigint)/alphaw; - - sig = MMath::SignificanceLiMaSigned(s, b); - - // ------------------------------------ - // This is always draw one update too late... why? - //fHist.SetTitle(Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ) \\omega=%.1f\\circ E=%d (\\chi_{b}^{2}/N=%.1f \\chi_{s}^{2}/N=%.1f)", - // sig, fSigInt, sigmagaus, (int)(s-b), chisq1, chisq2)); - - if (paint) - { - TLatex text(0.13, 0.94, Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ) \\omega=%.1f\\circ E=%d (\\chi_{b}^{2}/N=%.1f \\chi_{s}^{2}/N=%.1f)", - sig, fSigInt, sigmagaus, (int)(s-b), chisq1, chisq2)); - text.SetBit(TLatex::kTextNDC); - text.SetTextSize(0.04); - text.Paint(); - } - // ------------------------------------ - - return kTRUE; -} - -Bool_t MHAlpha::Finalize() -{ - Double_t sig = 0; - - if (!DoFit(sig)) - { - *fLog << warn << "Histogram empty." << endl; - return kFALSE; - } - - if (fResult) - fResult->SetVal(sig); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// You can use this function if you want to use a MHMatrix instead of -// MMcEvt. This function adds all necessary columns to the -// given matrix. Afterward you should fill the matrix with the corresponding -// data (eg from a file by using MHMatrix::Fill). If you now loop -// through the matrix (eg using MMatrixLoop) MHHadronness::Fill -// will take the values from the matrix instead of the containers. -// -void MHAlpha::InitMapping(MHMatrix *mat) -{ - if (fMatrix) - return; - - fMatrix = mat; - - fMap = fMatrix->AddColumn("MHillasSrc.fAlpha"); -} - -void MHAlpha::StopMapping() -{ - fMatrix = NULL; -} Index: trunk/MagicSoft/Mars/mhist/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHAlpha.h (revision 5065) +++ (revision ) @@ -1,55 +1,0 @@ -#ifndef MARS_MHAlpha -#define MARS_MHAlpha - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH1 -#include -#endif - -class MParList; -class MParameterD; -class MHMatrix; - -class MHAlpha : public MH -{ -private: - Float_t fAlphaCut; // Alpha cut - Float_t fBgMean; // Background mean - - Float_t fSigInt; - Float_t fSigMax; - Float_t fBgMin; - Float_t fBgMax; - Int_t fPolynom; - - TH1D fHist; // Alpha vs. x and y - - MParameterD *fResult; - - MHMatrix *fMatrix; //! - Int_t fMap; //! - - Bool_t DoFit(Double_t &, Bool_t); - Bool_t DoFit(Double_t &d) { return DoFit(d, kFALSE); } - void DoFit() { Double_t d; DoFit(d, kTRUE); } - -public: - MHAlpha(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pl); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - void Paint(Option_t *opt=""); - void Draw(Option_t *option=""); - - void InitMapping(MHMatrix *mat); - void StopMapping(); - - ClassDef(MHAlpha, 1) // Alpha-Plot which is fitted online -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc (revision 5065) +++ (revision ) @@ -1,397 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Wolfgang Wittek 5/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHEffOnTime // -// // -// calculates the effective on time for each bin in the variable Var // -// (Var may be time or Theta) // -// // -// Important : The sample of events used for this should be as close // -// to the sample of recorded events as possible. // -// This means that NO additional cuts (be it quality cuts or // -// gamma/hadron cuts) should be applied (see MAGIC-TDAS 02-02).// -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHEffOnTime.h" - -#include - -#include -#include - -#include -#include -#include - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MHTimeDiffTheta.h" - -ClassImp(MHEffOnTime); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets the name of the variable ("Time" or "Theta") -// and the units of the variable ("[s]" or "[\\circ])") -// -MHEffOnTime::MHEffOnTime(const char *varname, const char *unit) - : fHEffOn(), fHProb(), fHLambda(), fHRdead() -{ - fVarname = varname; - fUnit = unit; - - TString strg = fVarname + " " + fUnit; - - // - // set the name and title of this object - // - fName = TString("MHEffOnTime-")+fVarname; - fTitle = "1-D histogram of Eff On Time"; - - // effective on time versus Var - fHEffOn.SetName("EffOn"); - fHEffOn.SetTitle(TString("T_{on, eff} vs. ")+fVarname); - - fHEffOn.SetDirectory(NULL); - - fHEffOn.SetXTitle(strg); - fHEffOn.SetYTitle("T_{on, eff} [s]"); - - // chi2-probability versus Var - fHProb.SetName("Chi2-prob"); - TString strg3("\\chi^{2}-prob of OnTimeFit vs. "); - strg3 += fVarname; - fHProb.SetTitle(strg3); - - fHProb.SetDirectory(NULL); - - fHProb.SetXTitle(strg); - fHProb.SetYTitle("\\chi^{2}-probability"); - - // lambda versus Var - fHLambda.SetName("lambda"); - fHLambda.SetTitle(TString("\\lambda from OnTimeFit vs. ")+fVarname); - - fHLambda.SetDirectory(NULL); - - fHLambda.SetXTitle(strg); - fHLambda.SetYTitle("\\lambda [Hz]"); - - // Rdead versus Var - // Rdead is the fraction from real time lost by the dead time - fHRdead.SetName("Rdead"); - - fHRdead.SetTitle(TString("Rdead of OnTimeFit vs. ")+fVarname); - - fHRdead.SetDirectory(NULL); - - fHRdead.SetXTitle(strg); - fHRdead.SetYTitle("Rdead"); - -} - -// -------------------------------------------------------------------- -// -// determine range (yq[0], yq[1]) of time differences -// where fit should be performed; -// require a fraction >=min of all entries to lie below yq[0] -// and a fraction <=max of all entries to lie below yq[1]; -// within the range (yq[0], yq[1]) there must be no empty bin; -// -void MHEffOnTime::GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const -{ -#if ROOT_VERSION_CODE < ROOT_VERSION(3,02,07) - // WOrkaround for missing const qualifier of TH1::Integral - TH1 &h1 = (TH1&)h; - - // choose pedestrian approach as long as GetQuantiles is - // not available - const Int_t jbins = h1.GetNbinsX(); - const Double_t Nm = h1.Integral(); - - const Double_t xq[2] = { 0.15*Nm, 0.98*Nm }; - - yq[0] = yq[1] = h1.GetBinLowEdge(jbins+1); - - for (int j=1; j<=jbins; j++) - if (h1.Integral(2, j) >= xq[0]) - { - yq[0] = h1.GetBinLowEdge(j); - break; - } - - for (int j=1; j<=jbins; j++) - if (h1.Integral(1, j) >= xq[1] || h1.GetBinContent(j)==0) - { - yq[1] = h1.GetBinLowEdge(j); - break; - } -#else - // GetQuantiles doesn't seem to be available in root 3.01/06 - Double_t xq[2] = { min, max }; - ((TH1&)h).GetQuantiles(2, yq, xq); -#endif -} - -void MHEffOnTime::DrawBin(TH1 &h, Int_t i) const -{ - TString strg1 = fVarname+"-bin #"; - strg1 += i; - - new TCanvas(strg1, strg1); - - gPad->SetLogy(); - gStyle->SetOptFit(1011); - - TString name="Bin_"; - name += i; - - h.SetName(name); - h.SetXTitle("\\Delta t [s]"); - h.SetYTitle("Counts"); - h.DrawCopy(); -} - -Bool_t MHEffOnTime::CalcResults(const TF1 &func, Double_t Nm, Int_t i) -{ - const Double_t lambda = func.GetParameter(0); - const Double_t N0 = func.GetParameter(1); - const Double_t prob = func.GetProb(); - const Int_t NDF = func.GetNDF(); - - Double_t xmin, xmax; - ((TF1&)func).GetRange(xmin, xmax); - - *fLog << inf; - *fLog << "Fitted bin #" << i << " from " << xmin << " to " << xmax; - *fLog << ", got: lambda=" << lambda << "Hz N0=" << N0 << endl; - - if (prob<=0.01) - { - *fLog << warn << "WARN - Fit bin#" << i << " gives:"; - *fLog << " Chi^2-Probab(" << prob << ")<0.01"; - *fLog << " NoFitPts=" << func.GetNumberFitPoints(); - *fLog << " Chi^2=" << func.GetChisquare(); - } - - // was fit successful ? - if (NDF<=0 || /*prob<=0.001 ||*/ lambda<=0 || N0<=0) - { - *fLog << warn << dbginf << "Fit failed bin #" << i << ": "; - if (NDF<=0) - *fLog << " NDF(Number of Degrees of Freedom)=0"; - if (lambda<=0) - *fLog << " Parameter#0(lambda)=0"; - if (N0<=0) - *fLog << " Parameter#1(N0)=0"; - *fLog << endl; - - return kFALSE; - } - - // - // -------------- start error calculation ---------------- - // - Double_t emat[2][2]; - gMinuit->mnemat(&emat[0][0], 2); - - // - // Rdead : fraction of real time lost by the dead time - // kappa = number of observed events (Nm) divided by - // the number of genuine events (N0) - // Teff : effective on-time - // - const Double_t Teff = Nm/lambda; - const Double_t kappa = Nm/N0; - const Double_t Rdead = 1.0 - kappa; - - const Double_t dldl = emat[0][0]; - const Double_t dN0dN0 = emat[1][1]; - - const Double_t dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm); - const Double_t dl = sqrt(dldl); - const Double_t dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm); - // - // -------------- end error calculation ---------------- - // - - // the effective on time is Nm/lambda - fHEffOn.SetBinContent(i, Teff); - fHEffOn.SetBinError (i, dTeff); - - // plot chi2-probability of fit - fHProb.SetBinContent(i, prob); - - // lambda from fit - fHLambda.SetBinContent(i, lambda); - fHLambda.SetBinError (i, dl); - - // Rdead from fit - fHRdead.SetBinContent(i, Rdead); - fHRdead.SetBinError (i,dRdead); - - return kTRUE; -} - -void MHEffOnTime::ResetBin(Int_t i) -{ - fHEffOn.SetBinContent (i, 1.e-20); - fHProb.SetBinContent (i, 1.e-20); - fHLambda.SetBinContent(i, 1.e-20); - fHRdead.SetBinContent (i, 1.e-20); -} - -// ----------------------------------------------------------------------- -// -// Calculate the effective on time by fitting the distribution of -// time differences -// -void MHEffOnTime::Calc(const TH2D *hist, const Bool_t draw) -{ - // nbins = number of Var bins - const Int_t nbins = hist->GetNbinsY(); - - for (int i=1; i<=nbins; i++) - { - TH1D &h = *((TH2D*)hist)->ProjectionX(TString("Calc-")+fVarname, - i, i, "E"); - if (draw) - DrawBin(h, i); - - ResetBin(i); - - // Nmdel = Nm * binwidth, with Nm = number of observed events - const Double_t Nm = h.Integral(); - const Double_t Nmdel = h.Integral("width"); - if (Nm <= 0) - { - *fLog << warn << dbginf << "Nm<0 for bin #" << i << endl; - delete &h; - continue; - } - - Double_t yq[2]; - GetQuantiles(h, 0.15, 0.95, yq); - - // - // Setup Poisson function for the fit: - // lambda [Hz], N0 = ideal no of evts, del = bin width of dt - // - TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]); - func.SetParNames("lambda", "N0", "del"); - - func.SetParameter(0, 1); - func.SetParameter(1, Nm); - func.FixParameter(2, Nmdel/Nm); - - // For me (TB) it seems that setting parameter limits isn't necessary - - // For a description of the options see TH1::Fit - h.Fit("Poisson", "0IRQ"); - - // Calc and fill results of fit into the histograms - const Bool_t rc = CalcResults(func, Nm, i); - - // draw the distribution of time differences if requested - if (rc && draw) - func.DrawCopy("same"); - - delete &h; - } -} - -// ------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histograms -// -Bool_t MHEffOnTime::SetupFill(const MParList *plist) -{ - TString bn = "Binning"+fVarname; - - const MBinning* binsVar = (MBinning*)plist->FindObject(bn); - - if (!binsVar) - { - *fLog << err << dbginf << bn << " [MBinning] not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHEffOn, binsVar); - SetBinning(&fHProb, binsVar); - SetBinning(&fHLambda, binsVar); - SetBinning(&fHRdead, binsVar); - - fHEffOn.Sumw2(); - fHProb.Sumw2(); - fHLambda.Sumw2(); - fHRdead.Sumw2(); - - return kTRUE; -} - -// ------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHEffOnTime::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - pad->Divide(2,2); - - pad->cd(1); - gPad->SetBorderMode(0); - fHEffOn.Draw(opt); - - pad->cd(2); - gPad->SetBorderMode(0); - fHProb.Draw(opt); - - pad->cd(3); - gPad->SetBorderMode(0); - fHLambda.Draw(opt); - - pad->cd(4); - gPad->SetBorderMode(0); - fHRdead.Draw(opt); - - pad->Modified(); - pad->Update(); -} - -void MHEffOnTime::Calc(const MHTimeDiffTheta &hist, const Bool_t Draw) -{ - Calc(hist.GetHist(), Draw); -} - Index: trunk/MagicSoft/Mars/mhist/MHEffOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTime.h (revision 5065) +++ (revision ) @@ -1,48 +1,0 @@ -#ifndef MARS_MHEffOnTime -#define MARS_MHEffOnTime - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TH1 -#include -#endif - -class TH2D; -class MHTimeDiffTheta; -class MParList; - -class MHEffOnTime : public MH -{ -private: - TH1D fHEffOn; - TH1D fHProb; - TH1D fHLambda; - TH1D fHRdead; - - TString fVarname; - TString fUnit; - - void ResetBin(Int_t i); - Bool_t CalcResults(const TF1 &func, Double_t Nm, Int_t i); - void DrawBin(TH1 &h, Int_t i) const; - void GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const; - -public: - MHEffOnTime(const char *varname=NULL, const char *unit=NULL); - - Bool_t SetupFill(const MParList *pList); - - const TH1D *GetHist() { return &fHEffOn; } - const TH1D *GetHist() const { return &fHEffOn; } - - void Calc(const TH2D *hist, const Bool_t Draw); - void Calc(const MHTimeDiffTheta &hist, const Bool_t Draw); - - void Draw(Option_t *option=""); - - ClassDef(MHEffOnTime, 0) //1D-plot of Delta t vs. Var -}; - -#endif - Index: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.cc (revision 5065) +++ (revision ) @@ -1,286 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHEffOnTimeTheta // -// // -// calculates the effective on time for each bin in Theta // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHEffOnTimeTheta.h" - -#include - -#include -#include -#include - -#include "MTime.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHEffOnTimeTheta); - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histograms. -// -MHEffOnTimeTheta::MHEffOnTimeTheta(const char *name, const char *title) - : fHEffOn() -{ - // - // set the name and title of this object - // - fName = name ? name : "MHEffOnTimeTheta"; - fTitle = title ? title : "1-D histogram of Eff On Time"; - - // 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("chi^{2}/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; i<=nbins; i++) - { - // 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; - } -} - -// ------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histograms -// -Bool_t MHEffOnTimeTheta::SetupFill(const MParList *plist) -{ - const MBinning* bins = (MBinning*)plist->FindObject("BinningTheta"); - if (!bins) - { - *fLog << err << dbginf << "BinningTheta [MBinning] not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHEffOn, bins); - SetBinning(&fHChi2, bins); - SetBinning(&fHLambda, bins); - SetBinning(&fHN0del, bins); - - fHEffOn.Sumw2(); - fHChi2.Sumw2(); - fHLambda.Sumw2(); - fHN0del.Sumw2(); - - return kTRUE; -} - -// ------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHEffOnTimeTheta::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - pad->Divide(2,2); - - pad->cd(1); - gPad->SetBorderMode(0); - fHEffOn.Draw(opt); - - pad->cd(2); - gPad->SetBorderMode(0); - fHChi2.Draw(opt); - - pad->cd(3); - gPad->SetBorderMode(0); - fHLambda.Draw(opt); - - pad->cd(4); - gPad->SetBorderMode(0); - fHN0del.Draw(opt); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h (revision 5065) +++ (revision ) @@ -1,40 +1,0 @@ -#ifndef MARS_MHEffOnTimeTheta -#define MARS_MHEffOnTimeTheta - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TH1 -#include -#endif - -class MTime; -class TH2D; -class MParList; - -class MHEffOnTimeTheta : public MH -{ -private: - TH1D fHEffOn; - TH1D fHChi2; - TH1D fHN0del; - TH1D fHLambda; - -public: - MHEffOnTimeTheta(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - - const TH1D *GetHist() { return &fHEffOn; } - const TH1D *GetHist() const { return &fHEffOn; } - - void Calc(TH2D *hist); - - void Draw(Option_t *option=""); - - ClassDef(MHEffOnTimeTheta, 0) //1D-plot of Delta t vs. Theta -}; - -#endif - - Index: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.cc (revision 5065) +++ (revision ) @@ -1,298 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MHEffOnTimeTime // -// // -// calculates the effective on time for each bin in time // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MHEffOnTimeTime.h" - -#include - -#include -#include -#include - -#include "MTime.h" - -#include "MBinning.h" -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHEffOnTimeTime); - - -// -------------------------------------------------------------------------- -// -// Default Constructor. It sets name and title of the histograms. -// -MHEffOnTimeTime::MHEffOnTimeTime(const char *name, const char *title) - : fHEffOn() -{ - // - // set the name and title of this object - // - fName = name ? name : "MHEffOnTimeTime"; - fTitle = title ? title : "1-D histogram of Eff On Time"; - - // effective on time versus time - fHEffOn.SetName("EffOn"); - fHEffOn.SetTitle("Effective On Time vs. Time"); - - fHEffOn.SetDirectory(NULL); - - fHEffOn.SetXTitle("time [s]"); - fHEffOn.SetYTitle("t_{eff} [s]"); - - // chi2/NDF versus time - fHChi2.SetName("Chi2/NDF"); - fHChi2.SetTitle("Chi2/NDF of OnTimeFit vs. Time"); - - fHChi2.SetDirectory(NULL); - - fHChi2.SetXTitle("time [s]"); - fHChi2.SetYTitle("chi^{2}/NDF"); - - // lambda versus time - fHLambda.SetName("lambda"); - fHLambda.SetTitle("lambda of OnTimeFit vs. Time"); - - fHLambda.SetDirectory(NULL); - - fHLambda.SetXTitle("time [s]"); - fHLambda.SetYTitle("\\lambda [Hz]"); - - // N0del versus time - fHN0del.SetName("N0del"); - fHN0del.SetTitle("N0del of OnTimeFit vs. Time"); - - fHN0del.SetDirectory(NULL); - - fHN0del.SetXTitle("time [s]"); - fHN0del.SetYTitle("N0del"); -} - -// ----------------------------------------------------------------------- -// -// Calculate the effective on time by fitting the distribution of -// time differences -// -void MHEffOnTimeTime::Calc(TH2D *hist) -{ - // nbins = number of time bins - const Int_t nbins = hist->GetNbinsY(); - - for (int i=1; i<=nbins; i++) - { - - // TH1D &h = *hist->ProjectionX("Calc-time", i, i, "E"); - TH1D &h = *hist->ProjectionX("Calc-time", 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) - { - // char txt[100]; - // sprintf(txt, "time_bin:%d", i); - // new TCanvas(txt, txt); - - 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++) - { - const 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 (=zero) 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; - } -} - -// ------------------------------------------------------------------------- -// -// Set the binnings and prepare the filling of the histograms -// -Bool_t MHEffOnTimeTime::SetupFill(const MParList *plist) -{ - const MBinning* binstime = (MBinning*)plist->FindObject("BinningTime"); - if (!binstime) - { - *fLog << err << dbginf << "BinningTime [MBinning] not found... aborting." << endl; - return kFALSE; - } - - SetBinning(&fHEffOn, binstime); - SetBinning(&fHChi2, binstime); - SetBinning(&fHLambda, binstime); - SetBinning(&fHN0del, binstime); - - fHEffOn.Sumw2(); - fHChi2.Sumw2(); - fHLambda.Sumw2(); - fHN0del.Sumw2(); - - return kTRUE; -} - -// ------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHEffOnTimeTime::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - pad->Divide(2,2); - - pad->cd(1); - gPad->SetBorderMode(0); - fHEffOn.Draw(opt); - - pad->cd(2); - gPad->SetBorderMode(0); - fHChi2.Draw(opt); - - pad->cd(3); - gPad->SetBorderMode(0); - fHLambda.Draw(opt); - - pad->cd(4); - gPad->SetBorderMode(0); - fHN0del.Draw(opt); - - pad->Modified(); - pad->Update(); -} - - - - - Index: trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h (revision 5065) +++ (revision ) @@ -1,39 +1,0 @@ -#ifndef MARS_MHEffOnTimeTime -#define MARS_MHEffOnTimeTime - -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef ROOT_TH1 -#include -#endif - -class MTime; -class TH2D; -class MParList; - -class MHEffOnTimeTime : public MH -{ -private: - TH1D fHEffOn; - TH1D fHChi2; - TH1D fHLambda; - TH1D fHN0del; - -public: - MHEffOnTimeTime(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - - const TH1D *GetHist() { return &fHEffOn; } - const TH1D *GetHist() const { return &fHEffOn; } - - void Calc(TH2D *hist); - - void Draw(Option_t *option=""); - - ClassDef(MHEffOnTimeTime, 0) //1D-plot of Delta t vs. time -}; - -#endif - Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 5065) +++ (revision ) @@ -1,825 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -// -// Filling this you will get the effective on-time versus theta and -// observation time. -// -// 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 -// -// To determin the efective on time a poisson fit is done. For more details -// please have a look into the source code of FitH() it should be simple -// to understand. In this function a Delta-T distribution is fitted, while -// Delta-T is the time between two consecutive events. -// -// The fit is done for projections of a 2D histogram in Theta and Delta T. -// So you get the effective on time versus theta. -// -// To get the effective on-time versus time a histogram is filled with -// the Delta-T distribution of a number of events set by SetNumEvents(). -// The default is 12000 (roughly 1min at 200Hz) -// -// For each "time-bin" the histogram is fitted and the resulting effective -// 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. -// -// In addition the calculated effective on time is stored in a -// "MEffectiveOnTime [MParameterDerr]" and the corresponding time-stamp -// (the upper edge of the bin) "MTimeEffectiveOnTime [MTime]" -// -// The class takes two binnings from the Parameter list; if these binnings -// are not available the defaultbinning is used: -// MBinning("BinningDeltaT"); // Units of seconds -// MBinning("BinningTheta"); // Units of degrees -// -// -// Usage: -// ------ -// MFillH fill("MHEffectiveOnTime", "MTime"); -// tlist.AddToList(&fill); -// -// -// Input Container: -// MPointingPos -// -// Output Container: -// MEffectiveOnTime [MParameterDerr] -// MTimeEffectiveOnTime [MTime] -// -////////////////////////////////////////////////////////////////////////////// -#include "MHEffectiveOnTime.h" - -#include -#include -#include - -#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 initializes all histograms. -// -MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title) - : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), - fNumEvents(200*60), fNameProjDeltaT(Form("DeltaT_%p", this)), - fNameProjTheta(Form("Theta_%p", this)) -{ - // - // set the name and title of this object - // - fName = name ? name : "MHEffectiveOnTime"; - fTitle = title ? title : "Histogram to determin effective On-Time vs Time and Zenith Angle"; - - // Main histogram - fH2DeltaT.SetName("DeltaT"); - fH2DeltaT.SetXTitle("\\Delta t [s]"); - fH2DeltaT.SetYTitle("\\Theta [\\circ]"); - fH2DeltaT.SetZTitle("Count"); - fH2DeltaT.UseCurrentStyle(); - fH2DeltaT.SetDirectory(NULL); - - // Main histogram - fH1DeltaT.SetName("DeltaT"); - fH1DeltaT.SetXTitle("\\Delta t [s]"); - fH1DeltaT.SetYTitle("Counts"); - fH1DeltaT.UseCurrentStyle(); - fH1DeltaT.SetDirectory(NULL); - - // effective on time versus theta - 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 - 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 - 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 - 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 - 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 - /* - 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"); - btheta.SetEdgesCos(100, 0, 60); - - MBinning btime("BinningDeltaT"); - btime.SetEdges(50, 0, 0.1); - - MH::SetBinning(&fH2DeltaT, &btime, &btheta); - - btime.Apply(fH1DeltaT); - - btheta.Apply(fHThetaEffOn); - btheta.Apply(fHThetaLambda); - btheta.Apply(fHThetaNDF); - btheta.Apply(fHThetaProb); - //btheta.Apply(fHChi2); -} - -// -------------------------------------------------------------------------- -// -// 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("BinningDeltaT"); - const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); - if (binsdtime) - binsdtime->Apply(fH1DeltaT); - if (binstheta) - { - binstheta->Apply(fHThetaEffOn); - binstheta->Apply(fHThetaLambda); - binstheta->Apply(fHThetaNDF); - binstheta->Apply(fHThetaProb); - //binstheta->Apply(fHChi2); - } - if (binstheta && binsdtime) - SetBinning(&fH2DeltaT, binsdtime, binstheta); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fit a single Delta-T distribution. See source code for more details -// -Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res, Bool_t paint) 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 - // - static TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)"); - //func.SetParNames("lambda", "N0", "del"); - - func.SetParameter(0, 100); // Hz - func.SetParameter(1, Nm); - func.FixParameter(2, Nmdel/Nm); - - // 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, "NIQ", "", yq[0], yq[1]); - - const Double_t chi2 = func.GetChisquare(); - const Int_t NDF = func.GetNDF(); - - // was fit successful ? - const Bool_t ok = NDF>0 && chi2<2.5*NDF; - - if (paint) - { - func.SetLineWidth(2); - func.SetLineColor(ok ? kGreen : kRed); - func.Paint("same"); - } - - if (!ok) - 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[3] = lambda; - res[4] = dl; - - // NDoF of fit - res[5] = NDF; - - // Rdead (from fit) is the fraction from real time lost by the dead time - //fHRdead.SetBinContent(i, Rdead); - //fHRdead.SetBinError (i,dRdead); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fit a all bins of the distribution in theta. Store the result in the -// Theta-Histograms -// -void MHEffectiveOnTime::FitThetaBins() -{ - fHThetaEffOn.Reset(); - fHThetaProb.Reset(); - fHThetaLambda.Reset(); - fHThetaNDF.Reset(); - - const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999)); - - // nbins = number of Theta bins - const Int_t nbins = fH2DeltaT.GetNbinsY(); - - TH1D *h=0; - for (int i=1; i<=nbins; i++) - { - // TH1D &h = *hist->ProjectionX("Calc-theta", i, i); - h = fH2DeltaT.ProjectionX(name, i, i, "E"); - - Double_t res[6]; - if (!FitH(h, res)) - continue; - - // the effective on time is Nm/lambda - fHThetaEffOn.SetBinContent(i, res[0]); - fHThetaEffOn.SetBinError (i, res[1]); - - // plot chi2-probability of fit - fHThetaProb.SetBinContent(i, res[2]); - - // plot chi2/NDF of fit - //fHChi2.SetBinContent(i, res[3]); - - // lambda of fit - 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 - //fHRdead.SetBinContent(i, Rdead); - //fHRdead.SetBinError (i,dRdead); - } - - // Histogram is reused via gROOT->FindObject() - // Need to be deleted only once - if (h) - delete h; -} - -// -------------------------------------------------------------------------- -// -// Fit the single-time-bin histogram. Store the result in the -// Time-Histograms -// -void MHEffectiveOnTime::FitTimeBin() -{ - // - // Fit histogram - // - Double_t res[6]; - if (!FitH(&fH1DeltaT, res)) - return; - - // Reset Histogram - fH1DeltaT.Reset(); - - // - // Prepare Histogram - // - - // Get number of bins - const Int_t n = fHTimeEffOn.GetNbinsX(); - - // Enhance binning - MBinning bins; - bins.SetEdges(fHTimeEffOn, 'x'); - bins.AddEdge(fLastTime.GetAxisTime()); - bins.Apply(fHTimeEffOn); - bins.Apply(fHTimeProb); - bins.Apply(fHTimeLambda); - //bins.Apply(fHTimeNDF); - - // - // Fill histogram - // - 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]); - - // - // Now prepare output - // - fParam->SetVal(res[0], res[1]); - fParam->SetReadyToSave(); - - *fTime = fLastTime; - - // Include the current event - fTime->Plus1ns(); - - *fLog << fLastTime << ": Val=" << res[0] << " Err=" << res[1] << endl; -} - -// -------------------------------------------------------------------------- -// -// 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 - MHEffectiveOnTime::Fill without argument or container doesn't inherit from MTime... abort." << endl; - return kFALSE; - } - - // - // If this is the first call we have to initialize the time-histogram - // - if (fLastTime==MTime()) - { - MBinning bins; - bins.SetEdges(1, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime()); - bins.Apply(fHTimeEffOn); - bins.Apply(fHTimeProb); - bins.Apply(fHTimeLambda); - //bins.Apply(fHTimeNDF); - - fParam->SetVal(0, 0); - fParam->SetReadyToSave(); - - *fTime = *time; - - // Make this 1ns before the first event! - fTime->Minus1ns(); - } - - // - // Fill time difference into the histograms - // - const Double_t dt = *time-fLastTime; - - fH2DeltaT.Fill(dt, fPointPos->GetZd(), w); - fH1DeltaT.Fill(dt, w); - - fLastTime = *time; - - // - // If we reached the event number limit for the time-bins fit the histogram - // - if (fH1DeltaT.GetEntries()>=fNumEvents) - FitTimeBin(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fit the theta projections of the 2D histogram and the 1D Delta-T -// distribution -// -Bool_t MHEffectiveOnTime::Finalize() -{ - FitThetaBins(); - FitTimeBin(); - MH::RemoveFirstBin(fHTimeEffOn); - MH::RemoveFirstBin(fHTimeProb); - MH::RemoveFirstBin(fHTimeLambda); - //MH::RemoveFirstBin(fHTimeNDF); - - fIsFinalized = kTRUE; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Paint the integral and the error on top of the histogram -// -void MHEffectiveOnTime::PaintText(Double_t val, Double_t error) const -{ - TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", val, error)); - text.SetBit(TLatex::kTextNDC); - 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); -} - -// -------------------------------------------------------------------------- -// -// Prepare painting the histograms -// -void MHEffectiveOnTime::Paint(Option_t *opt) -{ - TH1D *h=0; - TPaveStats *st=0; - - TString o(opt); - if (o==(TString)"fit") - { - TVirtualPad *pad = gPad; - - for (int x=0; x<2; x++) - for (int y=0; y<3; y++) - { - TVirtualPad *p=gPad->GetPad(x+1)->GetPad(y+1); - if (!(st = (TPaveStats*)p->GetPrimitive("stats"))) - continue; - - if (st->GetOptStat()==11) - continue; - - const Double_t y1 = st->GetY1NDC(); - const Double_t y2 = st->GetY2NDC(); - const Double_t x1 = st->GetX1NDC(); - const Double_t x2 = st->GetX2NDC(); - - st->SetY1NDC((y2-y1)/3+y1); - st->SetX1NDC((x2-x1)/3+x1); - st->SetOptStat(11); - } - - pad->GetPad(1)->cd(1); - if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT))) - { - h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E"); - if (h->GetEntries()>0) - gPad->SetLogy(); - } - - pad->GetPad(2)->cd(1); - if ((h = (TH1D*)gPad->FindObject(fNameProjTheta))) - fH2DeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E"); - - if (!fIsFinalized) - FitThetaBins(); - return; - } - if (o==(TString)"paint") - { - if ((h = (TH1D*)gPad->FindObject(fNameProjDeltaT))) - { - Double_t res[6]; - if (FitH(h, res, kTRUE)) - PaintText(res[0], res[1]); - } - return; - } - - 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 = &fHThetaEffOn; - UpdateRightAxis(fHThetaLambda); - } - if (o==(TString)"time") - { - h = &fHTimeEffOn; - UpdateRightAxis(fHTimeLambda); - } - - if (!h) - return; - - Double_t error = 0; - for (int i=0; iGetXaxis()->GetNbins(); i++) - error += h->GetBinError(i); - - 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(); -} - -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHEffectiveOnTime::Draw(Option_t *opt) -{ - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad("fit"); - - pad->Divide(2, 1, 0, 0); - - TH1 *h; - - pad->cd(1); - gPad->SetBorderMode(0); - gPad->Divide(1, 3, 0, 0); - pad->GetPad(1)->cd(1); - gPad->SetBorderMode(0); - h = fH2DeltaT.ProjectionX(fNameProjDeltaT, -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(); - AppendPad("paint"); - - pad->GetPad(1)->cd(2); - gPad->SetBorderMode(0); - fHTimeProb.Draw(); - AppendPad("timendf"); - //fHTimeNDF.Draw("same"); - //DrawRightAxis("NDF"); - - pad->GetPad(1)->cd(3); - gPad->SetBorderMode(0); - fHTimeEffOn.Draw(); - AppendPad("time"); - fHTimeLambda.Draw("same"); - DrawRightAxis("Slope"); - - pad->cd(2); - gPad->SetBorderMode(0); - gPad->Divide(1, 3, 0, 0); - - pad->GetPad(2)->cd(1); - gPad->SetBorderMode(0); - h = fH2DeltaT.ProjectionY(fNameProjTheta, -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->GetPad(2)->cd(2); - gPad->SetBorderMode(0); - fHThetaProb.Draw(); - AppendPad("thetandf"); - fHThetaNDF.Draw("same"); - DrawRightAxis("NDF"); - - pad->GetPad(2)->cd(3); - gPad->SetBorderMode(0); - fHThetaEffOn.Draw(); - AppendPad("theta"); - fHThetaLambda.Draw("same"); - DrawRightAxis("Slope"); -} Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 5065) +++ (revision ) @@ -1,82 +1,0 @@ -#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 -#ifndef ROOT_TF1 -#include -#endif - -class MTime; -class MPointingPos; -class MParameterDerr; - -class MParList; - -class MHEffectiveOnTime : public MH -{ -private: - MPointingPos *fPointPos; //! Container to get the zenith angle from - MTime fLastTime; //! Time-Stamp of last event - - MTime *fTime; //! Time-Stamp of "effective on-time" event - MParameterDerr *fParam; //! Output container for effective on-time - - TH2D fH2DeltaT; // Counts vs Delta T and Theta - TH1D fH1DeltaT; //! Counts vs Delta T (for a time interval) - - 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 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 - - 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 - - const TString fNameProjDeltaT; //! This should make sure, that gROOT doen't confuse the projection with something else - const TString fNameProjTheta; //! This should make sure, that gROOT doen't confuse the projection with something else - - Bool_t FitH(TH1D *h, Double_t *res, Bool_t paint=kFALSE) const; - 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: - MHEffectiveOnTime(const char *name=NULL, const char *title=NULL); - - void SetNumEvents(Int_t i) { fNumEvents=i; } - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - const TH1D &GetHEffOnTheta() const { return fHThetaEffOn; } - const TH1D &GetHEffOnTime() const { return fHTimeEffOn; } - - void Draw(Option_t *option=""); - void Paint(Option_t *opt=""); - - ClassDef(MHEffectiveOnTime, 1) // Histogram to determin effective On-Time vs Time and Zenith Angle -}; - -#endif - Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 5065) +++ (revision ) @@ -1,88 +1,0 @@ -#ifndef MARS_MHFalseSource -#define MARS_MHFalseSource - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH3 -#include -#endif - -class TH2D; - -class MParList; -class MTime; -class MSrcPosCam; -class MPointingPos; -class MObservatory; - -class MHFalseSource : public MH -{ -private: - MTime *fTime; //! container to take the event time from - MPointingPos *fPointPos; //! container to take pointing position from - MSrcPosCam *fSrcPos; //! container for sopurce position in camera - MObservatory *fObservatory; //! conteiner to take observatory location from - - Float_t fMm2Deg; // conversion factor for display in degrees - - Float_t fAlphaCut; // Alpha cut - Float_t fBgMean; // Background mean - - Float_t fMinDist; // Min dist - Float_t fMaxDist; // Max dist - - Float_t fMinLD; // Minimum distance in percent of dist - Float_t fMaxLD; // Maximum distance in percent of dist - - TH3D fHist; // Alpha vs. x and y - - Double_t fRa; - Double_t fDec; - - Int_t DistancetoPrimitive(Int_t px, Int_t py); - void Modified(); - - void ProjectAll(TH2D *h); - void ProjectOff(TH2D *h, TH2D *all); - void ProjectOn(TH2D *h, TH2D *all); - - TObject *GetCatalog(); - - void MakeSymmetric(TH1 *h); - -public: - MHFalseSource(const char *name=NULL, const char *title=NULL); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - - TH1 *GetHistByName(const TString name) { return &fHist; } - - void FitSignificance(Float_t sigint=15, Float_t sigmax=70, Float_t bgmin=40, Float_t bgmax=70, Byte_t polynom=1); //*MENU* - void FitSignificanceStd() { FitSignificance(); } //*MENU* - - void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance - void SetMaxDist(Float_t dist) { fMaxDist = dist; } // Absolute maximum distance - void SetMinLD(Float_t ratio) { fMinLD = ratio; } // Minimum ratio between length/dist - void SetMaxLD(Float_t ratio) { fMaxLD = ratio; } // Maximum ratio between length/dist - - void SetAlphaCut(Float_t alpha); //*MENU* - void SetAlphaPlus5() { SetAlphaCut(fAlphaCut+5); } //*MENU* - void SetAlphaMinus5() { SetAlphaCut(fAlphaCut-5); } //*MENU* - - void SetBgMean(Float_t alpha); //*MENU* - void SetBgMeanPlus5() { SetBgMean(fBgMean+5); } //*MENU* - void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU* - - void Paint(Option_t *opt=""); - void Draw(Option_t *option=""); - - static Double_t Significance(Double_t s, Double_t b); - static Double_t SignificanceLiMa(Double_t s, Double_t b, Double_t alpha=1); - - ClassDef(MHFalseSource, 1) //3D-histogram in alpha, x and y -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 5065) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 5066) @@ -34,7 +34,4 @@ MHAlphaEnergyTime.cc \ MHAlphaEnergyTheta.cc \ - MHEffOnTime.cc \ - MHTimeDiffTime.cc \ - MHTimeDiffTheta.cc \ MHCompProb.cc \ MHHadronness.cc \