Index: trunk/MagicSoft/Mars/mhflux/FluxIncl.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/FluxIncl.h (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/FluxIncl.h (revision 4999) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 4999) @@ -0,0 +1,13 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MAlphaFitter+; + +#pragma link C++ class MHAlpha+; +#pragma link C++ class MHFalseSource+; +#pragma link C++ class MHEffectiveOnTime+; + +#endif Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 4999) @@ -0,0 +1,756 @@ +/* ======================================================================== *\ +! +! * +! * 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 + +#include "MGeomCam.h" +#include "MSrcPosCam.h" +#include "MHillasSrc.h" +#include "MEnergyEst.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); +ClassImp(MAlphaFitter); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MHAlpha::MHAlpha(const char *name, const char *title) + : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0), + fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHAlpha"; + fTitle = title ? title : "Alpha plot"; + + fHAlpha.SetName("Alpha"); + fHAlpha.SetTitle("Alpha"); + fHAlpha.SetXTitle("\\Theta [deg]"); + fHAlpha.SetYTitle("E_{est} [GeV]"); + fHAlpha.SetZTitle("|\\alpha| [\\circ]"); + fHAlpha.SetDirectory(NULL); + fHAlpha.UseCurrentStyle(); + + // Main histogram + fHAlphaTime.SetName("Alpha"); + fHAlphaTime.SetXTitle("|\\alpha| [\\circ]"); + fHAlphaTime.SetYTitle("Counts"); + fHAlphaTime.UseCurrentStyle(); + fHAlphaTime.SetDirectory(NULL); + + + fHEnergy.SetName("Energy"); + fHEnergy.SetTitle(" N_{exc} vs. E_{est} "); + fHEnergy.SetXTitle("E_{est} [GeV]"); + fHEnergy.SetYTitle("N_{exc}"); + fHEnergy.SetDirectory(NULL); + fHEnergy.UseCurrentStyle(); + + fHTheta.SetName("Theta"); + fHTheta.SetTitle(" N_{exc} vs. \\Theta "); + fHTheta.SetXTitle("\\Theta [\\circ]"); + fHTheta.SetYTitle("N_{exc}"); + fHTheta.SetDirectory(NULL); + fHTheta.UseCurrentStyle(); + + // effective on time versus time + fHTime.SetName("Time"); + fHTime.SetTitle(" N_{exc} vs. Time "); + fHTime.SetXTitle("Time"); + fHTime.SetYTitle("N_{exc} [s]"); + fHTime.UseCurrentStyle(); + fHTime.SetDirectory(NULL); + fHTime.GetYaxis()->SetTitleOffset(1.2); + fHTime.GetXaxis()->SetLabelSize(0.033); + fHTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fHTime.GetXaxis()->SetTimeDisplay(1); + fHTime.Sumw2(); + + MBinning binsa, binse, binst; + binsa.SetEdges(18, 0, 90); + binse.SetEdgesLog(25, 10, 100000); + binst.SetEdgesCos(50, 0, 60); + binse.Apply(fHEnergy); + binst.Apply(fHTheta); + binsa.Apply(fHAlphaTime); + + MH::SetBinning(&fHAlpha, &binst, &binse, &binsa); +} + +void MHAlpha::FitEnergySpec(Bool_t paint) +{ + TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]"); + f.SetParameter(0, -2.0); + f.SetParameter(1, 500); // 50 + f.SetParameter(2, 1); // 50 + f.SetLineWidth(2); + f.SetLineColor(kGreen); + + fHEnergy.Fit(&f, "0QI", "", 90, 1900); // Integral? + + if (paint) + { + f.Paint("same"); + + TLatex text(0.2, 0.94, Form("\\alpha=%.1f E_{0}=%.1fGeV N=%.1f", + f.GetParameter(0), + f.GetParameter(1), + 0/*f.GetParameter(2)*/)); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); + } +} + +void MHAlpha::FitEnergyBins(Bool_t paint) +{ + // Do not store the 'final' result in fFit + MAlphaFitter fit(fFit); + + const Int_t n = fHAlpha.GetNbinsY(); + + for (int i=1; i<=n; i++) + { + TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":""); + if (fit.Fit(*h)) + { + fHEnergy.SetBinContent(i, fit.GetExcessEvents()); + fHEnergy.SetBinError(i, fit.GetExcessEvents()*0.2); + } + delete h; + } +} + +void MHAlpha::FitThetaBins(Bool_t paint) +{ + // Do not store the 'final' result in fFit + MAlphaFitter fit(fFit); + + const Int_t n = fHAlpha.GetNbinsX(); + + for (int i=1; i<=n; i++) + { + TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":""); + if (fit.Fit(*h)) + { + fHTheta.SetBinContent(i, fit.GetExcessEvents()); + fHTheta.SetBinError(i, fit.GetExcessEvents()*0.2); + } + delete h; + } +} + +Bool_t MHAlpha::SetupFill(const MParList *pl) +{ + fHAlpha.Reset(); + fHEnergy.Reset(); + fHTheta.Reset(); + + fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst"); + if (!fEnergy) + *fLog << warn << "MEnergyEst not found... ignored." << endl; + + fPointPos = (MPointingPos*)pl->FindObject("MPointingPos"); + if (!fPointPos) + *fLog << warn << "MPointingPos not found... ignored." << endl; + + fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime"); + if (!fTimeEffOn) + *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl; + + fTime = (MTime*)pl->FindObject("MTime"); + if (!fTime) + *fLog << warn << "MTime not found... ignored." << endl; + + fResult = (MParameterD*)const_cast(pl)->FindCreateObj("MParameterD", "Significance"); + if (!fResult) + return kFALSE; + + //fExcess = (MParameterD*)const_cast(pl)->FindCreateObj("MParameterD", "MExcess"); + //if (!fExcess) + // return kFALSE; + + fLastTime = MTime(); + + MBinning binst, binse, binsa; + binst.SetEdges(fHAlpha, 'x'); + binse.SetEdges(fHAlpha, 'y'); + binsa.SetEdges(fHAlpha, 'z'); + + MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning"); + if (fPointPos && bins) + binst.SetEdges(*bins); + if (!fPointPos) + binst.SetEdges(1, 0, 90); + + bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning"); + if (fEnergy && bins) + binse.SetEdges(*bins); + if (!fEnergy) + binse.SetEdges(1, 10, 100000); + + bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning"); + if (bins) + binsa.SetEdges(*bins); + + binse.Apply(fHEnergy); + binst.Apply(fHTheta); + binsa.Apply(fHAlphaTime); + MH::SetBinning(&fHAlpha, &binst, &binse, &binsa); + + MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter"); + if (!fit) + *fLog << warn << "MAlphaFitter not found... ignored." << endl; + else + fit->Copy(fFit); + + return kTRUE; +} + +void MHAlpha::InitAlphaTime(const MTime &t) +{ + // + // If this is the first call we have to initialize the time-histogram + // + MBinning bins; + bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime()); + bins.Apply(fHTime); + + fLastTime=t; +} + +void MHAlpha::UpdateAlphaTime(Bool_t final) +{ + if (!fTimeEffOn) + return; + + const Int_t steps = 6; + + static int rebin = steps; + + if (!final) + { + if (fLastTime==MTime() && *fTimeEffOn!=MTime()) + { + InitAlphaTime(*fTimeEffOn); + return; + } + + if (fLastTime!=*fTimeEffOn) + { + fLastTime=*fTimeEffOn; + rebin--; + } + + if (rebin>0) + return; + } + + MAlphaFitter fit(fFit); + if (!fit.Fit(fHAlphaTime)) + return; + + // Reset Histogram + fHAlphaTime.Reset(); + + // Get number of bins + const Int_t n = fHTime.GetNbinsX(); + + // + // Prepare Histogram + // + + MTime *time = final ? fTime : fTimeEffOn; + if (final) + time->Plus1ns(); + + // Enhance binning + MBinning bins; + bins.SetEdges(fHTime, 'x'); + bins.AddEdge(time->GetAxisTime()); + bins.Apply(fHTime); + + // + // Fill histogram + // + fHTime.SetBinContent(n+1, fit.GetExcessEvents()); + fHTime.SetBinError(n+1, fit.GetExcessEvents()*0.1); + + *fLog << all << *fTimeEffOn << ": " << fit.GetExcessEvents() << endl; + + rebin = steps; +} + +// -------------------------------------------------------------------------- +// +// 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, energy, theta; + + if (fMatrix) + { + alpha = (*fMatrix)[fMap[0]]; + energy = 1000; //(*fMatrix)[fMap[1]]; + theta = 0; //(*fMatrix)[fMap[2]]; + } + else + { + const MHillasSrc *hil = dynamic_cast(par); + if (!par) + { + *fLog << err << dbginf << "MHillasSrc not found... abort." << endl; + return kFALSE; + } + + alpha = hil->GetAlpha(); + energy = fEnergy ? fEnergy->GetEnergy() : 1000; + theta = fPointPos ? fPointPos->GetZd() : 0; + } + + UpdateAlphaTime(); + + fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w); + fHAlphaTime.Fill(TMath::Abs(alpha), w); + + /* + // N_gamma vs Energy and Theta + const Double_t Ngam = fHUnfold.GetBinContent(m,n); + // A_eff vs Energy and Theta + const Double_t Aeff = fHAeff.GetBinContent(m,n); + // T_eff vs Theta + const Double_t Effon = teff.GetBinContent(n); + + const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam; + const Double_t c2 = teff.GetBinError(n) /Effon; + const Double_t c3 = fHAeff.GetBinError(m,n) /Aeff; + + const Double_t cont = Ngam / (DeltaE * Effon * Aeff); + const Double_t dcont = sqrt(c1*c1 + c2*c2 + c3*c3); + + // + // Out of Range + // + const Bool_t oor = Ngam<=0 || DeltaE<=0 || Effon<=0 || Aeff<=0; + + if (oor) + *fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") "; + + if (Ngam<=0) + *fLog << " Ngam=0"; + if (DeltaE<=0) + *fLog << " DeltaE=0"; + if (Effon<=0) + *fLog << " Effon=0"; + if (Aeff<=0) + *fLog << " Aeff=0"; + + if (oor) + *fLog << endl; + + fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont); + fHFlux.SetBinError(m,n, oor ? 1e-20 : dcont*cont); + */ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Paint the integral and the error on top of the histogram +// +void MHAlpha::PaintText(Double_t val, Double_t error) const +{ + TLatex text(0.45, 0.94, Form("N_{exc} = %.1fs \\pm %.1fs", val, error)); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); +} + +// -------------------------------------------------------------------------- +// +// Update the projections +// +void MHAlpha::Paint(Option_t *opt) +{ + // Note: Any object cannot be added to a pad twice! + // The object is searched by gROOT->FindObject only in + // gPad itself! + TVirtualPad *padsave = gPad; + + TH1D *h0=0; + + TString o(opt); + if (o==(TString)"proj") + { + TPaveStats *st=0; + for (int x=0; x<4; x++) + { + TVirtualPad *p=padsave->GetPad(x+1); + if (!p || !(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); + } + + padsave->cd(1); + fHAlpha.ProjectionZ(fNameProjAlpha); + FitEnergyBins(); + FitThetaBins(); + } + + if (o==(TString)"alpha") + if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha))) + { + // Do not store the 'final' result in fFit + MAlphaFitter fit(fFit); + fit.Fit(*h0, kTRUE); + fit.PaintResult(); + } + + if (o==(TString)"time") + PaintText(fHTime.Integral(), 0); + + if (o==(TString)"theta") + PaintText(fHTheta.Integral(), 0); + + if (o==(TString)"energy") + { + if (fHEnergy.GetEntries()>0) + { + gPad->SetLogx(); + gPad->SetLogy(); + } + FitEnergySpec(kTRUE); + } + + gPad = padsave; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHAlpha::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + // Do the projection before painting the histograms into + // the individual pads + AppendPad("proj"); + + pad->SetBorderMode(0); + pad->Divide(2,2); + + TH1D *h=0; + + pad->cd(1); + gPad->SetBorderMode(0); + h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E"); + h->SetBit(TH1::kNoTitle); + h->SetXTitle("\\alpha [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->Draw(); + // After the Alpha-projection has been drawn. Fit the histogram + // and paint the result into this pad + AppendPad("alpha"); + + if (fHEnergy.GetNbinsX()>1) + { + pad->cd(2); + gPad->SetBorderMode(0); + fHEnergy.Draw(); + AppendPad("energy"); + } + else + delete pad->GetPad(2); + + if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1) + { + pad->cd(3); + gPad->SetBorderMode(0); + fHTime.Draw(); + AppendPad("time"); + } + else + delete pad->GetPad(3); + + if (fHTheta.GetNbinsX()>1) + { + pad->cd(4); + gPad->SetBorderMode(0); + fHTheta.Draw(); + AppendPad("theta"); + } + else + delete pad->GetPad(4); + +} + +// -------------------------------------------------------------------------- +// +// 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 (alphaGetBinWidth(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; i0 && chi2<2.5*NDF; + + // ------------------------------------ + if (paint) + { + func.SetLineColor(kGreen); + func.SetLineWidth(2); + 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; + + fSignificance = MMath::SignificanceLiMaSigned(s, b); + //exc = s-b; + + const Double_t uin = 1.25*sigint; + const Int_t bin = h.GetXaxis()->FindFixBin(uin); + fIntegralMax = h.GetBinLowEdge(bin+1); + fExcessEvents = h.Integral(0, bin)-func.Integral(0, fIntegralMax)/alphaw; + + return kTRUE; +} + +void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const +{ + TLatex text(x, y, Form("\\sigma_{Li/Ma}=%.1f (\\alpha<%.1f\\circ) \\omega=%.1f\\circ E=%d (\\alpha<%.1f\\circ) (\\chi_{b}^{2}=%.1f \\chi_{s}^{2}=%.1f)", + fSignificance, fSigInt, fSigmaGaus, + (int)fExcessEvents, fIntegralMax, + fChiSqBg, fChiSqSignal)); + + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(size); + text.Paint(); +} + +Bool_t MHAlpha::Finalize() +{ + // Store the final result in fFit + TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E"); + Bool_t rc = fFit.Fit(*h); + delete h; + if (!rc) + { + *fLog << warn << "Histogram empty." << endl; + return kTRUE; + } + + if (fResult) + fResult->SetVal(fFit.GetSignificance()); + + FitEnergyBins(); + FitThetaBins(); + UpdateAlphaTime(kTRUE); + MH::RemoveFirstBin(fHTime); + + 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[0] = fMatrix->AddColumn("MHillasSrc.fAlpha"); + //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy"); +} + +void MHAlpha::StopMapping() +{ + fMatrix = NULL; +} Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 4999) @@ -0,0 +1,132 @@ +#ifndef MARS_MHAlpha +#define MARS_MHAlpha + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef MARS_MTime +#include "MTime.h" +#endif + +#ifndef ROOT_TH3 +#include +#endif + +class MParList; +class MParameterD; +class MEnergyEst; +class MHMatrix; +class MPointingPos; + + +class MAlphaFitter : public MParContainer +{ +private: + Float_t fSigInt; + Float_t fSigMax; + Float_t fBgMin; + Float_t fBgMax; + Int_t fPolynom; + + Double_t fSignificance; + Double_t fExcessEvents; + Double_t fChiSqSignal; + Double_t fChiSqBg; + Double_t fSigmaGaus; + Double_t fIntegralMax; + +public: + MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynom(2) + { + } + + MAlphaFitter(const MAlphaFitter &f) + { + f.Copy(*this); + } + + void Copy(TObject &o) const + { + MAlphaFitter &f = static_cast(o); + + f.fSigInt = fSigInt; + f.fSigMax = fSigMax; + f.fBgMin = fBgMin; + f.fBgMax = fBgMax; + f.fPolynom = fPolynom; + } + + void SetSignalIntegralMax(Float_t s) { fSigInt = s; } + void SetSignalFitMax(Float_t s) { fSigMax = s; } + void SetBackgroundFitMin(Float_t s) { fBgMin = s; } + void SetBackgroundFitMax(Float_t s) { fBgMax = s; } + void SetNumPolynom(Int_t s) { fPolynom = s; } + + Double_t GetExcessEvents() const { return fExcessEvents; } + Double_t GetSignificance() const { return fSignificance; } + Double_t GetChiSqSignal() const { return fChiSqSignal; } + Double_t GetChiSqBg() const { return fChiSqBg; } + Double_t GetSigmaGaus() const { return fSigmaGaus; } + + void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const; + + Bool_t Fit(TH1D &h, Bool_t paint=kFALSE); + ClassDef(MAlphaFitter, 1) +}; + + +class MHAlpha : public MH +{ +private: + MAlphaFitter fFit; + + TH3D fHAlpha; // Alpha vs. theta and energy + TH1D fHEnergy; // excess events vs energy + TH1D fHTheta; // excess events vs theta + TH1D fHTime; // excess events vs time; + TH1D fHAlphaTime; //! temporary histogram to get alpha vs. time + + MParameterD *fResult; //! + MEnergyEst *fEnergy; //! + MPointingPos *fPointPos; //! + + MTime *fTimeEffOn; //! Time to steer filling of fHTime + MTime *fTime; //! Event-Time to finalize fHTime + MTime fLastTime; //! Last fTimeEffOn + + const TString fNameProjAlpha; //! This should make sure, that gROOT doen't confuse the projection with something else + + MHMatrix *fMatrix; //! + Int_t fMap[2]; //! + + void FitEnergySpec(Bool_t paint=kFALSE); + void FitEnergyBins(Bool_t paint=kFALSE); + void FitThetaBins(Bool_t paint=kFALSE); + + void UpdateAlphaTime(Bool_t final=kFALSE); + void InitAlphaTime(const MTime &t); + void FinalAlphaTime(MBinning &bins); + + void PaintText(Double_t val, Double_t error) const; + +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(); + + const TH3D &GetHist() const { return fHAlpha; } + const MAlphaFitter &GetAlphaFitter() const { return fFit; } + + 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/mhflux/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.cc (revision 4999) @@ -0,0 +1,824 @@ +/* ======================================================================== *\ +! +! * +! * 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("\\labda"); + 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("\\lambda"); + 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; + + // 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; + fLastTime = *time; + + fH2DeltaT.Fill(dt, fPointPos->GetZd(), w); + fH1DeltaT.Fill(dt, w); + + // + // 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::PaintText(Double_t *res) const +{ + TLatex text(0.25, 0.94, Form("T_{eff}=%.1fs\\pm%.1fs \\labda=%.1f\\pm%.1f p=%.1f%% NDF=%d", + res[0], res[1], res[3], res[4], res[2], res[5])); + 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) +{ + 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; + + 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]; + FitH(h, res, kTRUE); + PaintText(res); + } + 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("\\lambda"); + + 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("\\lambda"); +} Index: trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.h (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHEffectiveOnTime.h (revision 4999) @@ -0,0 +1,84 @@ +#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 PaintText(Double_t *res) const; + void DrawRightAxis(const char *title); + void UpdateRightAxis(TH1 &h); + void PrintStatistics(); + +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/mhflux/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 4999) @@ -0,0 +1,1104 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHFalseSource +// +// Create a false source plot. For the Binning in x,y the object MBinning +// "BinningFalseSource" is taken from the paremeter list. +// +// The binning in alpha is currently fixed as 18bins from 0 to 90deg. +// +// If MTime, MObservatory and MPointingPos is available in the paremeter +// list each image is derotated. +// +// MHFalseSource fills a 3D histogram with alpha distribution for +// each (derotated) x and y position. +// +// Only a radius of 1.2deg around the camera center is taken into account. +// +// The displayed histogram is the projection of alpha90-fAlphaCut +// +// By using the context menu (find it in a small region between the single +// pads) you can change the AlphaCut 'online' +// +// Each Z-Projection (Alpha-histogram) is scaled such, that the number +// of entries fits the maximum number of entries in all Z-Projections. +// This should correct for the different acceptance in the center +// and at the vorder of the camera. This, however, produces more noise +// at the border. +// +// Here is a slightly simplified version of the algorithm: +// ------------------------------------------------------ +// MHillas hil; // Taken as second argument in MFillH +// +// MSrcPosCam src; +// MHillasSrc hsrc; +// hsrc.SetSrcPos(&src); +// +// for (int ix=0; ix +#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 "MAstroCatalog.h" +#include "MAstroSky2Local.h" +#include "MStatusDisplay.h" +#include "MMath.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHFalseSource); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MHFalseSource::MHFalseSource(const char *name, const char *title) + : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5), + fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHFalseSource"; + fTitle = title ? title : "3D-plot of Alpha vs x, y"; + + fHist.SetDirectory(NULL); + + fHist.SetName("Alpha"); + fHist.SetTitle("3D-plot of Alpha vs x, y"); + fHist.SetXTitle("x [\\circ]"); + fHist.SetYTitle("y [\\circ]"); + fHist.SetZTitle("\\alpha [\\circ]"); +} + +void MHFalseSource::MakeSymmetric(TH1 *h) +{ + const Float_t min = TMath::Abs(h->GetMinimum()); + const Float_t max = TMath::Abs(h->GetMaximum()); + + const Float_t absmax = TMath::Max(min, max)*1.002; + + h->SetMaximum( absmax); + h->SetMinimum(-absmax); +} + +// -------------------------------------------------------------------------- +// +// Set the alpha cut (|alpha|FindObject("MGeomCam"); + if (!geom) + { + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = geom->GetConvMm2Deg(); + + MBinning binsa; + binsa.SetEdges(18, 0, 90); + + const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource"); + if (!bins) + { + const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg; + + MBinning b; + b.SetEdges(20, -r, r); + SetBinning(&fHist, &b, &b, &binsa); + } + else + SetBinning(&fHist, bins, bins, &binsa); + + fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos")); + if (!fPointPos) + *fLog << warn << "MPointingPos not found... no derotation." << endl; + + fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime")); + if (!fTime) + *fLog << warn << "MTime not found... no derotation." << endl; + + fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + *fLog << warn << "MSrcPosCam not found... no translation." << endl; + + fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory")); + if (!fObservatory) + *fLog << warn << "MObservatory not found... no derotation." << endl; + + // FIXME: Because the pointing position could change we must check + // for the current pointing position and add a offset in the + // Fill function! + fRa = fPointPos->GetRa(); + fDec = fPointPos->GetDec(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram. For details see the code or the class description +// +Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w) +{ + const MHillas *hil = dynamic_cast(par); + if (!hil) + { + *fLog << err << "MHFalseSource::Fill: No container specified!" << endl; + return kFALSE; + } + + // Get max radius... + const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1)); + + // Get camera rotation angle + Double_t rho = 0; + if (fTime && fObservatory && fPointPos) + rho = fPointPos->RotationAngle(*fObservatory, *fTime); + //if (fPointPos) + // rho = fPointPos->RotationAngle(*fObservatory); + + // Create necessary containers for calculation + MSrcPosCam src; + MHillasSrc hsrc; + hsrc.SetSrcPos(&src); + + // Get number of bins and bin-centers + const Int_t nx = fHist.GetNbinsX(); + const Int_t ny = fHist.GetNbinsY(); + Axis_t cx[nx]; + Axis_t cy[ny]; + fHist.GetXaxis()->GetCenter(cx); + fHist.GetYaxis()->GetCenter(cy); + + for (int ix=0; ixmaxr) + continue; + + // rotate center of bin + // precalculation of sin/cos doesn't accelerate + TVector2 v(cx[ix], cy[iy]); + if (rho!=0) + v=v.Rotate(rho); + + // convert degrees to millimeters + v *= 1./fMm2Deg; + + if (fSrcPos) + v += fSrcPos->GetXY(); + + src.SetXY(v); + + // Source dependant hillas parameters + if (hsrc.Calc(*hil)>0) + { + *fLog << warn << "Calculation of MHillasSrc failed for x=" << cx[ix] << " y=" << cy[iy] << endl; + return kFALSE; + } + + // FIXME: This should be replaced by an external MFilter + // and/or MTaskList + // Source dependant distance cut + if (fMinDist>0 && hsrc.GetDist()*fMm2Deg0 && hsrc.GetDist()*fMm2Deg>fMaxDist) + continue; + + if (fMaxLD>0 && hil->GetLength()>fMaxLD*hsrc.GetDist()) + continue; + if (fMinLD>0 && hil->GetLength()SetTitle(Form("Distribution of %.1f\\circ<|\\alpha|<%.1f\\circ in x,y", cut1, cut2)); + + // Get projection for range + TH2D *p = (TH2D*)fHist.Project3D("yx_off"); + + // Reset range + axe.SetRange(0,9999); + + // Move contents from projection to h2 + h2->Reset(); + h2->Add(p, all->GetMaximum()); + h2->Divide(all); + + // Delete p + delete p; + + // Set Minimum as minimum value Greater Than 0 + h2->SetMinimum(GetMinimumGT(*h2)); +} + +// -------------------------------------------------------------------------- +// +// Create projection for on data, taking sum of bin contents of +// range (0, fAlphaCut) +// +void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all) +{ + TAxis &axe = *fHist.GetZaxis(); + + // Find range to cut + axe.SetRangeUser(0, fAlphaCut); + const Float_t cut = axe.GetBinUpEdge(axe.GetLast()); + h3->SetTitle(Form("Distribution of |\\alpha|<%.1f\\circ in x,y", cut)); + + // Get projection for range + TH2D *p = (TH2D*)fHist.Project3D("yx_on"); + + // Reset range + axe.SetRange(0,9999); + + // Move contents from projection to h3 + h3->Reset(); + h3->Add(p, all->GetMaximum()); + h3->Divide(all); + + // Delete p + delete p; + + // Set Minimum as minimum value Greater Than 0 + h3->SetMinimum(GetMinimumGT(*h3)); +} + +// -------------------------------------------------------------------------- +// +// Create projection for all data, taking sum of bin contents of +// range (0, 90) - corresponding to the number of entries in this slice. +// +void MHFalseSource::ProjectAll(TH2D *h3) +{ + h3->SetTitle("Number of entries"); + + // Get projection for range + TH2D *p = (TH2D*)fHist.Project3D("yx_all"); + + // Move contents from projection to h3 + h3->Reset(); + h3->Add(p); + delete p; + + // Set Minimum as minimum value Greater Than 0 + h3->SetMinimum(GetMinimumGT(*h3)); +} + +// -------------------------------------------------------------------------- +// +// Update the projections and paint them +// +void MHFalseSource::Paint(Option_t *opt) +{ + // Set pretty color palette + gStyle->SetPalette(1, 0); + + TVirtualPad *padsave = gPad; + + TH1D* h1; + TH2D* h0; + TH2D* h2; + TH2D* h3; + TH2D* h4; + TH2D* h5; + + /* + fHistProjAll = Form("All_%p", this); + fHistProjOn = Form("On_%p", this); + fHistProjOff = Form("Off_%p", this); + fHistProjDiff = Form("Diff_%p", this); + fHistProjAll = Form("All_%p", this); + */ + + // Update projection of all-events + padsave->GetPad(2)->cd(3); + if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all"))) + ProjectAll(h0); + + // Update projection of on-events + padsave->GetPad(1)->cd(1); + if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) + ProjectOn(h3, h0); + + // Update projection of off-events + padsave->GetPad(1)->cd(3); + if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) + ProjectOff(h2, h0); + + padsave->GetPad(2)->cd(2); + if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff"))) + { + h5->Add(h2, h3, -1); + MakeSymmetric(h5); + } + + // Update projection of significance + padsave->GetPad(1)->cd(2); + if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig"))) + { + const Int_t nx = h4->GetXaxis()->GetNbins(); + const Int_t ny = h4->GetYaxis()->GetNbins(); + //const Int_t nr = nx*nx + ny*ny; + + Int_t maxx=nx/2; + Int_t maxy=ny/2; + + Int_t max = h4->GetBin(nx, ny); + + for (int ix=1; ix<=nx; ix++) + for (int iy=1; iy<=ny; iy++) + { + const Int_t n = h4->GetBin(ix, iy); + + const Double_t s = h3->GetBinContent(n); + const Double_t b = h2->GetBinContent(n); + + const Double_t sig = SignificanceLiMa(s, b); + + h4->SetBinContent(n, sig); + + if (sig>h4->GetBinContent(max) && sig>0/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)GetPad(2)->cd(1); + if ((h1 = (TH1D*)gPad->FindObject("Alpha_z")) && max>0) + { + const Double_t x = h4->GetXaxis()->GetBinCenter(maxx); + const Double_t y = h4->GetYaxis()->GetBinCenter(maxy); + const Double_t s = h4->GetBinContent(max); + + // This is because a 3D histogram x and y are vice versa + // Than for their projections + TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy); + h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s)); + } + } + + gPad = padsave; +} + +// -------------------------------------------------------------------------- +// +// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude +// is set to 9, while the fov is equal to the current fov of the false +// source plot. +// +TObject *MHFalseSource::GetCatalog() +{ + const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1)); + + // Create catalog... + MAstroCatalog *stars = new MAstroCatalog; + stars->SetLimMag(9); + stars->SetGuiActive(kFALSE); + stars->SetRadiusFOV(maxr); + stars->SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad()); + stars->ReadBSC("bsc5.dat"); + + *fLog << err << "FIXME - The catalog will never be deleted, because this crashes!" << endl; + +// stars->SetBit(kCanDelete); + + return stars; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHFalseSource::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + pad->Divide(1, 2, 0, 0.03); + + *fLog << err << "FIXME - Plotting the catalog is broken!" << endl; + + TObject *catalog = GetCatalog(); + + // Initialize upper part + pad->cd(1); + gPad->SetBorderMode(0); + gPad->Divide(3, 1); + + // PAD #1 + pad->GetPad(1)->cd(1); + gPad->SetBorderMode(0); + fHist.GetZaxis()->SetRangeUser(0,fAlphaCut); + TH1 *h3 = fHist.Project3D("yx_on"); + fHist.GetZaxis()->SetRange(0,9999); + h3->SetDirectory(NULL); + h3->SetXTitle(fHist.GetXaxis()->GetTitle()); + h3->SetYTitle(fHist.GetYaxis()->GetTitle()); + h3->Draw("colz"); + h3->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + + // PAD #2 + pad->GetPad(1)->cd(2); + gPad->SetBorderMode(0); + fHist.GetZaxis()->SetRange(0,0); + TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning.... + fHist.GetZaxis()->SetRange(0,9999); + h4->SetTitle("Significance"); + h4->SetDirectory(NULL); + h4->SetXTitle(fHist.GetXaxis()->GetTitle()); + h4->SetYTitle(fHist.GetYaxis()->GetTitle()); + h4->Draw("colz"); + h4->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + + // PAD #3 + pad->GetPad(1)->cd(3); + gPad->SetBorderMode(0); + fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); + TH1 *h2 = fHist.Project3D("yx_off"); + h2->SetDirectory(NULL); + h2->SetXTitle(fHist.GetXaxis()->GetTitle()); + h2->SetYTitle(fHist.GetYaxis()->GetTitle()); + h2->Draw("colz"); + h2->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + + // Initialize lower part + pad->cd(2); + gPad->SetBorderMode(0); + gPad->Divide(3, 1); + + // PAD #4 + pad->GetPad(2)->cd(1); + gPad->SetBorderMode(0); + TH1 *h1 = fHist.ProjectionZ("Alpha_z"); + h1->SetDirectory(NULL); + h1->SetTitle("Distribution of \\alpha"); + h1->SetXTitle(fHist.GetZaxis()->GetTitle()); + h1->SetYTitle("Counts"); + h1->Draw(opt); + h1->SetBit(kCanDelete); + + // PAD #5 + pad->GetPad(2)->cd(2); + gPad->SetBorderMode(0); + TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff"); + h5->Add(h2, -1); + h5->SetTitle("Difference of on- and off-distribution"); + h5->SetDirectory(NULL); + h5->Draw("colz"); + h5->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + + // PAD #6 + pad->GetPad(2)->cd(3); + gPad->SetBorderMode(0); + TH1 *h0 = fHist.Project3D("yx_all"); + h0->SetDirectory(NULL); + h0->SetXTitle(fHist.GetXaxis()->GetTitle()); + h0->SetYTitle(fHist.GetYaxis()->GetTitle()); + h0->Draw("colz"); + h0->SetBit(kCanDelete); + catalog->Draw("mirror same *"); +} + +// -------------------------------------------------------------------------- +// +// Everything which is in the main pad belongs to this class! +// +Int_t MHFalseSource::DistancetoPrimitive(Int_t px, Int_t py) +{ + return 0; +} + +// -------------------------------------------------------------------------- +// +// Set all sub-pads to Modified() +// +void MHFalseSource::Modified() +{ + if (!gPad) + return; + + TVirtualPad *padsave = gPad; + padsave->Modified(); + padsave->GetPad(1)->cd(1); + gPad->Modified(); + padsave->GetPad(1)->cd(3); + gPad->Modified(); + padsave->GetPad(2)->cd(1); + gPad->Modified(); + padsave->GetPad(2)->cd(2); + gPad->Modified(); + padsave->GetPad(2)->cd(3); + gPad->Modified(); + gPad->cd(); +} + +// -------------------------------------------------------------------------- +// +// 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 (alphaSetDirectory(0); + TH1 *hists = fHist.Project3D("yx_fit"); + hists->SetDirectory(0); + TH1 *histb = fHist.Project3D("yx_fit"); + histb->SetDirectory(0); + hist->Reset(); + hists->Reset(); + histb->Reset(); + hist->SetNameTitle("Significance", + Form("Fit Region: Signal<%.1f\\circ, %.1f\\circSetName("Excess"); + histb->SetName("Background"); + hist->SetXTitle(fHist.GetXaxis()->GetTitle()); + hists->SetXTitle(fHist.GetXaxis()->GetTitle()); + histb->SetXTitle(fHist.GetXaxis()->GetTitle()); + hist->SetYTitle(fHist.GetYaxis()->GetTitle()); + hists->SetYTitle(fHist.GetYaxis()->GetTitle()); + histb->SetYTitle(fHist.GetYaxis()->GetTitle()); + + const Double_t w = fHist.GetZaxis()->GetBinWidth(1); + + // xmin, xmax, npar + //TF1 func("MyFunc", fcn, 0, 90, 6); + // Implementing the function yourself is only about 5% faster + TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90); + //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90); + TArrayD maxpar(func.GetNpar()); + + /* func.SetParName(0, "A"); + * func.SetParName(1, "mu"); + * func.SetParName(2, "sigma"); + */ + + func.FixParameter(1, 0); + func.FixParameter(4, 0); + func.SetParLimits(2, 0, 90); + func.SetParLimits(3, -1, 1); + + const Int_t nx = hist->GetXaxis()->GetNbins(); + const Int_t ny = hist->GetYaxis()->GetNbins(); + //const Int_t nr = nx*nx+ny*ny; + + Double_t maxalpha0=0; + Double_t maxs=3; + + Int_t maxx=0; + Int_t maxy=0; + + TStopwatch clk; + clk.Start(); + + *fLog << inf; + *fLog << "Signal fit: alpha < " << sigmax << endl; + *fLog << "Integration: alpha < " << sigint << endl; + *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl; + *fLog << "Polynom order: " << (int)polynom << endl; + *fLog << "Fitting False Source Plot..." << flush; + + TH1 *h0 = fHist.Project3D("yx_entries"); + Float_t entries = h0->GetMaximum(); + delete h0; + + TH1 *h=0; + for (int ix=1; ix<=nx; ix++) + for (int iy=1; iy<=ny; iy++) + { + // This is because a 3D histogram x and y are vice versa + // Than for their projections + h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy); + + const Double_t alpha0 = h->GetBinContent(1); + + // Check for the regios which is not filled... + if (alpha0==0) + continue; + + h->Scale(entries/h->GetEntries()); + + if (alpha0>maxalpha0) + maxalpha0=alpha0; + + // 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); + + h4a.Fill(func.GetChisquare()); + //h5a.Fill(func.GetProb()); + + //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1)); + //func.SetParLimits(2, 5, 90); + + func.ReleaseParameter(0); + //func.ReleaseParameter(1); + func.ReleaseParameter(2); + func.FixParameter(3, func.GetParameter(3)); + for (int i=5; iFit(&func, "N0Q", "", 0, sigmax); + + TArrayD p(func.GetNpar(), func.GetParameters()); + + // Fill results into some histograms + h0a.Fill(p[0]); + h0b.Fill(p[3]); + h1.Fill(p[1]); + h2.Fill(p[2]); + if (polynom>1) + h3.Fill(p[5]); + h4b.Fill(func.GetChisquare()); + //h5b.Fill(func.GetProb()); + + // Implementing the integral as analytical function + // gives the same result in the order of 10e-5 + // and it is not faster at all... + //const Bool_t ok = /*p[0]>=0 && /*p[0]1.75;// && p[2]<88.5; + /* + if (p[0]<-fabs(alpha0-p[3])*7 && p[2]GetBin(ix, iy); + hists->SetBinContent(n, s-b); + histb->SetBinContent(n, b); + + hist->SetBinContent(n, sig); + if (sig!=0) + h6.Fill(sig); + + if (sig>maxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)SetRangeUser(0, maxalpha0*1.5); + h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5); + + hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum())); + histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint)); + + //hists->SetMinimum(GetMinimumGT(*hists)); + histb->SetMinimum(GetMinimumGT(*histb)); + + MakeSymmetric(hists); + MakeSymmetric(hist); + + clk.Stop(); + clk.Print("m"); + + TCanvas *c=new TCanvas; + + gStyle->SetPalette(1, 0); + + c->Divide(3,2, 0, 0); + c->cd(1); + gPad->SetBorderMode(0); + hists->Draw("colz"); + hists->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + c->cd(2); + gPad->SetBorderMode(0); + hist->Draw("colz"); + hist->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + c->cd(3); + gPad->SetBorderMode(0); + histb->Draw("colz"); + histb->SetBit(kCanDelete); + catalog->Draw("mirror same *"); + c->cd(4); + gPad->Divide(1,3, 0, 0); + TVirtualPad *p=gPad; + p->SetBorderMode(0); + p->cd(1); + gPad->SetBorderMode(0); + h0b.DrawCopy(); + h0a.DrawCopy("same"); + p->cd(2); + gPad->SetBorderMode(0); + h3.DrawCopy(); + p->cd(3); + gPad->SetBorderMode(0); + h2.DrawCopy(); + c->cd(6); + gPad->Divide(1,2, 0, 0); + TVirtualPad *q=gPad; + q->SetBorderMode(0); + q->cd(1); + gPad->SetBorderMode(0); + gPad->SetBorderMode(0); + h4b.DrawCopy(); + h4a.DrawCopy("same"); + h6.DrawCopy("same"); + q->cd(2); + gPad->SetBorderMode(0); + //h5b.DrawCopy(); + //h5a.DrawCopy("same"); + c->cd(5); + gPad->SetBorderMode(0); + if (maxx>0 && maxy>0) + { + const char *title = Form(" \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f) ", + hist->GetXaxis()->GetBinCenter(maxx), + hist->GetYaxis()->GetBinCenter(maxy), maxs); + + TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); + result->Scale(entries/h->GetEntries()); + + result->SetDirectory(NULL); + result->SetNameTitle("Result \\alpha", title); + result->SetBit(kCanDelete); + result->SetXTitle("\\alpha [\\circ]"); + result->SetYTitle("Counts"); + result->Draw(); + + TF1 f1("", func.GetExpFormula(), 0, 90); + TF1 f2("", func.GetExpFormula(), 0, 90); + f1.SetParameters(maxpar.GetArray()); + f2.SetParameters(maxpar.GetArray()); + f2.FixParameter(0, 0); + f2.FixParameter(1, 0); + f2.FixParameter(2, 1); + f1.SetLineColor(kGreen); + f2.SetLineColor(kRed); + + f2.DrawCopy("same"); + f1.DrawCopy("same"); + + TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC"); + leg->SetBorderSize(1); + leg->SetTextSize(0.04); + leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22); + //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22); + leg->AddLine(0, 0.65, 0, 0.65); + leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11); + leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11); + leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11); + leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11); + leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11); + if (polynom>1) + leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11); + leg->SetBit(kCanDelete); + leg->Draw(); + + q->cd(2); + + TGraph *g = new TGraph; + g->SetPoint(0, 0, 0); + + Int_t max=0; + Float_t maxsig=0; + for (int i=1; i<89; i++) + { + const Double_t s = f1.Integral(0, (float)i)/w; + const Double_t b = f2.Integral(0, (float)i)/w; + + const Double_t sig = SignificanceLiMa(s, b); + + g->SetPoint(g->GetN(), i, sig); + + if (sig>maxsig) + { + max = i; + maxsig = sig; + } + } + + g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha"); + g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]"); + g->GetHistogram()->SetYTitle("Significance"); + g->SetBit(kCanDelete); + g->Draw("AP"); + + leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC"); + leg->SetBorderSize(1); + leg->SetTextSize(0.1); + leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max)); + leg->SetBit(kCanDelete); + leg->Draw(); + } +} Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 4999) @@ -0,0 +1,88 @@ +#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/mhflux/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhflux/Makefile (revision 4999) +++ trunk/MagicSoft/Mars/mhflux/Makefile (revision 4999) @@ -0,0 +1,36 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Flux + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhbase -I../mraw -I../manalysis \ + -I../mgui -I../mgeom -I../mdata -I../mfilter -I../mimage \ + -I../mmain -I../mmc -I../mreflector -I../mpointing \ + -I../mastro -I../mpedestal -I../msignal -I../mbadpixels + +SRCFILES = MHAlpha.cc \ + MHEffectiveOnTime.cc \ + MHFalseSource.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak