source: trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc@ 4329

Last change on this file since 4329 was 3847, checked in by wittek, 22 years ago
*** empty log message ***
File size: 25.5 KB
Line 
1
2/* ======================================================================== *\
3!
4! *
5! * This file is part of MARS, the MAGIC Analysis and Reconstruction
6! * Software. It is distributed to you in the hope that it can be a useful
7! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
8! * It is distributed WITHOUT ANY WARRANTY.
9! *
10! * Permission to use, copy, modify and distribute this software and its
11! * documentation for any purpose is hereby granted without fee,
12! * provided that the above copyright notice appear in all copies and
13! * that both that copyright notice and this permission notice appear
14! * in supporting documentation. It is provided "as is" without express
15! * or implied warranty.
16! *
17!
18! Author(s): Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MSourcePosfromStarPos
28//
29// This is a task which
30// - calculates the position of the source in the camera
31// from the position of a known star in the camera
32// - and puts the source position into the container MSrcPosCam
33//
34// Input :
35// ASCII file containing for each run
36// - the run number
37// - the direction (theta, phi) the telescope is pointing to in [deg]
38// - the position (xStar, yStar) of a known star in the camera in [mm]
39// - the error (dxStar, dyStar) of this position in [mm]
40//
41// Output Containers :
42// MSrcPosCam
43//
44/////////////////////////////////////////////////////////////////////////////
45#include <TList.h>
46#include <TSystem.h>
47#include <TMatrixD.h>
48
49#include <fstream>
50
51#include "MSourcePosfromStarPos.h"
52
53#include "MParList.h"
54#include "MRawRunHeader.h"
55#include "MGeomCam.h"
56#include "MSrcPosCam.h"
57#include "MPointingPos.h"
58
59#include "MLog.h"
60#include "MLogManip.h"
61#include "MObservatory.h"
62
63ClassImp(MSourcePosfromStarPos);
64
65using namespace std;
66
67// --------------------------------------------------------------------------
68//
69// Default constructor.
70//
71MSourcePosfromStarPos::MSourcePosfromStarPos(
72 const char *name, const char *title)
73 : fIn(NULL)
74{
75 fName = name ? name : "MSourcePosfromStarPos";
76 fTitle = title ? title : "Calculate source position from star position";
77
78 fFileNames = new TList;
79 fFileNames->SetOwner();
80
81 fRuns = 0;
82 fSize = 100;
83
84 fRunNr.Set(fSize);
85
86 fThetaTel.Set(fSize);
87 fPhiTel.Set(fSize);
88 fdThetaTel.Set(fSize);
89 fdPhiTel.Set(fSize);
90
91 Int_t fRows = 1;
92 fDecStar.Set(fRows);
93 fRaStar.Set(fRows);
94 fxStar.ResizeTo(fRows,fSize);
95 fyStar.ResizeTo(fRows,fSize);
96 fdxStar.ResizeTo(fRows,fSize);
97 fdyStar.ResizeTo(fRows,fSize);
98
99 fStars = 0;
100 fDecSource = 0.0;
101 fRaSource = 0.0;
102
103 // these are the default values when starting the excution;
104 // later these locations contain the values of the last event
105 fxSourceold = 25.0;
106 fySourceold = -40.0;
107 fThetaradold = 25.0/kRad2Deg;
108 fPhiradold = 180.0/kRad2Deg;
109}
110
111// --------------------------------------------------------------------------
112//
113// Delete the filename list and the input stream if one exists.
114//
115MSourcePosfromStarPos::~MSourcePosfromStarPos()
116{
117 delete fFileNames;
118 if (fIn)
119 delete fIn;
120}
121
122
123// --------------------------------------------------------------------------
124//
125// Set the sky coordinates of the source and of the star
126//
127// Input :
128// declination in units of (Deg, Min, Sec)
129// right ascension in units of (Hour, Min, Sec)
130//
131
132void MSourcePosfromStarPos::SetSourceAndStarPosition(
133 TString nameSource,
134 Double_t decSourceDeg, Double_t decSourceMin, Double_t decSourceSec,
135 Double_t raSourceHour, Double_t raSourceMin, Double_t raSourceSec,
136 TString nameStar,
137 Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec,
138 Double_t raStarHour, Double_t raStarMin, Double_t raStarSec )
139{
140 *fLog << "MSourcePosfromStarPos::SetSourceAndStarPosition :" << endl;
141 *fLog << " Source (dec) : " << nameSource << " " << decSourceDeg << ":"
142 << decSourceMin << ":" << decSourceSec << endl;
143 *fLog << " Source (ra) : " << nameSource << " " << raSourceHour << ":"
144 << raSourceMin << ":" << raSourceSec << endl;
145
146 *fLog << " Star (dec) : " << nameStar << " " << decStarDeg << ":"
147 << decStarMin << ":" << decStarSec << endl;
148 *fLog << " Star (ra) : " << nameStar << " " << raStarHour << ":"
149 << raStarMin << ":" << raStarSec << endl;
150
151 // convert into radians
152 fDecSource = (decSourceDeg + decSourceMin/60.0 + decSourceSec/3600.0)
153 / kRad2Deg;
154 fRaSource = (raSourceHour + raSourceMin/60.0 + raSourceSec/3600.0)
155 * 360.0 / (24.0 * kRad2Deg);
156
157 fStars += 1;
158 fDecStar.Set(fStars);
159 fRaStar.Set(fStars);
160 fxStar.ResizeTo(fStars,fSize);
161 fyStar.ResizeTo(fStars,fSize);
162 fdxStar.ResizeTo(fStars,fSize);
163 fdyStar.ResizeTo(fStars,fSize);
164
165 fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
166 / kRad2Deg;
167 fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0)
168 * 360.0 / (24.0 * kRad2Deg);
169
170 *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fStars = "
171 << fStars << endl;
172 *fLog << all << " fDecSource, fRaSource, fDecStar, fRaStar were set to : [radians] "
173 << fDecSource << ", " << fRaSource << ", "
174 << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl;
175}
176
177// --------------------------------------------------------------------------
178//
179// Set the sky coordinates of another star
180//
181// Input :
182// declination in units of (Deg, Min, Sec)
183// right ascension in units of (Hour, Min, Sec)
184//
185
186void MSourcePosfromStarPos::AddStar(
187 TString nameStar,
188 Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec,
189 Double_t raStarHour, Double_t raStarMin, Double_t raStarSec )
190{
191 *fLog << "MSourcePosfromStarPos::AddStar :" << endl;
192 *fLog << " Star (dec) : " << nameStar << " " << decStarDeg << ":"
193 << decStarMin << ":" << decStarSec << endl;
194 *fLog << " Star (ra) : " << nameStar << " " << raStarHour << ":"
195 << raStarMin << ":" << raStarSec << endl;
196
197 // convert into radians
198 fStars += 1;
199 fDecStar.Set(fStars);
200 fRaStar.Set(fStars);
201 fxStar.ResizeTo(fStars,fSize);
202 fyStar.ResizeTo(fStars,fSize);
203 fdxStar.ResizeTo(fStars,fSize);
204 fdyStar.ResizeTo(fStars,fSize);
205
206 fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
207 / kRad2Deg;
208 fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0)
209 * 360.0 / (24.0 * kRad2Deg);
210
211 *fLog << all << "MSourcePosfromStarPos::AddStar; fStars = " << fStars
212 << endl;
213 *fLog << all << " fDecStar, fRaStar were set to : [radians] "
214 << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl;
215}
216
217// --------------------------------------------------------------------------
218//
219//
220Int_t MSourcePosfromStarPos::PreProcess(MParList *pList)
221{
222 MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
223 if (!geom)
224 {
225 *fLog << err << "MSourcePosfromStarPos : MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
226 return kFALSE;
227 }
228 fMm2Deg = geom->GetConvMm2Deg();
229 // fDistCameraReflector is the distance of the camera from the reflector center in [mm]
230 fDistCameraReflector = kRad2Deg / fMm2Deg;
231 *fLog << all << "MSourcePosfromStarPos::PreProcess; fMm2Deg, fDistCameraReflector = "
232 << fMm2Deg << ", " << fDistCameraReflector << endl;
233
234 fObservatory = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory"));
235 if (!fObservatory)
236 {
237 *fLog << err << "MObservatory not found... aborting" << endl;
238 return kFALSE;
239 }
240
241 fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
242 if (!fRun)
243 {
244 *fLog << err << "MSourcePosfromStarPos::PreProcess; MRawRunHeader not found... aborting." << endl;
245 return kFALSE;
246 }
247
248
249
250 fPointPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
251 if (!fPointPos)
252 {
253 *fLog << err << "MSourcePosfromStarPos::PreProcess; MPointingPos not found... aborting." << endl;
254 return kFALSE;
255 }
256
257
258 fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber("MSrcPosCam"));
259 if (!fSrcPos)
260 {
261 *fLog << err << "MSourcePosfromStarPos::PreProcess; MSrcPosCam not found... aborting" << endl;
262 return kFALSE;
263 }
264
265 //---------------------------------------------------------------------
266 // read all files and call ReadData() to read and store the information
267 //
268
269 *fLog << all << "---------------------------------" << endl;
270 while(1)
271 {
272 if (!OpenNextFile())
273 {
274 *fLog << "there is no more file to open" << endl;
275 break;
276 }
277
278 *fLog << "read data" << endl;
279
280 // read "fStarsRead" = no.of (x,y) pairs to be read
281 *fIn >> fStarsRead;
282
283 while (1)
284 {
285 if (fIn->eof())
286 {
287 *fLog << "eof encountered; open next file" << endl;
288
289 if (!OpenNextFile()) break;
290 }
291
292 // FIXME! Set InputStreamID
293
294 ReadData();
295 }
296 }
297
298 *fLog << "all data were read" << endl;
299 FixSize();
300
301 if (fDecSource == 0.0 || fRaSource == 0.0 || fStars == 0)
302 {
303 *fLog << warn << "MSourcePosfromStarPos::PreProcess; there are no sky coordinates defined for the source or from stars; fStars, fStarsRead = "
304 << fStars << ", " << fStarsRead
305 << endl;
306 }
307 *fLog << all << "---------------------------------" << endl;
308
309 //-------------------------------------------------------------
310
311 return kTRUE;
312}
313
314//=========================================================================
315//
316// SourcefromStar
317//
318// this routine calculates the position of a source (for example Crab) in the camera
319// from the position of a star (for example ZetaTauri) in the camera. The latter
320// position may have been determined by analysing the DC currents in the different
321// pixels.
322//
323// Input : thetaTel, phiTel the direction the telescope is pointing to,
324// in local coordinates
325// f the distance between camera and reflector
326// decStar, raStar the position of the star in sky coordinates
327// decSource, raSource the position of the source in sky coordinates
328// xStar, yStar the position of the star in the camera
329// dxStar, dyStar error of the position of the star in the camera
330//
331// Output : xSource, ySource the calculated position of the source in the camera
332// dxSource, dySource error of the calculated position of the source in
333// the camera
334//
335// Useful formulas can be found in TDAS 00-11 and TDAS 01-05
336//
337
338void MSourcePosfromStarPos::SourcefromStar(Double_t &f,
339 TArrayD &decStar, TArrayD &raStar,
340 Double_t &decSource, Double_t &raSource,
341 Double_t &thetaTel, Double_t &phiTel,
342 TArrayD &xStar, TArrayD &yStar,
343 TArrayD &dxStar, TArrayD &dyStar,
344 Double_t &xSource, Double_t &ySource,
345 Double_t &dxSource, Double_t &dySource)
346{
347 /*
348 *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in degrees" << endl;
349 *fLog << " decStar, raStar = " << decStar[0]*kRad2Deg << ", "
350 << raStar[0]*kRad2Deg << endl;
351 *fLog << " decSource, raSource = " << decSource*kRad2Deg << ", "
352 << raSource*kRad2Deg << endl;
353 *fLog << " thetaTel, phiTel = " << thetaTel*kRad2Deg << ", "
354 << phiTel*kRad2Deg << endl;
355 *fLog << " xStar, yStar = " << xStar[0]*fMm2Deg << ", "
356 << yStar[0]*fMm2Deg << endl;
357
358 *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in radians and mm" << endl;
359 *fLog << " decStar, raStar = " << decStar[0] << ", "
360 << raStar[0] << endl;
361 *fLog << " decSource, raSource = " << decSource << ", "
362 << raSource << endl;
363 *fLog << " thetaTel, phiTel = " << thetaTel << ", "
364 << phiTel << endl;
365 *fLog << " xStar, yStar = " << xStar[0] << ", "
366 << yStar[0] << endl;
367 */
368
369 // the units are assumed to be radians for theta, phi, dec and ra
370 // and mm for f, x and y
371
372
373 // calculate rotation angle alpha of sky image in camera
374 // (see TDAS 00-11, eqs. (18) and (20))
375 // a1 = cos(Lat), a3 = -sin(Lat), where Lat is the geographical latitude of La Palma
376 Double_t a1 = 0.876627;
377 Double_t a3 = -0.481171;
378
379 Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
380 ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
381 ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) );
382 Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
383 Double_t sinal = a1 * sin(phiTel) * denom;
384
385 *fLog << "old thetaTel, phiTel, cosal, sinal = " << thetaTel << ", "
386 << phiTel << ", " << cosal << ", " << sinal << endl;
387
388
389 fObservatory->RotationAngle(thetaTel, phiTel, sinal, cosal);
390
391 *fLog << "new thetaTel, phiTel, cosal, sinal = " << thetaTel << ", "
392 << phiTel << ", " << cosal << ", " << sinal << endl;
393
394
395 // calculate coordinates of source in system B (see TDAS 00-11, eqs. (2))
396 // note that right ascension and hour angle go into opposite directions
397 Double_t xB0 = cos(decSource) * cos(-raSource);
398 Double_t yB0 = cos(decSource) * sin(-raSource);
399 Double_t zB0 = -sin(decSource);
400
401 //*fLog << "xB0, yB0, zB0 = " << xB0 << ", " << yB0 << ", "
402 // << zB0 << endl;
403
404 //-----------------------------------------------------
405 // loop over stars
406 Double_t sumx = 0.0;
407 Double_t sumy = 0.0;
408 Double_t sumwx = 0.0;
409 Double_t sumwy = 0.0;
410
411 for (Int_t i=0; i<decStar.GetSize(); i++)
412 {
413 // calculate weights
414 Double_t weightx = 1.0 / (dxStar[i]*dxStar[i]);
415 Double_t weighty = 1.0 / (dyStar[i]*dyStar[i]);
416 sumwx += weightx;
417 sumwy += weighty;
418
419 //*fLog << "weightx, weighty = " << weightx << ", " << weighty << endl;
420
421 // calculate coordinates of star in system B (see TDAS 00-11, eqs. (2))
422 // note that right ascension and hour angle go into opposite directions
423 Double_t xB = cos(decStar[i]) * cos(-raStar[i]);
424 Double_t yB = cos(decStar[i]) * sin(-raStar[i]);
425 Double_t zB = -sin(decStar[i]);
426
427
428 //*fLog << "xB, yB, zB = " << xB << ", " << yB << ", "
429 // << zB << endl;
430
431
432 // calculate coordinates of star in a system with the basis vectors e1, e2, e3
433 // where e1 is in the direction (r0 x a)
434 // e2 is in the direction (e1 x r0)
435 // and e3 is in the direction -r0;
436 // r0 is the direction to the source
437 // and a is the earth rotation axis (pointing to the celestial north pole)
438 //
439 Double_t x = (-xB*yB0 + xB0*yB) / sqrt( xB0*xB0 + yB0*yB0 );
440 Double_t y = ( xB*xB0*zB0 + yB*yB0*zB0 - zB*(xB0*xB0 + yB0*yB0) )
441 / sqrt( xB0*xB0 + yB0*yB0 );
442 Double_t z = -(xB*xB0 + yB*yB0 + zB*zB0);
443
444 //*fLog << "x, y, z = " << x << ", " << y << ", "
445 // << z << endl;
446
447
448 // calculate coordinates of star in camera
449 Double_t xtilde = -f/z * (cosal*x - sinal*y);
450 Double_t ytilde = -f/z * (sinal*x + cosal*y);
451
452 //*fLog << "i, xtilde, ytile = " << i << " : " << xtilde << ", "
453 // << ytilde << endl;
454
455
456 // calculate coordinates of source in camera
457 // note : in real camera signs are inverted (therefore s = -1.0)
458 Double_t s = -1.0;
459
460 Double_t xs = xStar[i] - s * xtilde;
461 Double_t ys = yStar[i] - s * ytilde;
462
463 *fLog << "i, xs, ys = " << i << " : " << xs << ", "
464 << ys << endl;
465
466 sumx += xs * weightx;
467 sumy += ys * weighty;
468 }
469 //-----------------------------------------------------
470
471 xSource = sumx / sumwx;
472 ySource = sumy / sumwy;
473 dxSource = 1.0 / sqrt(sumwx);
474 dySource = 1.0 / sqrt(sumwy);
475
476 /*
477 Int_t run = fRun->GetRunNumber();
478 *fLog << all << "MSourcePosfromStarPos::SourcefromStar; run, xSource, ySource = "
479 << run << " : "
480 << xSource << " +- " << dxSource << ", "
481 << ySource << " +- " << dySource << endl;
482 */
483}
484
485// --------------------------------------------------------------------------
486//
487// Get the source position and put it into MSrcPosCam
488//
489//
490Bool_t MSourcePosfromStarPos::ReInit(MParList *pList)
491{
492 //if (1 == 1) return kTRUE;
493
494
495 Int_t run = fRun->GetRunNumber();
496 *fLog << all << "MSourcePosfromStarPos::ReInit; run = " << run << endl;
497
498
499 //-------------------------------------------------------------------
500 // search this run in the list
501 for (Int_t i=0; i<fSize; i++)
502 {
503 if (run == fRunNr[i])
504 {
505 //-----------------------------------------
506 // put the zenith angle into MPointingPos
507
508 Double_t thetarad = fThetaTel[i];
509 Double_t phirad = fPhiTel[i];
510
511 if (fabs(thetarad*kRad2Deg+1.0) < 0.001)
512 thetarad = fThetaradold;
513 else
514 fThetaradold = thetarad;
515 if (fabs(phirad*kRad2Deg+1.0) < 0.001)
516 phirad = fPhiradold;
517 else
518 fPhiradold = phirad;
519
520 fPointPos->SetLocalPosition(thetarad*kRad2Deg, phirad*kRad2Deg);
521 fPointPos->SetReadyToSave();
522
523 *fLog << all << "theta, phi = " << thetarad*kRad2Deg << ", "
524 << phirad*kRad2Deg << " deg" << endl;
525
526 //-----------------------------------------
527 // Get source position and put it into MSrcPosCam
528
529
530 if (fStars > 0 && fStars == fStarsRead)
531 {
532 TArrayD xStar(fxStar.GetNrows());
533 TArrayD dxStar(fdxStar.GetNrows());
534 TArrayD yStar(fyStar.GetNrows());
535 TArrayD dyStar(fdyStar.GetNrows());
536 for (Int_t j=0; j<fxStar.GetNrows(); j++)
537 {
538 xStar[j] = fxStar(j, i);
539 dxStar[j] = fdxStar(j, i);
540 yStar[j] = fyStar(j, i);
541 dyStar[j] = fdyStar(j, i);
542 }
543
544 SourcefromStar( fDistCameraReflector,
545 fDecStar, fRaStar, fDecSource, fRaSource,
546 thetarad, phirad,
547 xStar, yStar,
548 dxStar, dyStar,
549 fxSource, fySource,
550 fdxSource, fdySource);
551
552 fSrcPos->SetXY(fxSource, fySource);
553
554 fxSourceold = fxSource;
555 fySourceold = fySource;
556
557 *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = "
558 << fRunNr[i] << ", " << fxSource << " +- " << fdxSource
559 << ", " << fySource << " +- " << fdySource << endl;
560
561 fSrcPos->SetReadyToSave();
562 }
563 else
564 {
565 // set default values
566 fxSource = fxSourceold;
567 fySource = fySourceold;
568 fSrcPos->SetXY(fxSource, fySource);
569 fSrcPos->SetReadyToSave();
570
571
572 *fLog << warn << "MSourcePosfromStarPos::ReInit; no information on source position for run number = "
573 << run << endl;
574 *fLog << warn << " set xSource, ySource = " << fxSource << ", "
575 << fySource << " mm" << endl;
576 }
577
578
579 return kTRUE;
580 }
581 }
582 //-------------------------------------------------------------------
583
584 // set default values
585 fxSource = fxSourceold;
586 fySource = fySourceold;
587 fSrcPos->SetXY(fxSource, fySource);
588 fSrcPos->SetReadyToSave();
589
590 Double_t thetarad = fThetaradold;
591 Double_t phirad = fPhiradold;
592 fPointPos->SetLocalPosition(thetarad*kRad2Deg, phirad*kRad2Deg);
593 fPointPos->SetReadyToSave();
594
595
596 *fLog << warn << "MSourcePosfromStarPos::ReInit; no information on theta, phi and source position for run number = "
597 << run << endl;
598 *fLog << warn << " set xSource, ySource = " << fxSource << ", "
599 << fySource << " mm" << endl;
600 *fLog << warn << " set theta, phi = " << thetarad*kRad2Deg << ", "
601 << phirad*kRad2Deg << " deg" << endl;
602
603
604 return kTRUE;
605}
606
607// --------------------------------------------------------------------------
608//
609//
610Int_t MSourcePosfromStarPos::Process()
611{
612 //Int_t run = fRun->GetRunNumber();
613 //*fLog << "MSourcePosfromStarPos::Process; run = " << run << endl;
614
615
616 return kTRUE;
617}
618
619// --------------------------------------------------------------------------
620//
621//
622Int_t MSourcePosfromStarPos::PostProcess()
623{
624
625 return kTRUE;
626}
627
628// --------------------------------------------------------------------------
629//
630// read the data from the ASCII file and store them
631//
632void MSourcePosfromStarPos::FixSize()
633{
634 fSize = fRuns;
635 if (fRuns <= 0)
636 fSize = 1;
637
638 fRunNr.Set(fSize);
639
640 fThetaTel.Set(fSize);
641 fPhiTel.Set(fSize);
642 fdThetaTel.Set(fSize);
643 fdPhiTel.Set(fSize);
644
645 Int_t fRows = fxStar.GetNrows();
646 fxStar.ResizeTo(fRows, fSize);
647 fyStar.ResizeTo(fRows, fSize);
648 fdxStar.ResizeTo(fRows, fSize);
649 fdyStar.ResizeTo(fRows, fSize);
650
651 *fLog << "MSourcePosfromStarPos::FixSize; fix size of arrays : fStars = "
652 << fStars << ", fRuns = " << fRuns << ", fRows = " << fRows
653 << ", fSize = " << fSize << endl;
654 *fLog << " first run : " << fRunNr[0] << ", last run : "
655 << fRunNr[fRuns-1] << endl;
656}
657
658// --------------------------------------------------------------------------
659//
660// read the data from the ASCII file and store them
661//
662void MSourcePosfromStarPos::ReadData()
663{
664 Float_t val;
665 Int_t ival;
666
667 // extend size of arrays if necessary
668 if ( fRuns >= fSize )
669 {
670 fSize += 100;
671
672 fRunNr.Set(fSize);
673
674 fThetaTel.Set(fSize);
675 fPhiTel.Set(fSize);
676 fdThetaTel.Set(fSize);
677 fdPhiTel.Set(fSize);
678
679 Int_t fRows = fxStar.GetNrows();
680 fxStar.ResizeTo(fRows, fSize);
681 fyStar.ResizeTo(fRows, fSize);
682 fdxStar.ResizeTo(fRows, fSize);
683 fdyStar.ResizeTo(fRows, fSize);
684
685 *fLog << "MSourcePosfromStarPos::ReadData(); size of arrays has been increased to (fRows, fSize) = ("
686 << fRows << ", " << fSize << ")" << endl;
687 }
688
689 //-------------------
690 // read header line
691 //*fIn >> val;
692
693 //*fLog << "val =" << val << endl;
694
695 //*fIn >> val;
696 //*fIn >> val;
697 //*fIn >> val;
698 //*fIn >> val;
699 //*fIn >> val;
700 //*fIn >> val;
701 //*fIn >> val;
702 //*fIn >> val;
703 //*fIn >> val;
704 //*fIn >> val;
705 //-------------------
706
707
708
709 while(1)
710 {
711 *fIn >> ival;
712
713 if (fIn->eof())
714 return;
715
716 // run number must be greater than 10000
717 if (TMath::Abs(ival) < 10000)
718 {
719 *fLog << err << "===========> Error when reading file with theta and phi <========="
720 << " ival = " << ival << endl;
721 }
722 else
723 break;
724 }
725
726 fRuns += 1;
727 *fLog << fRuns <<"th run : " << ival << endl;
728
729 fRunNr.AddAt(ival, fRuns-1);
730
731 //*fLog << "check : fRuns, fRunNr[fRuns-1], fRunNr[fRuns] = " << fRuns << ", "
732 // << fRunNr[fRuns-1] << ", " << fRunNr[fRuns] << endl;
733
734
735 // read mjdS, hmsS, mjdE, hmsE
736 // these data are present only for ON data (fStars > 0)
737 /*
738 if (fStars > 0)
739 {
740 *fIn >> val;
741 *fIn >> val;
742 *fIn >> val;
743 *fIn >> val;
744
745 *fIn >> val;
746 *fIn >> val;
747 *fIn >> val;
748 *fIn >> val;
749 }
750 */
751
752 *fIn >> val;
753 fThetaTel.AddAt(val/kRad2Deg, fRuns-1);
754 //*fLog << "val, fThetaTel[fRuns-1] = " << val << ", "
755 // << fThetaTel[fRuns-1] << endl;
756
757
758 *fIn >> val;
759 fPhiTel.AddAt(val/kRad2Deg, fRuns-1);
760 //*fLog << "val, fPhiTel[fRuns-1] = " << val << ", "
761 // << fPhiTel[fRuns-1] << endl;
762
763
764 //*fIn >> val;
765 //fdThetaTel.AddAt(val/kRad2Deg, fRuns-1);
766 //*fIn >> val;
767 //fdPhiTel.AddAt(val/kRad2Deg, fRuns-1);
768
769 // input is in [deg], convert to [mm]
770
771 //*fLog << "ReadData : fStarsRead = " << fStarsRead << endl;
772
773 for (Int_t i=0; i<fStarsRead; i++)
774 {
775 *fIn >> val;
776 if (i<fStars) fxStar(i, fRuns-1) = val;
777
778 *fIn >> val;
779 if (i<fStars) fyStar(i, fRuns-1) = val;
780
781 //*fIn >> val;
782 // if (i < fStars) fdxStar(i, fRuns-1) = val;
783 if (i < fStars) fdxStar(i, fRuns-1) = 1.0;
784
785 //*fIn >> val;
786 // if (i < fStars) fdyStar(i, fRuns-1) = val;
787 if (i < fStars) fdyStar(i, fRuns-1) = 1.0;
788 }
789}
790
791// --------------------------------------------------------------------------
792//
793// Add this file as the last entry in the chain
794//
795Int_t MSourcePosfromStarPos::AddFile(const char *txt, Int_t)
796{
797 TNamed *name = new TNamed(txt, "");
798 fFileNames->AddLast(name);
799
800 *fLog << "MSourcePosfromStarPos::AddFile; add file '" << txt << "'"
801 << endl;
802
803 return 1;
804}
805
806// --------------------------------------------------------------------------
807//
808// This opens the next file in the list and deletes its name from the list.
809//
810Bool_t MSourcePosfromStarPos::OpenNextFile()
811{
812 //
813 // open the input stream and check if it is really open (file exists?)
814 //
815 if (fIn)
816 delete fIn;
817 fIn = NULL;
818
819 //
820 // Check for the existence of a next file to read
821 //
822 TNamed *file = (TNamed*)fFileNames->First();
823 if (!file)
824 return kFALSE;
825
826 //
827 // open the file which is the first one in the chain
828 //
829 const char *name = file->GetName();
830
831 const char *expname = gSystem->ExpandPathName(name);
832 fIn = new ifstream(expname);
833 delete [] expname;
834
835 const Bool_t noexist = !(*fIn);
836
837
838 if (noexist)
839 *fLog << dbginf << "Cannot open file '" << name << "'" << endl;
840 else
841 *fLog << "Open file: '" << name << "'" << endl;
842
843 //
844 // Remove this file from the list of pending files
845 //
846 fFileNames->Remove(file);
847
848 return !noexist;
849}
850
851// --------------------------------------------------------------------------
852
853
854
855
856
857
858
859
860
861
862
863
864
865
Note: See TracBrowser for help on using the repository browser.