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

Last change on this file since 3819 was 3599, checked in by gaug, 21 years ago
*** empty log message ***
File size: 25.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17! Author(s): Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
18!
19! Copyright: MAGIC Software Development, 2000-2004
20!
21!
22\* ======================================================================== */
23
24/////////////////////////////////////////////////////////////////////////////
25//
26// MSourcePosfromStarPos
27//
28// This is a task which
29// - calculates the position of the source in the camera
30// from the position of a known star in the camera
31// - and puts the source position into the container MSrcPosCam
32//
33// Input :
34// ASCII file containing for each run
35// - the run number
36// - the direction (theta, phi) the telescope is pointing to in [deg]
37// - the position (xStar, yStar) of a known star in the camera in [mm]
38// - the error (dxStar, dyStar) of this position in [mm]
39//
40// Output Containers :
41// MSrcPosCam
42//
43/////////////////////////////////////////////////////////////////////////////
44#include <TList.h>
45#include <TSystem.h>
46#include <TMatrixD.h>
47
48#include <fstream>
49
50#include "MSourcePosfromStarPos.h"
51
52#include "MParList.h"
53#include "MRawRunHeader.h"
54#include "MGeomCam.h"
55#include "MSrcPosCam.h"
56#include "MMcEvt.hxx"
57
58#include "MLog.h"
59#include "MLogManip.h"
60#include "MObservatory.h"
61
62ClassImp(MSourcePosfromStarPos);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default constructor.
69//
70MSourcePosfromStarPos::MSourcePosfromStarPos(
71 const char *name, const char *title)
72 : fIn(NULL)
73{
74 fName = name ? name : "MSourcePosfromStarPos";
75 fTitle = title ? title : "Calculate source position from star position";
76
77 fFileNames = new TList;
78 fFileNames->SetOwner();
79
80 fRuns = 0;
81 fSize = 100;
82
83 fRunNr.Set(fSize);
84
85 fThetaTel.Set(fSize);
86 fPhiTel.Set(fSize);
87 fdThetaTel.Set(fSize);
88 fdPhiTel.Set(fSize);
89
90 Int_t fRows = 1;
91 fDecStar.Set(fRows);
92 fRaStar.Set(fRows);
93 fxStar.ResizeTo(fRows,fSize);
94 fyStar.ResizeTo(fRows,fSize);
95 fdxStar.ResizeTo(fRows,fSize);
96 fdyStar.ResizeTo(fRows,fSize);
97
98 fStars = 0;
99 fDecSource = 0.0;
100 fRaSource = 0.0;
101
102 // these are the default values when starting the excution;
103 // later these locations contain the values of the last event
104 fxSourceold = 25.0;
105 fySourceold = -40.0;
106 fThetaradold = 25.0/kRad2Deg;
107 fPhiradold = 180.0/kRad2Deg;
108}
109
110// --------------------------------------------------------------------------
111//
112// Delete the filename list and the input stream if one exists.
113//
114MSourcePosfromStarPos::~MSourcePosfromStarPos()
115{
116 delete fFileNames;
117 if (fIn)
118 delete fIn;
119}
120
121
122// --------------------------------------------------------------------------
123//
124// Set the sky coordinates of the source and of the star
125//
126// Input :
127// declination in units of (Deg, Min, Sec)
128// right ascension in units of (Hour, Min, Sec)
129//
130
131void MSourcePosfromStarPos::SetSourceAndStarPosition(
132 TString nameSource,
133 Double_t decSourceDeg, Double_t decSourceMin, Double_t decSourceSec,
134 Double_t raSourceHour, Double_t raSourceMin, Double_t raSourceSec,
135 TString nameStar,
136 Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec,
137 Double_t raStarHour, Double_t raStarMin, Double_t raStarSec )
138{
139 *fLog << "MSourcePosfromStarPos::SetSourceAndStarPosition :" << endl;
140 *fLog << " Source (dec) : " << nameSource << " " << decSourceDeg << ":"
141 << decSourceMin << ":" << decSourceSec << endl;
142 *fLog << " Source (ra) : " << nameSource << " " << raSourceHour << ":"
143 << raSourceMin << ":" << raSourceSec << endl;
144
145 *fLog << " Star (dec) : " << nameStar << " " << decStarDeg << ":"
146 << decStarMin << ":" << decStarSec << endl;
147 *fLog << " Star (ra) : " << nameStar << " " << raStarHour << ":"
148 << raStarMin << ":" << raStarSec << endl;
149
150 // convert into radians
151 fDecSource = (decSourceDeg + decSourceMin/60.0 + decSourceSec/3600.0)
152 / kRad2Deg;
153 fRaSource = (raSourceHour + raSourceMin/60.0 + raSourceSec/3600.0)
154 * 360.0 / (24.0 * kRad2Deg);
155
156 fStars += 1;
157 fDecStar.Set(fStars);
158 fRaStar.Set(fStars);
159 fxStar.ResizeTo(fStars,fSize);
160 fyStar.ResizeTo(fStars,fSize);
161 fdxStar.ResizeTo(fStars,fSize);
162 fdyStar.ResizeTo(fStars,fSize);
163
164 fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
165 / kRad2Deg;
166 fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0)
167 * 360.0 / (24.0 * kRad2Deg);
168
169 *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fStars = "
170 << fStars << endl;
171 *fLog << all << " fDecSource, fRaSource, fDecStar, fRaStar were set to : [radians] "
172 << fDecSource << ", " << fRaSource << ", "
173 << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl;
174}
175
176// --------------------------------------------------------------------------
177//
178// Set the sky coordinates of another star
179//
180// Input :
181// declination in units of (Deg, Min, Sec)
182// right ascension in units of (Hour, Min, Sec)
183//
184
185void MSourcePosfromStarPos::AddStar(
186 TString nameStar,
187 Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec,
188 Double_t raStarHour, Double_t raStarMin, Double_t raStarSec )
189{
190 *fLog << "MSourcePosfromStarPos::AddStar :" << endl;
191 *fLog << " Star (dec) : " << nameStar << " " << decStarDeg << ":"
192 << decStarMin << ":" << decStarSec << endl;
193 *fLog << " Star (ra) : " << nameStar << " " << raStarHour << ":"
194 << raStarMin << ":" << raStarSec << endl;
195
196 // convert into radians
197 fStars += 1;
198 fDecStar.Set(fStars);
199 fRaStar.Set(fStars);
200 fxStar.ResizeTo(fStars,fSize);
201 fyStar.ResizeTo(fStars,fSize);
202 fdxStar.ResizeTo(fStars,fSize);
203 fdyStar.ResizeTo(fStars,fSize);
204
205 fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0)
206 / kRad2Deg;
207 fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0)
208 * 360.0 / (24.0 * kRad2Deg);
209
210 *fLog << all << "MSourcePosfromStarPos::AddStar; fStars = " << fStars
211 << endl;
212 *fLog << all << " fDecStar, fRaStar were set to : [radians] "
213 << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl;
214}
215
216// --------------------------------------------------------------------------
217//
218//
219Int_t MSourcePosfromStarPos::PreProcess(MParList *pList)
220{
221 MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
222 if (!geom)
223 {
224 *fLog << err << "MSourcePosfromStarPos : MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
225 return kFALSE;
226 }
227 fMm2Deg = geom->GetConvMm2Deg();
228 // fDistCameraReflector is the distance of the camera from the reflector center in [mm]
229 fDistCameraReflector = kRad2Deg / fMm2Deg;
230 *fLog << all << "MSourcePosfromStarPos::PreProcess; fMm2Deg, fDistCameraReflector = "
231 << fMm2Deg << ", " << fDistCameraReflector << endl;
232
233 fObservatory = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory"));
234 if (!fObservatory)
235 {
236 *fLog << err << "MObservatory not found... aborting" << endl;
237 return kFALSE;
238 }
239
240 fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
241 if (!fRun)
242 {
243 *fLog << err << "MSourcePosfromStarPos::PreProcess; MRawRunHeader not found... aborting." << endl;
244 return kFALSE;
245 }
246
247
248 fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt");
249 if (!fMcEvt)
250 {
251 *fLog << err << "MSourcePosfromStarPos::PreProcess; MMcEvt not found... aborting." << endl;
252 return kFALSE;
253 }
254
255
256 fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber("MSrcPosCam"));
257 if (!fSrcPos)
258 {
259 *fLog << err << "MSourcePosfromStarPos::PreProcess; MSrcPosCam not found... aborting" << endl;
260 return kFALSE;
261 }
262
263 //---------------------------------------------------------------------
264 // read all files and call ReadData() to read and store the information
265 //
266
267 *fLog << all << "---------------------------------" << endl;
268 while(1)
269 {
270 if (!OpenNextFile())
271 {
272 *fLog << "there is no more file to open" << endl;
273 break;
274 }
275
276 *fLog << "read data" << endl;
277
278 // read "fStarsRead" = no.of (x,y) pairs to be read
279 *fIn >> fStarsRead;
280
281 while (1)
282 {
283 if (fIn->eof())
284 {
285 *fLog << "eof encountered; open next file" << endl;
286
287 if (!OpenNextFile()) break;
288 }
289
290 // FIXME! Set InputStreamID
291
292 ReadData();
293 }
294 }
295
296 *fLog << "all data were read" << endl;
297 FixSize();
298
299 if (fDecSource == 0.0 || fRaSource == 0.0 || fStars == 0)
300 {
301 *fLog << warn << "MSourcePosfromStarPos::PreProcess; there are no sky coordinates defined for the source or from stars; fStars, fStarsRead = "
302 << fStars << ", " << fStarsRead
303 << endl;
304 }
305 *fLog << all << "---------------------------------" << endl;
306
307 //-------------------------------------------------------------
308
309 return kTRUE;
310}
311
312//=========================================================================
313//
314// SourcefromStar
315//
316// this routine calculates the position of a source (for example Crab) in the camera
317// from the position of a star (for example ZetaTauri) in the camera. The latter
318// position may have been determined by analysing the DC currents in the different
319// pixels.
320//
321// Input : thetaTel, phiTel the direction the telescope is pointing to,
322// in local coordinates
323// f the distance between camera and reflector
324// decStar, raStar the position of the star in sky coordinates
325// decSource, raSource the position of the source in sky coordinates
326// xStar, yStar the position of the star in the camera
327// dxStar, dyStar error of the position of the star in the camera
328//
329// Output : xSource, ySource the calculated position of the source in the camera
330// dxSource, dySource error of the calculated position of the source in
331// the camera
332//
333// Useful formulas can be found in TDAS 00-11 and TDAS 01-05
334//
335
336void MSourcePosfromStarPos::SourcefromStar(Double_t &f,
337 TArrayD &decStar, TArrayD &raStar,
338 Double_t &decSource, Double_t &raSource,
339 Double_t &thetaTel, Double_t &phiTel,
340 TArrayD &xStar, TArrayD &yStar,
341 TArrayD &dxStar, TArrayD &dyStar,
342 Double_t &xSource, Double_t &ySource,
343 Double_t &dxSource, Double_t &dySource)
344{
345 /*
346 *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in degrees" << endl;
347 *fLog << " decStar, raStar = " << decStar[0]*kRad2Deg << ", "
348 << raStar[0]*kRad2Deg << endl;
349 *fLog << " decSource, raSource = " << decSource*kRad2Deg << ", "
350 << raSource*kRad2Deg << endl;
351 *fLog << " thetaTel, phiTel = " << thetaTel*kRad2Deg << ", "
352 << phiTel*kRad2Deg << endl;
353 *fLog << " xStar, yStar = " << xStar[0]*fMm2Deg << ", "
354 << yStar[0]*fMm2Deg << endl;
355
356 *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in radians and mm" << endl;
357 *fLog << " decStar, raStar = " << decStar[0] << ", "
358 << raStar[0] << endl;
359 *fLog << " decSource, raSource = " << decSource << ", "
360 << raSource << endl;
361 *fLog << " thetaTel, phiTel = " << thetaTel << ", "
362 << phiTel << endl;
363 *fLog << " xStar, yStar = " << xStar[0] << ", "
364 << yStar[0] << endl;
365 */
366
367 // the units are assumed to be radians for theta, phi, dec and ra
368 // and mm for f, x and y
369
370
371 // calculate rotation angle alpha of sky image in camera
372 // (see TDAS 00-11, eqs. (18) and (20))
373 // a1 = cos(Lat), a3 = -sin(Lat), where Lat is the geographical latitude of La Palma
374 Double_t a1 = 0.876627;
375 Double_t a3 = -0.481171;
376
377 Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
378 ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
379 ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) );
380 Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
381 Double_t sinal = a1 * sin(phiTel) * denom;
382
383 *fLog << "old thetaTel, phiTel, cosal, sinal = " << thetaTel << ", "
384 << phiTel << ", " << cosal << ", " << sinal << endl;
385
386
387 fObservatory->RotationAngle(thetaTel, phiTel, sinal, cosal);
388
389 *fLog << "new thetaTel, phiTel, cosal, sinal = " << thetaTel << ", "
390 << phiTel << ", " << cosal << ", " << sinal << endl;
391
392
393 // calculate coordinates of source in system B (see TDAS 00-11, eqs. (2))
394 // note that right ascension and hour angle go into opposite directions
395 Double_t xB0 = cos(decSource) * cos(-raSource);
396 Double_t yB0 = cos(decSource) * sin(-raSource);
397 Double_t zB0 = -sin(decSource);
398
399 //*fLog << "xB0, yB0, zB0 = " << xB0 << ", " << yB0 << ", "
400 // << zB0 << endl;
401
402 //-----------------------------------------------------
403 // loop over stars
404 Double_t sumx = 0.0;
405 Double_t sumy = 0.0;
406 Double_t sumwx = 0.0;
407 Double_t sumwy = 0.0;
408
409 for (Int_t i=0; i<decStar.GetSize(); i++)
410 {
411 // calculate weights
412 Double_t weightx = 1.0 / (dxStar[i]*dxStar[i]);
413 Double_t weighty = 1.0 / (dyStar[i]*dyStar[i]);
414 sumwx += weightx;
415 sumwy += weighty;
416
417 //*fLog << "weightx, weighty = " << weightx << ", " << weighty << endl;
418
419 // calculate coordinates of star in system B (see TDAS 00-11, eqs. (2))
420 // note that right ascension and hour angle go into opposite directions
421 Double_t xB = cos(decStar[i]) * cos(-raStar[i]);
422 Double_t yB = cos(decStar[i]) * sin(-raStar[i]);
423 Double_t zB = -sin(decStar[i]);
424
425
426 //*fLog << "xB, yB, zB = " << xB << ", " << yB << ", "
427 // << zB << endl;
428
429
430 // calculate coordinates of star in a system with the basis vectors e1, e2, e3
431 // where e1 is in the direction (r0 x a)
432 // e2 is in the direction (e1 x r0)
433 // and e3 is in the direction -r0;
434 // r0 is the direction to the source
435 // and a is the earth rotation axis (pointing to the celestial north pole)
436 //
437 Double_t x = (-xB*yB0 + xB0*yB) / sqrt( xB0*xB0 + yB0*yB0 );
438 Double_t y = ( xB*xB0*zB0 + yB*yB0*zB0 - zB*(xB0*xB0 + yB0*yB0) )
439 / sqrt( xB0*xB0 + yB0*yB0 );
440 Double_t z = -(xB*xB0 + yB*yB0 + zB*zB0);
441
442 //*fLog << "x, y, z = " << x << ", " << y << ", "
443 // << z << endl;
444
445
446 // calculate coordinates of star in camera
447 Double_t xtilde = -f/z * (cosal*x - sinal*y);
448 Double_t ytilde = -f/z * (sinal*x + cosal*y);
449
450 //*fLog << "i, xtilde, ytile = " << i << " : " << xtilde << ", "
451 // << ytilde << endl;
452
453
454 // calculate coordinates of source in camera
455 // note : in real camera signs are inverted (therefore s = -1.0)
456 Double_t s = -1.0;
457
458 Double_t xs = xStar[i] - s * xtilde;
459 Double_t ys = yStar[i] - s * ytilde;
460
461 *fLog << "i, xs, ys = " << i << " : " << xs << ", "
462 << ys << endl;
463
464 sumx += xs * weightx;
465 sumy += ys * weighty;
466 }
467 //-----------------------------------------------------
468
469 xSource = sumx / sumwx;
470 ySource = sumy / sumwy;
471 dxSource = 1.0 / sqrt(sumwx);
472 dySource = 1.0 / sqrt(sumwy);
473
474 /*
475 Int_t run = fRun->GetRunNumber();
476 *fLog << all << "MSourcePosfromStarPos::SourcefromStar; run, xSource, ySource = "
477 << run << " : "
478 << xSource << " +- " << dxSource << ", "
479 << ySource << " +- " << dySource << endl;
480 */
481}
482
483// --------------------------------------------------------------------------
484//
485// Get the source position and put it into MSrcPosCam
486//
487//
488Bool_t MSourcePosfromStarPos::ReInit(MParList *pList)
489{
490 //if (1 == 1) return kTRUE;
491
492
493 Int_t run = fRun->GetRunNumber();
494 *fLog << all << "MSourcePosfromStarPos::ReInit; run = " << run << endl;
495
496
497 //-------------------------------------------------------------------
498 // search this run in the list
499 for (Int_t i=0; i<fSize; i++)
500 {
501 if (run == fRunNr[i])
502 {
503 //-----------------------------------------
504 // put the zenith angle into MMcEvt
505
506 Double_t thetarad = fThetaTel[i];
507 Double_t phirad = fPhiTel[i];
508
509 if (fabs(thetarad*kRad2Deg+1.0) < 0.001)
510 thetarad = fThetaradold;
511 else
512 fThetaradold = thetarad;
513 if (fabs(phirad*kRad2Deg+1.0) < 0.001)
514 phirad = fPhiradold;
515 else
516 fPhiradold = phirad;
517
518 fMcEvt->SetTelescopeTheta(thetarad);
519 fMcEvt->SetTelescopePhi(phirad);
520 fMcEvt->SetReadyToSave();
521
522 *fLog << all << "theta, phi = " << thetarad*kRad2Deg << ", "
523 << phirad*kRad2Deg << " deg" << endl;
524
525 //-----------------------------------------
526 // Get source position and put it into MSrcPosCam
527
528
529 if (fStars > 0 && fStars == fStarsRead)
530 {
531 TArrayD xStar(fxStar.GetNrows());
532 TArrayD dxStar(fdxStar.GetNrows());
533 TArrayD yStar(fyStar.GetNrows());
534 TArrayD dyStar(fdyStar.GetNrows());
535 for (Int_t j=0; j<fxStar.GetNrows(); j++)
536 {
537 xStar[j] = fxStar(j, i);
538 dxStar[j] = fdxStar(j, i);
539 yStar[j] = fyStar(j, i);
540 dyStar[j] = fdyStar(j, i);
541 }
542
543 SourcefromStar( fDistCameraReflector,
544 fDecStar, fRaStar, fDecSource, fRaSource,
545 thetarad, phirad,
546 xStar, yStar,
547 dxStar, dyStar,
548 fxSource, fySource,
549 fdxSource, fdySource);
550
551 fSrcPos->SetXY(fxSource, fySource);
552
553 fxSourceold = fxSource;
554 fySourceold = fySource;
555
556 *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = "
557 << fRunNr[i] << ", " << fxSource << " +- " << fdxSource
558 << ", " << fySource << " +- " << fdySource << endl;
559
560 fSrcPos->SetReadyToSave();
561 }
562 else
563 {
564 // set default values
565 fxSource = fxSourceold;
566 fySource = fySourceold;
567 fSrcPos->SetXY(fxSource, fySource);
568 fSrcPos->SetReadyToSave();
569
570
571 *fLog << warn << "MSourcePosfromStarPos::ReInit; no information on source position for run number = "
572 << run << endl;
573 *fLog << warn << " set xSource, ySource = " << fxSource << ", "
574 << fySource << " mm" << endl;
575 }
576
577
578 return kTRUE;
579 }
580 }
581 //-------------------------------------------------------------------
582
583 // set default values
584 fxSource = fxSourceold;
585 fySource = fySourceold;
586 fSrcPos->SetXY(fxSource, fySource);
587 fSrcPos->SetReadyToSave();
588
589 Double_t thetarad = fThetaradold;
590 fMcEvt->SetTelescopeTheta(thetarad);
591
592 Double_t phirad = fPhiradold;
593 fMcEvt->SetTelescopePhi(phirad);
594 fMcEvt->SetReadyToSave();
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.