Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 7192) +++ /trunk/MagicSoft/Mars/Changelog (revision 7193) @@ -45,4 +45,7 @@ - comment have been updated. - some old code had been removed + + * mhflux/MHDisp.[h,cc]: + - updated with a more appropriate calculation of a background model Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 7192) +++ /trunk/MagicSoft/Mars/NEWS (revision 7193) @@ -115,4 +115,8 @@ more reliable in terms of gain-fluctuations and sudden VCSEL gain changes. + + - ganymed: The false source plot (MHDisp) is now based on Disp + and a background model determied in the first loop is + subtracted - sponde: the zenith angle distribution is now weighted instead of Index: /trunk/MagicSoft/Mars/mhflux/MHDisp.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7192) +++ /trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7193) @@ -46,4 +46,5 @@ #include "MMath.h" +#include "MBinning.h" #include "MParList.h" @@ -64,5 +65,5 @@ // MHDisp::MHDisp(const char *name, const char *title) - : fHilExt(0), fDisp(0)//, fIsWobble(kFALSE) + : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE) { // @@ -72,6 +73,4 @@ 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"); @@ -79,4 +78,11 @@ fHist.SetYTitle("y [\\circ]"); fHist.SetZTitle("Eq.cts"); + + fHist.SetDirectory(NULL); + fHistBg.SetDirectory(NULL); + fHistBg1.SetDirectory(NULL); + fHistBg2.SetDirectory(NULL); + + fHist.SetBit(TH1::kNoStats); } @@ -100,28 +106,4 @@ } - 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) @@ -131,18 +113,33 @@ } - //const Double_t xmax = fHist.GetXaxis()->GetXmax(); - - // Initialize all bins with a small (=0) value otherwise - // these bins are not displayed - for (int x=0; xGetBinCenter(x+1); - const Double_t cy = fHist.GetYaxis()->GetBinCenter(y+1); - //if (TMath::Hypot(cx, cy)FindObject("MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << "MHillasSrc not found... aborting." << endl; + return kFALSE; + } + + MBinning binsx, binsy; + binsx.SetEdges(fHist, 'x'); + binsy.SetEdges(fHist, 'y'); + MH::SetBinning(&fHistBg, &binsx, &binsy); + + if (!fHistOff) + { + MH::SetBinning(&fHistBg1, &binsx, &binsy); + MH::SetBinning(&fHistBg2, &binsx, &binsy); + } return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the delta angle between fSrcPos->GetXY() and v. +// Return result in deg. +// +Double_t MHDisp::DeltaPhiSrc(const TVector2 &v) const +{ + return TMath::Abs(fSrcPos->GetXY().DeltaPhi(v))*TMath::RadToDeg(); } @@ -165,79 +162,175 @@ rho = fPointPos->RotationAngle(*fObservatory, *fTime); - // FIXME: Do wobble-rotation when subtracting? - if (!fHistOff/* && fIsWobble*/) - rho += TMath::Pi(); - - // Get Disp from Parlist - const Double_t disp = fDisp->GetVal(); - - // Calculate both postiions where disp could point - TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta()); - TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta()); - pos1 *= disp; // Vector of length disp in direction of shower - pos2 *= -disp; // Vector of length -disp in direction of shower - - // Move origin of vector to center-of-gravity of shower - 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()*fMm2Deg) > 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(fMm2Deg, hsrc.GetCosDeltaAlpha()); - - gweight = m3l>fM3lCut ? 1 : 0; - } - - // Now we can safly derotate both position... - TVector2 srcp; - if (fSrcPos) - srcp = fSrcPos->GetXY(); - - // Derotate all position around camera center by -rho - if (rho!=0) - { + // Vector of length disp in direction of shower + const TVector2 p(hil->GetCosDelta(), hil->GetSinDelta()); + + // Move origin of vector to center-of-gravity of shower and derotate + TVector2 pos1 = hil->GetMean()*fMm2Deg + p*fDisp->GetVal(); + + Double_t w0 = 1; + if (fWobble) + { + const Double_t delta = DeltaPhiSrc(pos1); + + // Skip off-data not in the same half than the source (here: anti-source) + // Could be replaced by some kind of theta cut... + if (!fHistOff) + { + if (delta>165) + return kTRUE; + + // Because afterwards the plots are normalized with the number + // of entries the non-overlap region corresponds to half the + // contents, the overlap region to full contents. By taking + // the mean of both distributions we get the same result than + // we would have gotten using full off-events with a slightly + // increased uncertainty + // FIXME: The delta stuff could be replaced by a 2*antitheta cut... + w0 = delta>15 ? 1 : 2; + } + + // Define by the source position which histogram to fill + if (DeltaPhiSrc(fFormerSrc)>90) + fHalf = !fHalf; + fFormerSrc = fSrcPos->GetXY(); + } + + // If on-data: Derotate source and P. Translate P to center. + TVector2 m; + if (fHistOff) + { + // Derotate all position around camera center by -rho + // Move origin of vector to center-of-gravity of shower and derotate pos1=pos1.Rotate(-rho); - pos2=pos2.Rotate(-rho); - srcp=srcp.Rotate(-rho); - } - - // Shift the source position to 0/0 - pos1 -= srcp*fMm2Deg; - pos2 -= srcp*fMm2Deg; - - //const Double_t xmax = fHist.GetXaxis()->GetXmax(); - - // Fill histograms - //if (pos1.Mod()GetXY().Rotate(-rho)*fMm2Deg; + pos1 -= m; + } + } + + if (fSmearing<=0) + { + fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*w0); + return kTRUE; + } + + // ------------------------------------------------- + // The following algorithm may look complicated... + // It is optimized for speed + // ------------------------------------------------- + + // Define a vector used to calculate rotated x and y component + const TVector2 rot(TMath::Sin(rho), TMath::Cos(rho)); + + // Fold event position with the PSF and fill it into histogram + const TAxis &axex = *fHist.GetXaxis(); + const TAxis &axey = *fHist.GetYaxis(); + + const Int_t nx = axex.GetNbins(); + const Int_t ny = axey.GetNbins(); + + // normg: Normalization for Gauss [sqrt(2*pi)*sigma]^2 + const Double_t weight = axex.GetBinWidth(1)*axey.GetBinWidth(1)*w*w0; + const Double_t psf = 2*fSmearing*fSmearing; + const Double_t pi2 = fSmearing * TMath::Pi()/2; + const Double_t normg = pi2*pi2 * TMath::Sqrt(TMath::TwoPi()) / weight; + const Int_t bz = fHist.GetZaxis()->FindFixBin(0); + + TH1 &bg = fHalf ? fHistBg1 : fHistBg2; + + // To calculate significance map smear with 2*theta-cut and + // do not normalize gauss, for event map use theta-cut/2 instead + for (int x=1; x<=nx; x++) + { + const Double_t cx = axex.GetBinCenter(x); + const Double_t px = cx-pos1.X(); + + for (int y=1; y<=ny; y++) + { + const Double_t cy = axey.GetBinCenter(y); + const Double_t sp = Sq(px, cy-pos1.Y()); + const Double_t dp = sp/psf; + + // Values below 1e-3 (>3.5sigma) are not filled into the histogram + if (dp<4) + { + const Double_t rc = TMath::Exp(-dp)/normg; + if (!fHistOff) + bg.AddBinContent(bg.GetBin(x, y), rc); + else + fHist.AddBinContent(fHist.GetBin(x, y, bz), rc); + } + + if (!fHistOff) + continue; + + // If we are filling the signal already (fHistOff!=NULL) + // we also fill the background by projecting the + // background in the camera into the sky plot. + const TVector2 v = TVector2(cx+m.X(), cy+m.Y()); + + // Speed up further: Xmax-fWobble + if (v.Mod()>axex.GetXmax()) // Gains 10% speed + continue; + + const Int_t bx = axex.FindFixBin(v^rot); + const Int_t by = axey.FindFixBin(v*rot); + const Double_t bg = fHistOff->GetBinContent(bx, by, bz); + + fHistBg.AddBinContent(fHistBg.GetBin(x, y), bg); + } + } + fHist.SetEntries(fHist.GetEntries()+1); + + if (!fHistOff) + bg.SetEntries(fHistBg1.GetEntries()+1); return kTRUE; } +// -------------------------------------------------------------------------- +// +// Compile the background in camera coordinates from the two background +// histograms +// +Bool_t MHDisp::Finalize() +{ + if (fHistOff) + return kTRUE; + + const Int_t bz = fHist.GetZaxis()->FindFixBin(0); + + const Double_t n1 = fHistBg1.GetEntries(); + const Double_t n2 = fHistBg2.GetEntries(); + + for (int x=1; x<=fHist.GetNbinsX(); x++) + for (int y=1; y<=fHist.GetNbinsY(); y++) + { + const Int_t bin1 = fHistBg1.GetBin(x, y); + + const Double_t rc = + (n1==0?0:fHistBg1.GetBinContent(bin1)/n1)+ + (n2==0?0:fHistBg2.GetBinContent(bin1)/n2); + + fHist.SetBinContent(x, y, bz, rc/2); + } + + fHist.SetEntries(n1+n2); + + // Result corresponds to one smeared background event + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Return the mean signal in h around (0,0) in a distance between +// 0.33 and 0.47deg +// Double_t MHDisp::GetOffSignal(TH1 &h) const { @@ -245,17 +338,98 @@ const TAxis &axey = *h.GetYaxis(); + Int_t cnt = 0; Double_t sum = 0; for (int x=0; x0.35) + const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1)); + if (d>0.33 && d<0.47) + //if (d>0.4 && d<0.8) + { sum += h.GetBinContent(x+1,y+1); - } - - return sum; -} + cnt++; + } + } + + return sum/cnt; +} + +// -------------------------------------------------------------------------- +// +// Fill h2 with the radial profile of h1 around (x0, y0) +// +void MHDisp::Profile(TH1 &h2, const TH2 &h1, Axis_t x0, Axis_t y0) const +{ + const TAxis &axex = *h1.GetXaxis(); + const TAxis &axey = *h1.GetYaxis(); + + h2.Reset(); + + for (int x=1; x<=axex.GetNbins(); x++) + for (int y=1; y<=axey.GetNbins(); y++) + { + const Double_t dx = axex.GetBinCenter(x)-x0; + const Double_t dy = axey.GetBinCenter(y)-y0; + + const Double_t r = TMath::Hypot(dx, dy); + + h2.Fill(r, h1.GetBinContent(x, y)); + } +} + +// -------------------------------------------------------------------------- +// +// Remove contents of histogram around a circle. +// +void MHDisp::MakeDot(TH2 &h2) const +{ + const TAxis &axex = *h2.GetXaxis(); + const TAxis &axey = *h2.GetYaxis(); + + const Double_t rmax = fWobble ? axex.GetXmax()-0.4 : axex.GetXmax(); + + for (int x=1; x<=axex.GetNbins(); x++) + for (int y=1; y<=axey.GetNbins(); y++) + { + const Int_t bin = h2.GetBin(x,y); + + const Double_t px = h2.GetBinCenter(x); + const Double_t py = h2.GetBinCenter(y); + + if (rmaxSetRange(0,9999); - TH1 *h1=fHist.Project3D("yx_on"); + TH2 *h1=(TH2*)fHist.Project3D("yx_on"); // Set Glow-palette (PaintSurface doesn't allow more than 99 colors) @@ -275,68 +449,47 @@ // Project off data onto yx-plane and subtract it from on-data fHistOff->GetZaxis()->SetRange(0,9999); - TH1 *h=fHistOff->Project3D("yx_off"); + TH2 *h=(TH2*)fHistOff->Project3D("yx_off"); scale = -1; - //if (!fIsWobble) - { - const Double_t h1off = GetOffSignal(*h1); - const Double_t hoff = GetOffSignal(*h); - scale = hoff==0?-1:-h1off/hoff; - } - - h1->Add(h, scale); + const Double_t h1off = GetOffSignal(*h1); + const Double_t hoff = GetOffSignal(fHistBg); + + scale = hoff==0 ? -1 : -h1off/hoff; + + const Bool_t sig = kFALSE; + + if (!sig) + h1->Add(&fHistBg, scale); + else + MakeSignificance(*h1, *h1, fHistBg, scale); + + MakeDot(*h1); + delete h; - // Force calculation of minimum, maximum - h1->SetMinimum(); - h1->SetMaximum(); - - // Find maximum - const Double_t max = TMath::Max(TMath::Abs(h1->GetMinimum()), - TMath::Abs(h1->GetMaximum())); - - // Set new minimum, maximum - h1->SetMinimum(-max); - h1->SetMaximum( max); - } - - Int_t ix, iy, iz; - TF1 func("fcn", "gaus + [3]*x*x + [4]"); + MakeSymmetric(h1); + } pad->cd(3); TH1 *h2 = (TH1*)gPad->FindObject("RadProf"); - /* - if (!h2) - { - const Double_t maxr = TMath::Hypot(h1->GetXaxis()->GetXmax(), h1->GetYaxis()->GetXmax()); - const Int_t nbin = (h1->GetNbinsX()+h1->GetNbinsY())/2; - TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr); - //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr); - h->Sumw2(); - h->SetXTitle("\\vartheta [\\circ]"); - h->SetYTitle("/\\Delta R"); - h->SetBit(kCanDelete); - h->Draw(); - }*/ + + TString opt(o); + opt.ToLower(); + if (h1 && h2) { - h2->Reset(); - + Int_t ix, iy, iz; h1->GetMaximumBin(ix, iy, iz); const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix); const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy); - const Double_t w0 = h1->GetXaxis()->GetBinWidth(1); - - 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->GetYaxis()->GetBinCenter(y+1)-y0); - h2->Fill(r, h1->GetBinContent(x+1, y+1)); - } + //const Double_t w0 = h1->GetXaxis()->GetBinWidth(1); + + Profile(*h2, *h1, 0, 0); + //Profile(*h2, *h1, x0, y0); // Replace with MAlphaFitter? + TF1 func("fcn", "gaus + [3]*x*x + [4]"); func.SetLineWidth(1); func.SetLineColor(kBlue); @@ -347,49 +500,11 @@ func.FixParameter(1, 0); func.SetParameter(2, 0.15); - func.SetParameter(4, h2->GetBinContent(10)); - h2->Fit(&func, "IMQ", "", 0, 1.0); - - // No wintegrate the function f(x) per Delta Area - // which is f(x)/(pi*delta r*(2*r+delta r)) - TF1 func2("fcn2", Form("(gaus + [3]*x*x + [4])/(2*x+%.5f)", w0)); - for (int i=0; i<5; i++) - func2.SetParameter(i, func.GetParameter(i)); - - const Double_t r0 = 2*func2.GetParameter(2); - const Double_t e = func2.Integral(0, r0)/(w0*TMath::Pi()); - func2.SetParameter(0, 0); - const Double_t b = func2.Integral(0, r0)/(w0*TMath::Pi()); - const Double_t s = MMath::SignificanceLiMa(e, b); - - h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f f=%.2f", x0, y0, func.GetParameter(2), s, e-b, b, scale)); - } - /* - 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"); - }*/ + if (fHistOff) + func.FixParameter(3, 0); + func.SetParameter(4, h2->GetBinContent(15)); + h2->Fit(&func, "IQ", "", 0, 1.0); + + h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale))); + } } @@ -400,5 +515,5 @@ pad->SetBorderMode(0); - AppendPad(""); + AppendPad(o); TString name = Form("%s_1", pad->GetName()); @@ -416,5 +531,7 @@ h3->Draw("colz"); h3->SetBit(kCanDelete); -// catalog->Draw("mirror same *"); + + if (fHistOff) + GetCatalog()->Draw("white mirror same *"); pad->cd(); @@ -430,4 +547,5 @@ p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0); p->SetNumber(3); + p->SetGrid(); p->Draw(); p->cd(); @@ -449,33 +567,25 @@ // The following resources are available: // -// MHDisp.Wobble: on/off -// -/* -Int_t MHFalseSource::ReadEnv(const TEnv &env, TString prefix, Bool_t print) -{ - Bool_t rc = kFALSE; - if (IsEnvDefined(env, prefix, "DistMin", print)) +// MHDisp.Smearing: 0.1 +// MHDisp.Wobble: on/off +// +Int_t MHDisp::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Int_t rc = MHFalseSource::ReadEnv(env, prefix, print); + if (rc==kERROR) + return kERROR; + + if (IsEnvDefined(env, prefix, "Smearing", print)) { rc = kTRUE; - SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist)); - } - if (IsEnvDefined(env, prefix, "DistMax", print)) + SetSmearing(GetEnvValue(env, prefix, "Smearing", fSmearing)); + } + if (IsEnvDefined(env, prefix, "Wobble", print)) { rc = kTRUE; - SetMaxDist(GetEnvValue(env, prefix, "DistMax", fMaxDist)); - } - - if (IsEnvDefined(env, prefix, "DWMin", print)) - { - rc = kTRUE; - SetMinDW(GetEnvValue(env, prefix, "DWMin", fMinDist)); - } - if (IsEnvDefined(env, prefix, "DWMax", print)) - { - rc = kTRUE; - SetMaxDW(GetEnvValue(env, prefix, "DWMax", fMaxDist)); + SetWobble(GetEnvValue(env, prefix, "Wobble", fWobble)); } return rc; } -*/ + Index: /trunk/MagicSoft/Mars/mhflux/MHDisp.h =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7192) +++ /trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7193) @@ -5,7 +5,14 @@ #include "MHFalseSource.h" #endif +#ifndef ROOT_TH2 +#include +#endif +#ifndef ROOT_TVector2 +#include +#endif class MParList; class MHillasExt; +class MHillasSrc; class MParameterD; @@ -13,22 +20,58 @@ { private: - MHillasExt *fHilExt; //! - MParameterD *fDisp; //! + MHillasExt *fHilExt; //! + MHillasSrc *fHilSrc; //! + MParameterD *fDisp; //! - Double_t fM3lCut; //! - Double_t fXi; //! - Double_t fXiTheta; //! + TH2D fHistBg; - //Bool_t fIsWobble; + TH2D fHistBg1; + TH2D fHistBg2; + + TVector2 fFormerSrc; + Bool_t fHalf; + + Double_t fSmearing; + Bool_t fWobble; // MHDisp Double_t GetOffSignal(TH1 &h) const; + Double_t DeltaPhiSrc(const TVector2 &v) const; + + void Profile(TH1 &h2, const TH2 &h1, Axis_t x0=0, Axis_t y0=0) const; + void MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale=1) const; + void MakeDot(TH2 &h2) const; + + Double_t Sq(Double_t x, Double_t y) const { return x*x+y*y; } public: MHDisp(const char *name=NULL, const char *title=NULL); + // MHDisp + void SetSmearing(Double_t s=-1) { fSmearing=s; } + void SetWobble(Bool_t w=kTRUE) { fWobble=w; } + + // MHFalseSource (to be moved!) + void SetOffData(const MHFalseSource &fs) + { + MHFalseSource::SetOffData(fs); + if (dynamic_cast(&fs)) + { + const MHDisp &h = dynamic_cast(fs); + fWobble = h.fWobble; + fSmearing = h.fSmearing; + + h.fHistBg1.Copy(fHistBg1); + h.fHistBg2.Copy(fHistBg2); + + fHistBg1.SetDirectory(0); + fHistBg2.SetDirectory(0); + } + } + // MH Bool_t SetupFill(const MParList *pList); Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); // TObject @@ -37,7 +80,7 @@ // MParContainer - //Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); - ClassDef(MHDisp, 1) //3D-histogram in alpha, x and y + ClassDef(MHDisp, 2) //Class to provide a Disp map };