Coding convetion violations fixed.
[u/mrichter/AliRoot.git] / RALICE / AliAstrolab.cxx
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 //
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 //
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 // Also graphical facilities (e.g. DisplaySignals) are available to
51 // provide skymaps in various projections.
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);
93 // lab.SetSignal(&r,"gal","M",0,-1,"GRB890730");
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 //
122 // // Display all stored objects in a skymap (Hammer projection)
123 // lab.DisplaySignals("equ","M",&ts,"ham",1);
124 //
125 //--- Author: Nick van Eijndhoven 15-mar-2007 Utrecht University
126 //- Modified: NvE $Date$ Utrecht University
127 ///////////////////////////////////////////////////////////////////////////
128
129 #include <cstdlib>
130 #include "AliAstrolab.h"
131 #include "Riostream.h"
132  
133 ClassImp(AliAstrolab) // Class implementation to enable ROOT I/O
134  
135 AliAstrolab::AliAstrolab(const char* name,const char* title) : TTask(name,title),AliTimestamp()
136 {
137 // Default constructor
138
139  fToffset=0;
140  fXsig=0;
141  fRefs=0;
142  fBias=0;
143  fGal=0;
144  fIndices=0;
145  fMeridian=-999;
146  fProj="none";
147  fCanvas=0;
148  fHist=0;
149  fMarkers=0;
150 }
151 ///////////////////////////////////////////////////////////////////////////
152 AliAstrolab::~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  }
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  }
186 }
187 ///////////////////////////////////////////////////////////////////////////
188 AliAstrolab::AliAstrolab(const AliAstrolab& t) : TTask(t),AliTimestamp(t)
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;
210  fMeridian=-999;
211  fProj="none";
212  fCanvas=0;
213  fHist=0;
214  fMarkers=0;
215 }
216 ///////////////////////////////////////////////////////////////////////////
217 void 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 ///////////////////////////////////////////////////////////////////////////
250 void 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 ///////////////////////////////////////////////////////////////////////////
266 void 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 ///////////////////////////////////////////////////////////////////////////
299 AliPosition 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 ///////////////////////////////////////////////////////////////////////////
309 void 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 ///////////////////////////////////////////////////////////////////////////
334 Double_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 ///////////////////////////////////////////////////////////////////////////
345 Double_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 ///////////////////////////////////////////////////////////////////////////
359 Double_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 ///////////////////////////////////////////////////////////////////////////
373 void 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 ///////////////////////////////////////////////////////////////////////////
489 void 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.
507 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
508 //                    where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
509 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
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.
513 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
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 ///////////////////////////////////////////////////////////////////////////
706 void 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 ///////////////////////////////////////////////////////////////////////////
783 void 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 ///////////////////////////////////////////////////////////////////////////
835 AliSignal* 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.
848 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
849 //                    where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
850 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
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.
854 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
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 ///////////////////////////////////////////////////////////////////////////
968 AliSignal* 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.
982 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
983 //                    where the "sph" components of r correspond to theta=(pi/2)-b and phi=l.
984 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b),
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.
988 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b),
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 ///////////////////////////////////////////////////////////////////////////
1014 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1090 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1127 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1180 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1213 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1238 AliSignal* 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 ///////////////////////////////////////////////////////////////////////////
1248 void 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 ///////////////////////////////////////////////////////////////////////////
1283 void 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 ///////////////////////////////////////////////////////////////////////////
1299 Int_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 ///////////////////////////////////////////////////////////////////////////
1330 void 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 //
1350 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
1351 //
1352 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
1353 //
1354 //          "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1355 //
1356 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
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 //
1373 // The input parameter "mode" also determines which type of time and
1374 // local hour angle will appear in the printout.
1375 //
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)
1378 //
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
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
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);
1414   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
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); 
1425   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
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);
1436   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
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);
1447   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
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);
1465   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
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);
1475   cout << " " << slha.Data() << " : "; PrintAngle(lha,"deg","hms",ndig);
1476   return;
1477  }
1478 }
1479 ///////////////////////////////////////////////////////////////////////////
1480 void 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.
1484 // In case no such stored signal was available or one of the input arguments was
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 //
1501 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
1502 //
1503 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
1504 //
1505 //          "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1506 //
1507 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
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 //
1524 // The input parameter "mode" also determines which type of time and
1525 // local hour angle will appear in the printout.
1526 //
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)
1529 //
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 ///////////////////////////////////////////////////////////////////////////
1536 void 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 //
1551 // The default value is ndig=1.
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 //
1561 //          "gal" ==> Galactic coordinates with longitude (l) and latitude (b).
1562 //
1563 //          "ecl" ==> Ecliptic coordinates with longitude (l) and latitude (b).
1564 //
1565 //          "hor" ==> Horizontal azimuth and altitude coordinates at the AliAstrolab location.
1566 //
1567 //          "icr" ==> ICRS coordinates with longitude (l) and latitude (b).
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 //
1584 // The input parameter "mode" also determines which type of time and
1585 // local hour angle will appear in the listing.
1586 //
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)
1589
1590  Int_t iprint=0;
1591
1592  AliTimestamp* tx=0;
1593
1594  Int_t dform=1;
1595  if (mode=="T" || mode=="t") dform=-1;
1596
1597  if (fXsig)
1598  {
1599   tx=fXsig->GetTimestamp();
1600   if (!tx) tx=(AliTimestamp*)this;
1601   cout << " *AliAstrolab::ListSignals* List of all stored signals." << endl;
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;
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;
1621    cout << " Current timestamp of the laboratory" << endl;
1622    tx->Date(dform,fToffset);
1623    iprint=1;
1624   }
1625   if (iprint==1)
1626   {
1627    cout << " === All stored reference signals according to the above timestamp ===" << endl;
1628    iprint=2;
1629   }
1630   cout << " Index : " << i << " "; PrintSignal(frame,mode,tx,ndig,i); cout << endl;
1631  }
1632 }
1633 ///////////////////////////////////////////////////////////////////////////
1634 void 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 ///////////////////////////////////////////////////////////////////////////
1668 void 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 ///////////////////////////////////////////////////////////////////////////
1692 void 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 ///////////////////////////////////////////////////////////////////////////
1726 void 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 ///////////////////////////////////////////////////////////////////////////
1782 void 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 ///////////////////////////////////////////////////////////////////////////
1826 void 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 ///////////////////////////////////////////////////////////////////////////
1884 void 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 ///////////////////////////////////////////////////////////////////////////
1909 void 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 ///////////////////////////////////////////////////////////////////////////
1951 void 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
1956 // as the input arguments of TRotMatrix::SetAngles.
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 ///////////////////////////////////////////////////////////////////////////
1973 Double_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
1984 //      "hrs" : input angle provided in fractional hours
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
1990 //       "hrs" : output angle provided in fractional hours
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
2004  if (in=="hrs") b*=15.;
2005
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
2035  if (out=="hrs") b/=15.;
2036
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 ///////////////////////////////////////////////////////////////////////////
2107 Double_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
2121 // stored measurement (e.g. coincidence of the measurement with transient phenomena).
2122 //
2123 // jref = 0 --> Use the stored measurement
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;
2134  if (mode=="M" || mode=="m") GetSignal(a,d,"M",ts,jref);
2135  if (mode=="A" || mode=="a") GetSignal(a,d,"T",ts,jref);
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 ///////////////////////////////////////////////////////////////////////////
2150 void 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 ///////////////////////////////////////////////////////////////////////////
2171 void 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 ///////////////////////////////////////////////////////////////////////////
2197 Double_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 ///////////////////////////////////////////////////////////////////////////
2291 Double_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 ///////////////////////////////////////////////////////////////////////////
2321 TArrayI* 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 ///////////////////////////////////////////////////////////////////////////
2364 void 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
2457  if (fMeridian<-pi)
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 ///////////////////////////////////////////////////////////////////////////
2614 void 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 ///////////////////////////////////////////////////////////////////////////
2678 void 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 ///////////////////////////////////////////////////////////////////////////
2757 void AliAstrolab::SetCentralMeridian(Double_t phi,TString u)
2758 {
2759 // Set the central meridian for the sky display.
2760 // Setting a value smaller than -pi will induce automatic meridian setting
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 ///////////////////////////////////////////////////////////////////////////
2777 void 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 ///////////////////////////////////////////////////////////////////////////
2836 void 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 ///////////////////////////////////////////////////////////////////////////
2847 void 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 ///////////////////////////////////////////////////////////////////////////
2859 void 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 ///////////////////////////////////////////////////////////////////////////
2876 void 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 ///////////////////////////////////////////////////////////////////////////
2890 TObject* 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 ///////////////////////////////////////////////////////////////////////////