Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7216) +++ trunk/MagicSoft/Mars/Changelog (revision 7217) @@ -20,4 +20,21 @@ -*-*- END OF LINE -*-*- + + 2005/07/25 Thomas Bretz + + * mastro/MAstroSky2Local.[h,cc]: + - allow setting a pointing deviation when calculating the rotation angle + + * mhflux/MHDisp.[h,cc]: + - added pot for a significance map + - changed alignment range from 0.33/0.47 to 0.325/0.475 + - added fDeviation to calculation of Rotatio angle + - removed some obsolete code + + * mpointing/MPointingPos.[h,cc]: + - added option to calculate Rotation angle using MPointingDev + + + 2005/07/22 Daniela Dorner Index: trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc (revision 7216) +++ trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc (revision 7217) @@ -129,7 +129,11 @@ // seen with an Alt/Az telescope. // +// dzd and daz is a pointing offset, which is subtracted from the +// calculated local coordinates correspoding to time, observatory +// and ra/dec. +// // For more information see MAstro::RotationAngle // -Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec) const +Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec, Double_t dzd, Double_t daz) const { TVector3 v; @@ -137,4 +141,4 @@ v *= *this; - return MAstro::RotationAngle(ZZ(), XZ(), v.Theta(), v.Phi()); + return MAstro::RotationAngle(ZZ(), XZ(), v.Theta()-dzd, v.Phi()-daz); } Index: trunk/MagicSoft/Mars/mastro/MAstroSky2Local.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroSky2Local.h (revision 7216) +++ trunk/MagicSoft/Mars/mastro/MAstroSky2Local.h (revision 7217) @@ -18,5 +18,5 @@ MAstroSky2Local(const MTime &t, const MObservatory &obs); - Double_t RotationAngle(Double_t ra, Double_t dec) const; + Double_t RotationAngle(Double_t ra, Double_t dec, Double_t dzd=0, Double_t daz=0) const; ClassDef(MAstroSky2Local, 1) // Rotation Matrix to convert sky coordinates to ideal local coordinates Index: trunk/MagicSoft/Mars/mhflux/MHDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7216) +++ trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7217) @@ -51,8 +51,8 @@ #include "MParameters.h" -#include "MHillasExt.h" -#include "MHillasSrc.h" +#include "MHillas.h" #include "MSrcPosCam.h" #include "MPointingPos.h" +#include "MPointingDev.h" ClassImp(MHDisp); @@ -65,5 +65,5 @@ // MHDisp::MHDisp(const char *name, const char *title) - : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE) + : fDisp(0), fDeviation(0), fSmearing(-1), fWobble(kFALSE) { // @@ -106,17 +106,7 @@ } - fHilExt = (MHillasExt*)plist->FindObject("MHillasExt"); - if (!fHilExt) - { - *fLog << err << "MHillasExt not found... aborting." << endl; - return kFALSE; - } - - fHilSrc = (MHillasSrc*)plist->FindObject("MHillasSrc"); - if (!fHilSrc) - { - *fLog << err << "MHillasSrc not found... aborting." << endl; - return kFALSE; - } + fDeviation = (MPointingDev*)plist->FindObject("MPointingDev"); + if (!fDeviation) + *fLog << warn << "MPointingDev not found... ignored." << endl; MBinning binsx, binsy; @@ -160,5 +150,5 @@ Double_t rho = 0; if (fTime && fObservatory && fPointPos) - rho = fPointPos->RotationAngle(*fObservatory, *fTime); + rho = fPointPos->RotationAngle(*fObservatory, *fTime, fDeviation); // Vector of length disp in direction of shower @@ -347,5 +337,5 @@ // // Return the mean signal in h around (0,0) in a distance between -// 0.33 and 0.47deg +// 0.325 and 0.475deg // Double_t MHDisp::GetOffSignal(TH1 &h) const @@ -360,6 +350,5 @@ { 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) + if (d>0.325 && d<0.475) { sum += h.GetBinContent(x+1,y+1); @@ -367,5 +356,4 @@ } } - return sum/cnt; } @@ -415,4 +403,7 @@ if (rmaxn) + continue; + + const Int_t bin = s.GetBin(x+dx,y+dy); + + sig += h1.GetBinContent(bin); + bg += h2.GetBinContent(bin); + } + + const Double_t S = sig>0 ? MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale)) : 0; + + const Int_t bin = s.GetBin(x,y); s.SetBinContent(bin, S); } } +void MHDisp::Profile1D(const char *name, const TH2 &h) const +{ + if (!gPad) + return; + + TH1D *hrc = dynamic_cast(gPad->FindObject(name)); + if (!hrc) + return; + + hrc->Reset(); + + const Double_t max = h.GetMaximum(); + + MBinning(50, -max*1.1, max*1.1).Apply(*hrc); + + for (int x=1; x<=h.GetXaxis()->GetNbins(); x++) + for (int y=1; y<=h.GetXaxis()->GetNbins(); y++) + { + const Int_t bin = h.GetBin(x,y); + const Double_t sig = h.GetBinContent(bin); + if (sig!=0) + hrc->Fill(sig); + } + + gPad->SetLogy(hrc->GetMaximum()>0); + + if (!fHistOff || hrc->GetRMS()<0.1) + return; + + // ---------- Fix mean ---------- + // TF1 *g = (TF1*)gROOT->GetFunction("gaus"); + // g->FixParameter(1, 0); + // hrc->Fit("gaus", "BQI"); + + hrc->Fit("gaus", "QI"); + + TF1 *f = hrc->GetFunction("gaus"); + if (f) + { + f->SetLineWidth(1); + f->SetLineColor(kBlue); + } +} void MHDisp::Paint(Option_t *o) @@ -460,4 +507,6 @@ h1->SetContour(99); + TH2 *hx=0; + Double_t scale = 1; if (fHistOff) @@ -467,6 +516,4 @@ TH2 *h=(TH2*)fHistOff->Project3D("yx_off"); - scale = -1; - const Double_t h1off = GetOffSignal(*h1); const Double_t hoff = GetOffSignal(fHistBg); @@ -474,16 +521,19 @@ scale = hoff==0 ? -1 : -h1off/hoff; - const Bool_t sig = kFALSE; - - if (!sig) - h1->Add(&fHistBg, scale); - else - MakeSignificance(*h1, *h1, fHistBg, scale); + hx = (TH2*)pad->GetPad(4)->FindObject("Alpha_yx_sig"); + if (hx) + { + hx->SetContour(99); + MakeSignificance(*hx, *h1, fHistBg, scale); + MakeDot(*hx); + MakeSymmetric(hx); + } + + h1->Add(&fHistBg, scale); MakeDot(*h1); + MakeSymmetric(h1); delete h; - - MakeSymmetric(h1); } @@ -515,12 +565,20 @@ func.SetParameter(0, h2->GetBinContent(1)); func.FixParameter(1, 0); - func.SetParameter(2, 0.15); + func.SetParameter(2, 0.12); if (fHistOff) func.FixParameter(3, 0); - func.SetParameter(4, h2->GetBinContent(15)); + func.SetParameter(4, 0);//h2->GetBinContent(20)); 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))); } + + pad->cd(5); + if (h1) + Profile1D("Exc", *h1); + + pad->cd(6); + if (hx) + Profile1D("Sig", *hx); } @@ -533,11 +591,32 @@ AppendPad(o); - TString name = Form("%s_1", pad->GetName()); - TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0); + // ----- Pad number 4 ----- + TString name = Form("%s_4", pad->GetName()); + TPad *p = new TPad(name,name, 0.5025, 0.3355, 0.995, 0.995, col, 0, 0); + p->SetNumber(4); + p->Draw(); + p->cd(); + + TH1 *h3 = fHist.Project3D("yx_sig"); + h3->SetTitle("Significance Map"); + h3->SetDirectory(NULL); + h3->SetXTitle(fHist.GetXaxis()->GetTitle()); + h3->SetYTitle(fHist.GetYaxis()->GetTitle()); + h3->SetMinimum(0); + h3->Draw("colz"); + h3->SetBit(kCanDelete); + + if (fHistOff) + GetCatalog()->Draw("white mirror same *"); + + // ----- Pad number 1 ----- + pad->cd(); + name = Form("%s_1", pad->GetName()); + p = new TPad(name,name, 0.005, 0.3355, 0.4975, 0.995, col, 0, 0); p->SetNumber(1); p->Draw(); p->cd(); - TH1 *h3 = fHist.Project3D("yx_on"); + h3 = fHist.Project3D("yx_on"); h3->SetTitle("Distribution of equivalent events"); h3->SetDirectory(NULL); @@ -551,7 +630,8 @@ GetCatalog()->Draw("white mirror same *"); + // ----- Pad number 2 ----- 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 = new TPad(name,name, 0.005, 0.005, 0.2485, 0.3315, col, 0, 0); p->SetNumber(2); p->Draw(); @@ -559,7 +639,8 @@ h3->Draw("surf3"); + // ----- Pad number 3 ----- 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 = new TPad(name,name, 0.2525, 0.005, 0.4985, 0.3315, col, 0, 0); p->SetNumber(3); p->SetGrid(); @@ -577,4 +658,36 @@ h->SetBit(kCanDelete); h->Draw(); + + // ----- Pad number 5 ----- + pad->cd(); + name = Form("%s_5", pad->GetName()); + p = new TPad(name,name, 0.5025, 0.005, 0.7485, 0.3315, col, 0, 0); + p->SetNumber(5); + p->SetGrid(); + p->Draw(); + p->cd(); + + h3 = new TH1D("Exc", "Distribution of excess events", 10, -1, 1); + h3->SetDirectory(NULL); + h3->SetXTitle("N"); + h3->SetYTitle("Counts"); + h3->Draw(); + h3->SetBit(kCanDelete); + + // ----- Pad number 6 ----- + pad->cd(); + name = Form("%s_6", pad->GetName()); + p = new TPad(name,name, 0.7525, 0.005, 0.9995, 0.3315, col, 0, 0); + p->SetNumber(6); + p->SetGrid(); + p->Draw(); + p->cd(); + + h3 = new TH1D("Sig", "Distribution of significances", 10, -1, 1); + h3->SetDirectory(NULL); + h3->SetXTitle("N"); + h3->SetYTitle("Counts"); + h3->Draw(); + h3->SetBit(kCanDelete); } Index: trunk/MagicSoft/Mars/mhflux/MHDisp.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7216) +++ trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7217) @@ -16,14 +16,13 @@ class MHillasSrc; class MParameterD; +class MPointingDev; class MHDisp : public MHFalseSource { private: - MHillasExt *fHilExt; //! - MHillasSrc *fHilSrc; //! MParameterD *fDisp; //! + MPointingDev *fDeviation; //! TH2D fHistBg; - TH2D fHistBg1; TH2D fHistBg2; @@ -39,4 +38,5 @@ Double_t DeltaPhiSrc(const TVector2 &v) const; + void Profile1D(const char *name, const TH2 &h) 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; Index: trunk/MagicSoft/Mars/mpointing/MPointingPos.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MPointingPos.cc (revision 7216) +++ trunk/MagicSoft/Mars/mpointing/MPointingPos.cc (revision 7217) @@ -44,4 +44,5 @@ #include "MTime.h" #include "MObservatory.h" +#include "MPointingDev.h" #include "MAstroSky2Local.h" @@ -71,6 +72,8 @@ // For more information see MAstro::RotationAngle // -Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t) const +Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t, const MPointingDev *dev) const { - return MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad()); + return dev ? + MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad(), dev->GetDevZdRad(), dev->GetDevAzRad()): + MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad()); } Index: trunk/MagicSoft/Mars/mpointing/MPointingPos.h =================================================================== --- trunk/MagicSoft/Mars/mpointing/MPointingPos.h (revision 7216) +++ trunk/MagicSoft/Mars/mpointing/MPointingPos.h (revision 7217) @@ -13,4 +13,5 @@ class MTime; class MObservatory; +class MPointingDev; class MPointingPos : public MParContainer @@ -54,5 +55,5 @@ Double_t RotationAngle(const MObservatory &o) const; - Double_t RotationAngle(const MObservatory &o, const MTime &t) const; + Double_t RotationAngle(const MObservatory &o, const MTime &t, const MPointingDev *dev=0) const; Double_t RotationAngle(const MObservatory &o, const MTime *t) const {