]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/AliAstrolab.cxx
A few modifications in the macro + merging both macros into one
[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.
50//
51// Coding example :
52// ----------------
53// gSystem->Load("ralice");
54//
55// AliAstrolab lab("IceCube","The South Pole Neutrino Observatory");
56// lab.SetLabPosition(0,-90,"deg"); // South Pole
57//
58// lab.SetLocalFrame(90,180,90,270,0,0); // Local frame has X-axis to the North
59//
60// lab.Data(1,"dms"); // Print laboratory parameters
61//
62// // Enter some observed event to be investigated
63// AliTimestamp ts;
64// ts.SetUT(1989,7,30,8,14,23,738504,0);
65// Float_t vec[3]={1,23.8,118.65};
66// Ali3Vector r;
67// r.SetVector(vec,"sph","deg");
68// lab.SetSignal(&r,"loc","M",&ts,0,"Event10372");
69//
70// // Enter some reference signals
71// Float_t alpha=194818.0;
72// Float_t delta=84400.;
73// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"Altair");
74// alpha=124900.0;
75// delta=272400.;
76// lab.SetSignal(alpha,delta,"B",1950,"M",-1,"NGP");
77// alpha=64508.917;
78// delta=-164258.02;
79// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Sirius");
80// alpha=23149.08;
81// delta=891550.8;
82// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Polaris");
83// alpha=43600.;
84// delta=163100.;
85// lab.SetSignal(alpha,delta,"J",2000,"M",-1,"Aldebaran");
86// Float_t l=327.531;
87// Float_t b=-35.8903;
88// Float_t pos[3]={1,90.-b,l};
89// r.SetVector(pos,"sph","deg");
90// lab.SetUT(1989,7,30,8,14,16,0,0);
caa58e1a 91// lab.SetSignal(&r,"gal","M",0,-1,"GRB890730");
8bde545d 92//
93// // List all stored objects
94// lab.ListSignals("equ","M",5);
95//
96// // Determine minimal space and time differences with reference signals
97// Double_t da,dt;
98// Int_t ia,it;
99// da=lab.GetDifference(0,"deg",dt,"s",1,&ia,&it);
100// cout << " Minimal differences damin (deg) : " << da << " dtmin (s) : " << dt
101// << " damin-index : " << ia << " dtmin-index : " << it << endl;
102// cout << " damin for "; lab->PrintSignal("equ","T",&ts,5,ia); cout << endl;
103// cout << " dtmin for "; lab->PrintSignal("equ","T",&ts,5,it); cout << endl;
104//
105// // Search for space and time match with the reference signals
106// da=5;
107// dt=10;
108// TArrayI* arr=lab.MatchRefSignal(da,"deg",dt,"s");
109// Int_t index=0;
110// if (arr)
111// {
112// for (Int_t i=0; i<arr->GetSize(); i++)
113// {
114// index=arr->At(i);
115// cout << " Match found for index : " << index << endl;
116// cout << " Corresponding ref. object "; lab->PrintSignal("equ","T",&ts,5,index); cout << endl;
117// }
118// }
119//
120//
121//--- Author: Nick van Eijndhoven 15-mar-2007 Utrecht University
122//- Modified: NvE $Date$ Utrecht University
123///////////////////////////////////////////////////////////////////////////
124
125#include "AliAstrolab.h"
126#include "Riostream.h"
127
128ClassImp(AliAstrolab) // Class implementation to enable ROOT I/O
129
caa58e1a 130AliAstrolab::AliAstrolab(const char* name,const char* title) : TTask(name,title),AliTimestamp()
8bde545d 131{
132// Default constructor
133
134 fToffset=0;
135 fXsig=0;
136 fRefs=0;
137 fBias=0;
138 fGal=0;
139 fIndices=0;
140}
141///////////////////////////////////////////////////////////////////////////
142AliAstrolab::~AliAstrolab()
143{
144// Destructor to delete all allocated memory.
145
146 if (fXsig)
147 {
148 delete fXsig;
149 fXsig=0;
150 }
151 if (fRefs)
152 {
153 delete fRefs;
154 fRefs=0;
155 }
156 if (fIndices)
157 {
158 delete fIndices;
159 fIndices=0;
160 }
161}
162///////////////////////////////////////////////////////////////////////////
caa58e1a 163AliAstrolab::AliAstrolab(const AliAstrolab& t) : TTask(t),AliTimestamp(t)
8bde545d 164{
165// Copy constructor
166
167 fToffset=t.fToffset;
168 fLabPos=t.fLabPos;
169 fXsig=0;
170 if (t.fXsig) fXsig=new AliSignal(*(t.fXsig));
171 fRefs=0;
172 if (t.fRefs)
173 {
174 Int_t size=t.fRefs->GetSize();
175 fRefs=new TObjArray(size);
176 for (Int_t i=0; i<size; i++)
177 {
178 AliSignal* sx=(AliSignal*)t.fRefs->At(i);
179 if (sx) fRefs->AddAt(sx->Clone(),i);
180 }
181 }
182 fBias=0;
183 fGal=0;
184 fIndices=0;
185}
186///////////////////////////////////////////////////////////////////////////
187void AliAstrolab::Data(Int_t mode,TString u)
188{
189// Provide lab information.
190//
191// "mode" indicates the mode of the timestamp info (see AliTimestamp::Date).
192//
193// The string argument "u" allows to choose between different angular units
194// in case e.g. a spherical frame is selected.
195// u = "rad" : angles provided in radians
196// "deg" : angles provided in degrees
197// "dms" : angles provided in ddd:mm:ss.sss
198// "hms" : angles provided in hh:mm:ss.sss
199//
200// The defaults are mode=1 and u="deg".
201
202 const char* name=GetName();
203 const char* title=GetTitle();
204 cout << " *" << ClassName() << "::Data*";
205 if (strlen(name)) cout << " Name : " << GetName();
206 if (strlen(title)) cout << " Title : " << GetTitle();
207 cout << endl;
208
209 Double_t l,b;
210 GetLabPosition(l,b,"deg");
211 cout << " Lab position longitude : "; PrintAngle(l,"deg",u,2);
212 cout << " latitude : "; PrintAngle(b,"deg",u,2);
213 cout << endl;
214 cout << " Lab time offset w.r.t. UT : "; PrintTime(fToffset,12); cout << endl;
215
216 // UT and Local time info
217 Date(mode,fToffset);
218}
219///////////////////////////////////////////////////////////////////////////
220void AliAstrolab::SetLabPosition(Ali3Vector& p)
221{
222// Set the lab position in the terrestrial coordinates.
223// The right handed reference frame is defined such that the North Pole
224// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
225// to an azimuth angle phi=0, with phi increasing eastwards.
226
227 fLabPos.SetPosition(p);
228
229 // Determine local time offset in fractional hours w.r.t. UT.
230 Double_t vec[3];
231 p.GetVector(vec,"sph","deg");
232 Double_t l=vec[2];
233 fToffset=l/15.;
234}
235///////////////////////////////////////////////////////////////////////////
236void AliAstrolab::SetLabPosition(Double_t l,Double_t b,TString u)
237{
238// Set the lab position in the terrestrial longitude (l) and latitude (b).
239// Positions north of the equator have b>0, whereas b<0 indicates
240// locations south of the equator.
241// Positions east of Greenwich have l>0, whereas l<0 indicates
242// locations west of Greenwich.
243//
244// The string argument "u" allows to choose between different angular units
245// u = "rad" : angles provided in radians
246// "deg" : angles provided in degrees
247// "dms" : angles provided in dddmmss.sss
248// "hms" : angles provided in hhmmss.sss
249//
250// The default is u="deg".
251
252 Double_t r=1,theta=0,phi=0;
253
254 l=ConvertAngle(l,u,"deg");
255 b=ConvertAngle(b,u,"deg");
256
257 Double_t offset=90.;
258
259 theta=offset-b;
260 phi=l;
261
262 Double_t p[3]={r,theta,phi};
263 fLabPos.SetPosition(p,"sph","deg");
264
265 // Local time offset in fractional hours w.r.t. UT.
266 fToffset=l/15.;
267}
268///////////////////////////////////////////////////////////////////////////
269AliPosition AliAstrolab::GetLabPosition() const
270{
271// Provide the lab position in the terrestrial coordinates.
272// The right handed reference frame is defined such that the North Pole
273// corresponds to a polar angle theta=0 and the Greenwich meridian corresponds
274// to an azimuth angle phi=0, with phi increasing eastwards.
275
276 return fLabPos;
277}
278///////////////////////////////////////////////////////////////////////////
279void AliAstrolab::GetLabPosition(Double_t& l,Double_t& b,TString u) const
280{
281// Provide the lab position in the terrestrial longitude (l) and latitude (b).
282// Positions north of the equator have b>0, whereas b<0 indicates
283// locations south of the equator.
284// Positions east of Greenwich have l>0, whereas l<0 indicates
285// locations west of Greenwich.
286//
287// The string argument "u" allows to choose between different angular units
288// u = "rad" : angles provided in radians
289// "deg" : angles provided in degrees
290//
291// The default is u="deg".
292
293 Double_t pi=acos(-1.);
294
295 Double_t offset=90.;
296 if (u=="rad") offset=pi/2.;
297
298 Double_t p[3];
299 fLabPos.GetPosition(p,"sph",u);
300 b=offset-p[1];
301 l=p[2];
302}
303///////////////////////////////////////////////////////////////////////////
304Double_t AliAstrolab::GetLT()
305{
306// Provide the Lab's local time in fractional hours.
307// A mean solar day lasts 24h (i.e. 86400s).
308//
309// In case a hh:mm:ss format is needed, please use the Convert() facility.
310
311 Double_t h=GetLT(fToffset);
312 return h;
313}
314///////////////////////////////////////////////////////////////////////////
315Double_t AliAstrolab::GetLMST()
316{
317// Provide the Lab's Local Mean Sidereal Time (LMST) in fractional hours.
318// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
319// The definition of GMST is such that a sidereal clock corresponds with
320// 24 sidereal hours per revolution of the Earth.
321// As such, local time offsets w.r.t. UT and GMST can be treated similarly.
322//
323// In case a hh:mm:ss format is needed, please use the Convert() facility.
324
325 Double_t h=GetLMST(fToffset);
326 return h;
327}
328///////////////////////////////////////////////////////////////////////////
329Double_t AliAstrolab::GetLAST()
330{
331// Provide the Lab's Local Apparent Sidereal Time (LAST) in fractional hours.
332// A sidereal day corresponds to 23h 56m 04.09s (i.e. 86164.09s) mean solar time.
333// The definition of GMST and GAST is such that a sidereal clock corresponds with
334// 24 sidereal hours per revolution of the Earth.
335// As such, local time offsets w.r.t. UT, GMST and GAST can be treated similarly.
336//
337// In case a hh:mm:ss format is needed, please use the Convert() facility.
338
339 Double_t h=GetLAST(fToffset);
340 return h;
341}
342///////////////////////////////////////////////////////////////////////////
343void AliAstrolab::PrintAngle(Double_t a,TString in,TString out,Int_t ndig) const
344{
345// Printing of angles in various formats.
346//
347// The input argument "a" denotes the angle to be printed.
348// The string arguments "in" and "out" specify the angular I/O formats.
349//
350// in = "rad" : input angle provided in radians
351// "deg" : input angle provided in degrees
352// "dms" : input angle provided in dddmmss.sss
353// "hms" : input angle provided in hhmmss.sss
354//
355// out = "rad" : output angle provided in radians
356// "deg" : output angle provided in degrees
357// "dms" : output angle provided in dddmmss.sss
358// "hms" : output angle provided in hhmmss.sss
359//
360// The argument "ndig" specifies the number of digits for the fractional
361// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
362// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
363// will appear as 03.4 on the output.
364// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
365//
366// The default is ndig=1.
367//
368// Note : The angle info is printed without additional spaces or "endline".
369// This allows the print to be included in various composite output formats.
370
371 Double_t b=ConvertAngle(a,in,out);
372
373 if (out=="deg" || out=="rad")
374 {
375 cout.setf(ios::fixed,ios::floatfield);
376 cout << setprecision(ndig) << b << " " << out.Data();
377 cout.unsetf(ios::fixed);
378 return;
379 }
380
381 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
382 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
383 Double_t s;
384 ULong64_t sfrac=0;
385
386 if (out=="dms")
387 {
388 word=Int_t(b);
389 word=abs(word);
390 ddd=word/10000;
391 word=word%10000;
392 mm=word/100;
393 ss=word%100;
394 s=fabs(b)-Double_t(ddd*10000+mm*100+ss);
395 if (s>(1.-epsilon))
396 {
397 s=0.;
398 ss++;
399 }
400 while (ss>=60)
401 {
402 ss-=60;
403 mm++;
404 }
405 while (mm>=60)
406 {
407 mm-=60;
408 ddd++;
409 }
410 while (ddd>=360)
411 {
412 ddd-=360;
413 }
414 s*=pow(10.,ndig);
415 sfrac=ULong64_t(s);
416 if (b<0) cout << "-";
417 cout << ddd << "d " << mm << "' " << ss << "."
418 << setfill('0') << setw(ndig) << sfrac << "\"";
419 return;
420 }
421
422 if (out=="hms")
423 {
424 word=Int_t(b);
425 word=abs(word);
426 hh=word/10000;
427 word=word%10000;
428 mm=word/100;
429 ss=word%100;
430 s=fabs(b)-Double_t(hh*10000+mm*100+ss);
431 if (s>(1.-epsilon))
432 {
433 s=0.;
434 ss++;
435 }
436 while (ss>=60)
437 {
438 ss-=60;
439 mm++;
440 }
441 while (mm>=60)
442 {
443 mm-=60;
444 hh++;
445 }
446 while (hh>=24)
447 {
448 hh-=24;
449 }
450 s*=pow(10.,ndig);
451 sfrac=ULong64_t(s);
452 if (b<0) cout << "-";
453 cout << hh << "h " << mm << "m " << ss << "."
454 << setfill('0') << setw(ndig) << sfrac << "s";
455 return;
456 }
457}
458///////////////////////////////////////////////////////////////////////////
459void AliAstrolab::SetSignal(Ali3Vector* r,TString frame,TString mode,AliTimestamp* ts,Int_t jref,TString name)
460{
461// Store a signal as specified by the position r and the timestamp ts.
462// The position is stored in International Celestial Reference System (ICRS) coordinates.
463// The ICRS is a fixed, time independent frame and as such provides a unique reference
464// frame without the need of specifying any epoch etc...
465// The ICRS coordinate definitions match within 20 mas with the mean ones of the J2000.0
466// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
467// coordinate correction between J2000 and ICRS is performed here via the
468// so-called frame bias matrix.
469// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
470// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
471//
472// The input parameter "frame" allows the user to specify the frame to which
473// the components of r refer. Available options are :
474//
475// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
476// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
477// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
478// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
479// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
480// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
481// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
482// components of r correspond to theta=zenith angle and phi=pi-azimuth.
483// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
484// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
485// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
486// components of r correspond to the usual theta and phi angles.
487//
488// In case the coordinates are the equatorial right ascension and declination (a,d),
489// they can represent so-called "mean" and "true" values.
490// The distinction between these two representations is the following :
491//
492// mean values : (a,d) are only corrected for precession and not for nutation
493// true values : (a,d) are corrected for both precession and nutation
494//
495// The input parameter "mode" allows the user to specifiy either "mean" or "true"
496// values for the input in case of equatorial (a,d) coordinates.
497//
498// mode = "M" --> Input coordinates are the mean values
499// "T" --> Input coordinates are the true values
500//
501// The input parameter "jref" allows the user to store so-called "reference" signals.
502// These reference signals may be used to check space-time event coincidences with the
503// stored measurement (e.g. coincidence of the measurement with transient phenomena).
504//
505// jref = 0 --> Storage of the measurement
506// j --> Storage of a reference signal at the j-th position (j=1 is first)
507// < 0 --> Add a reference signal at the next available position
508//
509// Via the input argument "name" the user can give the stored signal also a name.
510//
511// The default values are jref=0 and name="".
512//
513// Note : In case ts=0 the current timestamp of the lab will be taken.
514
515 if (!r)
516 {
517 if (!jref && fXsig)
518 {
519 delete fXsig;
520 fXsig=0;
521 }
522 return;
523 }
524
525 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc")
526 {
527 if (!jref && fXsig)
528 {
529 delete fXsig;
530 fXsig=0;
531 }
532 return;
533 }
534
535 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t")
536 {
537 if (!jref && fXsig)
538 {
539 delete fXsig;
540 fXsig=0;
541 }
542 return;
543 }
544
545 if (!ts) ts=(AliTimestamp*)this;
546
547 Double_t vec[3]={1,0,0};
548 vec[1]=r->GetX(2,"sph","rad");
549 vec[2]=r->GetX(3,"sph","rad");
550 Ali3Vector q;
551 q.SetVector(vec,"sph","rad");
552
553 AliSignal* sxref=0;
554
555 if (!jref) // Storage of the measurement
556 {
557 if (!fXsig)
558 {
559 fXsig=new AliSignal();
560 }
561 else
562 {
563 fXsig->Reset(1);
564 }
565 if (name != "") fXsig->SetName(name);
566 fXsig->SetTitle("Event in ICRS coordinates");
567 fXsig->SetTimestamp(*ts);
568 }
569 else // Storage of a reference signal
570 {
571 if (!fRefs)
572 {
573 fRefs=new TObjArray();
574 fRefs->SetOwner();
575 }
576 // Expand array size if needed
577 if (jref>0 && jref>=fRefs->GetSize()) fRefs->Expand(jref+1);
578 sxref=new AliSignal();
579 if (name != "") sxref->SetName(name);
580 sxref->SetTitle("Reference event in ICRS coordinates");
581 sxref->SetTimestamp(*ts);
582 }
583
584 if (frame=="loc")
585 {
586 // Convert to horizontal coordinates
587 q=q.GetUnprimed(&fL);
588
589 // Store the signal
590 SetSignal(&q,"hor",mode,ts,jref);
591 return;
592 }
593
594 if (frame=="equ")
595 {
596 // Convert to "mean" values at specified epoch
597 if (mode=="T" || mode=="t")
598 {
599 SetNmatrix(ts);
600 q=q.GetUnprimed(&fN);
601 }
602
603 // Convert to "mean" values at J2000
604 SetPmatrix(ts);
605 q=q.GetUnprimed(&fP);
606
607 // Convert to ICRS values
608 if (!fBias) SetBmatrix();
609 q=q.GetUnprimed(&fB);
610 }
611
612 if (frame=="gal")
613 {
614 // Convert to J2000 equatorial mean coordinates
615 if (fGal != 2) SetGmatrix("J");
616 q=q.GetUnprimed(&fG);
617
618 // Convert to ICRS values
619 if (!fBias) SetBmatrix();
620 q=q.GetUnprimed(&fB);
621 }
622
623 if (frame=="ecl")
624 {
625 // Convert to mean equatorial values at specified epoch
626 SetEmatrix(ts);
627 q=q.GetUnprimed(&fE);
628
629 // Convert to "mean" values at J2000
630 SetPmatrix(ts);
631 q=q.GetUnprimed(&fP);
632
633 // Convert to ICRS values
634 if (!fBias) SetBmatrix();
635 q=q.GetUnprimed(&fB);
636 }
637
638 if (frame=="hor")
639 {
640 // Convert to "true" equatorial values at the specified timestamp
641 SetHmatrix(ts);
642 q=q.GetUnprimed(&fH);
643
644 // Convert to "mean" values at specified timestamp
645 SetNmatrix(ts);
646 q=q.GetUnprimed(&fN);
647
648 // Convert to "mean" values at J2000
649 SetPmatrix(ts);
650 q=q.GetUnprimed(&fP);
651
652 // Convert to ICRS values
653 if (!fBias) SetBmatrix();
654 q=q.GetUnprimed(&fB);
655 }
656
657 // Store the signal in ICRS coordinates
658 if (!jref) // Storage of a regular signal
659 {
660 fXsig->SetPosition(q);
661 }
662 else // Storage of a reference signal
663 {
664 sxref->SetPosition(q);
665 if (jref<0)
666 {
667 fRefs->Add(sxref);
668 }
669 else
670 {
671 fRefs->AddAt(sxref,jref-1);
672 }
673 }
674}
675///////////////////////////////////////////////////////////////////////////
676void AliAstrolab::SetSignal(Double_t a,Double_t d,TString s,Double_t e,TString mode,Int_t jref,TString name)
677{
678// Store a signal with right ascension (a) and declination (d) given for epoch e.
679// The position is stored in International Celestial Reference System (ICRS) coordinates.
680// The ICRS is a fixed, time independent frame and as such provides a unique reference
681// frame without the need of specifying any epoch etc...
682// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
683// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
684// coordinate correction between J2000 and ICRS is performed here via the
685// so-called frame bias matrix.
686// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
687// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
688//
689// The coordinates (a,d) can represent so-called "mean" and "true" values.
690// The distinction between these two representations is the following :
691//
692// mean values : (a,d) are only corrected for precession and not for nutation
693// true values : (a,d) are corrected for both precession and nutation
694//
695// The input parameter "mode" allows the user to specifiy either "mean" or "true"
696// values for the input (a,d) coordinates.
697//
698// a : Right ascension in hhmmss.sss
699// d : Declination in dddmmss.sss
700// s = "B" --> Besselian reference epoch.
701// "J" --> Julian reference epoch.
702// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
703// mode = "M" --> Input coordinates are the mean values
704// "T" --> Input coordinates are the true values
705//
706// The input parameter "jref" allows the user to store so-called "reference" signals.
707// These reference signals may be used to check space-time event coincidences with the
708// stored measurement (e.g. coincidence of the measurement with transient phenomena).
709//
710// jref = 0 --> Storage of the measurement
711// j --> Storage of a reference signal at the j-th position (j=1 is first)
712// < 0 --> Add a reference signal at the next available position
713//
714// Via the input argument "name" the user can give the stored signal also a name.
715//
716// The default values are jref=0 and name="".
717
718 if (s!="B" && s!="b" && s!="J" && s!="j")
719 {
720 if (!jref && fXsig)
721 {
722 delete fXsig;
723 fXsig=0;
724 }
725 return;
726 }
727
728 if (mode!="M" && mode!="m" && mode!="T" && mode!="t")
729 {
730 if (!jref && fXsig)
731 {
732 delete fXsig;
733 fXsig=0;
734 }
735 return;
736 }
737
738 // Convert coordinates to fractional degrees.
739 a=ConvertAngle(a,"hms","deg");
740 d=ConvertAngle(d,"dms","deg");
741
742
743 AliTimestamp tx;
744 tx.SetEpoch(e,s);
745
746 Ali3Vector r;
747 Double_t vec[3]={1.,90.-d,a};
748 r.SetVector(vec,"sph","deg");
749
750 SetSignal(&r,"equ",mode,&tx,jref,name);
751}
752///////////////////////////////////////////////////////////////////////////
753void AliAstrolab::SetSignal(Double_t a,Double_t d,TString mode,AliTimestamp* ts,Int_t jref,TString name)
754{
755// Store a signal with right ascension (a) and declination (d) given for timestamp ts.
756// The position is stored in International Celestial Reference System (ICRS) coordinates.
757// The ICRS is a fixed, time independent frame and as such provides a unique reference
758// frame without the need of specifying any epoch etc...
759// The ICRS coordinate definitions match within 20 mas the mean ones of the J2000.0
760// equatorial system. Nevertheless, to obtain the highest accuracy, the slight
761// coordinate correction between J2000 and ICRS is performed here via the
762// so-called frame bias matrix.
763// For further details see the U.S. Naval Observatory (USNO) circular 179 (2005),
764// which is available on http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
765//
766// The coordinates (a,d) can represent so-called "mean" and "true" values.
767// The distinction between these two representations is the following :
768//
769// mean values : (a,d) are only corrected for precession and not for nutation
770// true values : (a,d) are corrected for both precession and nutation
771//
772// The input parameter "mode" allows the user to specifiy either "mean" or "true"
773// values for the input (a,d) coordinates.
774//
775// a : Right ascension in hhmmss.sss
776// d : Declination in dddmmss.sss
777// mode = "M" --> Input coordinates are the mean values
778// "T" --> Input coordinates are the true values
779//
780// The input parameter "jref" allows the user to store so-called "reference" signals.
781// These reference signals may be used to check space-time event coincidences with the
782// stored measurement (e.g. coincidence of the measurement with transient phenomena).
783//
784// jref = 0 --> Storage of the measurement
785// j --> Storage of a reference signal at the j-th position (j=1 is first)
786// < 0 --> Add a reference signal at the next available position
787//
788// Via the input argument "name" the user can give the stored signal also a name.
789//
790// The default values are jref=0 and name="".
791//
792// Note : In case ts=0 the current timestamp of the lab will be taken.
793
794 // Convert coordinates to fractional degrees.
795 a=ConvertAngle(a,"hms","deg");
796 d=ConvertAngle(d,"dms","deg");
797
798 Ali3Vector r;
799 Double_t vec[3]={1.,90.-d,a};
800 r.SetVector(vec,"sph","deg");
801
802 SetSignal(&r,"equ",mode,ts,jref,name);
803}
804///////////////////////////////////////////////////////////////////////////
805AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,Int_t jref)
806{
807// Provide the user specified coordinates of a signal at the specific timestamp ts.
808// The coordinates are returned via the vector argument "r".
809// In addition also a pointer to the stored signal object is provided.
810// In case no stored signal was available or one of the input arguments was
811// invalid, the returned pointer will be 0.
812//
813// The input parameter "frame" allows the user to specify the frame to which
814// the components of r refer. Available options are :
815//
816// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
817// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
818// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
819// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
820// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
821// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
822// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
823// components of r correspond to theta=zenith angle and phi=pi-azimuth.
824// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
825// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
826// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
827// components of r correspond to the usual theta and phi angles.
828//
829// In case the coordinates are the equatorial right ascension and declination (a,d),
830// they can represent so-called "mean" and "true" values.
831// The distinction between these two representations is the following :
832//
833// mean values : (a,d) are only corrected for precession and not for nutation
834// true values : (a,d) are corrected for both precession and nutation
835//
836// The input parameter "mode" allows the user to specifiy either "mean" or "true"
837// values for the input in case of equatorial (a,d) coordinates.
838//
839// mode = "M" --> Input coordinates are the mean values
840// "T" --> Input coordinates are the true values
841//
842// The input parameter "jref" allows the user to access so-called "reference" signals.
843// These reference signals may be used to check space-time event coincidences with the
844// stored measurement (e.g. coincidence of the measurement with transient phenomena).
845//
846// jref = 0 --> Access to the measurement
847// j --> Access to the reference signal at the j-th position (j=1 is first)
848//
849// Default value is jref=0.
850//
851// Note : In case ts=0 the current timestamp of the lab will be taken.
852
853 r.SetZero();
854
855 if (frame!="equ" && frame!="gal" && frame!="ecl" && frame!="hor" && frame!="icr" && frame!="loc") return 0;
856
857 if (frame=="equ" && mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
858
859 AliSignal* sx=GetSignal(jref);
860
861 if (!sx) return 0;
862
863 if (!ts) ts=(AliTimestamp*)this;
864
865 Double_t vec[3];
866 sx->GetPosition(vec,"sph","rad");
867 Ali3Vector q;
868 q.SetVector(vec,"sph","rad");
869
870 if (frame=="icr")
871 {
872 r.Load(q);
873 return sx;
874 }
875
876 // Convert from ICRS to equatorial J2000 coordinates
877 if (!fBias) SetBmatrix();
878 q=q.GetPrimed(&fB);
879
880 if (frame=="equ")
881 {
882 // Precess to specified timestamp
883 AliTimestamp ts1;
884 ts1.SetEpoch(2000,"J");
885 Precess(q,&ts1,ts);
886
887 // Nutation correction if requested
888 if (mode=="T" || mode=="t") Nutate(q,ts);
889 }
890
891 if (frame=="gal")
892 {
893 // Convert from equatorial J2000 to galactic
894 if (fGal != 2) SetGmatrix("J");
895 q=q.GetPrimed(&fG);
896 }
897
898 if (frame=="ecl")
899 {
900 // Precess to specified timestamp
901 AliTimestamp ts1;
902 ts1.SetEpoch(2000,"J");
903 Precess(q,&ts1,ts);
904
905 // Convert from equatorial to ecliptic coordinates
906 SetEmatrix(ts);
907 q=q.GetPrimed(&fE);
908 }
909
910 if (frame=="hor")
911 {
912 // Precess to specified timestamp
913 AliTimestamp ts1;
914 ts1.SetEpoch(2000,"J");
915 Precess(q,&ts1,ts);
916
917 // Nutation correction
918 Nutate(q,ts);
919
920 // Convert from equatorial to horizontal coordinates
921 SetHmatrix(ts);
922 q=q.GetPrimed(&fH);
923 }
924
925 if (frame=="loc")
926 {
927 // Get the signal in horizontal coordinates
928 GetSignal(q,"hor",mode,ts);
929
930 // Convert from horizontal local-frame coordinates
931 q=q.GetPrimed(&fL);
932 }
933
934 r.Load(q);
935 return sx;
936}
937///////////////////////////////////////////////////////////////////////////
938AliSignal* AliAstrolab::GetSignal(Ali3Vector& r,TString frame,TString mode,AliTimestamp* ts,TString name)
939{
940// Provide the user specified coordinates of the signal with the specified
941// name at the specific timestamp ts.
942// The coordinates are returned via the vector argument "r".
943// In addition also a pointer to the stored signal object is provided.
944// In case no such stored signal was available or one of the input arguments was
945// invalid, the returned pointer will be 0.
946//
947// The input parameter "frame" allows the user to specify the frame to which
948// the components of r refer. Available options are :
949//
950// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d),
951// where the "sph" components of r correspond to theta=(pi/2)-d and phi=a.
952// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
953// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
954// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b),
955// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
956// "hor" ==> Horizontal coordinates at the AliAstrolab location, where the "sph"
957// components of r correspond to theta=zenith angle and phi=pi-azimuth.
958// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b),
959// where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
960// "loc" ==> Local coordinates at the AliAstrolab location, where the "sph"
961// components of r correspond to the usual theta and phi angles.
962//
963// In case the coordinates are the equatorial right ascension and declination (a,d),
964// they can represent so-called "mean" and "true" values.
965// The distinction between these two representations is the following :
966//
967// mean values : (a,d) are only corrected for precession and not for nutation
968// true values : (a,d) are corrected for both precession and nutation
969//
970// The input parameter "mode" allows the user to specifiy either "mean" or "true"
971// values for the input in case of equatorial (a,d) coordinates.
972//
973// mode = "M" --> Input coordinates are the mean values
974// "T" --> Input coordinates are the true values
975//
976// Note : In case ts=0 the current timestamp of the lab will be taken.
977
978 AliSignal* sx=0;
979 Int_t j=GetSignalIndex(name);
980 if (j>=0) sx=GetSignal(r,frame,mode,ts,j);
981 return sx;
982}
983///////////////////////////////////////////////////////////////////////////
984AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,Int_t jref)
985{
986// Provide precession (and nutation) corrected right ascension (a) and
987// declination (d) of the stored signal object at the specified timestamp.
988// In addition also a pointer to the stored signal object is provided.
989// In case no stored signal was available or one of the input arguments was
990// invalid, the returned pointer will be 0.
991//
992// The coordinates (a,d) can represent so-called "mean" and "true" values.
993// The distinction between these two representations is the following :
994//
995// mean values : (a,d) are only corrected for precession and not for nutation
996// true values : (a,d) are corrected for both precession and nutation
997//
998// The input parameter "mode" allows the user to select either
999// "mean" or "true" values for (a,d).
1000//
1001// The correction methods used are the new IAU 2000 ones as described in the
1002// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1003// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1004//
1005// a : Right ascension in hhmmss.sss
1006// d : Declination in dddmmss.sss
1007// mode = "M" --> Output coordinates are the mean values
1008// "T" --> Output coordinates are the true values
1009// ts : Timestamp at which the corrected coordinate values are requested.
1010//
1011// The input parameter "jref" allows the user to access so-called "reference" signals.
1012// These reference signals may be used to check space-time event coincidences with the
1013// stored regular signal (e.g. coincidence of the stored signal with transient phenomena).
1014//
1015// jref = 0 --> Access to the measurement
1016// j --> Access to the reference signal at the j-th position (j=1 is first)
1017//
1018// Default value is jref=0.
1019//
1020// Note : In case ts=0 the current timestamp of the lab will be taken.
1021
1022 a=0;
1023 d=0;
1024
1025 Ali3Vector r;
1026 AliSignal* sx=GetSignal(r,"equ",mode,ts,jref);
1027
1028 if (!sx) return 0;
1029
1030 // Retrieve the requested (a,d) values
1031 Double_t vec[3];
1032 r.GetVector(vec,"sph","deg");
1033 d=90.-vec[1];
1034 a=vec[2];
1035
1036 while (a<-360.)
1037 {
1038 a+=360.;
1039 }
1040 while (a>360.)
1041 {
1042 a-=360.;
1043 }
1044 while (d<-90.)
1045 {
1046 d+=90.;
1047 }
1048 while (d>90.)
1049 {
1050 d-=90.;
1051 }
1052
1053 // Convert coordinates to appropriate format
1054 a=ConvertAngle(a,"deg","hms");
1055 d=ConvertAngle(d,"deg","dms");
1056
1057 return sx;
1058}
1059///////////////////////////////////////////////////////////////////////////
1060AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString mode,AliTimestamp* ts,TString name)
1061{
1062// Provide precession (and nutation) corrected right ascension (a) and
1063// declination (d) of the stored signal object with the specified name
1064// at the specific timestamp ts.
1065// In addition also a pointer to the stored signal object is provided.
1066// In case no stored signal was available or one of the input arguments was
1067// invalid, the returned pointer will be 0.
1068//
1069// The coordinates (a,d) can represent so-called "mean" and "true" values.
1070// The distinction between these two representations is the following :
1071//
1072// mean values : (a,d) are only corrected for precession and not for nutation
1073// true values : (a,d) are corrected for both precession and nutation
1074//
1075// The input parameter "mode" allows the user to select either
1076// "mean" or "true" values for (a,d).
1077//
1078// The correction methods used are the new IAU 2000 ones as described in the
1079// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1080// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1081//
1082// a : Right ascension in hhmmss.sss
1083// d : Declination in dddmmss.sss
1084// mode = "M" --> Output coordinates are the mean values
1085// "T" --> Output coordinates are the true values
1086// ts : Timestamp at which the corrected coordinate values are requested.
1087// name : Name of the requested signal object
1088//
1089// Note : In case ts=0 the current timestamp of the lab will be taken.
1090
1091 AliSignal* sx=0;
1092 Int_t j=GetSignalIndex(name);
1093 if (j>=0) sx=GetSignal(a,d,mode,ts,j);
1094 return sx;
1095}
1096///////////////////////////////////////////////////////////////////////////
1097AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,Int_t jref)
1098{
1099// Provide precession (and nutation) corrected right ascension (a) and
1100// declination (d) of the stored signal object at the specified epoch e.
1101// In addition also a pointer to the stored signal object is provided.
1102// In case no stored signal was available or one of the input arguments was
1103// invalid, the returned pointer will be 0.
1104//
1105// The coordinates (a,d) can represent so-called "mean" and "true" values.
1106// The distinction between these two representations is the following :
1107//
1108// mean values : (a,d) are only corrected for precession and not for nutation
1109// true values : (a,d) are corrected for both precession and nutation
1110//
1111// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1112// values for the input (a,d) coordinates.
1113//
1114// a : Right ascension in hhmmss.sss
1115// d : Declination in dddmmss.sss
1116// s = "B" --> Besselian reference epoch.
1117// "J" --> Julian reference epoch.
1118// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1119// mode = "M" --> Input coordinates are the mean values
1120// "T" --> Input coordinates are the true values
1121//
1122// The input parameter "jref" allows the user to access so-called "reference" signals.
1123// These reference signals may be used to check space-time event coincidences with the
1124// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1125//
1126// jref = 0 --> Access to the measurement
1127// j --> Access to the reference signal at the j-th position (j=1 is first)
1128//
1129// Default value is jref=0.
1130
1131 a=0;
1132 d=0;
1133
1134 if (s!="B" && s!="b" && s!="J" && s!="j") return 0;
1135
1136 if (mode!="M" && mode!="m" && mode!="T" && mode!="t") return 0;
1137
1138 // Convert coordinates to fractional degrees.
1139 a=ConvertAngle(a,"hms","deg");
1140 d=ConvertAngle(d,"dms","deg");
1141
1142
1143 AliTimestamp tx;
1144 tx.SetEpoch(e,s);
1145
1146 AliSignal* sx=GetSignal(a,d,mode,&tx,jref);
1147 return sx;
1148}
1149///////////////////////////////////////////////////////////////////////////
1150AliSignal* AliAstrolab::GetSignal(Double_t& a,Double_t& d,TString s,Double_t e,TString mode,TString name)
1151{
1152// Provide precession (and nutation) corrected right ascension (a) and
1153// declination (d) of the stored signal object with the specified name
1154// at the specific epoch e.
1155// In addition also a pointer to the stored signal object is provided.
1156// In case no stored signal was available or one of the input arguments was
1157// invalid, the returned pointer will be 0.
1158//
1159// The coordinates (a,d) can represent so-called "mean" and "true" values.
1160// The distinction between these two representations is the following :
1161//
1162// mean values : (a,d) are only corrected for precession and not for nutation
1163// true values : (a,d) are corrected for both precession and nutation
1164//
1165// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1166// values for the input (a,d) coordinates.
1167//
1168// a : Right ascension in hhmmss.sss
1169// d : Declination in dddmmss.sss
1170// s = "B" --> Besselian reference epoch.
1171// "J" --> Julian reference epoch.
1172// e : Reference epoch for the input coordinates (e.g. 1900, 1950, 2000,...)
1173// mode = "M" --> Input coordinates are the mean values
1174// "T" --> Input coordinates are the true values
1175// name : Name of the requested signal object
1176
1177 AliSignal* sx=0;
1178 Int_t j=GetSignalIndex(name);
1179 if (j>=0) sx=GetSignal(a,d,s,e,mode,j);
1180 return sx;
1181}
1182///////////////////////////////////////////////////////////////////////////
1183AliSignal* AliAstrolab::GetSignal(Int_t jref)
1184{
1185// Provide the pointer to a stored signal object.
1186//
1187// The input parameter "jref" allows the user to access so-called "reference" signals.
1188// These reference signals may be used to check space-time event coincidences with the
1189// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1190//
1191// jref = 0 --> Access to the measurement
1192// j --> Access to the reference signal at the j-th position (j=1 is first)
1193//
1194// Default value is jref=0.
1195
1196 AliSignal* sx=0;
1197 if (!jref)
1198 {
1199 sx=fXsig;
1200 }
1201 else
1202 {
1203 if (jref>0 && jref<fRefs->GetSize()) sx=(AliSignal*)fRefs->At(jref-1);
1204 }
1205 return sx;
1206}
1207///////////////////////////////////////////////////////////////////////////
1208AliSignal* AliAstrolab::GetSignal(TString name)
1209{
1210// Provide the pointer to the stored signal object with the specified name.
1211
1212 AliSignal* sx=0;
1213 Int_t j=GetSignalIndex(name);
1214 if (j>=0) sx=GetSignal(j);
1215 return sx;
1216}
1217///////////////////////////////////////////////////////////////////////////
1218void AliAstrolab::RemoveRefSignal(Int_t j,Int_t compress)
1219{
1220// Remove the reference signal which was stored at the j-th position (j=1 is first).
1221// Note : j=0 means that all stored reference signals will be removed.
1222// j<0 allows array compression (see below) without removing any signals.
1223//
1224// The "compress" parameter allows compression of the ref. signal storage array.
1225//
1226// compress = 1 --> Array will be compressed
1227// 0 --> Array will not be compressed
1228//
1229// Note : Compression of the storage array means that the indices of the
1230// reference signals in the storage array will change.
1231
1232 if (!fRefs) return;
1233
1234 // Clearing of the complete storage
1235 if (!j)
1236 {
1237 delete fRefs;
1238 fRefs=0;
1239 return;
1240 }
1241
1242 // Removing a specific reference signal
1243 if (j>0 && j<fRefs->GetSize())
1244 {
1245 TObject* obj=fRefs->RemoveAt(j-1);
1246 if (obj) delete obj;
1247 }
1248
1249 // Compression of the storage array
1250 if (compress) fRefs->Compress();
1251}
1252///////////////////////////////////////////////////////////////////////////
1253void AliAstrolab::RemoveRefSignal(TString name,Int_t compress)
1254{
1255// Remove the reference signal with the specified name.
1256//
1257// The "compress" parameter allows compression of the ref. signal storage array.
1258//
1259// compress = 1 --> Array will be compressed
1260// 0 --> Array will not be compressed
1261//
1262// Note : Compression of the storage array means that the indices of the
1263// reference signals in the storage array will change.
1264
1265 Int_t j=GetSignalIndex(name);
1266 if (j>0) RemoveRefSignal(j,compress);
1267}
1268///////////////////////////////////////////////////////////////////////////
1269Int_t AliAstrolab::GetSignalIndex(TString name)
1270{
1271// Provide storage index of the signal with the specified name.
1272// In case the name matches with the stored measurement,
1273// the value 0 is returned.
1274// In case no signal with the specified name was found, the value -1 is returned.
1275
1276 Int_t index=-1;
1277
1278 if (fXsig)
1279 {
1280 if (name==fXsig->GetName()) return 0;
1281 }
1282
1283 if (!fRefs) return -1;
1284
1285 for (Int_t i=0; i<fRefs->GetSize(); i++)
1286 {
1287 AliSignal* sx=(AliSignal*)fRefs->At(i);
1288 if (!sx) continue;
1289
1290 if (name==sx->GetName())
1291 {
1292 index=i+1;
1293 break;
1294 }
1295 }
1296
1297 return index;
1298}
1299///////////////////////////////////////////////////////////////////////////
1300void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,Int_t jref)
1301{
1302// Print data of a stored signal in user specified coordinates at the specific timestamp ts.
1303// In case no stored signal was available or one of the input arguments was
1304// invalid, no printout is produced.
1305//
1306// The argument "ndig" specifies the number of digits for the fractional
1307// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1308// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1309// will appear as 03.4 on the output.
1310// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1311//
1312// Note : The angle info is printed without additional spaces or "endline".
1313// This allows the print to be included in various composite output formats.
1314//
1315// The input parameter "frame" allows the user to specify the frame to which
1316// the coordinates refer. Available options are :
1317//
1318// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1319//
1320// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1321//
1322// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1323//
1324// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1325//
1326// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1327//
1328// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1329//
1330// In case the coordinates are the equatorial right ascension and declination (a,d),
1331// they can represent so-called "mean" and "true" values.
1332// The distinction between these two representations is the following :
1333//
1334// mean values : (a,d) are only corrected for precession and not for nutation
1335// true values : (a,d) are corrected for both precession and nutation
1336//
1337// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1338// values for the input in case of equatorial (a,d) coordinates.
1339//
1340// mode = "M" --> Input coordinates are the mean values
1341// "T" --> Input coordinates are the true values
1342//
661385a7 1343// The input parameter "mode" also determines which type of time and
1344// local hour angle will appear in the printout.
1345//
caa58e1a 1346// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1347// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
661385a7 1348//
8bde545d 1349// The input parameter "jref" allows printing of a so-called "reference" signal.
1350// These reference signals may serve to check space-time event coincidences with the
1351// stored measurement (e.g. coincidence of the measurement with transient phenomena).
1352//
1353// jref = 0 --> Printing of the measurement
1354// j --> Printing of the j-th reference signal
1355//
1356// Default value is jref=0.
1357//
1358// Note : In case ts=0 the current timestamp of the lab will be taken.
1359
1360 Ali3Vector r;
1361 AliSignal* sx=GetSignal(r,frame,mode,ts,jref);
1362
1363 if (!sx) return;
1364
661385a7 1365 // Local Hour Angle of the signal
1366 Double_t lha=GetHourAngle("A",ts,jref);
1367 TString slha="LAHA";
1368 if (mode=="M" || mode=="m")
1369 {
1370 lha=GetHourAngle("M",ts,jref);
1371 slha="LMHA";
1372 }
1373
8bde545d 1374 TString name=sx->GetName();
1375 if (name != "") cout << name.Data() << " ";
1376
1377 if (frame=="equ")
1378 {
1379 Double_t a,d;
1380 d=90.-r.GetX(2,"sph","deg");
1381 a=r.GetX(3,"sph","rad");
1382 cout << "Equatorial (" << mode.Data() <<") a : "; PrintAngle(a,"rad","hms",ndig);
1383 cout << " d : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1384 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1385 return;
1386 }
1387
1388 if (frame=="gal")
1389 {
1390 Double_t l,b;
1391 b=90.-r.GetX(2,"sph","deg");
1392 l=r.GetX(3,"sph","deg");
1393 cout << "Galactic l : "; PrintAngle(l,"deg","deg",ndig);
1394 cout << " b : "; PrintAngle(b,"deg","deg",ndig);
661385a7 1395 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1396 return;
1397 }
1398
1399 if (frame=="icr")
1400 {
1401 Double_t a,d;
1402 d=90.-r.GetX(2,"sph","deg");
1403 a=r.GetX(3,"sph","rad");
1404 cout << "ICRS l : "; PrintAngle(a,"rad","hms",ndig);
1405 cout << " b : "; PrintAngle(d,"deg","dms",ndig);
661385a7 1406 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1407 return;
1408 }
1409
1410 if (frame=="ecl")
1411 {
1412 Double_t a,d;
1413 d=90.-r.GetX(2,"sph","deg");
1414 a=r.GetX(3,"sph","deg");
1415 cout << "Ecliptic l : "; PrintAngle(a,"deg","deg",ndig);
1416 cout << " b : "; PrintAngle(d,"deg","deg",ndig);
661385a7 1417 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1418 return;
1419 }
1420
1421 if (frame=="hor")
1422 {
1423 Double_t alt=90.-r.GetX(2,"sph","deg");
1424 Double_t azi=180.-r.GetX(3,"sph","deg");
1425 while (azi>360)
1426 {
1427 azi-=360.;
1428 }
1429 while (azi<0)
1430 {
1431 azi+=360.;
1432 }
1433 cout << "Horizontal azi : "; PrintAngle(azi,"deg","deg",ndig);
1434 cout << " alt : "; PrintAngle(alt,"deg","deg",ndig);
661385a7 1435 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1436 return;
1437 }
1438
1439 if (frame=="loc")
1440 {
1441 Double_t theta=r.GetX(2,"sph","deg");
1442 Double_t phi=r.GetX(3,"sph","deg");
1443 cout << "Local-frame phi : "; PrintAngle(phi,"deg","deg",ndig);
1444 cout << " theta : "; PrintAngle(theta,"deg","deg",ndig);
661385a7 1445 cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
8bde545d 1446 return;
1447 }
1448}
1449///////////////////////////////////////////////////////////////////////////
1450void AliAstrolab::PrintSignal(TString frame,TString mode,AliTimestamp* ts,Int_t ndig,TString name)
1451{
1452// Print data of the stored signal with the specified name in user specified coordinates
1453// at the specific timestamp ts.
1454// In case such stored signal was available or one of the input arguments was
1455// invalid, no printout is produced.
1456//
1457// The argument "ndig" specifies the number of digits for the fractional
1458// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1459// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1460// will appear as 03.4 on the output.
1461// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1462//
1463// Note : The angle info is printed without additional spaces or "endline".
1464// This allows the print to be included in various composite output formats.
1465//
1466// The input parameter "frame" allows the user to specify the frame to which
1467// the coordinates refer. Available options are :
1468//
1469// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1470//
1471// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1472//
1473// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1474//
1475// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1476//
1477// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1478//
1479// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1480//
1481// In case the coordinates are the equatorial right ascension and declination (a,d),
1482// they can represent so-called "mean" and "true" values.
1483// The distinction between these two representations is the following :
1484//
1485// mean values : (a,d) are only corrected for precession and not for nutation
1486// true values : (a,d) are corrected for both precession and nutation
1487//
1488// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1489// values for the input in case of equatorial (a,d) coordinates.
1490//
1491// mode = "M" --> Input coordinates are the mean values
1492// "T" --> Input coordinates are the true values
1493//
661385a7 1494// The input parameter "mode" also determines which type of time and
1495// local hour angle will appear in the printout.
1496//
caa58e1a 1497// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1498// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
661385a7 1499//
8bde545d 1500// Note : In case ts=0 the current timestamp of the lab will be taken.
1501
1502 Int_t j=GetSignalIndex(name);
1503 if (j>=0) PrintSignal(frame,mode,ts,ndig,j);
1504}
1505///////////////////////////////////////////////////////////////////////////
1506void AliAstrolab::ListSignals(TString frame,TString mode,Int_t ndig)
1507{
1508// List all stored signals in user specified coordinates at the timestamp
1509// of the actual recording of the stored measurement under investigation.
1510// In case no (timestamp of the) actual measurement is available,
1511// the current timestamp of the lab will be taken.
1512// In case no stored signal is available or one of the input arguments is
1513// invalid, no printout is produced.
1514//
1515// The argument "ndig" specifies the number of digits for the fractional
1516// part (e.g. ndig=6 for "dms" corresponds to micro-arcsecond precision).
1517// No rounding will be performed, so an arcsecond count of 3.473 with ndig=1
1518// will appear as 03.4 on the output.
1519// Due to computer accuracy, precision on the pico-arcsecond level may get lost.
1520//
661385a7 1521// The default value is ndig=1.
8bde545d 1522//
1523// Note : The angle info is printed without additional spaces or "endline".
1524// This allows the print to be included in various composite output formats.
1525//
1526// The input parameter "frame" allows the user to specify the frame to which
1527// the coordinates refer. Available options are :
1528//
1529// frame = "equ" ==> Equatorial coordinates with right ascension (a) and declination (d).
1530//
1531// "gal" ==> Galactic coordinates with longitude (l) and lattitude (b).
1532//
1533// "ecl" ==> Ecliptic coordinates with longitude (l) and lattitude (b).
1534//
1535// "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1536//
1537// "icr" ==> ICRS coordinates with longitude (l) and lattitude (b).
1538//
1539// "loc" ==> Local spherical angles theta and phi at the AliAstrolab location.
1540//
1541// In case the coordinates are the equatorial right ascension and declination (a,d),
1542// they can represent so-called "mean" and "true" values.
1543// The distinction between these two representations is the following :
1544//
1545// mean values : (a,d) are only corrected for precession and not for nutation
1546// true values : (a,d) are corrected for both precession and nutation
1547//
1548// The input parameter "mode" allows the user to specifiy either "mean" or "true"
1549// values for the input in case of equatorial (a,d) coordinates.
1550//
1551// mode = "M" --> Input coordinates are the mean values
1552// "T" --> Input coordinates are the true values
1553//
661385a7 1554// The input parameter "mode" also determines which type of time and
1555// local hour angle will appear in the listing.
1556//
caa58e1a 1557// mode = "M" --> Mean Sidereal Time (MST) and Local Mean Hour Angle (LMHA)
1558// "T" --> Apparent Sidereal Time (AST) and Local Apparent Hour Angle (LAHA)
8bde545d 1559
1560 Int_t iprint=0;
1561
1562 AliTimestamp* tx=0;
1563
661385a7 1564 Int_t dform=1;
1565 if (mode=="T" || mode=="t") dform=-1;
1566
8bde545d 1567 if (fXsig)
1568 {
661385a7 1569 tx=fXsig->GetTimestamp();
8bde545d 1570 if (!tx) tx=(AliTimestamp*)this;
1571 cout << " *AliAstrolab::ListSignals* List of all stored signals." << endl;
661385a7 1572 cout << " === The measurement under investigation ===" << endl;
1573 cout << " Timestamp of the actual observation" << endl;
1574 tx->Date(dform,fToffset);
1575 cout << " Location of the actual observation" << endl;
1576 cout << " "; PrintSignal(frame,mode,tx,ndig); cout << endl;
8bde545d 1577 iprint=1;
1578 }
1579
1580 if (!fRefs) return;
1581
1582 for (Int_t i=1; i<=fRefs->GetSize(); i++)
1583 {
1584 AliSignal* sx=GetSignal(i);
1585 if (!sx) continue;
1586
1587 if (!iprint)
1588 {
1589 cout << " *AliAstrolab::ListRefSignals* List of all stored signals." << endl;
1590 tx=(AliTimestamp*)this;
661385a7 1591 cout << " Current timestamp of the laboratory" << endl;
1592 tx->Date(dform,fToffset);
8bde545d 1593 iprint=1;
1594 }
1595 if (iprint==1)
1596 {
661385a7 1597 cout << " === All stored reference signals according to the above timestamp ===" << endl;
8bde545d 1598 iprint=2;
1599 }
661385a7 1600 cout << " Index : " << i << " "; PrintSignal(frame,mode,tx,ndig,i); cout << endl;
8bde545d 1601 }
1602}
1603///////////////////////////////////////////////////////////////////////////
1604void AliAstrolab::Precess(Ali3Vector& r,AliTimestamp* ts1,AliTimestamp* ts2)
1605{
1606// Correct mean right ascension and declination given for Julian date "jd"
1607// for the earth's precession corresponding to the specified timestamp.
1608// The results are the so-called "mean" (i.e. precession corrected) values,
1609// corresponding to the specified timestamp.
1610// The method used is the new IAU 2000 one as described in the
1611// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1612// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1613// Since the standard reference epoch is J2000, this implies that all
1614// input (a,d) coordinates will be first internally converted to the
1615// corresponding J2000 values before the precession correction w.r.t. the
1616// specified lab timestamp will be applied.
1617//
1618// r : Input vector containing the right ascension and declination information
1619// in the form of standard Ali3Vector coordinates.
1620// In spherical coordinates the phi corresponds to the right ascension,
1621// whereas the declination corresponds to (pi/2)-theta.
1622// jd : Julian date corresponding to the input coordinate values.
1623// ts : Timestamp corresponding to the requested corrected coordinate values.
1624//
1625// Note : In case ts=0 the current timestamp of the lab will be taken.
1626
1627 // Convert back to J2000 values
1628 Ali3Vector r0;
1629 SetPmatrix(ts1);
1630 r0=r.GetUnprimed(&fP);
1631
1632 // Precess to the specified timestamp
1633 if (!ts2) ts2=(AliTimestamp*)this;
1634 SetPmatrix(ts2);
1635 r=r0.GetPrimed(&fP);
1636}
1637///////////////////////////////////////////////////////////////////////////
1638void AliAstrolab::Nutate(Ali3Vector& r,AliTimestamp* ts)
1639{
1640// Correct mean right ascension and declination for the earth's nutation
1641// corresponding to the specified timestamp.
1642// The results are the so-called "true" (i.e. nutation corrected) values,
1643// corresponding to the specified timestamp.
1644// The method used is the new IAU 2000 one as described in the
1645// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1646// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
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// ts : Timestamp for which the corrected coordinate values are requested.
1653//
1654// Note : In case ts=0 the current timestamp of the lab will be taken.
1655
1656 // Nutation correction for the specified timestamp
1657 if (!ts) ts=(AliTimestamp*)this;
1658 SetNmatrix(ts);
1659 r=r.GetPrimed(&fN);
1660}
1661///////////////////////////////////////////////////////////////////////////
1662void AliAstrolab::SetBmatrix()
1663{
1664// Set the frame bias matrix elements.
1665// The formulas and numerical constants used are the ones from the
1666// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1667// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1668
1669 Double_t pi=acos(-1.);
1670
1671 // Parameters in mas
1672 Double_t a=-14.6;
1673 Double_t x=-16.6170;
1674 Double_t e=-6.8192;
1675
1676 // Convert to radians
1677 a*=pi/(180.*3600.*1000.);
1678 x*=pi/(180.*3600.*1000.);
1679 e*=pi/(180.*3600.*1000.);
1680
1681 Double_t mat[9];
1682 mat[0]=1.-0.5*(a*a+x*x);
1683 mat[1]=a;
1684 mat[2]=-x;
1685 mat[3]=-a-e*x;
1686 mat[4]=1.-0.5*(a*a+e*e);
1687 mat[5]=-e;
1688 mat[6]=x-e*a;
1689 mat[7]=e+x*a;
1690 mat[8]=1.-0.5*(e*e+x*x);
1691
1692 fB.SetMatrix(mat);
1693 fBias=1;
1694}
1695///////////////////////////////////////////////////////////////////////////
1696void AliAstrolab::SetPmatrix(AliTimestamp* ts)
1697{
1698// Set precession matrix elements for Julian date jd w.r.t. J2000.
1699// The formulas and numerical constants used are the ones from the
1700// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1701// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1702// All numerical constants refer to the standard reference epoch J2000.
1703
1704 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1705 if (!ts)
1706 {
1707 fP.SetMatrix(mat);
1708 return;
1709 }
1710
1711 Double_t pi=acos(-1.);
1712
1713 Double_t t=(ts->GetJD()-2451545.0)/36525.; // Julian centuries since J2000.0
1714
1715 // Parameters for the precession matrix in arcseconds
1716 Double_t eps0=84381.406; // Mean ecliptic obliquity at J2000.0
1717 Double_t psi=5038.481507*t-1.0790069*pow(t,2)-0.00114045*pow(t,3)+0.000132851*pow(t,4)
1718 -0.0000000951*pow(t,4);
1719 Double_t om=eps0-0.025754*t+0.0512623*pow(t,2)-0.00772503*pow(t,3)-0.000000467*pow(t,4)
1720 +0.0000003337*pow(t,5);
1721 Double_t chi=10.556403*t-2.3814292*pow(t,2)-0.00121197*pow(t,3)+0.000170663*pow(t,4)
1722 -0.0000000560*pow(t,5);
1723
1724 // Convert to radians
1725 eps0*=pi/(180.*3600.);
1726 psi*=pi/(180.*3600.);
1727 om*=pi/(180.*3600.);
1728 chi*=pi/(180.*3600.);
1729
1730 Double_t s1=sin(eps0);
1731 Double_t s2=sin(-psi);
1732 Double_t s3=sin(-om);
1733 Double_t s4=sin(chi);
1734 Double_t c1=cos(eps0);
1735 Double_t c2=cos(-psi);
1736 Double_t c3=cos(-om);
1737 Double_t c4=cos(chi);
1738
1739 mat[0]=c4*c2-s2*s4*c3;
1740 mat[1]=c4*s2*c1+s4*c3*c2*c1-s1*s4*s3;
1741 mat[2]=c4*s2*s1+s4*c3*c2*s1+c1*s4*s3;
1742 mat[3]=-s4*c2-s2*c4*c3;
1743 mat[4]=-s4*s2*c1+c4*c3*c2*c1-s1*c4*s3;
1744 mat[5]=-s4*s2*s1+c4*c3*c2*s1+c1*c4*s3;
1745 mat[6]=s2*s3;
1746 mat[7]=-s3*c2*c1-s1*c3;
1747 mat[8]=-s3*c2*s1+c3*c1;
1748
1749 fP.SetMatrix(mat);
1750}
1751///////////////////////////////////////////////////////////////////////////
1752void AliAstrolab::SetNmatrix(AliTimestamp* ts)
1753{
1754// Set nutation matrix elements for the specified Julian date jd.
1755// The formulas and numerical constants used are the ones from the
1756// U.S. Naval Observatory (USNO) circular 179 (2005), which is available on
1757// http://aa.usno,navy.mil/publications/docs/Circular_179.pdf.
1758
1759 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
1760 if (!ts)
1761 {
1762 fN.SetMatrix(mat);
1763 return;
1764 }
1765
1766 Double_t pi=acos(-1.);
1767
1768 Double_t dpsi,deps,eps;
1769 ts->Almanac(&dpsi,&deps,&eps);
1770
1771 // Convert to radians
1772 dpsi*=pi/(180.*3600.);
1773 deps*=pi/(180.*3600.);
1774 eps*=pi/(180.*3600.);
1775
1776 Double_t s1=sin(eps);
1777 Double_t s2=sin(-dpsi);
1778 Double_t s3=sin(-(eps+deps));
1779 Double_t c1=cos(eps);
1780 Double_t c2=cos(-dpsi);
1781 Double_t c3=cos(-(eps+deps));
1782
1783 mat[0]=c2;
1784 mat[1]=s2*c1;
1785 mat[2]=s2*s1;
1786 mat[3]=-s2*c3;
1787 mat[4]=c3*c2*c1-s1*s3;
1788 mat[5]=c3*c2*s1+c1*s3;
1789 mat[6]=s2*s3;
1790 mat[7]=-s3*c2*c1-s1*c3;
1791 mat[8]=-s3*c2*s1+c3*c1;
1792
1793 fN.SetMatrix(mat);
1794}
1795///////////////////////////////////////////////////////////////////////////
1796void AliAstrolab::SetGmatrix(TString mode)
1797{
1798// Set the mean equatorial to galactic coordinate conversion matrix.
1799// The B1950 parameters were taken from section 22.3 of the book
1800// "An Introduction to Modern Astrophysics" by Carrol and Ostlie (1996).
1801// The J2000 parameters are obtained by precession of the B1950 values.
1802//
1803// Via the input argument "mode" the required epoch can be selected
1804// mode = "B" ==> B1950
1805// "J" ==> J2000
1806
1807 Ali3Vector x; // The Galactic x-axis in the equatorial frame
1808 Ali3Vector y; // The Galactic y-axis in the equatorial frame
1809 Ali3Vector z; // The Galactic z-axis in the equatorial frame
1810
1811 Double_t a,d;
1812 Double_t vec[3]={1,0,0};
1813
1814 fGal=1; // Set flag to indicate B1950 matrix values
1815
1816 // B1950 equatorial coordinates of the North Galactic Pole (NGP)
1817 a=124900.;
1818 d=272400.;
1819 a=ConvertAngle(a,"hms","deg");
1820 d=ConvertAngle(d,"dms","deg");
1821 vec[1]=90.-d;
1822 vec[2]=a;
1823 z.SetVector(vec,"sph","deg");
1824
1825 // B1950 equatorial coordinates of the Galactic l=b=0 point
1826 a=174224.;
1827 d=-285500.;
1828 a=ConvertAngle(a,"hms","deg");
1829 d=ConvertAngle(d,"dms","deg");
1830 vec[1]=90.-d;
1831 vec[2]=a;
1832 x.SetVector(vec,"sph","deg");
1833
1834 // Precess to the corresponding J2000 values if requested
1835 if (mode=="J")
1836 {
1837 fGal=2; // Set flag to indicate J2000 matrix values
1838 AliTimestamp t1;
1839 t1.SetEpoch(1950,"B");
1840 AliTimestamp t2;
1841 t2.SetEpoch(2000,"J");
1842 Precess(z,&t1,&t2);
1843 Precess(x,&t1,&t2);
1844 }
1845
1846 // The Galactic y-axis is determined for the right handed frame
1847 y=z.Cross(x);
1848
1849 fG.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1850 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1851 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1852}
1853///////////////////////////////////////////////////////////////////////////
1854void AliAstrolab::SetEmatrix(AliTimestamp* ts)
1855{
1856// Set the mean equatorial to ecliptic coordinate conversion matrix
1857// for the specified timestamp.
1858// A nice sketch and explanation of the two frames can be found
1859// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1860
1861 Double_t dpsi,deps,eps;
1862 ts->Almanac(&dpsi,&deps,&eps);
1863
1864 // Convert to degrees
1865 eps/=3600.;
1866
1867 // Positions of the ecliptic axes w.r.t. the equatorial ones
1868 // at the moment of the specified timestamp
1869 Double_t theta1=90; // Ecliptic x-axis
1870 Double_t phi1=0;
1871 Double_t theta2=90.-eps; //Ecliptic y-axis
1872 Double_t phi2=90;
1873 Double_t theta3=eps; // Ecliptic z-axis
1874 Double_t phi3=270;
1875
1876 fE.SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);
1877}
1878///////////////////////////////////////////////////////////////////////////
1879void AliAstrolab::SetHmatrix(AliTimestamp* ts)
1880{
1881// Set the mean equatorial to horizontal coordinate conversion matrix
1882// for the specified timestamp.
1883// A nice sketch and explanation of the two frames can be found
1884// in chapter 3 of the book "Astronomy Methods" by Hale Bradt (2004).
1885//
1886// Note : In order to simplify the calculations, we use here a
1887// right-handed horizontal frame.
1888
1889 Ali3Vector x; // The (South pointing) horizontal x-axis in the equatorial frame
1890 Ali3Vector y; // The (East pointing) horizontal y-axis in the equatorial frame
1891 Ali3Vector z; // The (Zenith pointing) horizontal z-axis in the equatorial frame
1892
1893 Double_t l,b;
1894 GetLabPosition(l,b,"deg");
1895
1896 Double_t a;
1897 Double_t vec[3]={1,0,0};
1898
1899 // Equatorial coordinates of the horizontal z-axis
1900 // at the moment of the specified timestamp
1901 a=ts->GetLAST(fToffset);
1902 a*=15.; // Convert fractional hours to degrees
1903 vec[1]=90.-b;
1904 vec[2]=a;
1905 z.SetVector(vec,"sph","deg");
1906
1907 // Equatorial coordinates of the horizontal x-axis
1908 // at the moment of the specified timestamp
1909 vec[1]=180.-b;
1910 vec[2]=a;
1911 x.SetVector(vec,"sph","deg");
1912
1913 // The horizontal y-axis is determined for the right handed frame
1914 y=z.Cross(x);
1915
1916 fH.SetAngles(x.GetX(2,"sph","deg"),x.GetX(3,"sph","deg"),
1917 y.GetX(2,"sph","deg"),y.GetX(3,"sph","deg"),
1918 z.GetX(2,"sph","deg"),z.GetX(3,"sph","deg"));
1919}
1920///////////////////////////////////////////////////////////////////////////
1921void AliAstrolab::SetLocalFrame(Double_t t1,Double_t p1,Double_t t2,Double_t p2,Double_t t3,Double_t p3)
1922{
1923// Specification of the orientations of the local-frame axes.
1924// The input arguments represent the angles (in degrees) of the local-frame axes
1925// w.r.t. a so called Master Reference Frame (MRF), with the same convention
1926// as the input arguments of TRrotMatix::SetAngles.
1927//
1928// The right handed Master Reference Frame is defined as follows :
1929// Z-axis : Points to the local Zenith
1930// X-axis : Makes an angle of 90 degrees with the Z-axis and points South
1931// Y-axis : Makes an angle of 90 degrees with the Z-axis and points East
1932//
1933// Once the user has specified the local reference frame, any observed event
1934// can be related to astronomical space-time locations via the SetSignal
1935// and GetSignal memberfunctions.
1936
1937 // Set the matrix for the conversion of our reference frame coordinates
1938 // into the local-frame ones.
1939
1940 fL.SetAngles(t1,p1,t2,p2,t3,p3);
1941}
1942///////////////////////////////////////////////////////////////////////////
1943Double_t AliAstrolab::ConvertAngle(Double_t a,TString in,TString out) const
1944{
1945// Conversion of various angular formats.
1946//
1947// The input argument "a" denotes the angle to be converted.
1948// The string arguments "in" and "out" specify the angular I/O formats.
1949//
1950// in = "rad" : input angle provided in radians
1951// "deg" : input angle provided in degrees
1952// "dms" : input angle provided in dddmmss.sss
1953// "hms" : input angle provided in hhmmss.sss
1954//
1955// out = "rad" : output angle provided in radians
1956// "deg" : output angle provided in degrees
1957// "dms" : output angle provided in dddmmss.sss
1958// "hms" : output angle provided in hhmmss.sss
1959
1960 if (in==out) return a;
1961
1962 // Convert input to its absolute value in (fractional) degrees.
1963 Double_t pi=acos(-1.);
1964 Double_t epsilon=1.e-12; // Accuracy in (arc)seconds
1965 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
1966 Double_t s=0;
1967
1968 Double_t b=fabs(a);
1969
1970 if (in=="rad") b*=180./pi;
1971
1972 if (in=="dms")
1973 {
1974 word=Int_t(b);
1975 ddd=word/10000;
1976 word=word%10000;
1977 mm=word/100;
1978 ss=word%100;
1979 s=b-Double_t(ddd*10000+mm*100+ss);
1980 b=Double_t(ddd)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.;
1981 }
1982
1983 if (in=="hms")
1984 {
1985 word=Int_t(b);
1986 hh=word/10000;
1987 word=word%10000;
1988 mm=word/100;
1989 ss=word%100;
1990 s=b-Double_t(hh*10000+mm*100+ss);
1991 b=15.*(Double_t(hh)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.);
1992 }
1993
1994 while (b>360)
1995 {
1996 b-=360.;
1997 }
1998
1999 if (out=="rad") b*=pi/180.;
2000
2001 if (out=="dms")
2002 {
2003 ddd=Int_t(b);
2004 b=b-Double_t(ddd);
2005 b*=60.;
2006 mm=Int_t(b);
2007 b=b-Double_t(mm);
2008 b*=60.;
2009 ss=Int_t(b);
2010 s=b-Double_t(ss);
2011 if (s>(1.-epsilon))
2012 {
2013 s=0.;
2014 ss++;
2015 }
2016 while (ss>=60)
2017 {
2018 ss-=60;
2019 mm++;
2020 }
2021 while (mm>=60)
2022 {
2023 mm-=60;
2024 ddd++;
2025 }
2026 while (ddd>=360)
2027 {
2028 ddd-=360;
2029 }
2030 b=Double_t(10000*ddd+100*mm+ss)+s;
2031 }
2032
2033 if (out=="hms")
2034 {
2035 b/=15.;
2036 hh=Int_t(b);
2037 b=b-Double_t(hh);
2038 b*=60.;
2039 mm=Int_t(b);
2040 b=b-Double_t(mm);
2041 b*=60.;
2042 ss=Int_t(b);
2043 s=b-Double_t(ss);
2044 if (s>(1.-epsilon))
2045 {
2046 s=0.;
2047 ss++;
2048 }
2049 while (ss>=60)
2050 {
2051 ss-=60;
2052 mm++;
2053 }
2054 while (mm>=60)
2055 {
2056 mm-=60;
2057 hh++;
2058 }
2059 while (hh>=24)
2060 {
2061 hh-=24;
2062 }
2063 b=Double_t(10000*hh+100*mm+ss)+s;
2064 }
2065
2066 if (a<0) b=-b;
2067
2068 return b;
2069}
2070///////////////////////////////////////////////////////////////////////////
2071Double_t AliAstrolab::GetHourAngle(TString mode,AliTimestamp* ts,Int_t jref)
2072{
2073// Provide the Local Hour Angle (in fractional degrees) of a stored signal
2074// object at the specified timestamp.
2075//
2076// The input parameter "mode" allows the user to select either the
2077// "mean" or "apparent" value for the returned Hour Angle.
2078//
2079// mode = "M" --> Output is the Mean Hour Angle
2080// "A" --> Output is the Apparent Hour Angle
2081// ts : Timestamp at which the hour angle is requested.
2082//
2083// The input parameter "jref" allows the user to specify so-called "reference" signals.
2084// These reference signals may be used to check space-time event coincidences with the
661385a7 2085// stored measurement (e.g. coincidence of the measurement with transient phenomena).
8bde545d 2086//
661385a7 2087// jref = 0 --> Use the stored measurement
8bde545d 2088// j --> Use the reference signal at the j-th position (j=1 is first)
2089//
2090// Default value is jref=0.
2091//
2092// Note : In case ts=0 the current timestamp of the lab will be taken.
2093
2094 if (!ts) ts=(AliTimestamp*)this;
2095
2096 // Get corrected right ascension and declination for the specified timestamp.
2097 Double_t a,d;
661385a7 2098 if (mode=="M" || mode=="m") GetSignal(a,d,"M",ts,jref);
2099 if (mode=="A" || mode=="a") GetSignal(a,d,"T",ts,jref);
8bde545d 2100
2101 // Convert coordinates to fractional degrees.
2102 a=ConvertAngle(a,"hms","deg");
2103 d=ConvertAngle(d,"dms","deg");
2104
2105 a/=15.; // Convert a to fractional hours
2106 Double_t ha=0;
2107 if (mode=="M" || mode=="m") ha=ts->GetLMST(fToffset)-a;
2108 if (mode=="A" || mode=="a") ha=ts->GetLAST(fToffset)-a;
2109 ha*=15.; // Convert to (fractional) degrees
2110
2111 return ha;
2112}
2113///////////////////////////////////////////////////////////////////////////
2114void 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)
2115{
2116// Set the AliTimestamp parameters corresponding to the LT date and time
2117// in the Gregorian calendar as specified by the input arguments.
2118//
2119// The input arguments represent the following :
2120// y : year in LT (e.g. 1952, 2003 etc...)
2121// m : month in LT (1=jan 2=feb etc...)
2122// d : day in LT (1-31)
2123// hh : elapsed hours in LT (0-23)
2124// mm : elapsed minutes in LT (0-59)
2125// ss : elapsed seconds in LT (0-59)
2126// ns : remaining fractional elapsed second of LT in nanosecond
2127// ps : remaining fractional elapsed nanosecond of LT in picosecond
2128//
2129// Note : ns=0 and ps=0 are the default values.
2130//
2131
2132 SetLT(fToffset,y,m,d,hh,mm,ss,ns,ps);
2133}
2134///////////////////////////////////////////////////////////////////////////
2135void AliAstrolab::SetLT(Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps)
2136{
2137// Set the AliTimestamp parameters corresponding to the specified elapsed
2138// timespan since the beginning of the new LT year.
2139//
2140// The LT year and elapsed time span is entered via the following input arguments :
2141//
2142// y : year in LT (e.g. 1952, 2003 etc...)
2143// d : elapsed number of days
2144// s : (remaining) elapsed number of seconds
2145// ns : (remaining) elapsed number of nanoseconds
2146// ps : (remaining) elapsed number of picoseconds
2147//
2148// The specified d, s, ns and ps values will be used in an additive
2149// way to determine the elapsed timespan.
2150// So, specification of d=1, s=100, ns=0, ps=0 will result in the
2151// same elapsed time span as d=0, s=24*3600+100, ns=0, ps=0.
2152// However, by making use of the latter the user should take care
2153// of possible integer overflow problems in the input arguments,
2154// which obviously will provide incorrect results.
2155//
2156// Note : ns=0 and ps=0 are the default values.
2157
2158 SetLT(fToffset,y,d,s,ns,ps);
2159}
2160///////////////////////////////////////////////////////////////////////////
2161Double_t AliAstrolab::GetDifference(Int_t j,TString au,Double_t& dt,TString tu,Int_t mode,Int_t* ia,Int_t* it)
2162{
2163// Provide space and time difference between the j-th reference signal
2164// (j=1 indicates first) and the stored measurement.
2165//
2166// The return value of this memberfunction provides the positional angular
2167// difference, whereas the output argument "dt" provides the time difference.
2168//
2169// The units of the angular difference can be specified via the the "au"
2170// input argument, where
2171//
2172// au = "rad" --> Angular difference in (fractional) radians
2173// "deg" --> Angular difference in (fractional) degrees
2174//
2175// The units of the time difference can be specified via the "tu" and "mode"
2176// input arguments. For details please refer to AliTimestamp::GetDifference().
2177// Also here mode=1 is the default value.
2178//
2179// For the time difference the reference signal is used as the standard.
2180// This means that in case of a positive time difference, the stored
2181// measurement occurred later than the reference signal.
2182//
2183// In case j=0, the stored measurement will be compared with each
2184// reference signal and the returned angular and time differences are
2185// the minimal differences which were encountered.
2186// In this case the user may obtain the indices of the two stored reference signals
2187// which had the minimal angular and minimal time difference via the output
2188// arguments "ia" and "it" as follows :
2189//
2190// ia = Index of the stored reference signal with minimial angular difference
2191// it = Index of the stored reference signal with minimial time difference
2192//
2193// In case these indices are the same, there obviously was 1 single reference signal
2194// which showed both the minimal angular and time difference.
2195//
2196// The default values are mode=1, ia=0 and it=0;
2197
2198 Double_t dang=999;
2199 dt=1.e20;
2200 if (ia) *ia=0;
2201 if (it) *it=0;
2202
2203 if (!fRefs) return dang;
2204
2205 Ali3Vector rx; // Position of the measurement
2206 Ali3Vector r0; // Position of the reference signal
2207
2208 AliSignal* sx=GetSignal(rx,"icr","M",0);
2209 if (!sx) return dang;
2210
2211 AliTimestamp* tx=sx->GetTimestamp();
2212 if (!tx) return dang;
2213
2214 // Space and time difference w.r.t. a specific reference signal
2215 if (j>0)
2216 {
2217 AliSignal* s0=GetSignal(r0,"icr","M",0,j);
2218 if (!s0) return dang;
2219 AliTimestamp* t0=s0->GetTimestamp();
2220
2221 if (!t0) return dang;
2222
2223 dang=r0.GetOpeningAngle(rx,au);
2224 dt=t0->GetDifference(tx,tu,mode);
2225 return dang;
2226 }
2227
2228 // Minimal space and time difference encountered over all reference signals
2229 Double_t dangmin=dang;
2230 Double_t dtmin=dt;
2231 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2232 {
2233 AliSignal* s0=GetSignal(r0,"icr","M",0,i);
2234 if (!s0) continue;
2235 AliTimestamp* t0=s0->GetTimestamp();
2236 if (!t0) continue;
2237 dang=r0.GetOpeningAngle(rx,au);
2238 dt=t0->GetDifference(tx,tu,mode);
2239 if (fabs(dang)<dangmin)
2240 {
2241 dangmin=fabs(dang);
2242 *ia=i;
2243 }
2244 if (fabs(dt)<dtmin)
2245 {
2246 dtmin=fabs(dt);
2247 *it=i;
2248 }
2249 }
2250
2251 dt=dtmin;
2252 return dangmin;
2253}
2254///////////////////////////////////////////////////////////////////////////
2255Double_t AliAstrolab::GetDifference(TString name,TString au,Double_t& dt,TString tu,Int_t mode)
2256{
2257// Provide space and time difference between the reference signal
2258// with the specified name and the stored measurement.
2259//
2260// The return value of this memberfunction provides the positional angular
2261// difference, whereas the output argument "dt" provides the time difference.
2262//
2263// The units of the angular difference can be specified via the the "au"
2264// input argument, where
2265//
2266// au = "rad" --> Angular difference in (fractional) radians
2267// "deg" --> Angular difference in (fractional) degrees
2268//
2269// The units of the time difference can be specified via the "tu" and "mode"
2270// input arguments. For details please refer to AliTimestamp::GetDifference().
2271// Also here mode=1 is the default value.
2272//
2273// For the time difference the reference signal is used as the standard.
2274// This means that in case of a positive time difference, the stored
2275// measurement occurred later than the reference signal.
2276
2277 Double_t dang=999;
2278 dt=1.e20;
2279
2280 Int_t j=GetSignalIndex(name);
2281 if (j>0) dang=GetDifference(j,au,dt,tu,mode);
2282 return dang;
2283}
2284///////////////////////////////////////////////////////////////////////////
2285TArrayI* AliAstrolab::MatchRefSignal(Double_t da,TString au,Double_t dt,TString tu,Int_t mode)
2286{
2287// Provide the storage indices of the reference signals which match in space
2288// and time with the stored measurement.
2289// The indices are returned via a pointer to a TArrayI object.
2290// In case no matches were found, the null pointer is returned.
2291// A reference signal is regarded as matching with the stored measurement
2292// if the positional angular difference doesn't exceed "da" and the absolute
2293// value of the time difference doesn't exceed "dt".
2294//
2295// The units of the angular difference "da" can be specified via the the "au"
2296// input argument, where
2297//
2298// au = "rad" --> Angular difference in (fractional) radians
2299// "deg" --> Angular difference in (fractional) degrees
2300//
2301// The units of the time difference "dt" can be specified via the "tu" and "mode"
2302// input arguments. For details please refer to AliTimestamp::GetDifference().
2303// Also here mode=1 is the default value.
2304
2305 if (!fXsig || !fRefs) return 0;
2306
2307 Int_t nrefs=fRefs->GetEntries();
2308
2309 if (fIndices) delete fIndices;
2310 fIndices=new TArrayI(nrefs);
2311
2312 Double_t dang,dtime;
2313 Int_t jfill=0;
2314 for (Int_t i=1; i<=fRefs->GetSize(); i++)
2315 {
2316 dang=GetDifference(i,au,dtime,tu,mode);
2317 if (fabs(dang)<=da && fabs(dtime)<=dt)
2318 {
2319 fIndices->AddAt(i,jfill);
2320 jfill++;
2321 }
2322 }
2323
2324 fIndices->Set(jfill);
2325 return fIndices;
2326}
2327///////////////////////////////////////////////////////////////////////////
2328TObject* AliAstrolab::Clone(const char* name) const
2329{
2330// Make a deep copy of the current object and provide the pointer to the copy.
2331// This memberfunction enables automatic creation of new objects of the
2332// correct type depending on the object type, a feature which may be very useful
2333// for containers when adding objects in case the container owns the objects.
2334
2335 AliAstrolab* lab=new AliAstrolab(*this);
2336 if (name)
2337 {
2338 if (strlen(name)) lab->SetName(name);
2339 }
2340 return lab;
2341}
2342///////////////////////////////////////////////////////////////////////////