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

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