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