source: branches/Corsika7405Compatibility/manalysis/MSourcePosfromStarPos.cc@ 18815

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