Index: trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 7094) +++ trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h (revision 7109) @@ -7,4 +7,5 @@ #pragma link C++ class MAlphaFitter+; +#pragma link C++ class MHDisp+; #pragma link C++ class MHAlpha+; #pragma link C++ class MHThetaSq+; @@ -13,4 +14,9 @@ #pragma link C++ class MHEffectiveOnTime+; #pragma link C++ class MHCollectionArea+; +#pragma link C++ class MMcSpectrumWeight+; + +//#pragma link C++ class MMatrixHist+; +//#pragma link C++ class MMatrixUnfold+; +//#pragma link C++ class MUnfoldBert+; #endif Index: trunk/MagicSoft/Mars/mhflux/MHDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7109) +++ trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7109) @@ -0,0 +1,351 @@ +/* ======================================================================== *\ +! +! * +! * 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, 5/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHDisp +// +// Create a false source plot using disp. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHDisp.h" + +#include +#include + +#include +#include +#include + +#include "MParList.h" +#include "MParameters.h" + +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MSrcPosCam.h" +#include "MPointingPos.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHDisp); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MHDisp::MHDisp(const char *name, const char *title) + : fHilExt(0), fDisp(0) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHDisp"; + fTitle = title ? title : "3D-plot using Disp vs x, y"; + + fHist.SetDirectory(NULL); + + fHist.SetName("Alpha"); + fHist.SetTitle("3D-plot of ThetaSq vs x, y"); + fHist.SetXTitle("x [\\circ]"); + fHist.SetYTitle("y [\\circ]"); + fHist.SetZTitle("\\vartheta [deg^{2}]"); +} + +// -------------------------------------------------------------------------- +// +// Set binnings (takes BinningFalseSource) and prepare filling of the +// histogram. +// +// Also search for MTime, MObservatory and MPointingPos +// +//MHillasExt *ext=0; +Bool_t MHDisp::SetupFill(const MParList *plist) +{ + if (!MHFalseSource::SetupFill(plist)) + return kFALSE; + + fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD"); + if (!fDisp) + { + *fLog << err << "Disp [MParameterD] not found... abort." << endl; + return kFALSE; + } + + MParameterD *p = (MParameterD*)plist->FindObject("M3lCut", "MParameterD"); + if (!p) + { + *fLog << err << "M3lCut [MParameterD] not found... abort." << endl; + return kFALSE; + } + fM3lCut = TMath::Abs(p->GetVal()); + + p = (MParameterD*)plist->FindObject("DispXi", "MParameterD"); + if (!p) + { + *fLog << err << "DispXi [MParameterD] not found... abort." << endl; + return kFALSE; + } + fXi = p->GetVal(); + + p = (MParameterD*)plist->FindObject("DispXiTheta", "MParameterD"); + if (!p) + { + *fLog << err << "DispXiTheta [MParameterD] not found... abort." << endl; + return kFALSE; + } + fXiTheta = p->GetVal(); + + fHilExt = (MHillasExt*)plist->FindObject("MHillasExt"); + if (!fHilExt) + { + *fLog << err << "MHillasExt not found... aborting." << endl; + return kFALSE; + } + + // Initialize all bins with a small (=0) value otherwise + // these bins are not displayed + for (int x=0; xGetBinCenter(x+1), + fHist.GetYaxis()->GetBinCenter(y+1), + 0.0, 1e-10); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram. For details see the code or the class description +// +Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w) +{ + const MHillas *hil = dynamic_cast(par); + if (!hil) + { + *fLog << err << "MHDisp::Fill: No container specified!" << endl; + return kFALSE; + } + + // Get camera rotation angle + Double_t rho = 0; + if (fTime && fObservatory && fPointPos) + rho = fPointPos->RotationAngle(*fObservatory, *fTime); + + // Get Disp from Parlist + const Double_t disp = fDisp->GetVal(); + + // Calculate where disp is pointing + TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta()); + TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta()); + pos1 *= disp; + pos2 *= -disp; + + pos1 += hil->GetMean()*fMm2Deg; + pos2 += hil->GetMean()*fMm2Deg; + + // gammaweight: If we couldn't decide which position makes the + // event a gamma, both position are assigned 'half of a gamma' + Double_t gweight = 0.5; + + // Check whether our g/h-separation allows to asign definitly + // to one unique position. Therefor we requeire that the event + // would be considered a gamma for one, but not for the other + // position. This can only be the case if the third moment + // has a value higher than the absolute cut value. + if (TMath::Abs(fHilExt->GetM3Long()) > fM3lCut) + { + // Because at one position the event is considered a gamma + // we have to find out which position it is... + MSrcPosCam src; + MHillasSrc hsrc; + hsrc.SetSrcPos(&src); + + // Calculate the sign for one of the desired source positions + // The other position must have the opposite sign + src.SetXY(pos1/fMm2Deg); + if (hsrc.Calc(*hil)>0) + { + *fLog << warn << "Calculation of MHillasSrc failed." << endl; + return kFALSE; + } + const Double_t m3l = fHilExt->GetM3Long()*TMath::Sign(1.0f, hsrc.GetCosDeltaAlpha())*fMm2Deg; + + gweight = m3l>fM3lCut ? 1 : 0; + } + + // Now we can safly derotate both position... + //TVector2 srcp(fSrcPos->GetXY()); + if (rho!=0) + { + pos1=pos1.Rotate(-rho); + pos2=pos2.Rotate(-rho); + //srcp=srcp.Rotate(-rho); + } + + /*if (fSrcPos) + { + pos1 -= srcp*fMm2Deg; + pos2 -= srcp*fMm2Deg; + }*/ + + // Workaround: Number-of-entries + if (gweight>0.25) + fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight); + if (gweight<0.75) + fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight)); + + return kTRUE; +} +/* +static Double_t FcnGauss2d(Double_t *x, Double_t *par) +{ + TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad()); + + const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]); + const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]); + + //Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE); + return par[0]*g0*g1 + par[6]*(v.X()*v.X() + v.Y()*v.Y()) + par[7]; +}*/ + +void MHDisp::Paint(Option_t *o) +{ + TVirtualPad *pad = gPad; + + pad->cd(1); + + fHist.GetZaxis()->SetRange(0,9999); + TH1 *h1=fHist.Project3D("yx_on"); + gStyle->SetPalette(1, 0); + + Int_t ix, iy, iz; + TF1 func("fcn", "gaus + [3]*x*x + [4]"); + + pad->cd(3); + TH1 *h2 = (TH1*)gPad->FindObject("Radial"); + if (h1 && h2) + { + h2->Reset(); + + h1->GetMaximumBin(ix, iy, iz); + + const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix); + const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy); + + h2->SetTitle(Form("Profile @ (%.2f\\circ, %.2f\\circ)", x0, y0)); + + for (int x=0; xGetNbinsX(); x++) + for (int y=0; yGetNbinsY(); y++) + { + const Double_t r = TMath::Hypot(h1->GetXaxis()->GetBinCenter(x+1)-x0, + h1->GetXaxis()->GetBinCenter(y+1)-y0); + h2->Fill(r, h1->GetBinContent(x+1, y+1)); + } + + func.SetLineWidth(1); + func.SetLineColor(kBlue); + func.SetParameter(0, h2->GetBinContent(1)); + func.FixParameter(1, 0); + func.SetParameter(2, 0.05); + func.SetParameter(4, h2->GetBinContent(10)); + h2->Fit(&func, "IMQ"); + } + /* + if (h1) + { + const Double_t maxr = 0.9*TMath::Abs(fHist.GetXaxis()->GetXmax()); + + TF2 f2d("Gaus2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 8); + f2d.SetLineWidth(1); + + f2d.SetParameter(0, h1->GetMaximum()*5); // A + f2d.SetParLimits(1, h1->GetXaxis()->GetBinCenter(ix)-h1->GetXaxis()->GetBinWidth(ix)*5, h1->GetXaxis()->GetBinCenter(ix)+h1->GetXaxis()->GetBinWidth(ix)); // mu_1 + f2d.SetParLimits(3, h1->GetYaxis()->GetBinCenter(iy)-h1->GetYaxis()->GetBinWidth(iy)*5, h1->GetYaxis()->GetBinCenter(iy)+h1->GetYaxis()->GetBinWidth(iy)); // mu_2 + f2d.SetParLimits(2, 0, func.GetParameter(2)*5); // sigma_1 + f2d.SetParLimits(4, 0, func.GetParameter(2)*5); // sigma_2 + f2d.SetParLimits(5, 0, 45); // phi + f2d.SetParLimits(6, 0, func.GetParameter(3)*5); + f2d.SetParLimits(7, 0, func.GetParameter(4)*5); + + f2d.SetParameter(0, h1->GetMaximum()); // A + f2d.SetParameter(1, h1->GetXaxis()->GetBinCenter(ix)); // mu_1 + f2d.SetParameter(2, func.GetParameter(2)); // sigma_1 + f2d.SetParameter(3, h1->GetYaxis()->GetBinCenter(iy)); // mu_2 + f2d.SetParameter(4, func.GetParameter(2)); // sigma_2 + f2d.FixParameter(5, 0); // phi + f2d.SetParameter(6, func.GetParameter(3)); + f2d.SetParameter(7, func.GetParameter(4)); + h1->Fit(&f2d, "Q", "cont2"); + //f2d.DrawCopy("cont2same"); + }*/ +} + +void MHDisp::Draw(Option_t *o) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + const Int_t col = pad->GetFillColor(); + pad->SetBorderMode(0); + + AppendPad(""); + + TString name = Form("%s_1", pad->GetName()); + TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0); + p->SetNumber(1); + p->Draw(); + p->cd(); + + TH1 *h3 = fHist.Project3D("yx_on"); + h3->SetTitle("Distribution of equivalent events"); + h3->SetDirectory(NULL); + h3->SetXTitle(fHist.GetXaxis()->GetTitle()); + h3->SetYTitle(fHist.GetYaxis()->GetTitle()); + h3->SetMinimum(0); + h3->Draw("colz"); + h3->SetBit(kCanDelete); +// catalog->Draw("mirror same *"); + + pad->cd(); + name = Form("%s_2", pad->GetName()); + p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0); + p->SetNumber(2); + p->Draw(); + p->cd(); + h3->Draw("surf3"); + + pad->cd(); + name = Form("%s_3", pad->GetName()); + p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0); + p->SetNumber(3); + p->Draw(); + p->cd(); + + const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax()); + const Int_t nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2; + TProfile *h = new TProfile("Radial", "Radial Profile", nbin, 0, maxr, "s"); + h->SetBit(kCanDelete); + h->Draw(); +} Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 7094) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 7109) @@ -246,12 +246,14 @@ fMm2Deg = geom->GetConvMm2Deg(); - if (fName!=(TString)"MHFalseSourceOff" && fHistOff==NULL) - { - MHFalseSource *hoff = (MHFalseSource*)plist->FindObject("MHFalseSourceOff", "MHFalseSource"); + const TString off(Form("%sOff", fName.Data())); + if (fName!=off && fHistOff==NULL) + { + const TString desc(Form("%s [%s] found... using ", off.Data(), ClassName())); + MHFalseSource *hoff = (MHFalseSource*)plist->FindObject(off, ClassName()); if (!hoff) - *fLog << inf << "No MHFalseSourceOff [MHFalseSource] found... using current data only!" << endl; + *fLog << inf << "No " << desc << "current data only!" << endl; else { - *fLog << inf << "MHFalseSource [MHFalseSource] found... using on-off mode!" << endl; + *fLog << inf << desc << "on-off mode!" << endl; SetOffData(*hoff); } @@ -294,9 +296,13 @@ *fLog << warn << "MObservatory not found... no derotation." << endl; + MPointingPos *point = (MPointingPos*)plist->FindObject("MSourcePos", "MPointingPos"); + if (!point) + point = fPointPos; + // 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 ? fPointPos->GetRa() : 0; - fDec = fPointPos ? fPointPos->GetDec() : 90; + fRa = point ? point->GetRa() : 0; + fDec = point ? point->GetDec() : 90; return kTRUE; @@ -615,4 +621,6 @@ const Int_t l = h0->GetXaxis()->FindFixBin(fAlphaCut*3)+f-1; h0->Scale(h1->Integral(f, l)/h0->Integral(f, l)); + + //h0->Scale(h1->GetEntries()/h0->GetEntries()); @@ -854,5 +862,5 @@ } -Double_t FcnGauss2d(Double_t *x, Double_t *par) +static Double_t FcnGauss2d(Double_t *x, Double_t *par) { TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad()); Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 7094) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 7109) @@ -20,5 +20,5 @@ class MHFalseSource : public MH { -private: +protected: MTime *fTime; //! container to take the event time from MPointingPos *fPointPos; //! container to take pointing position from @@ -28,4 +28,5 @@ Float_t fMm2Deg; // conversion factor for display in degrees +private: Float_t fAlphaCut; // Alpha cut Float_t fBgMean; // Background mean @@ -37,4 +38,5 @@ Float_t fMaxDW; // Maximum distance in percent of dist +protected: TH3D fHist; // Alpha vs. x and y @@ -44,4 +46,5 @@ Double_t fDec; +private: Int_t DistancetoPrimitive(Int_t px, Int_t py); void Modified();