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