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