Coding convetion violations fixed.
[u/mrichter/AliRoot.git] / RALICE / AliAstrolab.cxx
CommitLineData
8bde545d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// $Id$
17
18///////////////////////////////////////////////////////////////////////////
19// Class AliAstrolab
20// Virtual lab to correlate measurements with astrophysical phenomena.
21//
caa58e1a 22// This class is derived from TTask, but the only reason for this
23// is to enable this class to serve as a base class for other TTask
24// derived classes (e.g. AliEventSelector) without the need for
25// multiple virtual inheritance.
26// So, this AliAstrolab class itself does not provide any TTask
27// related functionality.
28//
8bde545d 29// The lab can be given a terrestrial location via the usual longitude
30// and latitude specifications.
31// Since this class is derived from AliTimestamp, a lab can also be
32// given a specific timestamp. Together with the terrestrial location
33// this provides access to local (sidereal) times etc...
34// In addition to the usual astronomical reference frames, a local
35// lab reference frame can also be specified. Together with the lab's
36// timestamp this uniquely defines all the coordinate transformations
37// between the various reference frames.
38// These lab characteristics allow a space-time correlation of lab
39// observations with external (astrophysical) phenomena.
40//
41// Observations are entered as generic signals containing a position,
42// reference frame specification and a timestamp.
43// These observations can then be analysed in various reference frames
44// via the available GET functions.
45//
46// Various external (astrophysical) phenomena may be entered as
47// so-called reference signals.
48// This class provides facilities (e.g. MatchRefSignal) to check
49// correlations of the stored measurement with these reference signals.
b867da20 50// Also graphical facilities (e.g. DisplaySignals) are available to
51// provide skymaps in various projections.
8bde545d 52//
53// Coding example :
54// ----------------
55// gSystem->Load("ralice");
56//
57// AliAstrolab lab("IceCube","The South Pole Neutrino Observatory");
58// lab.SetLabPosition(0,-90,"deg"); // South Pole
59//
60// lab.SetLocalFrame(90,180,90,270,0,0); // Local frame has X-axis to the North
61//
62// lab.Data(1,"dms"); // Print laboratory parameters
63//
64// // Enter some observed event to be investigated
65// AliTimestamp ts;
66// ts.SetUT(1989,7,30,8,14,23,738504,0);
67// Float_t vec[3]={1,23.8,118.65};
68// Ali3Vector r;
69// r.SetVector(vec,"sph","deg");
70// lab.SetSignal(&r,"loc","M",&ts,0,"Event10372");
71//
72// // Enter some reference signals
73// Float_t alpha=194818.0;
74// Float_t delta=84400.;
75// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"Altair");
76// alpha=124900.0;
77// delta=272400.;
78// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"NGP");
79// alpha=64508.917;
80// delta=-164258.02;
81// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Sirius");
82// alpha=23149.08;
83// delta=891550.8;
84// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Polaris");
85// alpha=43600.;
86// delta=163100.;
87// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Aldebaran");
88// Float_t l=327.531;
89// Float_t b=-35.8903;
90// Float_t pos[3]={1,90.-b,l};
91// r.SetVector(pos,"sph","deg");
92// lab.SetUT(1989,7,30,8,14,16,0,0);
caa58e1a 93// lab.SetSignal(&r,"gal","M",0,-1,"GRB890730");
8bde545d 94//
95// // List all stored objects
96// lab.ListSignals("equ","M",5);
97//
98// // Determine minimal space and time differences with reference signals
99// Double_t da,dt;
100// Int_t ia,it;
101// da=lab.GetDifference(0,"deg",dt,"s",1,&ia,&it);
102// cout << " Minimal differences damin (deg) : " << da << " dtmin (s) : " << dt
103// << " damin-index : " << ia << " dtmin-index : " << it << endl;
104// cout << " damin for "; lab->PrintSignal("equ","T",&ts,5,ia); cout << endl;
105// cout << " dtmin for "; lab->PrintSignal("equ","T",&ts,5,it); cout << endl;
106//
107// // Search for space and time match with the reference signals
108// da=5;
109// dt=10;
110// TArrayI* arr=lab.MatchRefSignal(da,"deg",dt,"s");
111// Int_t index=0;
112// if (arr)
113// {
114// for (Int_t i=0; i<arr->GetSize(); i++)
115// {
116// index=arr->At(i);
117// cout << " Match found for index : " << index << endl;
118// cout << " Corresponding ref. object "; lab->PrintSignal("equ","T",&ts,5,index); cout << endl;
119// }
120// }
121//
b867da20 122// // Display all stored objects in a skymap (Hammer projection)
123// lab.DisplaySignals("equ","M",&ts,"ham",1);
8bde545d 124//
125//--- Author: Nick van Eijndhoven 15-mar-2007 Utrecht University
126//- Modified: NvE $Date$ Utrecht University
127///////////////////////////////////////////////////////////////////////////
128
b09247a2 129#include <cstdlib>
8bde545d 130#include "AliAstrolab.h"
131#include "Riostream.h"
132
133ClassImp(AliAstrolab) // Class implementation to enable ROOT I/O
134
caa58e1a 135AliAstrolab::AliAstrolab(const char* name,const char* title) : TTask(name,title),AliTimestamp()
8bde545d 136{
137// Default constructor
138
139 fToffset=0;
140 fXsig=0;
141 fRefs=0;
142 fBias=0;
143 fGal=0;
144 fIndices=0;
b867da20 145 fMeridian=-999;
146 fProj="none";
147 fCanvas=0;
148 fHist=0;
149 fMarkers=0;
8bde545d 150}
151///////////////////////////////////////////////////////////////////////////
152AliAstrolab::~AliAstrolab()
153{
154// Destructor to delete all allocated memory.
155
156 if (fXsig)
157 {
158 delete fXsig;
159 fXsig=0;
160 }
161 if (fRefs)
162 {
163 delete fRefs;
164 fRefs=0;
165 }
166 if (fIndices)
167 {
168 delete fIndices;
169 fIndices=0;
170 }
b867da20 171 if (fHist)
172 {
173 delete fHist;
174 fHist=0;
175 }
176 if (fMarkers)
177 {
178 delete fMarkers;
179 fMarkers=0;
180 }
181 if (fCanvas)
182 {
183 delete fCanvas;
184 fCanvas=0;
185 }
8bde545d 186}
187///////////////////////////////////////////////////////////////////////////
caa58e1a 188AliAstrolab::AliAstrolab(const AliAstrolab& t) : TTask(t),AliTimestamp(t)
8bde545d 189{
190// Copy constructor
191
192 fToffset=t.fToffset;
193 fLabPos=t.fLabPos;
194 fXsig=0;
195 if (t.fXsig) fXsig=new AliSignal(*(t.fXsig));
196 fRefs=0;
197 if (t.fRefs)
198 {
199 Int_t size=t.fRefs->GetSize();
200 fRefs=new TObjArray(size);
201 for (Int_t i=0; i<size; i++)
202 {
203 AliSignal* sx=(AliSignal*)t.fRefs->At(i);
204 if (sx) fRefs->AddAt(sx->Clone(),i);
205 }
206 }
207 fBias=0;
208 fGal=0;
209 fIndices=0;
b867da20 210 fMeridian=-999;
211 fProj="none";
212 fCanvas=0;
213 fHist=0;
214 fMarkers=0;
8bde545d 215}
216///////////////////////////////////////////////////////////////////////////
217void AliAstrolab::Data(Int_t mode,TString u)
218{
219// Provide lab information.
220//
221// "mode" indicates the mode of the timestamp info (see AliTimestamp::Date).
222//
223// The string argument "u" allows to choose between different angular units
224// in case e.g. a spherical frame is selected.
225// u = "rad" : angles provided in radians
226// "deg" : angles provided in degrees
227// "dms" : angles provided in ddd:mm:ss.sss
228// "hms" : angles provided in hh:mm:ss.sss
229//
230// The defaults are mode=1 and u="deg".
231
232 const char* name=GetName();
233 const char* title=GetTitle();
234 cout << " *" << ClassName() << "::Data*";
235 if (strlen(name)) cout << " Name : " << GetName();
236 if (strlen(title)) cout << " Title : " << GetTitle();
237 cout << endl;
238
239 Double_t l,b;
240 GetLabPosition(l,b,"deg");
241 cout << " Lab position longitude : "; PrintAngle(l,"deg",u,2);
242 cout << " latitude : "; PrintAngle(b,"deg",u,2);
243 cout << endl;
244 cout << " Lab time offset w.r.t. UT : "; PrintTime(fToffset,12); cout << endl;
245
246 // UT and Local time info
247 Date(mode,fToffset);
248}
249///////////////////////////////////////////////////////////////////////////
250void AliAstrolab::SetLabPosition(Ali3Vector& p)
251{
252// Set the lab position in the terrestrial coordinates.
253// The right handed reference frame is defined such that the North Pole
254// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
255// to an azimuth angle phi=0, with phi increasing eastwards.
256
257 fLabPos.SetPosition(p);
258
259 // Determine local time offset in fractional hours w.r.t. UT.
260 Double_t vec[3];
261 p.GetVector(vec,"sph","deg");
262 Double_t l=vec[2];
263 fToffset=l/15.;
264}
265///////////////////////////////////////////////////////////////////////////
266void AliAstrolab::SetLabPosition(Double_t l,Double_t b,TString u)
267{
268// Set the lab position in the terrestrial longitude (l) and latitude (b).
269// Positions north of the equator have b>0, whereas b<0 indicates
270// locations south of the equator.
271// Positions east of Greenwich have l>0, whereas l<0 indicates
272// locations west of Greenwich.
273//
274// The string argument "u" allows to choose between different angular units
275// u = "rad" : angles provided in radians
276// "deg" : angles provided in degrees
277// "dms" : angles provided in dddmmss.sss
278// "hms" : angles provided in hhmmss.sss
279//
280// The default is u="deg".
281
282 Double_t r=1,theta=0,phi=0;
283
284 l=ConvertAngle(l,u,"deg");
285 b=ConvertAngle(b,u,"deg");
286
287 Double_t offset=90.;
288
289 theta=offset-b;
290 phi=l;
291
292 Double_t p[3]={r,theta,phi};
293 fLabPos.SetPosition(p,"sph","deg");
294
295 // Local time offset in fractional hours w.r.t. UT.
296 fToffset=l/15.;
297}
298///////////////////////////////////////////////////////////////////////////
299AliPosition AliAstrolab::GetLabPosition() const
300{
301// Provide the lab position in the terrestrial coordinates.
302// The right handed reference frame is defined such that the North Pole
303// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
304// to an azimuth angle phi=0, with phi increasing eastwards.
305
306 return fLabPos;
307}
308///////////////////////////////////////////////////////////////////////////
309void AliAstrolab::GetLabPosition(Double_t& l,Double_t& b,TString u) const
310{
311// Provide the lab position in the terrestrial longitude (l) and latitude (b).
312// Positions north of the equator have b>0, whereas b<0 indicates
313// locations south of the equator.
314// Positions east of Greenwich have l>0, whereas l<0 indicates
315// locations west of Greenwich.
316//
317// The string argument "u" allows to choose between different angular units
318// u = "rad" : angles provided in radians
319// "deg" : angles provided in degrees
320//
321// The default is u="deg".
322
323 Double_t pi=acos(-1.);
324
325 Double_t offset=90.;
326 if (u=="rad") offset=pi/2.;
327
328 Double_t p[3];
329 fLabPos.GetPosition(p,"sph",u);
330 b=offset-p[1];
331 l=p[2];
332}
333///////////////////////////////////////////////////////////////////////////
334Double_t AliAstrolab::GetLT()
335{
336// Provide the Lab's local time in fractional hours.
337// A mean solar day lasts 24h (i.e. 86400s).
338//
339// In case a hh:mm:ss format is needed, please use the Convert() facility.
340
341 Double_t h=GetLT(fToffset);
342 return h;
343}
344///////////////////////////////////////////////////////////////////////////
345Double_t AliAstrolab::GetLMST()
346{
347// Provide the Lab's Local Mean Sidereal Time (LMST) in fractional hours.
348// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
349// The definition of GMST is such that a sidereal clock corresponds with
350// 24 sidereal hours per revolution of the Earth.
351// As such, local time offsets w.r.t. UT and GMST can be treated similarly.
352//
353// In case a hh:mm:ss format is needed, please use the Convert() facility.
354
355 Double_t h=GetLMST(fToffset);
356 return h;
357}
358///////////////////////////////////////////////////////////////////////////
359Double_t AliAstrolab::GetLAST()
360{
361// Provide the Lab's Local Apparent Sidereal Time (LAST) in fractional hours.
362// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
363// The definition of GMST and GAST is such that a sidereal clock corresponds with
364// 24 sidereal hours per revolution of the Earth.
365// As such, local time offsets w.r.t. UT, GMST and GAST can be treated similarly.
366//
367// In case a hh:mm:ss format is needed, please use the Convert() facility.
368
369 Double_t h=GetLAST(fToffset);
370 return h;
371}
372///////////////////////////////////////////////////////////////////////////
373void AliAstrolab::PrintAngle(Double_t a,TString in,TString out,Int_t ndig) const
374{
375// Printing of angles in various formats.
376//
377// The input argument "a" denotes the angle to be printed.
378// The string arguments "in" and "out" specify the angular I/O formats.
379//
380// in = "rad" : input angle provided in radians
381// "deg" : input angle provided in degrees
382// "dms" : input angle provided in dddmmss.sss
383// "hms" : input angle provided in hhmmss.sss
384//
385// out = "rad" : output angle provided in radians
386// "deg" : output angle provided in degrees
387// "dms" : output angle provided in dddmmss.sss
388// "hms" : output angle provided in hhmmss.sss
389//
390// The argument "ndig" specifies the number of digits for the fractional
391// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
392// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
393// will appear as 03.4 on the output.
394// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
395//
396// The default is ndig=1.
397//
398// Note : The angle info is printed without additional spaces or "endline".
399// This allows the print to be included in various composite output formats.
400
401 Double_t b=ConvertAngle(a,in,out);
402
403 if (out=="deg" || out=="rad")
404 {
405 cout.setf(ios::fixed,ios::floatfield);
406 cout << setprecision(ndig) << b << " " << out.Data();
407 cout.unsetf(ios::fixed);
408 return;
409 }
410
411 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
412 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
413 Double_t s;
414 ULong64_t sfrac=0;
415
416 if (out=="dms")
417 {
418 word=Int_t(b);
419 word=abs(word);
420 ddd=word/10000;
421 word=word%10000;
422 mm=word/100;
423 ss=word%100;
424 s=fabs(b)-Double_t(ddd*10000+mm*100+ss);
425 if (s>(1.-epsilon))
426 {
427 s=0.;
428 ss++;
429 }
430 while (ss>=60)
431 {
432 ss-=60;
433 mm++;
434 }
435 while (mm>=60)
436 {
437 mm-=60;
438 ddd++;
439 }
440 while (ddd>=360)
441 {
442 ddd-=360;
443 }
444 s*=pow(10.,ndig);
445 sfrac=ULong64_t(s);
446 if (b<0) cout << "-";
447 cout << ddd << "d " << mm << "' " << ss << "."
448 << setfill('0') << setw(ndig) << sfrac << "\"";
449 return;
450 }
451
452 if (out=="hms")
453 {
454 word=Int_t(b);
455 word=abs(word);
456 hh=word/10000;
457 word=word%10000;
458 mm=word/100;
459 ss=word%100;
460 s=fabs(b)-Double_t(hh*10000+mm*100+ss);
461 if (s>(1.-epsilon))
462 {
463 s=0.;
464 ss++;
465 }
466 while (ss>=60)
467 {
468 ss-=60;
469 mm++;
470 }
471 while (mm>=60)
472 {
473 mm-=60;
474 hh++;
475 }
476 while (hh>=24)
477 {
478 hh-=24;
479 }
480 s*=pow(10.,ndig);
481 sfrac=ULong64_t(s);
482 if (b<0) cout << "-";
483 cout << hh << "h " << mm << "m " << ss << "."
484 << setfill('0') << setw(ndig) << sfrac << "s";
485 return;
486 }
487}
488///////////////////////////////////////////////////////////////////////////
489void AliAstrolab::SetSignal(Ali3Vector* r,TString frame,TString mode,AliTimestamp* ts,Int_t jref,TString name)
490{
491// Store a signal as specified by the position r and the timestamp ts.
492// The position is stored in International Celestial Reference System (ICRS) coordinates.
493// The ICRS is a fixed, time independent frame and as such provides a unique reference
494// frame without the need of specifying any epoch etc...
495// The ICRS coordinate definitions match within 20 mas with the mean ones of the J2000.0
496// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
497// coordinate correction between J2000 and ICRS is performed here via the
498// so-called frame bias matrix.
499// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
500// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
501//
502// The input parameter "frame" allows the user to specify the frame to which
503// the components of r refer. Available options are :
504//
505// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
506// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
b867da20 507// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 508// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
b867da20 509// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
8bde545d 510// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
511// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
512// components of r correspond to theta=zenith angle and phi=pi-azimuth.
b867da20 513// "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
8bde545d 514// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
515// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
516// components of r correspond to the usual theta and phi angles.
517//
518// In case the coordinates are the equatorial right ascension and declination (a,d),
519// they can represent so-called "mean" and "true" values.
520// The distinction between these two representations is the following :
521//
522// mean values : (a,d) are only corrected for precession and not for nutation
523// true values : (a,d) are corrected for both precession and nutation
524//
525// The input parameter "mode" allows the user to specifiy either "mean" or "true"
526// values for the input in case of equatorial (a,d) coordinates.
527//
528// mode = "M" --> Input coordinates are the mean values
529// "T" --> Input coordinates are the true values
530//
531// The input parameter "jref" allows the user to store so-called "reference" signals.
532// These reference signals may be used to check space-time event coincidences with the
533// stored measurement (e.g. coincidence of the measurement with transient phenomena).
534//
535// jref = 0 --> Storage of the measurement
536// j --> Storage of a reference signal at the j-th position (j=1 is first)
537// < 0 --> Add a reference signal at the next available position
538//
539// Via the input argument "name" the user can give the stored signal also a name.
540//
541// The default values are jref=0 and name="".
542//
543// Note : In case ts=0 the current timestamp of the lab will be taken.
544
545 if (!r)
546 {
547 if (!jref && fXsig)
548 {
549 delete fXsig;
550 fXsig=0;
551 }
552 return;
553 }
554
555 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc")
556 {
557 if (!jref && fXsig)
558 {
559 delete fXsig;
560 fXsig=0;
561 }
562 return;
563 }
564
565 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t")
566 {
567 if (!jref && fXsig)
568 {
569 delete fXsig;
570 fXsig=0;
571 }
572 return;
573 }
574
575 if (!ts) ts=(AliTimestamp*)this;
576
577 Double_t vec[3]={1,0,0};
578 vec[1]=r->GetX(2,"sph","rad");
579 vec[2]=r->GetX(3,"sph","rad");
580 Ali3Vector q;
581 q.SetVector(vec,"sph","rad");
582
583 AliSignal* sxref=0;
584
585 if (!jref) // Storage of the measurement
586 {
587 if (!fXsig)
588 {
589 fXsig=new AliSignal();
590 }
591 else
592 {
593 fXsig->Reset(1);
594 }
595 if (name != "") fXsig->SetName(name);
596 fXsig->SetTitle("Event in ICRS coordinates");
597 fXsig->SetTimestamp(*ts);
598 }
599 else // Storage of a reference signal
600 {
601 if (!fRefs)
602 {
603 fRefs=new TObjArray();
604 fRefs->SetOwner();
605 }
606 // Expand array size if needed
607 if (jref>0 && jref>=fRefs->GetSize()) fRefs->Expand(jref+1);
608 sxref=new AliSignal();
609 if (name != "") sxref->SetName(name);
610 sxref->SetTitle("Reference event in ICRS coordinates");
611 sxref->SetTimestamp(*ts);
612 }
613
614 if (frame=="loc")
615 {
616 // Convert to horizontal coordinates
617 q=q.GetUnprimed(&fL);
618
619 // Store the signal
620 SetSignal(&q,"hor",mode,ts,jref);
621 return;
622 }
623
624 if (frame=="equ")
625 {
626 // Convert to "mean" values at specified epoch
627 if (mode=="T" || mode=="t")
628 {
629 SetNmatrix(ts);
630 q=q.GetUnprimed(&fN);
631 }
632
633 // Convert to "mean" values at J2000
634 SetPmatrix(ts);
635 q=q.GetUnprimed(&fP);
636
637 // Convert to ICRS values
638 if (!fBias) SetBmatrix();
639 q=q.GetUnprimed(&fB);
640 }
641
642 if (frame=="gal")
643 {
644 // Convert to J2000 equatorial mean coordinates
645 if (fGal != 2) SetGmatrix("J");
646 q=q.GetUnprimed(&fG);
647
648 // Convert to ICRS values
649 if (!fBias) SetBmatrix();
650 q=q.GetUnprimed(&fB);
651 }
652
653 if (frame=="ecl")
654 {
655 // Convert to mean equatorial values at specified epoch
656 SetEmatrix(ts);
657 q=q.GetUnprimed(&fE);
658
659 // Convert to "mean" values at J2000
660 SetPmatrix(ts);
661 q=q.GetUnprimed(&fP);
662
663 // Convert to ICRS values
664 if (!fBias) SetBmatrix();
665 q=q.GetUnprimed(&fB);
666 }
667
668 if (frame=="hor")
669 {
670 // Convert to "true" equatorial values at the specified timestamp
671 SetHmatrix(ts);
672 q=q.GetUnprimed(&fH);
673
674 // Convert to "mean" values at specified timestamp
675 SetNmatrix(ts);
676 q=q.GetUnprimed(&fN);
677
678 // Convert to "mean" values at J2000
679 SetPmatrix(ts);
680 q=q.GetUnprimed(&fP);
681
682 // Convert to ICRS values
683 if (!fBias) SetBmatrix();
684 q=q.GetUnprimed(&fB);
685 }
686
687 // Store the signal in ICRS coordinates
688 if (!jref) // Storage of a regular signal
689 {
690 fXsig->SetPosition(q);
691 }
692 else // Storage of a reference signal
693 {
694 sxref->SetPosition(q);
695 if (jref<0)
696 {
697 fRefs->Add(sxref);
698 }
699 else
700 {
701 fRefs->AddAt(sxref,jref-1);
702 }
703 }
704}
705///////////////////////////////////////////////////////////////////////////
706void AliAstrolab::SetSignal(Double_t a,Double_t d,TString s,Double_t e,TString mode,Int_t jref,TString name)
707{
708// Store a signal with right ascension (a) and declination (d) given for epoch e.
709// The position is stored in International Celestial Reference System (ICRS) coordinates.
710// The ICRS is a fixed, time independent frame and as such provides a unique reference
711// frame without the need of specifying any epoch etc...
712// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
713// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
714// coordinate correction between J2000 and ICRS is performed here via the
715// so-called frame bias matrix.
716// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
717// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
718//
719// The coordinates (a,d) can represent so-called "mean" and "true" values.
720// The distinction between these two representations is the following :
721//
722// mean values : (a,d) are only corrected for precession and not for nutation
723// true values : (a,d) are corrected for both precession and nutation
724//
725// The input parameter "mode" allows the user to specifiy either "mean" or "true"
726// values for the input (a,d) coordinates.
727//
728// a : Right ascension in hhmmss.sss
729// d : Declination in dddmmss.sss
730// s = "B" --> Besselian reference epoch.
731// "J" --> Julian reference epoch.
732// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
733// mode = "M" --> Input coordinates are the mean values
734// "T" --> Input coordinates are the true values
735//
736// The input parameter "jref" allows the user to store so-called "reference" signals.
737// These reference signals may be used to check space-time event coincidences with the
738// stored measurement (e.g. coincidence of the measurement with transient phenomena).
739//
740// jref = 0 --> Storage of the measurement
741// j --> Storage of a reference signal at the j-th position (j=1 is first)
742// < 0 --> Add a reference signal at the next available position
743//
744// Via the input argument "name" the user can give the stored signal also a name.
745//
746// The default values are jref=0 and name="".
747
748 if (s!="B" && s!="b" && s!="J" && s!="j")
749 {
750 if (!jref && fXsig)
751 {
752 delete fXsig;
753 fXsig=0;
754 }
755 return;
756 }
757
758 if (mode!="M" && mode!="m" && mode!="T" && mode!="t")
759 {
760 if (!jref && fXsig)
761 {
762 delete fXsig;
763 fXsig=0;
764 }
765 return;
766 }
767
768 // Convert coordinates to fractional degrees.
769 a=ConvertAngle(a,"hms","deg");
770 d=ConvertAngle(d,"dms","deg");
771
772
773 AliTimestamp tx;
774 tx.SetEpoch(e,s);
775
776 Ali3Vector r;
777 Double_t vec[3]={1.,90.-d,a};
778 r.SetVector(vec,"sph","deg");
779
780 SetSignal(&r,"equ",mode,&tx,jref,name);
781}
782///////////////////////////////////////////////////////////////////////////
783void AliAstrolab::SetSignal(Double_t a,Double_t d,TString mode,AliTimestamp* ts,Int_t jref,TString name)
784{
785// Store a signal with right ascension (a) and declination (d) given for timestamp ts.
786// The position is stored in International Celestial Reference System (ICRS) coordinates.
787// The ICRS is a fixed, time independent frame and as such provides a unique reference
788// frame without the need of specifying any epoch etc...
789// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
790// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
791// coordinate correction between J2000 and ICRS is performed here via the
792// so-called frame bias matrix.
793// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
794// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
795//
796// The coordinates (a,d) can represent so-called "mean" and "true" values.
797// The distinction between these two representations is the following :
798//
799// mean values : (a,d) are only corrected for precession and not for nutation
800// true values : (a,d) are corrected for both precession and nutation
801//
802// The input parameter "mode" allows the user to specifiy either "mean" or "true"
803// values for the input (a,d) coordinates.
804//
805// a : Right ascension in hhmmss.sss
806// d : Declination in dddmmss.sss
807// mode = "M" --> Input coordinates are the mean values
808// "T" --> Input coordinates are the true values
809//
810// The input parameter "jref" allows the user to store so-called "reference" signals.
811// These reference signals may be used to check space-time event coincidences with the
812// stored measurement (e.g. coincidence of the measurement with transient phenomena).
813//
814// jref = 0 --> Storage of the measurement
815// j --> Storage of a reference signal at the j-th position (j=1 is first)
816// < 0 --> Add a reference signal at the next available position
817//
818// Via the input argument "name" the user can give the stored signal also a name.
819//
820// The default values are jref=0 and name="".
821//
822// Note : In case ts=0 the current timestamp of the lab will be taken.
823
824 // Convert coordinates to fractional degrees.
825 a=ConvertAngle(a,"hms","deg");
826 d=ConvertAngle(d,"dms","deg");
827
828 Ali3Vector r;
829 Double_t vec[3]={1.,90.-d,a};
830 r.SetVector(vec,"sph","deg");
831
832 SetSignal(&r,"equ",mode,ts,jref,name);
833}
834///////////////////////////////////////////////////////////////////////////
835AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,Int_t jref)
836{
837// Provide the user specified coordinates of a signal at the specific timestamp ts.
838// The coordinates are returned via the vector argument "r".
839// In addition also a pointer to the stored signal object is provided.
840// In case no stored signal was available or one of the input arguments was
841// invalid, the returned pointer will be 0.
842//
843// The input parameter "frame" allows the user to specify the frame to which
844// the components of r refer. Available options are :
845//
846// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
847// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
b867da20 848// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 849// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
b867da20 850// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
8bde545d 851// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
852// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
853// components of r correspond to theta=zenith angle and phi=pi-azimuth.
b867da20 854// "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
8bde545d 855// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
856// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
857// components of r correspond to the usual theta and phi angles.
858//
859// In case the coordinates are the equatorial right ascension and declination (a,d),
860// they can represent so-called "mean" and "true" values.
861// The distinction between these two representations is the following :
862//
863// mean values : (a,d) are only corrected for precession and not for nutation
864// true values : (a,d) are corrected for both precession and nutation
865//
866// The input parameter "mode" allows the user to specifiy either "mean" or "true"
867// values for the input in case of equatorial (a,d) coordinates.
868//
869// mode = "M" --> Input coordinates are the mean values
870// "T" --> Input coordinates are the true values
871//
872// The input parameter "jref" allows the user to access so-called "reference" signals.
873// These reference signals may be used to check space-time event coincidences with the
874// stored measurement (e.g. coincidence of the measurement with transient phenomena).
875//
876// jref = 0 --> Access to the measurement
877// j --> Access to the reference signal at the j-th position (j=1 is first)
878//
879// Default value is jref=0.
880//
881// Note : In case ts=0 the current timestamp of the lab will be taken.
882
883 r.SetZero();
884
885 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc") return 0;
886
887 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
888
889 AliSignal* sx=GetSignal(jref);
890
891 if (!sx) return 0;
892
893 if (!ts) ts=(AliTimestamp*)this;
894
895 Double_t vec[3];
896 sx->GetPosition(vec,"sph","rad");
897 Ali3Vector q;
898 q.SetVector(vec,"sph","rad");
899
900 if (frame=="icr")
901 {
902 r.Load(q);
903 return sx;
904 }
905
906 // Convert from ICRS to equatorial J2000 coordinates
907 if (!fBias) SetBmatrix();
908 q=q.GetPrimed(&fB);
909
910 if (frame=="equ")
911 {
912 // Precess to specified timestamp
913 AliTimestamp ts1;
914 ts1.SetEpoch(2000,"J");
915 Precess(q,&ts1,ts);
916
917 // Nutation correction if requested
918 if (mode=="T" || mode=="t") Nutate(q,ts);
919 }
920
921 if (frame=="gal")
922 {
923 // Convert from equatorial J2000 to galactic
924 if (fGal != 2) SetGmatrix("J");
925 q=q.GetPrimed(&fG);
926 }
927
928 if (frame=="ecl")
929 {
930 // Precess to specified timestamp
931 AliTimestamp ts1;
932 ts1.SetEpoch(2000,"J");
933 Precess(q,&ts1,ts);
934
935 // Convert from equatorial to ecliptic coordinates
936 SetEmatrix(ts);
937 q=q.GetPrimed(&fE);
938 }
939
940 if (frame=="hor")
941 {
942 // Precess to specified timestamp
943 AliTimestamp ts1;
944 ts1.SetEpoch(2000,"J");
945 Precess(q,&ts1,ts);
946
947 // Nutation correction
948 Nutate(q,ts);
949
950 // Convert from equatorial to horizontal coordinates
951 SetHmatrix(ts);
952 q=q.GetPrimed(&fH);
953 }
954
955 if (frame=="loc")
956 {
957 // Get the signal in horizontal coordinates
958 GetSignal(q,"hor",mode,ts);
959
960 // Convert from horizontal local-frame coordinates
961 q=q.GetPrimed(&fL);
962 }
963
964 r.Load(q);
965 return sx;
966}
967///////////////////////////////////////////////////////////////////////////
968AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,TString name)
969{
970// Provide the user specified coordinates of the signal with the specified
971// name at the specific timestamp ts.
972// The coordinates are returned via the vector argument "r".
973// In addition also a pointer to the stored signal object is provided.
974// In case no such stored signal was available or one of the input arguments was
975// invalid, the returned pointer will be 0.
976//
977// The input parameter "frame" allows the user to specify the frame to which
978// the components of r refer. Available options are :
979//
980// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
981// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
b867da20 982// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 983// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
b867da20 984// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
8bde545d 985// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
986// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
987// components of r correspond to theta=zenith angle and phi=pi-azimuth.
b867da20 988// "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
8bde545d 989// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
990// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
991// components of r correspond to the usual theta and phi angles.
992//
993// In case the coordinates are the equatorial right ascension and declination (a,d),
994// they can represent so-called "mean" and "true" values.
995// The distinction between these two representations is the following :
996//
997// mean values : (a,d) are only corrected for precession and not for nutation
998// true values : (a,d) are corrected for both precession and nutation
999//
1000// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1001// values for the input in case of equatorial (a,d) coordinates.
1002//
1003// mode = "M" --> Input coordinates are the mean values
1004// "T" --> Input coordinates are the true values
1005//
1006// Note : In case ts=0 the current timestamp of the lab will be taken.
1007
1008 AliSignal* sx=0;
1009 Int_t j=GetSignalIndex(name);
1010 if (j>=0) sx=GetSignal(r,frame,mode,ts,j);
1011 return sx;
1012}
1013///////////////////////////////////////////////////////////////////////////
1014AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,Int_t jref)
1015{
1016// Provide precession (and nutation) corrected right ascension (a) and
1017// declination (d) of the stored signal object at the specified timestamp.
1018// In addition also a pointer to the stored signal object is provided.
1019// In case no stored signal was available or one of the input arguments was
1020// invalid, the returned pointer will be 0.
1021//
1022// The coordinates (a,d) can represent so-called "mean" and "true" values.
1023// The distinction between these two representations is the following :
1024//
1025// mean values : (a,d) are only corrected for precession and not for nutation
1026// true values : (a,d) are corrected for both precession and nutation
1027//
1028// The input parameter "mode" allows the user to select either
1029// "mean" or "true" values for (a,d).
1030//
1031// The correction methods used are the new IAU 2000 ones as described in the
1032// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1033// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1034//
1035// a : Right ascension in hhmmss.sss
1036// d : Declination in dddmmss.sss
1037// mode = "M" --> Output coordinates are the mean values
1038// "T" --> Output coordinates are the true values
1039// ts : Timestamp at which the corrected coordinate values are requested.
1040//
1041// The input parameter "jref" allows the user to access so-called "reference" signals.
1042// These reference signals may be used to check space-time event coincidences with the
1043// stored regular signal (e.g. coincidence of the stored signal with transient phenomena).
1044//
1045// jref = 0 --> Access to the measurement
1046// j --> Access to the reference signal at the j-th position (j=1 is first)
1047//
1048// Default value is jref=0.
1049//
1050// Note : In case ts=0 the current timestamp of the lab will be taken.
1051
1052 a=0;
1053 d=0;
1054
1055 Ali3Vector r;
1056 AliSignal* sx=GetSignal(r,"equ",mode,ts,jref);
1057
1058 if (!sx) return 0;
1059
1060 // Retrieve the requested (a,d) values
1061 Double_t vec[3];
1062 r.GetVector(vec,"sph","deg");
1063 d=90.-vec[1];
1064 a=vec[2];
1065
1066 while (a<-360.)
1067 {
1068 a+=360.;
1069 }
1070 while (a>360.)
1071 {
1072 a-=360.;
1073 }
1074 while (d<-90.)
1075 {
1076 d+=90.;
1077 }
1078 while (d>90.)
1079 {
1080 d-=90.;
1081 }
1082
1083 // Convert coordinates to appropriate format
1084 a=ConvertAngle(a,"deg","hms");
1085 d=ConvertAngle(d,"deg","dms");
1086
1087 return sx;
1088}
1089///////////////////////////////////////////////////////////////////////////
1090AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,TString name)
1091{
1092// Provide precession (and nutation) corrected right ascension (a) and
1093// declination (d) of the stored signal object with the specified name
1094// at the specific timestamp ts.
1095// In addition also a pointer to the stored signal object is provided.
1096// In case no stored signal was available or one of the input arguments was
1097// invalid, the returned pointer will be 0.
1098//
1099// The coordinates (a,d) can represent so-called "mean" and "true" values.
1100// The distinction between these two representations is the following :
1101//
1102// mean values : (a,d) are only corrected for precession and not for nutation
1103// true values : (a,d) are corrected for both precession and nutation
1104//
1105// The input parameter "mode" allows the user to select either
1106// "mean" or "true" values for (a,d).
1107//
1108// The correction methods used are the new IAU 2000 ones as described in the
1109// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1110// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1111//
1112// a : Right ascension in hhmmss.sss
1113// d : Declination in dddmmss.sss
1114// mode = "M" --> Output coordinates are the mean values
1115// "T" --> Output coordinates are the true values
1116// ts : Timestamp at which the corrected coordinate values are requested.
1117// name : Name of the requested signal object
1118//
1119// Note : In case ts=0 the current timestamp of the lab will be taken.
1120
1121 AliSignal* sx=0;
1122 Int_t j=GetSignalIndex(name);
1123 if (j>=0) sx=GetSignal(a,d,mode,ts,j);
1124 return sx;
1125}
1126///////////////////////////////////////////////////////////////////////////
1127AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,Int_t jref)
1128{
1129// Provide precession (and nutation) corrected right ascension (a) and
1130// declination (d) of the stored signal object at the specified epoch e.
1131// In addition also a pointer to the stored signal object is provided.
1132// In case no stored signal was available or one of the input arguments was
1133// invalid, the returned pointer will be 0.
1134//
1135// The coordinates (a,d) can represent so-called "mean" and "true" values.
1136// The distinction between these two representations is the following :
1137//
1138// mean values : (a,d) are only corrected for precession and not for nutation
1139// true values : (a,d) are corrected for both precession and nutation
1140//
1141// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1142// values for the input (a,d) coordinates.
1143//
1144// a : Right ascension in hhmmss.sss
1145// d : Declination in dddmmss.sss
1146// s = "B" --> Besselian reference epoch.
1147// "J" --> Julian reference epoch.
1148// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1149// mode = "M" --> Input coordinates are the mean values
1150// "T" --> Input coordinates are the true values
1151//
1152// The input parameter "jref" allows the user to access so-called "reference" signals.
1153// These reference signals may be used to check space-time event coincidences with the
1154// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1155//
1156// jref = 0 --> Access to the measurement
1157// j --> Access to the reference signal at the j-th position (j=1 is first)
1158//
1159// Default value is jref=0.
1160
1161 a=0;
1162 d=0;
1163
1164 if (s!="B" && s!="b" && s!="J" && s!="j") return 0;
1165
1166 if (mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
1167
1168 // Convert coordinates to fractional degrees.
1169 a=ConvertAngle(a,"hms","deg");
1170 d=ConvertAngle(d,"dms","deg");
1171
1172
1173 AliTimestamp tx;
1174 tx.SetEpoch(e,s);
1175
1176 AliSignal* sx=GetSignal(a,d,mode,&tx,jref);
1177 return sx;
1178}
1179///////////////////////////////////////////////////////////////////////////
1180AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,TString name)
1181{
1182// Provide precession (and nutation) corrected right ascension (a) and
1183// declination (d) of the stored signal object with the specified name
1184// at the specific epoch e.
1185// In addition also a pointer to the stored signal object is provided.
1186// In case no stored signal was available or one of the input arguments was
1187// invalid, the returned pointer will be 0.
1188//
1189// The coordinates (a,d) can represent so-called "mean" and "true" values.
1190// The distinction between these two representations is the following :
1191//
1192// mean values : (a,d) are only corrected for precession and not for nutation
1193// true values : (a,d) are corrected for both precession and nutation
1194//
1195// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1196// values for the input (a,d) coordinates.
1197//
1198// a : Right ascension in hhmmss.sss
1199// d : Declination in dddmmss.sss
1200// s = "B" --> Besselian reference epoch.
1201// "J" --> Julian reference epoch.
1202// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1203// mode = "M" --> Input coordinates are the mean values
1204// "T" --> Input coordinates are the true values
1205// name : Name of the requested signal object
1206
1207 AliSignal* sx=0;
1208 Int_t j=GetSignalIndex(name);
1209 if (j>=0) sx=GetSignal(a,d,s,e,mode,j);
1210 return sx;
1211}
1212///////////////////////////////////////////////////////////////////////////
1213AliSignal* AliAstrolab::GetSignal(Int_t jref)
1214{
1215// Provide the pointer to a stored signal object.
1216//
1217// The input parameter "jref" allows the user to access so-called "reference" signals.
1218// These reference signals may be used to check space-time event coincidences with the
1219// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1220//
1221// jref = 0 --> Access to the measurement
1222// j --> Access to the reference signal at the j-th position (j=1 is first)
1223//
1224// Default value is jref=0.
1225
1226 AliSignal* sx=0;
1227 if (!jref)
1228 {
1229 sx=fXsig;
1230 }
1231 else
1232 {
1233 if (jref>0 && jref<fRefs->GetSize()) sx=(AliSignal*)fRefs->At(jref-1);
1234 }
1235 return sx;
1236}
1237///////////////////////////////////////////////////////////////////////////
1238AliSignal* AliAstrolab::GetSignal(TString name)
1239{
1240// Provide the pointer to the stored signal object with the specified name.
1241
1242 AliSignal* sx=0;
1243 Int_t j=GetSignalIndex(name);
1244 if (j>=0) sx=GetSignal(j);
1245 return sx;
1246}
1247///////////////////////////////////////////////////////////////////////////
1248void AliAstrolab::RemoveRefSignal(Int_t j,Int_t compress)
1249{
1250// Remove the reference signal which was stored at the j-th position (j=1 is first).
1251// Note : j=0 means that all stored reference signals will be removed.
1252// j<0 allows array compression (see below) without removing any signals.
1253//
1254// The "compress" parameter allows compression of the ref. signal storage array.
1255//
1256// compress = 1 --> Array will be compressed
1257// 0 --> Array will not be compressed
1258//
1259// Note : Compression of the storage array means that the indices of the
1260// reference signals in the storage array will change.
1261
1262 if (!fRefs) return;
1263
1264 // Clearing of the complete storage
1265 if (!j)
1266 {
1267 delete fRefs;
1268 fRefs=0;
1269 return;
1270 }
1271
1272 // Removing a specific reference signal
1273 if (j>0 && j<fRefs->GetSize())
1274 {
1275 TObject* obj=fRefs->RemoveAt(j-1);
1276 if (obj) delete obj;
1277 }
1278
1279 // Compression of the storage array
1280 if (compress) fRefs->Compress();
1281}
1282///////////////////////////////////////////////////////////////////////////
1283void AliAstrolab::RemoveRefSignal(TString name,Int_t compress)
1284{
1285// Remove the reference signal with the specified name.
1286//
1287// The "compress" parameter allows compression of the ref. signal storage array.
1288//
1289// compress = 1 --> Array will be compressed
1290// 0 --> Array will not be compressed
1291//
1292// Note : Compression of the storage array means that the indices of the
1293// reference signals in the storage array will change.
1294
1295 Int_t j=GetSignalIndex(name);
1296 if (j>0) RemoveRefSignal(j,compress);
1297}
1298///////////////////////////////////////////////////////////////////////////
1299Int_t AliAstrolab::GetSignalIndex(TString name)
1300{
1301// Provide storage index of the signal with the specified name.
1302// In case the name matches with the stored measurement,
1303// the value 0 is returned.
1304// In case no signal with the specified name was found, the value -1 is returned.
1305
1306 Int_t index=-1;
1307
1308 if (fXsig)
1309 {
1310 if (name==fXsig->GetName()) return 0;
1311 }
1312
1313 if (!fRefs) return -1;
1314
1315 for (Int_t i=0; i<fRefs->GetSize(); i++)
1316 {
1317 AliSignal* sx=(AliSignal*)fRefs->At(i);
1318 if (!sx) continue;
1319
1320 if (name==sx->GetName())
1321 {
1322 index=i+1;
1323 break;
1324 }
1325 }
1326
1327 return index;
1328}
1329///////////////////////////////////////////////////////////////////////////
1330void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,Int_t jref)
1331{
1332// Print data of a stored signal in user specified coordinates at the specific timestamp ts.
1333// In case no stored signal was available or one of the input arguments was
1334// invalid, no printout is produced.
1335//
1336// The argument "ndig" specifies the number of digits for the fractional
1337// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1338// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1339// will appear as 03.4 on the output.
1340// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1341//
1342// Note : The angle info is printed without additional spaces or "endline".
1343// This allows the print to be included in various composite output formats.
1344//
1345// The input parameter "frame" allows the user to specify the frame to which
1346// the coordinates refer. Available options are :
1347//
1348// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1349//
b867da20 1350// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 1351//
b867da20 1352// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
8bde545d 1353//
1354// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1355//
b867da20 1356// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
8bde545d 1357//
1358// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1359//
1360// In case the coordinates are the equatorial right ascension and declination (a,d),
1361// they can represent so-called "mean" and "true" values.
1362// The distinction between these two representations is the following :
1363//
1364// mean values : (a,d) are only corrected for precession and not for nutation
1365// true values : (a,d) are corrected for both precession and nutation
1366//
1367// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1368// values for the input in case of equatorial (a,d) coordinates.
1369//
1370// mode = "M" --> Input coordinates are the mean values
1371// "T" --> Input coordinates are the true values
1372//
661385a7 1373// The input parameter "mode" also determines which type of time and
1374// local hour angle will appear in the printout.
1375//
caa58e1a 1376// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1377// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
661385a7 1378//
8bde545d 1379// The input parameter "jref" allows printing of a so-called "reference" signal.
1380// These reference signals may serve to check space-time event coincidences with the
1381// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1382//
1383// jref = 0 --> Printing of the measurement
1384// j --> Printing of the j-th reference signal
1385//
1386// Default value is jref=0.
1387//
1388// Note : In case ts=0 the current timestamp of the lab will be taken.
1389
1390 Ali3Vector r;
1391 AliSignal* sx=GetSignal(r,frame,mode,ts,jref);
1392
1393 if (!sx) return;
1394
661385a7 1395 // Local Hour Angle of the signal
1396 Double_t lha=GetHourAngle("A",ts,jref);
1397 TString slha="LAHA";
1398 if (mode=="M" || mode=="m")
1399 {
1400 lha=GetHourAngle("M",ts,jref);
1401 slha="LMHA";
1402 }
1403
8bde545d 1404 TString name=sx->GetName();
1405 if (name != "") cout << name.Data() << " ";
1406
1407 if (frame=="equ")
1408 {
1409 Double_t a,d;
1410 d=90.-r.GetX(2,"sph","deg");
1411 a=r.GetX(3,"sph","rad");
1412 cout << "Equatorial (" << mode.Data() <<") a : "; PrintAngle(a,"rad","hms",ndig);
1413 cout << " d : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1414 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1415 return;
1416 }
1417
1418 if (frame=="gal")
1419 {
1420 Double_t l,b;
1421 b=90.-r.GetX(2,"sph","deg");
1422 l=r.GetX(3,"sph","deg");
1423 cout << "Galactic l : "; PrintAngle(l,"deg","deg",ndig);
1424 cout << " b : "; PrintAngle(b,"deg","deg",ndig);
661385a7 1425 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1426 return;
1427 }
1428
1429 if (frame=="icr")
1430 {
1431 Double_t a,d;
1432 d=90.-r.GetX(2,"sph","deg");
1433 a=r.GetX(3,"sph","rad");
1434 cout << "ICRS l : "; PrintAngle(a,"rad","hms",ndig);
1435 cout << " b : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1436 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1437 return;
1438 }
1439
1440 if (frame=="ecl")
1441 {
1442 Double_t a,d;
1443 d=90.-r.GetX(2,"sph","deg");
1444 a=r.GetX(3,"sph","deg");
1445 cout << "Ecliptic l : "; PrintAngle(a,"deg","deg",ndig);
1446 cout << " b : "; PrintAngle(d,"deg","deg",ndig);
661385a7 1447 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1448 return;
1449 }
1450
1451 if (frame=="hor")
1452 {
1453 Double_t alt=90.-r.GetX(2,"sph","deg");
1454 Double_t azi=180.-r.GetX(3,"sph","deg");
1455 while (azi>360)
1456 {
1457 azi-=360.;
1458 }
1459 while (azi<0)
1460 {
1461 azi+=360.;
1462 }
1463 cout << "Horizontal azi : "; PrintAngle(azi,"deg","deg",ndig);
1464 cout << " alt : "; PrintAngle(alt,"deg","deg",ndig);
661385a7 1465 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1466 return;
1467 }
1468
1469 if (frame=="loc")
1470 {
1471 Double_t theta=r.GetX(2,"sph","deg");
1472 Double_t phi=r.GetX(3,"sph","deg");
1473 cout << "Local-frame phi : "; PrintAngle(phi,"deg","deg",ndig);
1474 cout << " theta : "; PrintAngle(theta,"deg","deg",ndig);
661385a7 1475 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1476 return;
1477 }
1478}
1479///////////////////////////////////////////////////////////////////////////
1480void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,TString name)
1481{
1482// Print data of the stored signal with the specified name in user specified coordinates
1483// at the specific timestamp ts.
b867da20 1484// In case no such stored signal was available or one of the input arguments was
8bde545d 1485// invalid, no printout is produced.
1486//
1487// The argument "ndig" specifies the number of digits for the fractional
1488// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1489// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1490// will appear as 03.4 on the output.
1491// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1492//
1493// Note : The angle info is printed without additional spaces or "endline".
1494// This allows the print to be included in various composite output formats.
1495//
1496// The input parameter "frame" allows the user to specify the frame to which
1497// the coordinates refer. Available options are :
1498//
1499// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1500//
b867da20 1501// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 1502//
b867da20 1503// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
8bde545d 1504//
1505// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1506//
b867da20 1507// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
8bde545d 1508//
1509// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1510//
1511// In case the coordinates are the equatorial right ascension and declination (a,d),
1512// they can represent so-called "mean" and "true" values.
1513// The distinction between these two representations is the following :
1514//
1515// mean values : (a,d) are only corrected for precession and not for nutation
1516// true values : (a,d) are corrected for both precession and nutation
1517//
1518// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1519// values for the input in case of equatorial (a,d) coordinates.
1520//
1521// mode = "M" --> Input coordinates are the mean values
1522// "T" --> Input coordinates are the true values
1523//
661385a7 1524// The input parameter "mode" also determines which type of time and
1525// local hour angle will appear in the printout.
1526//
caa58e1a 1527// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1528// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
661385a7 1529//
8bde545d 1530// Note : In case ts=0 the current timestamp of the lab will be taken.
1531
1532 Int_t j=GetSignalIndex(name);
1533 if (j>=0) PrintSignal(frame,mode,ts,ndig,j);
1534}
1535///////////////////////////////////////////////////////////////////////////
1536void AliAstrolab::ListSignals(TString frame,TString mode,Int_t ndig)
1537{
1538// List all stored signals in user specified coordinates at the timestamp
1539// of the actual recording of the stored measurement under investigation.
1540// In case no (timestamp of the) actual measurement is available,
1541// the current timestamp of the lab will be taken.
1542// In case no stored signal is available or one of the input arguments is
1543// invalid, no printout is produced.
1544//
1545// The argument "ndig" specifies the number of digits for the fractional
1546// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1547// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1548// will appear as 03.4 on the output.
1549// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1550//
661385a7 1551// The default value is ndig=1.
8bde545d 1552//
1553// Note : The angle info is printed without additional spaces or "endline".
1554// This allows the print to be included in various composite output formats.
1555//
1556// The input parameter "frame" allows the user to specify the frame to which
1557// the coordinates refer. Available options are :
1558//
1559// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1560//
b867da20 1561// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
8bde545d 1562//
b867da20 1563// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
8bde545d 1564//
1565// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1566//
b867da20 1567// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
8bde545d 1568//
1569// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1570//
1571// In case the coordinates are the equatorial right ascension and declination (a,d),
1572// they can represent so-called "mean" and "true" values.
1573// The distinction between these two representations is the following :
1574//
1575// mean values : (a,d) are only corrected for precession and not for nutation
1576// true values : (a,d) are corrected for both precession and nutation
1577//
1578// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1579// values for the input in case of equatorial (a,d) coordinates.
1580//
1581// mode = "M" --> Input coordinates are the mean values
1582// "T" --> Input coordinates are the true values
1583//
661385a7 1584// The input parameter "mode" also determines which type of time and
1585// local hour angle will appear in the listing.
1586//
caa58e1a 1587// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1588// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
8bde545d 1589
1590 Int_t iprint=0;
1591
1592 AliTimestamp* tx=0;
1593
661385a7 1594 Int_t dform=1;
1595 if (mode=="T" || mode=="t") dform=-1;
1596
8bde545d 1597 if (fXsig)
1598 {
661385a7 1599 tx=fXsig->GetTimestamp();
8bde545d 1600 if (!tx) tx=(AliTimestamp*)this;
1601 cout << " *AliAstrolab::ListSignals* List of all stored signals." << endl;
661385a7 1602 cout << " === The measurement under investigation ===" << endl;
1603 cout << " Timestamp of the actual observation" << endl;
1604 tx->Date(dform,fToffset);
1605 cout << " Location of the actual observation" << endl;
1606 cout << " "; PrintSignal(frame,mode,tx,ndig); cout << endl;
8bde545d 1607 iprint=1;
1608 }
1609
1610 if (!fRefs) return;
1611
1612 for (Int_t i=1; i<=fRefs->GetSize(); i++)
1613 {
1614 AliSignal* sx=GetSignal(i);
1615 if (!sx) continue;
1616
1617 if (!iprint)
1618 {
1619 cout << " *AliAstrolab::ListRefSignals* List of all stored signals." << endl;
1620 tx=(AliTimestamp*)this;
661385a7 1621 cout << " Current timestamp of the laboratory" << endl;
1622 tx->Date(dform,fToffset);
8bde545d 1623 iprint=1;
1624 }
1625 if (iprint==1)
1626 {
661385a7 1627 cout << " === All stored reference signals according to the above timestamp ===" << endl;
8bde545d 1628 iprint=2;
1629 }
661385a7 1630 cout << " Index : " << i << " "; PrintSignal(frame,mode,tx,ndig,i); cout << endl;
8bde545d 1631 }
1632}
1633///////////////////////////////////////////////////////////////////////////
1634void AliAstrolab::Precess(Ali3Vector& r,AliTimestamp* ts1,AliTimestamp* ts2)
1635{
1636// Correct mean right ascension and declination given for Julian date "jd"
1637// for the earth's precession corresponding to the specified timestamp.
1638// The results are the so-called "mean" (i.e. precession corrected) values,
1639// corresponding to the specified timestamp.
1640// The method used is the new IAU 2000 one as described in the
1641// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1642// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1643// Since the standard reference epoch is J2000, this implies that all
1644// input (a,d) coordinates will be first internally converted to the
1645// corresponding J2000 values before the precession correction w.r.t. the
1646// specified lab timestamp will be applied.
1647//
1648// r : Input vector containing the right ascension and declination information
1649// in the form of standard Ali3Vector coordinates.
1650// In spherical coordinates the phi corresponds to the right ascension,
1651// whereas the declination corresponds to (pi/2)-theta.
1652// jd : Julian date corresponding to the input coordinate values.
1653// ts : Timestamp corresponding to the requested corrected coordinate values.
1654//
1655// Note : In case ts=0 the current timestamp of the lab will be taken.
1656
1657 // Convert back to J2000 values
1658 Ali3Vector r0;
1659 SetPmatrix(ts1);
1660 r0=r.GetUnprimed(&fP);
1661
1662 // Precess to the specified timestamp
1663 if (!ts2) ts2=(AliTimestamp*)this;
1664 SetPmatrix(ts2);
1665 r=r0.GetPrimed(&fP);
1666}
1667///////////////////////////////////////////////////////////////////////////
1668void AliAstrolab::Nutate(Ali3Vector& r,AliTimestamp* ts)
1669{
1670// Correct mean right ascension and declination for the earth's nutation
1671// corresponding to the specified timestamp.
1672// The results are the so-called "true" (i.e. nutation corrected) values,
1673// corresponding to the specified timestamp.
1674// The method used is the new IAU 2000 one as described in the
1675// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1676// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1677//
1678// r : Input vector containing the right ascension and declination information
1679// in the form of standard Ali3Vector coordinates.
1680// In spherical coordinates the phi corresponds to the right ascension,
1681// whereas the declination corresponds to (pi/2)-theta.
1682// ts : Timestamp for which the corrected coordinate values are requested.
1683//
1684// Note : In case ts=0 the current timestamp of the lab will be taken.
1685
1686 // Nutation correction for the specified timestamp
1687 if (!ts) ts=(AliTimestamp*)this;
1688 SetNmatrix(ts);
1689 r=r.GetPrimed(&fN);
1690}
1691///////////////////////////////////////////////////////////////////////////
1692void AliAstrolab::SetBmatrix()
1693{
1694// Set the frame bias matrix elements.
1695// The formulas and numerical constants used are the ones from the
1696// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1697// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1698
1699 Double_t pi=acos(-1.);
1700
1701 // Parameters in mas
1702 Double_t a=-14.6;
1703 Double_t x=-16.6170;
1704 Double_t e=-6.8192;
1705
1706 // Convert to radians
1707 a*=pi/(180.*3600.*1000.);
1708 x*=pi/(180.*3600.*1000.);
1709 e*=pi/(180.*3600.*1000.);
1710
1711 Double_t mat[9];
1712 mat[0]=1.-0.5*(a*a+x*x);
1713 mat[1]=a;
1714 mat[2]=-x;
1715 mat[3]=-a-e*x;
1716 mat[4]=1.-0.5*(a*a+e*e);
1717 mat[5]=-e;
1718 mat[6]=x-e*a;
1719 mat[7]=e+x*a;
1720 mat[8]=1.-0.5*(e*e+x*x);
1721
1722 fB.SetMatrix(mat);
1723 fBias=1;
1724}
1725///////////////////////////////////////////////////////////////////////////
1726void AliAstrolab::SetPmatrix(AliTimestamp* ts)
1727{
1728// Set precession matrix elements for Julian date jd w.r.t. J2000.
1729// The formulas and numerical constants used are the ones from the
1730// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1731// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1732// All numerical constants refer to the standard reference epoch J2000.
1733
1734 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1735 if (!ts)
1736 {
1737 fP.SetMatrix(mat);
1738 return;
1739 }
1740
1741 Double_t pi=acos(-1.);
1742
1743 Double_t t=(ts->GetJD()-2451545.0)/36525.; // Julian centuries since J2000.0
1744
1745 // Parameters for the precession matrix in arcseconds
1746 Double_t eps0=84381.406; // Mean ecliptic obliquity at J2000.0
1747 Double_t psi=5038.481507*t-1.0790069*pow(t,2)-0.00114045*pow(t,3)+0.000132851*pow(t,4)
1748 -0.0000000951*pow(t,4);
1749 Double_t om=eps0-0.025754*t+0.0512623*pow(t,2)-0.00772503*pow(t,3)-0.000000467*pow(t,4)
1750 +0.0000003337*pow(t,5);
1751 Double_t chi=10.556403*t-2.3814292*pow(t,2)-0.00121197*pow(t,3)+0.000170663*pow(t,4)
1752 -0.0000000560*pow(t,5);
1753
1754 // Convert to radians
1755 eps0*=pi/(180.*3600.);
1756 psi*=pi/(180.*3600.);
1757 om*=pi/(180.*3600.);
1758 chi*=pi/(180.*3600.);
1759
1760 Double_t s1=sin(eps0);
1761 Double_t s2=sin(-psi);
1762 Double_t s3=sin(-om);
1763 Double_t s4=sin(chi);
1764 Double_t c1=cos(eps0);
1765 Double_t c2=cos(-psi);
1766 Double_t c3=cos(-om);
1767 Double_t c4=cos(chi);
1768
1769 mat[0]=c4*c2-s2*s4*c3;
1770 mat[1]=c4*s2*c1+s4*c3*c2*c1-s1*s4*s3;
1771 mat[2]=c4*s2*s1+s4*c3*c2*s1+c1*s4*s3;
1772 mat[3]=-s4*c2-s2*c4*c3;
1773 mat[4]=-s4*s2*c1+c4*c3*c2*c1-s1*c4*s3;
1774 mat[5]=-s4*s2*s1+c4*c3*c2*s1+c1*c4*s3;
1775 mat[6]=s2*s3;
1776 mat[7]=-s3*c2*c1-s1*c3;
1777 mat[8]=-s3*c2*s1+c3*c1;
1778
1779 fP.SetMatrix(mat);
1780}
1781///////////////////////////////////////////////////////////////////////////
1782void AliAstrolab::SetNmatrix(AliTimestamp* ts)
1783{
1784// Set nutation matrix elements for the specified Julian date jd.
1785// The formulas and numerical constants used are the ones from the
1786// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1787// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1788
1789 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1790 if (!ts)
1791 {
1792 fN.SetMatrix(mat);
1793 return;
1794 }
1795
1796 Double_t pi=acos(-1.);
1797
1798 Double_t dpsi,deps,eps;
1799 ts->Almanac(&dpsi,&deps,&eps);
1800
1801 // Convert to radians
1802 dpsi*=pi/(180.*3600.);
1803 deps*=pi/(180.*3600.);
1804 eps*=pi/(180.*3600.);
1805
1806 Double_t s1=sin(eps);
1807 Double_t s2=sin(-dpsi);
1808 Double_t s3=sin(-(eps+deps));
1809 Double_t c1=cos(eps);
1810 Double_t c2=cos(-dpsi);
1811 Double_t c3=cos(-(eps+deps));
1812
1813 mat[0]=c2;
1814 mat[1]=s2*c1;
1815 mat[2]=s2*s1;
1816 mat[3]=-s2*c3;
1817 mat[4]=c3*c2*c1-s1*s3;
1818 mat[5]=c3*c2*s1+c1*s3;
1819 mat[6]=s2*s3;
1820 mat[7]=-s3*c2*c1-s1*c3;
1821 mat[8]=-s3*c2*s1+c3*c1;
1822
1823 fN.SetMatrix(mat);
1824}
1825///////////////////////////////////////////////////////////////////////////
1826void AliAstrolab::SetGmatrix(TString mode)
1827{
1828// Set the mean equatorial to galactic coordinate conversion matrix.
1829// The B1950 parameters were taken from section 22.3 of the book
1830// "An Introduction to Modern Astrophysics" by Carrol and Ostlie (1996).
1831// The J2000 parameters are obtained by precession of the B1950 values.
1832//
1833// Via the input argument "mode" the required epoch can be selected
1834// mode = "B" ==> B1950
1835// "J" ==> J2000
1836
1837 Ali3Vector x; // The Galactic x-axis in the equatorial frame
1838 Ali3Vector y; // The Galactic y-axis in the equatorial frame
1839 Ali3Vector z; // The Galactic z-axis in the equatorial frame
1840
1841 Double_t a,d;
1842 Double_t vec[3]={1,0,0};
1843
1844 fGal=1; // Set flag to indicate B1950 matrix values
1845
1846 // B1950 equatorial coordinates of the North Galactic Pole (NGP)
1847 a=124900.;
1848 d=272400.;
1849 a=ConvertAngle(a,"hms","deg");
1850 d=ConvertAngle(d,"dms","deg");
1851 vec[1]=90.-d;
1852 vec[2]=a;
1853 z.SetVector(vec,"sph","deg");
1854
1855 // B1950 equatorial coordinates of the Galactic l=b=0 point
1856 a=174224.;
1857 d=-285500.;
1858 a=ConvertAngle(a,"hms","deg");
1859 d=ConvertAngle(d,"dms","deg");
1860 vec[1]=90.-d;
1861 vec[2]=a;
1862 x.SetVector(vec,"sph","deg");
1863
1864 // Precess to the corresponding J2000 values if requested
1865 if (mode=="J")
1866 {
1867 fGal=2; // Set flag to indicate J2000 matrix values
1868 AliTimestamp t1;
1869 t1.SetEpoch(1950,"B");
1870 AliTimestamp t2;
1871 t2.SetEpoch(2000,"J");
1872 Precess(z,&t1,&t2);
1873 Precess(x,&t1,&t2);
1874 }
1875
1876 // The Galactic y-axis is determined for the right handed frame
1877 y=z.Cross(x);
1878
1879 fG.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1880 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1881 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1882}
1883///////////////////////////////////////////////////////////////////////////
1884void AliAstrolab::SetEmatrix(AliTimestamp* ts)
1885{
1886// Set the mean equatorial to ecliptic coordinate conversion matrix
1887// for the specified timestamp.
1888// A nice sketch and explanation of the two frames can be found
1889// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1890
1891 Double_t dpsi,deps,eps;
1892 ts->Almanac(&dpsi,&deps,&eps);
1893
1894 // Convert to degrees
1895 eps/=3600.;
1896
1897 // Positions of the ecliptic axes w.r.t. the equatorial ones
1898 // at the moment of the specified timestamp
1899 Double_t theta1=90; // Ecliptic x-axis
1900 Double_t phi1=0;
1901 Double_t theta2=90.-eps; //Ecliptic y-axis
1902 Double_t phi2=90;
1903 Double_t theta3=eps; // Ecliptic z-axis
1904 Double_t phi3=270;
1905
1906 fE.SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);
1907}
1908///////////////////////////////////////////////////////////////////////////
1909void AliAstrolab::SetHmatrix(AliTimestamp* ts)
1910{
1911// Set the mean equatorial to horizontal coordinate conversion matrix
1912// for the specified timestamp.
1913// A nice sketch and explanation of the two frames can be found
1914// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1915//
1916// Note : In order to simplify the calculations, we use here a
1917// right-handed horizontal frame.
1918
1919 Ali3Vector x; // The (South pointing) horizontal x-axis in the equatorial frame
1920 Ali3Vector y; // The (East pointing) horizontal y-axis in the equatorial frame
1921 Ali3Vector z; // The (Zenith pointing) horizontal z-axis in the equatorial frame
1922
1923 Double_t l,b;
1924 GetLabPosition(l,b,"deg");
1925
1926 Double_t a;
1927 Double_t vec[3]={1,0,0};
1928
1929 // Equatorial coordinates of the horizontal z-axis
1930 // at the moment of the specified timestamp
1931 a=ts->GetLAST(fToffset);
1932 a*=15.; // Convert fractional hours to degrees
1933 vec[1]=90.-b;
1934 vec[2]=a;
1935 z.SetVector(vec,"sph","deg");
1936
1937 // Equatorial coordinates of the horizontal x-axis
1938 // at the moment of the specified timestamp
1939 vec[1]=180.-b;
1940 vec[2]=a;
1941 x.SetVector(vec,"sph","deg");
1942
1943 // The horizontal y-axis is determined for the right handed frame
1944 y=z.Cross(x);
1945
1946 fH.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1947 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1948 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1949}
1950///////////////////////////////////////////////////////////////////////////
1951void AliAstrolab::SetLocalFrame(Double_t t1,Double_t p1,Double_t t2,Double_t p2,Double_t t3,Double_t p3)
1952{
1953// Specification of the orientations of the local-frame axes.
1954// The input arguments represent the angles (in degrees) of the local-frame axes
1955// w.r.t. a so called Master Reference Frame (MRF), with the same convention
f535bc3b 1956// as the input arguments of TRotMatrix::SetAngles.
8bde545d 1957//
1958// The right handed Master Reference Frame is defined as follows :
1959// Z-axis : Points to the local Zenith
1960// X-axis : Makes an angle of 90 degrees with the Z-axis and points South
1961// Y-axis : Makes an angle of 90 degrees with the Z-axis and points East
1962//
1963// Once the user has specified the local reference frame, any observed event
1964// can be related to astronomical space-time locations via the SetSignal
1965// and GetSignal memberfunctions.
1966
1967 // Set the matrix for the conversion of our reference frame coordinates
1968 // into the local-frame ones.
1969
1970 fL.SetAngles(t1,p1,t2,p2,t3,p3);
1971}
1972///////////////////////////////////////////////////////////////////////////
1973Double_t AliAstrolab::ConvertAngle(Double_t a,TString in,TString out) const
1974{
1975// Conversion of various angular formats.
1976//
1977// The input argument "a" denotes the angle to be converted.
1978// The string arguments "in" and "out" specify the angular I/O formats.
1979//
1980// in = "rad" : input angle provided in radians
1981// "deg" : input angle provided in degrees
1982// "dms" : input angle provided in dddmmss.sss
1983// "hms" : input angle provided in hhmmss.sss
b867da20 1984// "hrs" : input angle provided in fractional hours
8bde545d 1985//
1986// out = "rad" : output angle provided in radians
1987// "deg" : output angle provided in degrees
1988// "dms" : output angle provided in dddmmss.sss
1989// "hms" : output angle provided in hhmmss.sss
b867da20 1990// "hrs" : output angle provided in fractional hours
8bde545d 1991
1992 if (in==out) return a;
1993
1994 // Convert input to its absolute value in (fractional) degrees.
1995 Double_t pi=acos(-1.);
1996 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
1997 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
1998 Double_t s=0;
1999
2000 Double_t b=fabs(a);
2001
2002 if (in=="rad") b*=180./pi;
2003
b867da20 2004 if (in=="hrs") b*=15.;
2005
8bde545d 2006 if (in=="dms")
2007 {
2008 word=Int_t(b);
2009 ddd=word/10000;
2010 word=word%10000;
2011 mm=word/100;
2012 ss=word%100;
2013 s=b-Double_t(ddd*10000+mm*100+ss);
2014 b=Double_t(ddd)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.;
2015 }
2016
2017 if (in=="hms")
2018 {
2019 word=Int_t(b);
2020 hh=word/10000;
2021 word=word%10000;
2022 mm=word/100;
2023 ss=word%100;
2024 s=b-Double_t(hh*10000+mm*100+ss);
2025 b=15.*(Double_t(hh)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.);
2026 }
2027
2028 while (b>360)
2029 {
2030 b-=360.;
2031 }
2032
2033 if (out=="rad") b*=pi/180.;
2034
b867da20 2035 if (out=="hrs") b/=15.;
2036
8bde545d 2037 if (out=="dms")
2038 {
2039 ddd=Int_t(b);
2040 b=b-Double_t(ddd);
2041 b*=60.;
2042 mm=Int_t(b);
2043 b=b-Double_t(mm);
2044 b*=60.;
2045 ss=Int_t(b);
2046 s=b-Double_t(ss);
2047 if (s>(1.-epsilon))
2048 {
2049 s=0.;
2050 ss++;
2051 }
2052 while (ss>=60)
2053 {
2054 ss-=60;
2055 mm++;
2056 }
2057 while (mm>=60)
2058 {
2059 mm-=60;
2060 ddd++;
2061 }
2062 while (ddd>=360)
2063 {
2064 ddd-=360;
2065 }
2066 b=Double_t(10000*ddd+100*mm+ss)+s;
2067 }
2068
2069 if (out=="hms")
2070 {
2071 b/=15.;
2072 hh=Int_t(b);
2073 b=b-Double_t(hh);
2074 b*=60.;
2075 mm=Int_t(b);
2076 b=b-Double_t(mm);
2077 b*=60.;
2078 ss=Int_t(b);
2079 s=b-Double_t(ss);
2080 if (s>(1.-epsilon))
2081 {
2082 s=0.;
2083 ss++;
2084 }
2085 while (ss>=60)
2086 {
2087 ss-=60;
2088 mm++;
2089 }
2090 while (mm>=60)
2091 {
2092 mm-=60;
2093 hh++;
2094 }
2095 while (hh>=24)
2096 {
2097 hh-=24;
2098 }
2099 b=Double_t(10000*hh+100*mm+ss)+s;
2100 }
2101
2102 if (a<0) b=-b;
2103
2104 return b;
2105}
2106///////////////////////////////////////////////////////////////////////////
2107Double_t AliAstrolab::GetHourAngle(TString mode,AliTimestamp* ts,Int_t jref)
2108{
2109// Provide the Local Hour Angle (in fractional degrees) of a stored signal
2110// object at the specified timestamp.
2111//
2112// The input parameter "mode" allows the user to select either the
2113// "mean" or "apparent" value for the returned Hour Angle.
2114//
2115// mode = "M" --> Output is the Mean Hour Angle
2116// "A" --> Output is the Apparent Hour Angle
2117// ts : Timestamp at which the hour angle is requested.
2118//
2119// The input parameter "jref" allows the user to specify so-called "reference" signals.
2120// These reference signals may be used to check space-time event coincidences with the
661385a7 2121// stored measurement (e.g. coincidence of the measurement with transient phenomena).
8bde545d 2122//
661385a7 2123// jref = 0 --> Use the stored measurement
8bde545d 2124// j --> Use the reference signal at the j-th position (j=1 is first)
2125//
2126// Default value is jref=0.
2127//
2128// Note : In case ts=0 the current timestamp of the lab will be taken.
2129
2130 if (!ts) ts=(AliTimestamp*)this;
2131
2132 // Get corrected right ascension and declination for the specified timestamp.
2133 Double_t a,d;
661385a7 2134 if (mode=="M" || mode=="m") GetSignal(a,d,"M",ts,jref);
2135 if (mode=="A" || mode=="a") GetSignal(a,d,"T",ts,jref);
8bde545d 2136
2137 // Convert coordinates to fractional degrees.
2138 a=ConvertAngle(a,"hms","deg");
2139 d=ConvertAngle(d,"dms","deg");
2140
2141 a/=15.; // Convert a to fractional hours
2142 Double_t ha=0;
2143 if (mode=="M" || mode=="m") ha=ts->GetLMST(fToffset)-a;
2144 if (mode=="A" || mode=="a") ha=ts->GetLAST(fToffset)-a;
2145 ha*=15.; // Convert to (fractional) degrees
2146
2147 return ha;
2148}
2149///////////////////////////////////////////////////////////////////////////
2150void AliAstrolab::SetLT(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns,Int_t ps)
2151{
2152// Set the AliTimestamp parameters corresponding to the LT date and time
2153// in the Gregorian calendar as specified by the input arguments.
2154//
2155// The input arguments represent the following :
2156// y : year in LT (e.g. 1952, 2003 etc...)
2157// m : month in LT (1=jan 2=feb etc...)
2158// d : day in LT (1-31)
2159// hh : elapsed hours in LT (0-23)
2160// mm : elapsed minutes in LT (0-59)
2161// ss : elapsed seconds in LT (0-59)
2162// ns : remaining fractional elapsed second of LT in nanosecond
2163// ps : remaining fractional elapsed nanosecond of LT in picosecond
2164//
2165// Note : ns=0 and ps=0 are the default values.
2166//
2167
2168 SetLT(fToffset,y,m,d,hh,mm,ss,ns,ps);
2169}
2170///////////////////////////////////////////////////////////////////////////
2171void AliAstrolab::SetLT(Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps)
2172{
2173// Set the AliTimestamp parameters corresponding to the specified elapsed
2174// timespan since the beginning of the new LT year.
2175//
2176// The LT year and elapsed time span is entered via the following input arguments :
2177//
2178// y : year in LT (e.g. 1952, 2003 etc...)
2179// d : elapsed number of days
2180// s : (remaining) elapsed number of seconds
2181// ns : (remaining) elapsed number of nanoseconds
2182// ps : (remaining) elapsed number of picoseconds
2183//
2184// The specified d, s, ns and ps values will be used in an additive
2185// way to determine the elapsed timespan.
2186// So, specification of d=1, s=100, ns=0, ps=0 will result in the
2187// same elapsed time span as d=0, s=24*3600+100, ns=0, ps=0.
2188// However, by making use of the latter the user should take care
2189// of possible integer overflow problems in the input arguments,
2190// which obviously will provide incorrect results.
2191//
2192// Note : ns=0 and ps=0 are the default values.
2193
2194 SetLT(fToffset,y,d,s,ns,ps);
2195}
2196///////////////////////////////////////////////////////////////////////////
2197Double_t AliAstrolab::GetDifference(Int_t j,TString au,Double_t& dt,TString tu,Int_t mode,Int_t* ia,Int_t* it)
2198{
2199// Provide space and time difference between the j-th reference signal
2200// (j=1 indicates first) and the stored measurement.
2201//
2202// The return value of this memberfunction provides the positional angular
2203// difference, whereas the output argument "dt" provides the time difference.
2204//
2205// The units of the angular difference can be specified via the the "au"
2206// input argument, where
2207//
2208// au = "rad" --> Angular difference in (fractional) radians
2209// "deg" --> Angular difference in (fractional) degrees
2210//
2211// The units of the time difference can be specified via the "tu" and "mode"
2212// input arguments. For details please refer to AliTimestamp::GetDifference().
2213// Also here mode=1 is the default value.
2214//
2215// For the time difference the reference signal is used as the standard.
2216// This means that in case of a positive time difference, the stored
2217// measurement occurred later than the reference signal.
2218//
2219// In case j=0, the stored measurement will be compared with each
2220// reference signal and the returned angular and time differences are
2221// the minimal differences which were encountered.
2222// In this case the user may obtain the indices of the two stored reference signals
2223// which had the minimal angular and minimal time difference via the output
2224// arguments "ia" and "it" as follows :
2225//
2226// ia = Index of the stored reference signal with minimial angular difference
2227// it = Index of the stored reference signal with minimial time difference
2228//
2229// In case these indices are the same, there obviously was 1 single reference signal
2230// which showed both the minimal angular and time difference.
2231//
2232// The default values are mode=1, ia=0 and it=0;
2233
2234 Double_t dang=999;
2235 dt=1.e20;
2236 if (ia) *ia=0;
2237 if (it) *it=0;
2238
2239 if (!fRefs) return dang;
2240
2241 Ali3Vector rx; // Position of the measurement
2242 Ali3Vector r0; // Position of the reference signal
2243
2244 AliSignal* sx=GetSignal(rx,"icr","M",0);
2245 if (!sx) return dang;
2246
2247 AliTimestamp* tx=sx->GetTimestamp();
2248 if (!tx) return dang;
2249
2250 // Space and time difference w.r.t. a specific reference signal
2251 if (j>0)
2252 {
2253 AliSignal* s0=GetSignal(r0,"icr","M",0,j);
2254 if (!s0) return dang;
2255 AliTimestamp* t0=s0->GetTimestamp();
2256
2257 if (!t0) return dang;
2258
2259 dang=r0.GetOpeningAngle(rx,au);
2260 dt=t0->GetDifference(tx,tu,mode);
2261 return dang;
2262 }
2263
2264 // Minimal space and time difference encountered over all reference signals
2265 Double_t dangmin=dang;
2266 Double_t dtmin=dt;
2267 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2268 {
2269 AliSignal* s0=GetSignal(r0,"icr","M",0,i);
2270 if (!s0) continue;
2271 AliTimestamp* t0=s0->GetTimestamp();
2272 if (!t0) continue;
2273 dang=r0.GetOpeningAngle(rx,au);
2274 dt=t0->GetDifference(tx,tu,mode);
2275 if (fabs(dang)<dangmin)
2276 {
2277 dangmin=fabs(dang);
2278 *ia=i;
2279 }
2280 if (fabs(dt)<dtmin)
2281 {
2282 dtmin=fabs(dt);
2283 *it=i;
2284 }
2285 }
2286
2287 dt=dtmin;
2288 return dangmin;
2289}
2290///////////////////////////////////////////////////////////////////////////
2291Double_t AliAstrolab::GetDifference(TString name,TString au,Double_t& dt,TString tu,Int_t mode)
2292{
2293// Provide space and time difference between the reference signal
2294// with the specified name and the stored measurement.
2295//
2296// The return value of this memberfunction provides the positional angular
2297// difference, whereas the output argument "dt" provides the time difference.
2298//
2299// The units of the angular difference can be specified via the the "au"
2300// input argument, where
2301//
2302// au = "rad" --> Angular difference in (fractional) radians
2303// "deg" --> Angular difference in (fractional) degrees
2304//
2305// The units of the time difference can be specified via the "tu" and "mode"
2306// input arguments. For details please refer to AliTimestamp::GetDifference().
2307// Also here mode=1 is the default value.
2308//
2309// For the time difference the reference signal is used as the standard.
2310// This means that in case of a positive time difference, the stored
2311// measurement occurred later than the reference signal.
2312
2313 Double_t dang=999;
2314 dt=1.e20;
2315
2316 Int_t j=GetSignalIndex(name);
2317 if (j>0) dang=GetDifference(j,au,dt,tu,mode);
2318 return dang;
2319}
2320///////////////////////////////////////////////////////////////////////////
2321TArrayI* AliAstrolab::MatchRefSignal(Double_t da,TString au,Double_t dt,TString tu,Int_t mode)
2322{
2323// Provide the storage indices of the reference signals which match in space
2324// and time with the stored measurement.
2325// The indices are returned via a pointer to a TArrayI object.
2326// In case no matches were found, the null pointer is returned.
2327// A reference signal is regarded as matching with the stored measurement
2328// if the positional angular difference doesn't exceed "da" and the absolute
2329// value of the time difference doesn't exceed "dt".
2330//
2331// The units of the angular difference "da" can be specified via the the "au"
2332// input argument, where
2333//
2334// au = "rad" --> Angular difference in (fractional) radians
2335// "deg" --> Angular difference in (fractional) degrees
2336//
2337// The units of the time difference "dt" can be specified via the "tu" and "mode"
2338// input arguments. For details please refer to AliTimestamp::GetDifference().
2339// Also here mode=1 is the default value.
2340
2341 if (!fXsig || !fRefs) return 0;
2342
2343 Int_t nrefs=fRefs->GetEntries();
2344
2345 if (fIndices) delete fIndices;
2346 fIndices=new TArrayI(nrefs);
2347
2348 Double_t dang,dtime;
2349 Int_t jfill=0;
2350 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2351 {
2352 dang=GetDifference(i,au,dtime,tu,mode);
2353 if (fabs(dang)<=da && fabs(dtime)<=dt)
2354 {
2355 fIndices->AddAt(i,jfill);
2356 jfill++;
2357 }
2358 }
2359
2360 fIndices->Set(jfill);
2361 return fIndices;
2362}
2363///////////////////////////////////////////////////////////////////////////
b867da20 2364void AliAstrolab::DisplaySignal(TString frame,TString mode,AliTimestamp* ts,Int_t jref,TString proj,Int_t clr)
2365{
2366// Display a stored signal in a user specified coordinate projection
2367// at the specific timestamp ts.
2368//
2369// In case no stored signal was available or one of the input arguments was
2370// invalid, no display is produced.
2371//
2372// Note : In case ts=0 the current timestamp of the lab will be taken.
2373//
2374// The input parameter "frame" allows the user to specify the frame to which
2375// the coordinates refer. Available options are :
2376//
2377// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
2378//
2379// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
2380//
2381// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
2382//
2383// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
2384//
2385// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
2386//
2387// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
2388//
2389// In case the coordinates are the equatorial right ascension and declination (a,d),
2390// they can represent so-called "mean" and "true" values.
2391// The distinction between these two representations is the following :
2392//
2393// mean values : (a,d) are only corrected for precession and not for nutation
2394// true values : (a,d) are corrected for both precession and nutation
2395//
2396// The input parameter "mode" allows the user to specifiy either "mean" or "true"
2397// values for the input in case of equatorial (a,d) coordinates.
2398//
2399// mode = "M" --> Input coordinates are the mean values
2400// "T" --> Input coordinates are the true values
2401//
2402// The input parameter "jref" allows display of a so-called "reference" signal.
2403// These reference signals may serve to check space-time event coincidences with the
2404// stored measurement (e.g. coincidence of the measurement with transient phenomena).
2405//
2406// jref = 0 --> Display of the measurement
2407// j --> Display of the j-th reference signal
2408//
2409// Default value is jref=0.
2410//
2411// The input parameter "proj" allows the user to specify the desired projection.
2412// The available projection modes are :
2413//
2414// cyl : Cylindrical projection plotted with colored markers
2415// cylh : Cylindrical projection plotted in a 2-D histogram
2416// ham : Hammer equal area projection plotted with colored markers
2417// hamh : Hammer equal area projection plotted in a 2-D histogram
2418// ait : Aitoff projection plotted with colored markers
2419// aith : Aitoff projection plotted in a 2-D histogram
2420// mer : Mercator projection plotted with colored markers
2421// merh : Mercator projection plotted in a 2-D histogram
2422// ang : Straight cos(b) vs. l plot with colored markers
2423// angh : Straight cos(b) vs. l plot in a 2-D histogram
2424//
2425// Note : The ang(h) plot allows for easy identification of a homogeneous distribution.
2426//
2427// The input argument "clr" allows to clear (1) the display before drawing or not (0).
2428//
2429// The default values are : jref=0, proj="ham" and clr=0.
2430//
2431// This routine is based on the work by Garmt de Vries-Uiterweerd.
2432
2433 Ali3Vector r;
2434 AliSignal* sx=GetSignal(r,frame,mode,ts,jref);
2435
2436 if (!sx) return;
2437
2438 // The generic input angles (in rad) for the projections
2439 Double_t theta=0;
2440 Double_t phi=0;
2441
2442 Double_t pi=acos(-1.);
2443
2444 if (frame=="equ" || frame=="gal" || frame=="icr" || frame=="ecl" || frame=="loc")
2445 {
2446 theta=(pi/2.)-r.GetX(2,"sph","rad");
2447 phi=r.GetX(3,"sph","rad");
2448 }
2449
2450 if (frame=="hor")
2451 {
2452 theta=(pi/2.)-r.GetX(2,"sph","rad");
2453 phi=pi-r.GetX(3,"sph","rad");
2454 }
2455
2456 // Automatic choice of central meridian if not selected by the user
7f4551ca 2457 if (fMeridian<-pi)
b867da20 2458 {
2459 if (frame=="gal")
2460 {
2461 fMeridian=0;
2462 }
2463 else
2464 {
2465 fMeridian=pi;
2466 }
2467 }
2468
2469 // Obtain the projected (x,y) position
2470 Double_t x=0;
2471 Double_t y=0;
2472 Project(phi,theta,proj,x,y);
2473
2474 Int_t hist=0;
2475 if (proj=="hamh" || proj=="aith" || proj=="merh" || proj=="cylh" || proj=="angh") hist=1;
2476
2477 // Update the display for this signal position
2478
2479 // Create a new canvas if needed
2480 if (!fCanvas) fCanvas=new TCanvas("AliAstrolab","Skymap");
2481
2482 // Construct the title string for this map
2483 TString title="Projection : ";
2484 title+=frame;
2485 title+=" ";
2486 title+=proj;
2487 title+=" Center : ";
2488 Int_t ang,h,m,s,d;
2489 if (frame=="equ" || frame=="icr")
2490 {
2491 ang=int(ConvertAngle(fMeridian,"rad","hms"));
2492 h=ang/10000;
2493 ang=ang%10000;
2494 m=ang/100;
2495 s=ang%100;
2496 title+=h;
2497 title+="h ";
2498 title+=m;
2499 title+="m ";
2500 title+=s;
2501 title+="s";
2502 }
2503 else
2504 {
2505 ang=int(ConvertAngle(fMeridian,"rad","dms"));
2506 d=ang/10000;
2507 ang=ang%10000;
2508 m=ang/100;
2509 s=ang%100;
2510 title+=d;
2511 title+="d ";
2512 title+=m;
2513 title+="' ";
2514 title+=s;
2515 title+="\"";
2516 }
2517
2518 if (!hist) // 2-D Marker display (i.e. not a histogram)
2519 {
2520 // Remove existing markers, grid and outline from display if needed
2521 if (clr==1 || proj!=fProj)
2522 {
2523 if (fMarkers)
2524 {
2525 delete fMarkers;
2526 fMarkers=0;
2527 }
2528 fCanvas->Clear();
2529 fProj=proj;
2530 }
2531
2532 if (!fMarkers)
2533 {
2534 fMarkers=new TObjArray();
2535 fMarkers->SetOwner();
2536 // Set canvas range, header and axes
2537 if (proj=="mer")
2538 {
2539 fCanvas->Range(-2.2,-6.,2.2,6.);
2540 TText* header=new TText;
2541 fMarkers->Add(header);
2542 header->DrawText(-1.4,5.5,title.Data());
2543 TLine* line=new TLine(-2,0,2,0);
2544 fMarkers->Add(line);
2545 line->Draw();
2546 line=new TLine(0,5.5,0,-5.5);
2547 fMarkers->Add(line);
2548 line->Draw();
2549 }
2550 else
2551 {
2552 fCanvas->Range(-2.2,-1.2,2.2,1.2);
2553 TText* header=new TText;
2554 fMarkers->Add(header);
2555 header->DrawText(-1.4,1.1,title.Data());
2556 TLine* line=new TLine(-2,0,2,0);
2557 fMarkers->Add(line);
2558 line->Draw();
2559 line=new TLine(0,1,0,-1);
2560 fMarkers->Add(line);
2561 line->Draw();
2562 }
2563 // Draw the outline
2564 if (proj=="ham" || proj=="ait")
2565 {
2566 // Draw ellips outline with axes
2567 TEllipse* outline=new TEllipse(0,0,2,1);
2568 fMarkers->Add(outline);
2569 outline->Draw();
2570 }
2571 }
2572 TMarker* marker=0;
2573 if (!jref) // The measurement
2574 {
2575 marker=new TMarker(x,y,29);
2576 marker->SetMarkerColor(kRed);
2577 }
2578 else // The reference signal(s)
2579 {
2580 marker=new TMarker(x,y,20);
2581 marker->SetMarkerColor(kBlue);
2582 }
2583 marker->SetMarkerSize(1);
2584 fMarkers->Add(marker);
2585 marker->Draw();
2586 }
2587 else if (hist==1) // 2-D display via histogram
2588 {
2589 // Reset the histogram if needed
2590 if (clr==1 || proj!=fProj)
2591 {
2592 fCanvas->Clear();
2593 if (!fHist) fHist=new TH2F();
2594 fHist->Reset();
2595 fHist->SetMarkerStyle(20);
2596 fHist->SetMarkerSize(1);
2597 fHist->SetMarkerColor(kBlue);
2598 fHist->SetNameTitle("SkyMap",title.Data());
2599 if (proj=="merh")
2600 {
2601 fHist->SetBins(100,-2.2,2.2,100,-5.5,5.5);
2602 }
2603 else
2604 {
2605 fHist->SetBins(100,-2.2,2.2,100,-1.1,1.1);
2606 }
2607 fProj=proj;
2608 }
2609 fHist->Fill(x,y);
2610 fHist->Draw();
2611 }
2612}
2613///////////////////////////////////////////////////////////////////////////
2614void AliAstrolab::DisplaySignal(TString frame,TString mode,AliTimestamp* ts,TString name,TString proj,Int_t clr)
2615{
2616// Display the stored signal with the specified name in a user specified
2617// coordinate projection at the specific timestamp ts.
2618//
2619// Note : In case ts=0 the current timestamp of the lab will be taken.
2620//
2621// In case no such stored signal was available or one of the input arguments was
2622// invalid, no display is produced.
2623//
2624// The input parameter "frame" allows the user to specify the frame to which
2625// the coordinates refer. Available options are :
2626//
2627// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
2628//
2629// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
2630//
2631// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
2632//
2633// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
2634//
2635// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
2636//
2637// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
2638//
2639// In case the coordinates are the equatorial right ascension and declination (a,d),
2640// they can represent so-called "mean" and "true" values.
2641// The distinction between these two representations is the following :
2642//
2643// mean values : (a,d) are only corrected for precession and not for nutation
2644// true values : (a,d) are corrected for both precession and nutation
2645//
2646// The input parameter "mode" allows the user to specifiy either "mean" or "true"
2647// values for the input in case of equatorial (a,d) coordinates.
2648//
2649// mode = "M" --> Input coordinates are the mean values
2650// "T" --> Input coordinates are the true values
2651//
2652// The input parameter "proj" allows the user to specify the desired projection.
2653// The available projection modes are :
2654//
2655// cyl : Cylindrical projection plotted with colored markers
2656// cylh : Cylindrical projection plotted in a 2-D histogram
2657// ham : Hammer equal area projection plotted with colored markers
2658// hamh : Hammer equal area projection plotted in a 2-D histogram
2659// ait : Aitoff projection plotted with colored markers
2660// aith : Aitoff projection plotted in a 2-D histogram
2661// mer : Mercator projection plotted with colored markers
2662// merh : Mercator projection plotted in a 2-D histogram
2663// ang : Straight cos(b) vs. l plot with colored markers
2664// angh : Straight cos(b) vs. l plot in a 2-D histogram
2665//
2666// Note : The ang(h) plot allows for easy identification of a homogeneous distribution.
2667//
2668// The input argument "clr" allows to clear (1) the display before drawing or not (0).
2669//
2670// The default values are : proj="ham" and clr=0.
2671//
2672// This routine is based on the work by Garmt de Vries-Uiterweerd.
2673
2674 Int_t j=GetSignalIndex(name);
2675 if (j>=0) DisplaySignal(frame,mode,ts,j,proj,clr);
2676}
2677///////////////////////////////////////////////////////////////////////////
2678void AliAstrolab::DisplaySignals(TString frame,TString mode,AliTimestamp* ts,TString proj,Int_t clr)
2679{
2680// Display of all stored signals in user specified coordinate projection
2681// at the specific timestamp ts.
2682// In case ts=0, the timestamp of the actual recording of the stored measurement
2683// under investigation will be used.
2684// In case no (timestamp of the) actual measurement is available,
2685// the current timestamp of the lab will be taken.
2686// In case no stored signal is available or one of the input arguments is
2687// invalid, no display is produced.
2688//
2689// The input parameter "frame" allows the user to specify the frame to which
2690// the coordinates refer. Available options are :
2691//
2692// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
2693//
2694// "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
2695//
2696// "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
2697//
2698// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
2699//
2700// "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
2701//
2702// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
2703//
2704// In case the coordinates are the equatorial right ascension and declination (a,d),
2705// they can represent so-called "mean" and "true" values.
2706// The distinction between these two representations is the following :
2707//
2708// mean values : (a,d) are only corrected for precession and not for nutation
2709// true values : (a,d) are corrected for both precession and nutation
2710//
2711// The input parameter "mode" allows the user to specifiy either "mean" or "true"
2712// values for the input in case of equatorial (a,d) coordinates.
2713//
2714// mode = "M" --> Input coordinates are the mean values
2715// "T" --> Input coordinates are the true values
2716//
2717// The input parameter "proj" allows the user to specify the desired projection.
2718// The available projection modes are :
2719//
2720// cyl : Cylindrical projection plotted with colored markers
2721// cylh : Cylindrical projection plotted in a 2-D histogram
2722// ham : Hammer equal area projection plotted with colored markers
2723// hamh : Hammer equal area projection plotted in a 2-D histogram
2724// ait : Aitoff projection plotted with colored markers
2725// aith : Aitoff projection plotted in a 2-D histogram
2726// mer : Mercator projection plotted with colored markers
2727// merh : Mercator projection plotted in a 2-D histogram
2728// ang : Straight cos(b) vs. l plot with colored markers
2729// angh : Straight cos(b) vs. l plot in a 2-D histogram
2730//
2731// Note : The ang(h) plot allows for easy identification of a homogeneous distribution.
2732//
2733// The input argument "clr" allows to clear (1) the display before drawing or not (0).
2734//
2735// The default values are : proj="ham" and clr=0.
2736//
2737// This routine is based on the work by Garmt de Vries-Uiterweerd.
2738
2739 AliTimestamp* tx=ts;
2740 if (fXsig && !tx) tx=fXsig->GetTimestamp();
2741 if (!tx) tx=(AliTimestamp*)this;
2742
2743 // Display all stored reference signals
2744 if (fRefs)
2745 {
2746 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2747 {
2748 DisplaySignal(frame,mode,tx,i,proj,clr);
2749 clr=0; // No display clear for subsequent signals
2750 }
2751 }
2752
2753 // Display the measurement
2754 DisplaySignal(frame,mode,tx,0,proj,clr);
2755}
2756///////////////////////////////////////////////////////////////////////////
2757void AliAstrolab::SetCentralMeridian(Double_t phi,TString u)
2758{
2759// Set the central meridian for the sky display.
7f4551ca 2760// Setting a value smaller than -pi will induce automatic meridian setting
b867da20 2761// in the display.
2762// By default the central meridian is set at -999 in the constructor.
2763//
2764// The string argument "u" allows to choose between different angular units
2765// u = "rad" : angle provided in radians
2766// "deg" : angle provided in degrees
2767// "dms" : angle provided in dddmmss.sss
2768// "hms" : angle provided in hhmmss.sss
2769//
2770// The default is u="deg".
2771//
2772// This routine is based on the work by Garmt de Vries-Uiterweerd.
2773
2774 fMeridian=ConvertAngle(phi,u,"rad");
2775}
2776///////////////////////////////////////////////////////////////////////////
2777void AliAstrolab::Project(Double_t l,Double_t b,TString proj,Double_t& x,Double_t& y)
2778{
2779// Generic interface for projection of a (long,lat) pair onto a (x,y) pair.
2780//
2781// Meaning of the arguments :
2782// --------------------------
2783// l : Longitude (e.g. right ascension)
2784// b : Latitude (e.g. declination)
2785// proj : Projection mode (e.g. "ham")
2786// x : The resulting x coordinate for the selected projection
2787// y : The resulting x coordinate for the selected projection
2788//
2789// The available projection modes are :
2790//
2791// cyl : Cylindrical projection plotted with colored markers
2792// cylh : Cylindrical projection plotted in a 2-D histogram
2793// ham : Hammer equal area projection plotted with colored markers
2794// hamh : Hammer equal area projection plotted in a 2-D histogram
2795// ait : Aitoff projection plotted with colored markers
2796// aith : Aitoff projection plotted in a 2-D histogram
2797// mer : Mercator projection plotted with colored markers
2798// merh : Mercator projection plotted in a 2-D histogram
2799// ang : Straight cos(b) vs. l plot with colored markers
2800// angh : Straight cos(b) vs. l plot in a 2-D histogram
2801//
2802// Note : The ang(h) plot allows for easy identification of a homogeneous distribution.
2803//
2804// This routine is based on the work by Garmt de Vries-Uiterweerd.
2805
2806 Double_t pi=acos(-1.);
2807
2808 // Subtract central meridian from longitude
2809 l-=fMeridian;
2810
2811 // Take l between -180 and 180 degrees
2812 while (l>pi)
2813 {
2814 l-=2*pi;
2815 }
2816 while (l<-pi)
2817 {
2818 l+=2*pi;
2819 }
2820
2821 x=0;
2822 y=0;
2823
2824 // Convert (l,b) to (x,y) with -2 < x <= 2
2825 if (proj=="cyl" || proj=="cylh") ProjectCylindrical(l,b,x,y);
2826 if (proj=="ham" || proj=="hamh") ProjectHammer(l,b,x,y);
2827 if (proj=="ait" || proj=="aith") ProjectAitoff(l,b,x,y);
2828 if (proj=="mer" || proj=="merh") ProjectMercator(l,b,x,y);
2829 if (proj=="ang" || proj=="angh")
2830 {
2831 x=2*l/pi;
2832 y=cos(b);
2833 }
2834}
2835///////////////////////////////////////////////////////////////////////////
2836void AliAstrolab::ProjectCylindrical(Double_t l,Double_t b,Double_t& x, Double_t& y)
2837{
2838// Cylindrical projection of a (long,lat) coordinate pair
2839//
2840// This routine is based on the work by Garmt de Vries-Uiterweerd.
2841
2842 Double_t pi=acos(-1.);
2843 x=2*l/pi;
2844 y=2*b/pi;
2845}
2846///////////////////////////////////////////////////////////////////////////
2847void AliAstrolab::ProjectHammer(Double_t l,Double_t b,Double_t& x,Double_t& y)
2848{
2849// Hammer-Aitoff projection of a (long,lat) coordinate pair.
2850// This is an equal-area projection.
2851//
2852// This routine is based on the work by Garmt de Vries-Uiterweerd.
2853
2854 Double_t k=1./sqrt(1.+cos(b)*cos(l/2.));
2855 x=2.*k*cos(b)*sin(l/2.);
2856 y=k*sin(b);
2857}
2858///////////////////////////////////////////////////////////////////////////
2859void AliAstrolab::ProjectAitoff(Double_t l,Double_t b,Double_t& x,Double_t& y)
2860{
2861// Aitoff projection of a (long,lat) coordinate pair.
2862//
2863// This routine is based on the work by Garmt de Vries-Uiterweerd.
2864
2865 Double_t pi=acos(-1.);
2866 x=0;
2867 y=0;
2868 Double_t k=acos(cos(b)*cos(l/2.));
2869 if(sin(k)!=0)
2870 {
2871 x=4.*k*cos(b)*sin(l/2.)/(pi*sin(k));
2872 y=2.*k*sin(b)/(pi*sin(k));
2873 }
2874}
2875///////////////////////////////////////////////////////////////////////////
2876void AliAstrolab::ProjectMercator(Double_t l,Double_t b,Double_t& x,Double_t& y)
2877{
2878// Mercator projection of a (long,lat) coordinate pair
2879//
2880// This routine is based on the work by Garmt de Vries-Uiterweerd.
2881
2882 Double_t pi=acos(-1.);
2883 x=2.*l/pi;
2884 y=0;
2885 if (b >= 89.*pi/180.) b=89.*pi/180.;
2886 if (b <= -89.*pi/180.) y=-89.*pi/180.;
2887 if (fabs(y) < 1) y=log((1.+sin(b))/(1.-sin(b)))/2.;
2888}
2889///////////////////////////////////////////////////////////////////////////
8bde545d 2890TObject* AliAstrolab::Clone(const char* name) const
2891{
2892// Make a deep copy of the current object and provide the pointer to the copy.
2893// This memberfunction enables automatic creation of new objects of the
2894// correct type depending on the object type, a feature which may be very useful
2895// for containers when adding objects in case the container owns the objects.
2896
2897 AliAstrolab* lab=new AliAstrolab(*this);
2898 if (name)
2899 {
2900 if (strlen(name)) lab->SetName(name);
2901 }
2902 return lab;
2903}
2904///////////////////////////////////////////////////////////////////////////