/* ======================================================================== *\ ! ! * ! * 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, 03/2004 ! ! Copyright: MAGIC Software Development, 2002-2004 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MAstroCatalog // // THIS IMPLEMENTATION IS PRELIMINARY AND WILL BE MERGED WITH // SOME PARTS OF THE DRIVE SOFTWARE SOON! // ////////////////////////////////////////////////////////////////////////////// #include "MAstroCatalog.h" #include #include #include #include // INT_MAX (Suse 7.3/gcc 2.95) #include #include // TPad::GetMaxPickDistance #include #include #include #include #include #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MTime.h" #include "MAstro.h" #include "MAstroSky2Local.h" #include "MObservatory.h" ClassImp(MVector3); ClassImp(MAstroCatalog); using namespace std; /* class MRotation : public TRotation { public: MRotation(Double_t gmst, const MObservatory &obs) : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1) { RotateZ(gmst + obs.GetElong()); RotateY(obs.GetPhi()-TMath::Pi()/2); RotateZ(TMath::Pi()); } MRotation(const MTime &t, const MObservatory &obs) : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1) { RotateZ(t.GetGmst() + obs.GetElong()); RotateY(obs.GetPhi()-TMath::Pi()/2); RotateZ(TMath::Pi()); } }; */ /* MVector3 MVector3::GetZdAz(const MObservatory &obs, Double_t gmst) const { if (!fType==kIsRaDec) return MVector3(); MVector3 v(*this); v *= MAstroSky2Local(gmst, obs); return v; // ------ Using vectors ------- // v(1) = -v(1); // phi --> -phi // v.RotateZ(gmst + obs.GetElong()); // -phi --> alpha-phi // v.RotateY(obs.GetPhi()-TMath::Pi()/2); // v.RotateZ(TMath::Pi()); // ------ The same using slalib, tested in the drive system ------- // const Double_t alpha = slaGmst(mjd) + obs.GetElong(); // Double_t el; // slaDe2h(fAlpha-ra, dec, obs.GetPhi(), &az, &el); // zd = TMath::Pi()/2-el; } MVector3 MVector3::GetZdAz(const MTime &time, MObservatory &obs) const { return GetZdAz(obs, time.GetGmst()); } MVector3 MVector3::GetRaDec(const MObservatory &obs, Double_t gmst) const { if (!fType==kIsZdAz) return MVector3(); MVector3 v(*this); v *= MAstroSky2Local(gmst, obs).Inverse(); return v; // ------ Using vectors ------- // v.RotateZ(-TMath::Pi()); // v.RotateY(TMath::Pi()/2-obs.GetPhi()); // v.RotateZ(-gmst - obs.GetElong()); // alpha-phi --> -phi // v(1) = -v(1); // -phi --> phi // ------ The same using slalib, tested in the drive system ------- // const Double_t alpha = slaGmst(mjd) + obs.GetElong(); // Double_t el; // slaDe2h(fAlpha-ra, dec, obs.GetPhi(), &az, &el); // zd = TMath::Pi()/2-el; } MVector3 MVector3::GetRaDec(const MTime &time, MObservatory &obs) const { return GetRaDec(obs, time.GetGmst()); } */ // -------------------------------------------------------------------------- MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(90), fToolTip(0), fObservatory(0), fTime(0) { fList.SetOwner(); fToolTip = new TGToolTip(0, "", 0); } // -------------------------------------------------------------------------- MAstroCatalog::~MAstroCatalog() { if (fTime) delete fTime; if (fObservatory) delete fObservatory; fToolTip->Hide(); delete fToolTip; DeleteMap(); // FIXME: There must be an easier way! TIter Next(gROOT->GetListOfCanvases()); TCanvas *c; while ((c=(TCanvas*)Next())) c->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this, "EventInfo(Int_t,Int_t,Int_t,TObject*)"); } // -------------------------------------------------------------------------- TString MAstroCatalog::FindToken(TString &line, Char_t tok) { Ssiz_t token = line.First(tok); if (token<0) { const TString copy(line); line = ""; return copy; } const TString res = line(0, token); line.Remove(0, token+1); return res; } // -------------------------------------------------------------------------- Int_t MAstroCatalog::atoi(const TSubString &sub) { return atoi(TString(sub)); } // -------------------------------------------------------------------------- Float_t MAstroCatalog::atof(const TSubString &sub) { return atof(TString(sub)); } // -------------------------------------------------------------------------- Int_t MAstroCatalog::atoi(const TString &s) { return std::atoi(s); } // -------------------------------------------------------------------------- Float_t MAstroCatalog::atof(const TString &s) { return std::atof(s); } // -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadXephem(TString catalog) { SetBit(kHasChanged); gLog << inf << "Reading Xephem catalog: " << catalog << endl; ifstream fin(catalog); if (!fin) { gLog << err << "Cannot open file " << catalog << ": "; gLog << strerror(errno) << endl; return 0; } Int_t add =0; Int_t cnt =0; Int_t line=0; Double_t maxmag=0; while (1) { TString row; row.ReadLine(fin); if (!fin) break; line++; if (row[0]=='#') continue; TString line(row); TString name = FindToken(line); TString dummy = FindToken(line); TString r = FindToken(line); TString d = FindToken(line); TString m = FindToken(line); TString epoch = FindToken(line); if (name.IsNull() || r.IsNull() || d.IsNull() || m.IsNull() || epoch.IsNull()) { gLog << warn << "Invalid Entry Line #" << line << ": " << row << endl; continue; } cnt++; const Double_t mag = atof(m); maxmag = TMath::Max(maxmag, mag); if (mag>fLimMag) continue; if (epoch!="2000") { gLog << warn << "Epoch != 2000... skipped." << endl; continue; } Double_t ra0, dec0; MAstro::Coordinate2Angle(r, ra0); MAstro::Coordinate2Angle(d, dec0); ra0 *= TMath::Pi()/12; dec0 *= TMath::Pi()/180; MVector3 *star=new MVector3; star->SetRaDec(ra0, dec0, mag); star->SetName(name); if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV) { delete star; continue; } fList.Add(star); add++; } gLog << inf << "Read " << add << " out of " << cnt << " (Total max mag=" << maxmag << ")" << endl; return add; } // -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadNGC2000(TString catalog) { SetBit(kHasChanged); gLog << inf << "Reading NGC2000 catalog: " << catalog << endl; ifstream fin(catalog); if (!fin) { gLog << err << "Cannot open file " << catalog << ": "; gLog << strerror(errno) << endl; return 0; } Int_t add=0; Int_t cnt=0; Int_t n =0; Double_t maxmag=0; while (1) { TString row; row.ReadLine(fin); if (!fin) break; cnt++; const Int_t rah = atoi(row(13, 2)); const Float_t ram = atof(row(16, 4)); const Char_t decs = row(22); const Int_t decd = atoi(row(23, 2)); const Int_t decm = atoi(row(26, 2)); const TString m = row(43, 4); if (m.Strip().IsNull()) continue; n++; const Double_t mag = atof(m); maxmag = TMath::Max(maxmag, mag); if (mag>fLimMag) continue; const Double_t ra = MAstro::Hms2Rad(rah, (int)ram, fmod(ram, 1)*60); const Double_t dec = MAstro::Dms2Rad(decd, decm, 0, decs); MVector3 *star=new MVector3; star->SetRaDec(ra, dec, mag); star->SetName(row(0, 8)); if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV) { delete star; continue; } fList.Add(star); add++; } gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl; return add; } // -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadBSC(TString catalog) { SetBit(kHasChanged); gLog << inf << "Reading Bright Star Catalog (BSC5) catalog: " << catalog << endl; ifstream fin(catalog); if (!fin) { gLog << err << "Cannot open file " << catalog << ": "; gLog << strerror(errno) << endl; return 0; } Int_t add=0; Int_t cnt=0; Int_t n =0; Double_t maxmag=0; while (1) { TString row; row.ReadLine(fin); if (!fin) break; cnt++; const Int_t rah = atoi(row(75, 2)); const Int_t ram = atoi(row(77, 2)); const Float_t ras = atof(row(79, 4)); const Char_t decsgn = row(83); const Int_t decd = atoi(row(84, 2)); const Int_t decm = atoi(row(86, 2)); const Int_t decs = atoi(row(88, 2)); const TString m = row(102, 5); if (m.Strip().IsNull()) continue; n++; const Double_t mag = atof(m.Data()); maxmag = TMath::Max(maxmag, mag); if (mag>fLimMag) continue; const Double_t ra = MAstro::Hms2Rad(rah, ram, ras); const Double_t dec = MAstro::Dms2Rad(decd, decm, decs, decsgn); MVector3 *star=new MVector3; star->SetRaDec(ra, dec, mag); star->SetName(row(4,9)); if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV) { delete star; continue; } fList.Add(star); add++; } gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl; return add; } // -------------------------------------------------------------------------- void MAstroCatalog::Paint(Option_t *o) { SetRangePad(); if (TestBit(kHasChanged)) DrawPrimitives(o); } // -------------------------------------------------------------------------- void MAstroCatalog::DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t transparent, const char *txt) { const Double_t ra = v.Phi()*TMath::RadToDeg()/15; const Double_t dec = (TMath::Pi()/2-v.Theta())*TMath::RadToDeg(); TString str = v.GetName(); str += Form(": Ra=%.2fh", ra); str += Form(" Dec=%.1fd", dec); str += Form(" Mag=%.1f", -2.5*log10(v.Mag())); if (txt) str += Form(" (%s)", txt); // draw star on the camera display TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotLarge);; tip->SetMarkerColor(kBlack); AddMap(tip, new TString(str)); } // -------------------------------------------------------------------------- void MAstroCatalog::Update(Bool_t upd) { if (gPad && TestBit(kMustCleanup)) { SetBit(kHasChanged); gPad->Modified(); if (upd) gPad->Update(); } } // -------------------------------------------------------------------------- // // Set the observation time. Necessary to use local coordinate // system. // void MAstroCatalog::SetTime(const MTime &time) { if (fTime) delete fTime; fTime=(MTime*)time.Clone(); } // -------------------------------------------------------------------------- // // Set the observatory location. Necessary to use local coordinate // system. // void MAstroCatalog::SetObservatory(const MObservatory &obs) { if (fObservatory) delete fObservatory; fObservatory=(MObservatory*)obs.Clone(); } // -------------------------------------------------------------------------- // // Convert the vector to pad coordinates. After conversion // the x- coordinate of the vector must be the x coordinate // of the pad - the same for y. If the coordinate is inside // the current draw area return kTRUE, otherwise kFALSE. // If it is an invalid coordinate return kERROR // Int_t MAstroCatalog::ConvertToPad(const TVector3 &w0, TVector2 &v) const { TVector3 w(w0); // Stretch such, that the Z-component is alwas the same. Now // X and Y contains the intersection point between the star-light // and the plain of a virtual plain screen (ccd...) if (TestBit(kPlainScreen)) w *= 1./w(2); w *= TMath::RadToDeg(); v.Set(w(0), w(1)); if (w(2)<0) return kERROR; return w(0)>gPad->GetX1() && w(1)>gPad->GetY1() && w(0)GetX2() && w(1)GetY2(); } // -------------------------------------------------------------------------- Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v) const { MVector3 w; w.SetMagThetaPhi(1, v.Y(), v.X()); w *= rot; return ConvertToPad(w, v); } // -------------------------------------------------------------------------- Bool_t MAstroCatalog::DrawLine(const TVector2 &v, Double_t dx, Double_t dy, const TRotation &rot, Int_t type) { const TVector2 add(dx*TMath::DegToRad(), dy*TMath::DegToRad()); TVector2 v0 = v; TVector2 v1 = v+add; const Int_t rc0 = Convert(rot, v0); const Int_t rc1 = Convert(rot, v1); // Both are kFALSE or both are kERROR if ((rc0|rc1)==kFALSE || (rc0&rc1)==kERROR) return kFALSE; TLine *line = new TLine(v0.X(), v0.Y(), v1.X(), v1.Y()); line->SetLineStyle(kDashDotted); //kDashed, kDotted, kDashDotted line->SetLineColor(kWhite+type*2); AddMap(line); const TVector2 deg = v*TMath::RadToDeg(); TString txt = type==1 ? Form("Ra=%.2fh Dec=%.1fd", fmod(deg.X()/15+48, 24), fmod(90-deg.Y()+270,180)-90) : Form("Zd=%.1fd Az=%.1fd", fmod(deg.Y()+270,180)-90, fmod(deg.X()+720, 360)); TMarker *tip=new TMarker(v0.X(), v0.Y(), kDot); tip->SetMarkerColor(kWhite+type*2); AddMap(tip, new TString(txt)); return kTRUE; } // -------------------------------------------------------------------------- // // Use "local" draw option to align the display to the local // coordinate system instead of the sky coordinate system. // void MAstroCatalog::Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type) { const TVector2 v1 = v0 + TVector2(dx[0]*TMath::DegToRad(), dy[0]*TMath::DegToRad()); Int_t idx[] = {1, 1, 1, 1}; Int_t dirs[4][2] = { {0, stepy}, {stepx, 0}, {0, -stepy}, {-stepx, 0} }; for (int i=0; i180-fRadiusFOV) stepx=36; else { // This is a rough estimate how many degrees are visible const Float_t m = log(fRadiusFOV/180.)/log(90./(fRadiusFOV+1)+1); const Float_t t = log(180.)-m*log(fRadiusFOV); const Float_t f = m*log(90-fabs(90-v.Y()))+t; const Int_t nx = (Int_t)(exp(f)+0.5); stepx = nx<4 ? 1 : nx/4; if (stepx>36) stepx=36; } const Int_t ny = (Int_t)(fRadiusFOV+1); Int_t stepy = ny<4 ? 1 : ny/4; // align stepsizes to be devisor or 180 and 90 while (180%stepx) stepx++; while (90%stepy) stepy++; // align to step-size boundary (search for the nearest one) Int_t dv = 1; while ((int)(v.X())%stepx) { v.Set(v.X()+dv, v.Y()); dv = -TMath::Sign(TMath::Abs(dv)+1, dv); } dv = 1; while ((int)(v.Y())%stepy) { v.Set(v.X(), v.Y()+dv); dv = -TMath::Sign(TMath::Abs(dv)+1, dv); } // draw... v *= TMath::DegToRad(); Draw(v, rot, dx, dy, stepx, stepy, type); } // -------------------------------------------------------------------------- // // Get a rotation matrix which aligns the pointing position // to the center of the x,y plain // TRotation MAstroCatalog::AlignCoordinates(const TVector3 &v) const { TRotation trans; trans.RotateZ(-v.Phi()); trans.RotateY(-v.Theta()); trans.RotateZ(-TMath::Pi()/2); return trans; } // -------------------------------------------------------------------------- // // Return the rotation matrix which converts either sky or // local coordinates to coordinates which pole is the current // pointing direction. // TRotation MAstroCatalog::GetGrid(Bool_t local) { const Bool_t enable = fTime && fObservatory; if (!local) { const TRotation trans(AlignCoordinates(fRaDec)); DrawGrid(fRaDec, trans, 1); if (enable) { const MAstroSky2Local rot(*fTime, *fObservatory); DrawGrid(rot*fRaDec, trans*rot.Inverse(), 2); } return trans; } if (local && enable) { const MAstroSky2Local rot(*fTime, *fObservatory); const TRotation trans(AlignCoordinates(rot*fRaDec)); DrawGrid(fRaDec, trans*rot, 1); DrawGrid(rot*fRaDec, trans, 2); return trans*rot; } return TRotation(); } // -------------------------------------------------------------------------- TString MAstroCatalog::GetPadTitle() const { const Double_t ra = fRaDec.Phi()*TMath::RadToDeg(); const Double_t dec = (TMath::Pi()/2-fRaDec.Theta())*TMath::RadToDeg(); TString txt; txt += Form("\\alpha=%.2fh ", fmod(ra/15+48, 24)); txt += Form("\\delta=%.1f\\circ ", fmod(dec+270,180)-90); txt += Form("/ FOV=%.1f\\circ", fRadiusFOV); if (!fTime || !fObservatory) return txt; const MAstroSky2Local rot(*fTime, *fObservatory); const TVector3 loc = rot*fRaDec; const Double_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta()); const Double_t zd = TMath::RadToDeg()*loc.Theta(); const Double_t az = TMath::RadToDeg()*loc.Phi(); txt.Prepend("#splitline{"); txt += Form(" \\theta=%.1fh ", fmod(zd+270,180)-90); txt += Form("\\phi=%.1f\\circ ", fmod(az+720, 360)); txt += Form(" / \\rho=%.1f\\circ", rho*TMath::RadToDeg()); txt += "}{<"; txt += fTime->GetSqlDateTime(); txt += ">}"; return txt; } // -------------------------------------------------------------------------- void MAstroCatalog::AddPrimitives(Option_t *o) { const Bool_t local = TString(o).Contains("local", TString::kIgnoreCase); cout << "Opt: " << o << endl; const TRotation rot(GetGrid(local)); TIter Next(&fList); TVector3 *v=0; while ((v=(TVector3*)Next())) { // FIXME: Check Magnitude! TVector2 s(v->Phi(), v->Theta()); if (Convert(rot, s)==kTRUE) DrawStar(s.X(), TMath::Pi()/2-s.Y(), *v, kFALSE); } TPaveText *pv = new TPaveText(0.01, 0.90, 0.63, 0.99, "brNDC"); pv->AddText(GetPadTitle()); AddMap(pv); TMarker *mk=new TMarker(0, 0, kMultiply); mk->SetMarkerColor(kBlack); mk->SetMarkerSize(1.5); AddMap(mk); } // -------------------------------------------------------------------------- void MAstroCatalog::SetRangePad() { const Double_t edge = fRadiusFOV/TMath::Sqrt(2.); gPad->Range(-edge, -edge, edge, edge); const Float_t w = gPad->GetWw(); const Float_t h = gPad->GetWh(); if (wRange(-edge, -edge*h/w, edge, edge*h/w); else gPad->Range(-edge*w/h, -edge, edge*w/h, edge); } // -------------------------------------------------------------------------- void MAstroCatalog::DrawPrimitives(Option_t *o) { DeleteMap(); SetRangePad(); TStopwatch clk; clk.Start(); AddPrimitives(o); clk.Stop(); clk.Print(); // Append to a possible second pad if (!gPad->GetListOfPrimitives()->FindObject(this)) AppendPad(o); // Append all objects to pad DrawMap(); ResetBit(kHasChanged); } // -------------------------------------------------------------------------- void MAstroCatalog::DrawMap() { Long_t key, val; TExMapIter map(&fMapG); while (map.Next(key, val)) ((TObject*)key)->Draw(); } // -------------------------------------------------------------------------- void MAstroCatalog::Draw(Option_t *o) { // Append to first pad AppendPad(o); // If contents have not previously changed make sure that // all primitives are recreated. if (!TestBit(kHasChanged)) DrawPrimitives(o); // Connect all TCanvas::ProcessedEvent to this->EventInfo // This means, that after TCanvas has processed an event // EventInfo of this class is called, see TCanvas::HandleInput gPad->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "MAstroCatalog", this, "EventInfo(Int_t,Int_t,Int_t,TObject*)"); // Do this instead of fListG.Draw, because // TCollection overwrites Draw // Would be nice, but doesn't work because the single // graphical object are not handled by TPad anymore... // fListG.AppendPad(); } // -------------------------------------------------------------------------- // // This function was connected to all created canvases. It is used // to redirect GetObjectInfo into our own status bar. // // The 'connection' is done in AddTab // void MAstroCatalog::EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected) { TCanvas *c = (TCanvas*)gTQSender; gPad = c ? c->GetSelectedPad() : NULL; if (!gPad) return; // Try to find a corresponding object with kCannotPick set and // an available TString (for a tool tip) TString *str=0; if (!selected || selected==this) { Long_t key, val; TExMapIter map(&fMapG); while (map.Next(key, val)) { if (!val) continue; TObject *o=(TObject*)key; if (o->DistancetoPrimitive(px, py)>TPad::GetMaxPickDistance()) continue; selected = o; str = (TString*)val; break; } } if (!selected) return; switch (event) { case kMouseMotion: if (!fToolTip->IsMapped() && str) ShowToolTip(px, py, *str); break; case kMouseLeave: if (fToolTip->IsMapped()) fToolTip->Hide(); break; case kKeyPress: ExecuteEvent(kKeyPress, px, py); break; } } // -------------------------------------------------------------------------- void MAstroCatalog::ExecuteEventKbd(Int_t keycode, Int_t keysym) { Double_t dra =0; Double_t ddec=0; switch (keysym) { case kKey_Left: dra = -TMath::DegToRad(); break; case kKey_Right: dra = +TMath::DegToRad(); break; case kKey_Up: ddec = +TMath::DegToRad(); break; case kKey_Down: ddec = -TMath::DegToRad(); break; case kKey_Plus: SetRadiusFOV(fRadiusFOV+1); break; case kKey_Minus: SetRadiusFOV(fRadiusFOV-1); break; default: return; } const Double_t r = fRaDec.Phi(); const Double_t d = TMath::Pi()/2-fRaDec.Theta(); SetRaDec(r+dra, d+ddec); gPad->Update(); } // ------------------------------------------------------------------------ // // Execute a gui event on the camera // void MAstroCatalog::ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2) { if (!TestBit(kGuiActive)) return; if (event==kKeyPress) ExecuteEventKbd(mp1, mp2); } // -------------------------------------------------------------------------- Int_t MAstroCatalog::DistancetoPrimitive(Int_t px, Int_t py) { Int_t min = INT_MAX; Long_t key, val; TExMapIter map(&fMapG); while (map.Next(key, val)) { TObject *o=(TObject*)key; const Int_t d = o->DistancetoPrimitive(px, py); if (dGetWindowID(gPad->GetCanvasID()); const Window_t id2 = fToolTip->GetParent()->GetId(); Window_t id3; gVirtualX->TranslateCoordinates(id1, id2, px, py, x, y, id3); // Show tool tip fToolTip->SetText(txt); fToolTip->Show(x+4, y+4); }