Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6987) +++ trunk/MagicSoft/Mars/Changelog (revision 6988) @@ -21,4 +21,18 @@ -*-*- END OF LINE -*-*- + + 2005/05/02 Thomas Bretz + + * ganymed.rc: + - updated to the latest changes + + * mhflux/MHAlpha.[h,cc]: + - made a base class with the necessary interface to derived classes + to support other variables like dca or theta + + * mjobs/MJCut.[h,cc]: + - added support for MHAlpha derived classes + + 2005/04/29 Thomas Bretz Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 6987) +++ trunk/MagicSoft/Mars/NEWS (revision 6988) @@ -46,4 +46,9 @@ - MEventRateCalc handles the calculation of the event rate more accurate now in case of the start of a new run inside a sequence + + - ganymed got support for using other variables than Alpha, eg. Theta. + Therefor you need a class deriving from MHAlpha which supports + this variable (one is already existing: MHTheta). It is setup + through ganymed.rc Index: trunk/MagicSoft/Mars/ganymed.rc =================================================================== --- trunk/MagicSoft/Mars/ganymed.rc (revision 6987) +++ trunk/MagicSoft/Mars/ganymed.rc (revision 6988) @@ -47,5 +47,12 @@ #MAlphaFitter.ScaleMax: 80 #MAlphaFitter.PolynomOrder: 2 -#/*MISSING*/ ScaleMode +#MAlphaFitter.ScaleMode: kSignificance + +# ------------------------------------------------------------------------- +# Define here which histogram class to use to determin the signal. +# Currently availyble: MHAlpha, MHTheta +# ------------------------------------------------------------------------- +#MJCut.NameHist: MHAlpha + # ------------------------------------------------------------------------- @@ -66,5 +73,5 @@ # Get rid of triangular events... Cut0.Condition: {0} -Cut0.0: 0.33 * log10(MHillas.fSize) < 1.05 - log10(MNewImagePar.fConc) +Cut0.0: log10(MNewImagePar.fConc1) < (-0.467)*log10(MHillas.fSize) +0.75 # If you could setup MFEventSelector by resource file you could use it here @@ -92,22 +99,8 @@ # Setup cuts in false source plot # ------------------------------------------------------------------------- -MHFalseSource.DistMin: 0.55 +#MHFalseSource.DistMin: 0.55 #MHFalseSource.DistMax: 0.55 #MHFalseSource.DWMin: 0.55 #MHFalseSource.DWMax: 0.55 - -# ------------------------------------------------------------------------- -# Energy Estimation -# ------------------------------------------------------------------------- - -#EstimateEnergy: MEnergyEstimate -EstimateEnergy.Rule: {0} + {1} + {2} + {3} + {4} + {5} -EstimateEnergy.0: 6.3 -EstimateEnergy.1: 0.04*MHillasSrc.fDist*MGeomCam.fConvMm2Deg -EstimateEnergy.2: -0.13*MHillas.fLength*MGeomCam.fConvMm2Deg -EstimateEnergy.3: 0.15*MHillas.fSize -EstimateEnergy.4: 0.0000519*MHillas.fSize*MHillasSrc.fDist*MGeomCam.fConvMm2Deg -EstimateEnergy.5: 0.0000519*MHillas.fLength*MGeomCam.fConvMm2Deg*MHillasSrc.fDist*MGeomCam.fConvMm2Deg - # ------------------------------------------------------------------------- Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 6987) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 6988) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 3/2004 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2005 ! ! @@ -81,5 +81,6 @@ // MHAlpha::MHAlpha(const char *name, const char *title) - : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0), + : fNameParameter("MHillasSrc"), fParameter(0), + fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0), fPointPos(0), fTimeEffOn(0), fTime(0), fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE), @@ -93,18 +94,18 @@ 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(); + fHist.SetName("Alpha"); + fHist.SetTitle("Alpha"); + fHist.SetXTitle("\\Theta [deg]"); + //fHist.SetYTitle("E_{est} [GeV]"); + fHist.SetZTitle("|\\alpha| [\\circ]"); + fHist.SetDirectory(NULL); + fHist.UseCurrentStyle(); // Main histogram - fHAlphaTime.SetName("Alpha"); - fHAlphaTime.SetXTitle("|\\alpha| [\\circ]"); - fHAlphaTime.SetYTitle("Counts"); - fHAlphaTime.UseCurrentStyle(); - fHAlphaTime.SetDirectory(NULL); + fHistTime.SetName("Alpha"); + fHistTime.SetXTitle("|\\alpha| [\\circ]"); + fHistTime.SetYTitle("Counts"); + fHistTime.UseCurrentStyle(); + fHistTime.SetDirectory(NULL); @@ -142,44 +143,12 @@ 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]*x");// - f.SetParameter(0, -2.1); - f.SetParameter(1, 1400); // 50 - f.SetParameter(2, fHEnergy.Integral()); // 50 - f.SetLineWidth(2); - - fHEnergy.Fit(&f, "0QI", "", 100, 2400); // Integral? 90, 2800 - - const Double_t chi2 = f.GetChisquare(); - const Int_t NDF = f.GetNDF(); - - // was fit successful ? - const Bool_t ok = NDF>0 && chi2<2.5*NDF; - f.SetLineColor(ok ? kGreen : kRed); - - if (paint) - { - f.Paint("same"); - - if (ok) - { - TString s; - s += Form("\\alpha=%.1f ", f.GetParameter(0)); - s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000); - s += Form("N=%.1e", f.GetParameter(2)); - - TLatex text(0.25, 0.94, s.Data()); - text.SetBit(TLatex::kTextNDC); - text.SetTextSize(0.04); - text.Paint(); - } - }*/ + binsa.Apply(fHistTime); + + MH::SetBinning(&fHist, &binst, &binse, &binsa); +} + +Double_t MHAlpha::GetVal() const +{ + return static_cast(fParameter)->GetAlpha(); } @@ -189,5 +158,5 @@ MAlphaFitter fit(fFit); - const Int_t n = fHAlpha.GetNbinsY(); + const Int_t n = fHist.GetNbinsY(); fHEnergy.SetEntries(0); @@ -197,5 +166,5 @@ for (int i=1; i<=n; i++) { - if (fit.FitEnergy(fHAlpha, fOffData, i)) + if (fit.FitEnergy(fHist, fOffData, i)) { fHEnergy.SetBinContent(i, fit.GetEventsExcess()); @@ -220,5 +189,5 @@ MAlphaFitter fit(fFit); - const Int_t n = fHAlpha.GetNbinsX(); + const Int_t n = fHist.GetNbinsX(); fHTheta.SetEntries(0); @@ -226,5 +195,5 @@ for (int i=1; i<=n; i++) { - if (fit.FitTheta(fHAlpha, fOffData, i)) + if (fit.FitTheta(fHist, fOffData, i)) { fHTheta.SetBinContent(i, fit.GetEventsExcess()); @@ -237,22 +206,37 @@ } +Bool_t MHAlpha::GetParameter(const MParList &pl) +{ + fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc"); + if (fParameter) + return kTRUE; + + *fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl; + return kFALSE; +} + Bool_t MHAlpha::SetupFill(const MParList *pl) { - fHAlpha.Reset(); + fHist.Reset(); fHEnergy.Reset(); fHTheta.Reset(); fHTime.Reset(); - if (fName!=(TString)"MHAlphaOff" && fOffData==NULL) - { - MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha"); + const TString off(Form("%sOff", ClassName())); + if (fName!=off && fOffData==NULL) + { + const TString desc(Form("%s [%s] found... using ", off.Data(), ClassName())); + MHAlpha *hoff = (MHAlpha*)pl->FindObject(off, ClassName()); if (!hoff) - *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl; + *fLog << inf << "No " << desc << "current data only!" << endl; else { - *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl; + *fLog << inf << desc << "on-off mode!" << endl; SetOffData(*hoff); } } + + if (!GetParameter(*pl)) + return kFALSE; fHillas = 0; @@ -271,5 +255,5 @@ fHEnergy.SetTitle(" N_{exc} vs. Size "); fHEnergy.SetXTitle("Size [phe]"); - fHAlpha.SetYTitle("Size [phe]"); + fHist.SetYTitle("Size [phe]"); } else @@ -277,5 +261,5 @@ fHEnergy.SetTitle(" N_{exc} vs. E_{est} "); fHEnergy.SetXTitle("E_{est} [GeV]"); - fHAlpha.SetYTitle("E_{est} [GeV]"); + fHist.SetYTitle("E_{est} [GeV]"); } @@ -306,7 +290,7 @@ MBinning binst, binse, binsa; - binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x'); - binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y'); - binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z'); + binst.SetEdges(fOffData ? *fOffData : fHist, 'x'); + binse.SetEdges(fOffData ? *fOffData : fHist, 'y'); + binsa.SetEdges(fOffData ? *fOffData : fHist, 'z'); if (!fOffData) { @@ -324,5 +308,5 @@ binse.SetEdges(*pl, "BinningSize"); - binsa.SetEdges(*pl, "BinningAlpha"); + binsa.SetEdges(*pl, Form("Binning%s", ClassName()+2)); } else @@ -330,11 +314,11 @@ fHEnergy.SetTitle(fOffData->GetTitle()); fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle()); - fHAlpha.SetYTitle(fOffData->GetYaxis()->GetTitle()); + fHist.SetYTitle(fOffData->GetYaxis()->GetTitle()); } binse.Apply(fHEnergy); binst.Apply(fHTheta); - binsa.Apply(fHAlphaTime); - MH::SetBinning(&fHAlpha, &binst, &binse, &binsa); + binsa.Apply(fHistTime); + MH::SetBinning(&fHist, &binst, &binse, &binsa); MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter"); @@ -402,5 +386,5 @@ TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0; - const Bool_t rc = fit.ScaleAndFit(fHAlphaTime, h); + const Bool_t rc = fit.ScaleAndFit(fHistTime, h); if (h) delete h; @@ -410,6 +394,6 @@ // Reset Histogram - fHAlphaTime.Reset(); - fHAlphaTime.SetEntries(0); + fHistTime.Reset(); + fHistTime.SetEntries(0); // @@ -469,12 +453,5 @@ else { - const MHillasSrc *hil = dynamic_cast(par); - if (!par) - { - *fLog << err << dbginf << "MHillasSrc not found... abort." << endl; - return kFALSE; - } - - alpha = hil->GetAlpha(); + alpha = GetVal(); if (fHillas) @@ -491,5 +468,5 @@ // Fill histograms - fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w); + fHist.Fill(theta, energy, TMath::Abs(alpha), w); // Check cuts @@ -503,5 +480,5 @@ if (!fSkipHistTime) - fHAlphaTime.Fill(TMath::Abs(alpha), w); + fHistTime.Fill(TMath::Abs(alpha), w); return kTRUE; @@ -557,8 +534,8 @@ padsave->cd(1); - TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha"); + TH1D *hon = fHist.ProjectionZ("Proj"); if (fOffData) { - TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff"); + TH1D *hoff = fOffData->ProjectionZ("ProjOff"); const Double_t alpha = fFit.Scale(*hoff, *hon); @@ -569,5 +546,5 @@ hoff->SetMaximum(alpha); - if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff"))) + if ((h0=(TH1D*)gPad->FindObject("ProjOnOff"))) { h0->Reset(); @@ -583,9 +560,9 @@ } - if (o==(TString)"alpha") - if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha"))) + if (o==(TString)"variable") + if ((h0 = (TH1D*)gPad->FindObject("Proj"))) { // Check whether we are dealing with on-off analysis - TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff"); + TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff"); // BE CARFEULL: This is a really weird workaround! @@ -603,8 +580,8 @@ if (o==(TString)"theta") { - TH1 *h = (TH1*)gPad->FindObject("Alpha_x"); + TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName())); if (h) { - TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x"); + TH1D *h2 = (TH1D*)fHist.Project3D("dum_x"); h2->SetDirectory(0); h2->Scale(fHTheta.Integral()/h2->Integral()); @@ -618,8 +595,8 @@ if (o==(TString)"energy") { - TH1 *h = (TH1*)gPad->FindObject("Alpha_y"); + TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName())); if (h) { - TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y"); + TH1D *h2 = (TH1D*)fHist.Project3D("dum_y"); h2->SetDirectory(0); h2->Scale(fHEnergy.Integral()/h2->Integral()); @@ -634,5 +611,4 @@ gPad->SetLogy(); } - FitEnergySpec(kTRUE); } @@ -660,8 +636,8 @@ gPad->SetBorderMode(0); - h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E"); + h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E"); h->SetBit(TH1::kNoTitle); h->SetStats(kTRUE); - h->SetXTitle("\\alpha [\\circ]"); + h->SetXTitle(fHist.GetZaxis()->GetTitle()); h->SetYTitle("Counts"); h->SetDirectory(NULL); @@ -674,7 +650,7 @@ h->SetMarkerColor(kGreen); - h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E"); + h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E"); h->SetBit(TH1::kNoTitle); - h->SetXTitle("\\alpha [\\circ]"); + h->SetXTitle(fHist.GetZaxis()->GetTitle()); h->SetYTitle("Counts"); h->SetDirectory(NULL); @@ -684,7 +660,7 @@ h->Draw("same"); - h = (TH1D*)h->Clone("ProjAlphaOnOff"); + h = (TH1D*)h->Clone("ProjOnOff"); h->SetBit(TH1::kNoTitle); - h->SetXTitle("\\alpha [\\circ]"); + h->SetXTitle(fHist.GetZaxis()->GetTitle()); h->SetYTitle("Counts"); h->SetDirectory(NULL); @@ -696,10 +672,10 @@ TLine lin; lin.SetLineStyle(kDashed); - lin.DrawLine(0, 0, 90, 0); + lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0); } // After the Alpha-projection has been drawn. Fit the histogram // and paint the result into this pad - AppendPad("alpha"); + AppendPad("variable"); if (fHEnergy.GetNbinsX()>1) @@ -713,5 +689,5 @@ AppendPad("energy"); - h = (TH1D*)fHAlpha.Project3D("y"); + h = (TH1D*)fHist.Project3D("y"); h->SetBit(TH1::kNoTitle|TH1::kNoStats); h->SetXTitle("E [GeV]"); @@ -749,5 +725,5 @@ AppendPad("theta"); - h = (TH1D*)fHAlpha.Project3D("x"); + h = (TH1D*)fHist.Project3D("x"); h->SetBit(TH1::kNoTitle|TH1::kNoStats); h->SetXTitle("\\theta [\\circ]"); @@ -768,5 +744,5 @@ // FIXME: Do in Paint if special option given! TCanvas *c = new TCanvas; - Int_t n = fHAlpha.GetNbinsY(); + Int_t n = fHist.GetNbinsY(); Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1); c->Divide(nc, nc, 1e-10, 1e-10); @@ -775,12 +751,12 @@ MAlphaFitter fit(fFit); - for (int i=1; i<=fHAlpha.GetNbinsY(); i++) + for (int i=1; i<=fHist.GetNbinsY(); i++) { c->cd(i); - TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E"); + TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E"); hon->SetBit(TH1::kNoTitle); hon->SetStats(kTRUE); - hon->SetXTitle("\\alpha [\\circ]"); + hon->SetXTitle(fHist.GetZaxis()->GetTitle()); hon->SetYTitle("Counts"); hon->SetDirectory(NULL); @@ -796,7 +772,7 @@ hon->SetMarkerColor(kGreen); - hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E"); + hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E"); hof->SetBit(TH1::kNoTitle|TH1::kNoStats); - hof->SetXTitle("\\alpha [\\circ]"); + hof->SetXTitle(fHist.GetZaxis()->GetTitle()); hof->SetYTitle("Counts"); hof->SetDirectory(NULL); @@ -814,5 +790,5 @@ diff->Add(hof, -1); diff->SetBit(TH1::kNoTitle); - diff->SetXTitle("\\alpha [\\circ]"); + diff->SetXTitle(fHist.GetZaxis()->GetTitle()); diff->SetYTitle("Counts"); diff->SetDirectory(NULL); @@ -824,5 +800,5 @@ TLine lin; lin.SetLineStyle(kDashed); - lin.DrawLine(0, 0, 90, 0); + lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0); const Float_t min = diff->GetMinimum()*1.05; @@ -833,5 +809,5 @@ *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl; /* - if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE)) + if (fit.FitEnergy(fHist, fOffData, i, kTRUE)) { fHEnergy.SetBinContent(i, fit.GetEventsExcess()); @@ -845,10 +821,10 @@ Bool_t MHAlpha::Finalize() { - //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E"); + //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E"); //h->SetDirectory(0); //Bool_t rc = fFit.Fit(*h); //delete h; - if (!fFit.FitAlpha(fHAlpha, fOffData)) + if (!fFit.FitAlpha(fHist, fOffData)) { *fLog << warn << "MAlphaFitter - Fit failed..." << endl; @@ -900,5 +876,5 @@ fMatrix = mat; - fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha"); + fMap[0] = fMatrix->AddColumn(GetParameterRule()); fMap[1] = -1; fMap[2] = -1; Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 6987) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.h (revision 6988) @@ -24,63 +24,14 @@ 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 { +protected: + TH3D fHist; // Alpha vs. theta and energy + TH1D fHistTime; //! temporary histogram to get alpha vs. time + + TString fNameParameter; + + MParContainer *fParameter; //! + private: const TH3D *fOffData; @@ -88,14 +39,12 @@ MAlphaFitter fFit; // SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT) - 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; //! - MParameterD *fEnergy; //! - MHillas *fHillas; //! - MPointingPos *fPointPos; //! + MParameterD *fResult; //! + MParameterD *fEnergy; //! + MHillas *fHillas; //! + MPointingPos *fPointPos; //! MTime *fTimeEffOn; //! Time to steer filling of fHTime @@ -116,5 +65,4 @@ Int_t fMap[5]; //! - void FitEnergySpec(Bool_t paint=kFALSE); Float_t FitEnergyBins(Bool_t paint=kFALSE); void FitThetaBins(Bool_t paint=kFALSE); @@ -128,19 +76,30 @@ Int_t DistancetoPrimitive(Int_t px, Int_t py) { return 0; } + virtual Bool_t GetParameter(const MParList &pl); + virtual Double_t GetVal() const; + virtual const char *GetParameterRule() const + { + return "MHillasSrc.fAlpha"; + } + public: MHAlpha(const char *name=NULL, const char *title=NULL); + // MH 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; } + // MHAlpha + const TH3D &GetHist() const { return fHist; } const MAlphaFitter &GetAlphaFitter() const { return fFit; } const TH1D &GetHEnergy() const { return fHEnergy; } + // Setter + void SetNameParameter(const char *name) { fNameParameter=name; } void SetOffData(const MHAlpha &h) { - fOffData = &h.fHAlpha; + fOffData = &h.fHist; fForceUsingSize = h.fForceUsingSize; fNumTimeBins = h.fNumTimeBins; Index: trunk/MagicSoft/Mars/mhflux/MHTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHTheta.cc (revision 6988) +++ trunk/MagicSoft/Mars/mhflux/MHTheta.cc (revision 6988) @@ -0,0 +1,103 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHTheta +// +// This is a MHAlpha using "ThetaSquared [MParameterD]" as 'alpha' +// +// For more detailes see MHAlpha. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHTheta.h" + +#include "MParameters.h" +#include "MHMatrix.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHTheta); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MHTheta::MHTheta(const char *name, const char *title) + : MHAlpha(name, title), fThetaSq(0) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHTheta"; + fTitle = title ? title : "Theta Squared plot"; + + fNameParameter = "ThetaSquared"; + + fHist.SetName("Theta"); + fHist.SetTitle("Theta"); + fHist.SetZTitle("\\theta^{2} [deg^{2}]"); + fHist.SetDirectory(NULL); + + // Main histogram + fHistTime.SetName("Theta"); + fHistTime.SetXTitle("\\theta^{2} [deg^{2}]"); + fHistTime.SetDirectory(NULL); + + //fHistTime.SetYTitle("\\theta^{2] [deg^{2}]"); + /* + TArrayD arr(50); + for (int i=1; i<50; i++) + arr[i] = sqrt((arr[i-1]+1)*(arr[i-1]+1) + 1.1)-1; + */ + MBinning binsa, binse, binst; + binsa.SetEdges(75, 0, 1.5); + //binsa.SetEdges(arr); + binse.SetEdgesLog(15, 10, 100000); + binst.SetEdgesCos(50, 0, 60); + binsa.Apply(fHistTime); + + MH::SetBinning(&fHist, &binst, &binse, &binsa); +} + +Bool_t MHTheta::GetParameter(const MParList &pl) +{ + fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MParameterD"); + if (fParameter) + return kTRUE; + + *fLog << err << fNameParameter << " [MParameterD] not found... abort." << endl; + return kFALSE; +} + +Double_t MHTheta::GetVal() const +{ + return static_cast(fParameter)->GetVal(); +} Index: trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6987) +++ trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6988) @@ -77,6 +77,6 @@ MJCut::MJCut(const char *name, const char *title) : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE), - fIsWobble(kFALSE), fIsMonteCarlo(kFALSE), fFullDisplay(kFALSE), /*fSubstraction(kFALSE),*/ - /*fEstimateEnergy(0),*/ fCalcHadronness(0) + fIsWobble(kFALSE), fIsMonteCarlo(kFALSE), fFullDisplay(kFALSE), + fNameHist("MHTheta"), fCalcHadronness(0) { fName = name ? name : "MJCut"; @@ -258,5 +258,5 @@ // Save also the result, not only the setup - const MHAlpha *halpha = (MHAlpha*)plist.FindObject("MHAlpha"); + const MHAlpha *halpha = (MHAlpha*)plist.FindObject(fNameHist, "MHAlpha"); if (halpha) arr.Add((TObject*)(&halpha->GetAlphaFitter())); @@ -285,4 +285,5 @@ // MJCut.WriteResult: yes, no // MJCut.ResultFile: filename +// MJCut.HistName: MHAlpha // Bool_t MJCut::CheckEnvLocal() @@ -301,4 +302,6 @@ EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay)); //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction)); + + SetNameHist(GetEnv("NameHist", fNameHist)); return kTRUE; @@ -341,15 +344,33 @@ } -void MJCut::DisplayResult(const MParList &plist) const -{ - /* - TObject *h1 = plist.FindObject("MHHillasOffPre", "MHHillas"); - TObject *h2 = plist.FindObject("MHHillasOffPost", "MHHillas"); - TObject *h3 = plist.FindObject("MHVsSizeOffPost", "MHVsSize"); - TObject *h4 = plist.FindObject("MHHilExtOffPost", "MHHillasExt"); - TObject *h5 = plist.FindObject("MHHilSrcOffPost", "MHHillasSrc"); - TObject *h6 = plist.FindObject("MHImgParOffPost", "MHImagePar"); - TObject *h7 = plist.FindObject("MHNewParOffPost", "MHNewImagePar"); - */ +// -------------------------------------------------------------------------- +// +// Create a new instance of an object with name name of class +// type fNameHist in parlist. It must derive from MHAlpha. +// Call ForceUsingSize for it and return its pointer. +// If something fails NULL is returned. +// +MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const +{ + TClass *cls = gROOT->GetClass(fNameHist); + if (!cls) + { + *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl; + return NULL; + } + if (!cls->InheritsFrom(MHAlpha::Class())) + { + *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl; + return NULL; + } + + const TString objname(Form("%s%s", fNameHist.Data(), name)); + MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname); + if (!h) + return NULL; + + h->ForceUsingSize(); + + return h; } @@ -424,9 +445,6 @@ fit.SetScaleMode(MAlphaFitter::kNone); - MHAlpha halphaoff("MHAlphaOff"); - halphaoff.ForceUsingSize(); - plist.AddToList(&halphaoff); - - MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha"); + MHAlpha *halphaoff = CreateNewHist(plist, "Off"); + MFillH falpha(halphaoff, "MHillasSrc", "FillAlpha"); MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS"); @@ -655,9 +673,6 @@ } */ - MHAlpha halphaon("MHAlpha"); - halphaon.ForceUsingSize(); - plist.AddToList(&halphaon); - - MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha"); + MHAlpha *halphaon=CreateNewHist(plist); + MFillH falpha2(halphaon, "MHillasSrc", "FillAlpha"); MFillH ffs2("MHFalseSource", "MHillas", "FillFS"); @@ -706,6 +721,4 @@ tlist.PrintStatistics(); - DisplayResult(plist); - // FIXME: Perform fit and plot energy dependant alpha plots // and fit result to new tabs! Index: trunk/MagicSoft/Mars/mjobs/MJCut.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6987) +++ trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6988) @@ -9,4 +9,5 @@ class MDataSet; class MParList; +class MHAlpha; class MWriteRootFile; @@ -26,16 +27,18 @@ TString fNameOutput; + TString fNameHist; + //MTask *fEstimateEnergy; MTask *fCalcHadronness; - TString GetOutputFile(UInt_t num) const; - Bool_t CheckEnvLocal(); - void SetupWriter(MWriteRootFile *write, const char *name) const; - Bool_t WriteTasks(UInt_t num, TObjArray &cont) const; - Bool_t WriteResult(const MParList &plist, UInt_t num) const; - void DisplayResult(const MParList &plist) const; + TString GetOutputFile(UInt_t num) const; + Bool_t CheckEnvLocal(); + void SetupWriter(MWriteRootFile *write, const char *name) const; + Bool_t WriteTasks(UInt_t num, TObjArray &cont) const; + Bool_t WriteResult(const MParList &plist, UInt_t num) const; + MHAlpha *CreateNewHist(MParList &plist, const char *name="") const; - Bool_t CanStoreSummary() const { return !fPathOut.IsNull() && fStoreSummary; } - Bool_t CanStoreResult() const { return !fPathOut.IsNull() && fStoreResult; } + Bool_t CanStoreSummary() const { return !fPathOut.IsNull() && fStoreSummary; } + Bool_t CanStoreResult() const { return !fPathOut.IsNull() && fStoreResult; } public: @@ -57,4 +60,6 @@ void SetNameOutFile(const char *name="") { fNameOutput=name; } + void SetNameHist(const char *name) { fNameHist=name; } + //void SetEnergyEstimator(const MTask *task=0); void SetHadronnessCalculator(const MTask *task=0);