source: trunk/MagicSoft/Mars/mpointing/MPointingPosInterpolate.cc@ 8664

Last change on this file since 8664 was 6508, checked in by mazin, 21 years ago
*** empty log message ***
File size: 10.8 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!
18! Author(s): Marcos Lopez 03/2004 <mailto:marcos@gae.ucm.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MExtrapolatePointingPos
28//
29// In the PreProcess, read the drive report and store it in an TSpline.
30// In the Process, use the TSpline to calculate the PointingPos for the
31// time of the current event.
32//
33// Input Containers:
34// MRawEvtData
35// MReportDrive
36// MTimeDrive
37// MTime
38//
39// Output Containers:
40// MPointingPos
41//
42//
43/////////////////////////////////////////////////////////////////////////////
44
45#include "MPointingPosInterpolate.h"
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50#include <TSpline.h>
51
52#include "MTaskList.h"
53#include "MParList.h"
54#include "MEvtLoop.h"
55#include "MReadReports.h"
56#include "MReportDrive.h"
57#include "MPointingPos.h"
58#include "MTime.h"
59#include "MRawRunHeader.h"
60#include "MDirIter.h"
61
62#include <TCanvas.h>
63
64ClassImp(MPointingPosInterpolate);
65
66using namespace std;
67
68const Int_t MPointingPosInterpolate::fgNumStartEvents = 1;
69// --------------------------------------------------------------------------
70//
71// Constructor
72//
73MPointingPosInterpolate::MPointingPosInterpolate(const char *name, const char *title)
74 : fEvtTime(NULL), fPointingPos(NULL), fRunHeader(NULL), fDirIter(NULL),
75 fSplineZd(NULL), fSplineAz(NULL), fRa(0.), fDec(0.),
76 fTimeMode(MPointingPosInterpolate::kEventTime)
77
78{
79 fName = name ? name : "MPointingPosInterpolate";
80 fTitle = title ? title : "Task to interpolate the pointing positons from drive reports";
81
82 fFilename = "";
83 fDebug = kFALSE;
84
85 SetNumStartEvents();
86}
87
88
89MPointingPosInterpolate::~MPointingPosInterpolate()
90{
91 Clear();
92}
93
94void MPointingPosInterpolate::Clear(Option_t *o)
95{
96 if(fSplineZd)
97 delete fSplineZd;
98 if(fSplineAz)
99 delete fSplineAz;
100}
101
102// ---------------------------------------------------------------------------
103//
104// Read the drive report file for the whole night, a build from it the splines
105//
106Bool_t MPointingPosInterpolate::ReadDriveReport()
107{
108
109 *fLog << inf << "Loading report file \"" << fFilename << "\" into TSpline..." << endl;
110
111 if (!fDebug)
112 gLog.SetNullOutput();
113
114 //
115 // ParList
116 //
117 MParList plist;
118 MTaskList tlist;
119 plist.AddToList(&tlist);
120
121
122 //
123 // TaskList
124 //
125 MReadReports read;
126 read.AddTree("Drive");
127 if (fDirIter)
128 read.AddFiles(*fDirIter);
129 else
130 read.AddFile(fFilename); // after the reading of the trees!!!
131
132 read.AddToBranchList("MReportDrive.*");
133
134 tlist.AddToList(&read);
135
136 //
137 // EventLoop
138 //
139 MEvtLoop evtloop;
140 evtloop.SetParList(&plist);
141
142 if (!evtloop.PreProcess())
143 {
144 gLog.SetNullOutput(kFALSE);
145 return kFALSE;
146 }
147
148 TArrayD reportTime(fNumStartEvents);
149 TArrayD currentZd(fNumStartEvents);
150 TArrayD currentAz(fNumStartEvents);
151 TArrayD nominalZd(fNumStartEvents);
152 TArrayD nominalAz(fNumStartEvents);
153
154 Int_t n=1;
155 while (tlist.Process())
156 {
157 MReportDrive* report = (MReportDrive*)plist.FindObject("MReportDrive");
158 MTime* reporttime = (MTime*)plist.FindObject("MTimeDrive");
159
160 if(n==1)
161 fFirstDriveTime = *reporttime;
162 else
163 fLastDriveTime = *reporttime;
164
165 //
166 // Update the number of entries
167 //
168 if (n>fNumStartEvents)
169 {
170 reportTime.Set(n);
171 currentZd.Set(n);
172 currentAz.Set(n);
173 nominalZd.Set(n);
174 nominalAz.Set(n);
175 }
176 //
177 // Sometimes there are two reports with the same time
178 //
179 if (n>1)
180 if (reporttime->GetTime() == reportTime[n-2])
181 {
182 *fLog << warn <<"["<< GetName()
183 << "]: Warning: this report has the same time that the previous one...skipping it " << endl;
184 continue;
185 }
186
187 reportTime[n-1] = reporttime->GetTime();
188 currentZd [n-1] = report->GetCurrentZd();
189 currentAz [n-1] = report->GetCurrentAz();
190 nominalZd [n-1] = report->GetNominalZd();
191 nominalAz [n-1] = report->GetNominalAz();
192 fRa = report->GetRa();
193 fDec = report->GetDec();
194
195 if (fDebug)
196 {
197 *fLog << " GetTime(): " << reporttime->GetTime() << endl;
198 *fLog << " GetCurrentZd(): " << report->GetCurrentZd() << endl;
199 }
200 n++;
201 }
202
203 tlist.PrintStatistics();
204
205 gLog.SetNullOutput(kFALSE);
206
207 *fLog << inf << GetDescriptor() << ": loaded " << n-1 << " ReportDrive from "
208 << fFirstDriveTime << " to " << fLastDriveTime << endl;
209
210 if (fDebug)
211 {
212 for (int i=0;i<reportTime.GetSize();i++)
213 *fLog << i << " " << reportTime[i] << " "
214 << currentZd[i] << " " << currentAz[i] << " "
215 << nominalZd[i] << " " << nominalAz[i] << endl;
216 }
217
218 fSplineZd = new TSpline3("zenith",
219 reportTime.GetArray(), nominalZd.GetArray(), n-1);
220 fSplineAz = new TSpline3("azimuth",
221 reportTime.GetArray(), nominalAz.GetArray(), n-1);
222
223
224 if (fDebug)
225 {
226 *fLog << n-1 << " " << reportTime.GetSize() << endl;
227 *fLog << n-1 << " " << nominalZd.GetSize() << endl;
228 *fLog << fFirstDriveTime.GetTime() << " " << fSplineZd->Eval((fFirstDriveTime.GetTime()+fLastDriveTime.GetTime())/2.) << endl;
229 TCanvas* c = new TCanvas();
230 c->Divide(2,1);
231 c->cd(1);
232 fSplineZd->Draw();
233 c->cd(2);
234 fSplineAz->Draw();
235 c->Modified();
236 c->Update();
237 c->SaveAs("pointing.root");
238 }
239
240 return kTRUE;
241}
242
243
244// --------------------------------------------------------------------------
245//
246// Input:
247// - MTime
248//
249// Output:
250// - MPointingPos
251//
252Int_t MPointingPosInterpolate::PreProcess( MParList *pList )
253{
254
255 Clear();
256
257 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
258 if (!fRunHeader)
259 {
260 *fLog << err << "MRunHeader not found... aborting." << endl;
261 return kFALSE;
262 }
263
264 fEvtTime = (MTime*)pList->FindObject("MTime");
265 if (!fEvtTime)
266 {
267 *fLog << err << "MTime not found... aborting." << endl;
268 return kFALSE;
269 }
270
271 fPointingPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
272 if (!fPointingPos)
273 return kFALSE;
274
275 if (fFilename.IsNull() && !fDirIter)
276 {
277 *fLog << err << "File name is empty or no MDirIter has been handed over... aborting" << endl;
278 return kFALSE;
279 }
280
281 if( !ReadDriveReport() )
282 return kFALSE;
283
284 return kTRUE;
285}
286
287
288// --------------------------------------------------------------------------
289//
290// Get the run start time, and get the pointing position for that time
291//
292Int_t MPointingPosInterpolate::Process()
293{
294
295 //const Int_t run = fRunHeader->GetRunNumber();
296 const MTime &StartRunTime = fRunHeader->GetRunStart();
297 const MTime &EndRunTime = fRunHeader->GetRunEnd();
298 Int_t time = StartRunTime.GetTime();
299
300 switch(fTimeMode)
301 {
302 case kRunTime:
303 //
304 // Check that we have drive report for this time
305 //
306
307 time = (EndRunTime.GetTime() + StartRunTime.GetTime())/2;
308
309 if( StartRunTime<fFirstDriveTime || StartRunTime>fLastDriveTime)
310 {
311 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
312 << " outside range of drive reports (" << fFirstDriveTime
313 << ", " << fLastDriveTime << ")" << endl;
314
315 if ( *fEvtTime<fFirstDriveTime ) time = fFirstDriveTime.GetTime();
316 if ( *fEvtTime>fLastDriveTime ) time = fLastDriveTime.GetTime();
317
318 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
319 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
320 }
321 break;
322
323 case kEventTime:
324
325 time = fEvtTime->GetTime();
326
327 //
328 // Check that we have drive report for this time
329 //
330 if( *fEvtTime<fFirstDriveTime || *fEvtTime>fLastDriveTime)
331 {
332 if (fDebug)
333 {
334 *fLog << err << GetDescriptor() << ": Run time = "
335 << *fEvtTime << " outside range of drive reports ("
336 << fFirstDriveTime << ", "<< fLastDriveTime << ")" << endl;
337 }
338
339 if ( *fEvtTime < (fFirstDriveTime) ) time = fFirstDriveTime.GetTime();
340 else time = fLastDriveTime.GetTime();
341
342 }
343 break;
344 }
345
346 if (fDebug)
347 {
348 *fLog << " real time : " << time
349 << " first time: " << fFirstDriveTime.GetTime()
350 << " last time: " << fLastDriveTime.GetTime() << endl;
351 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
352 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
353 }
354 //
355 // Check that we have drive report for this time
356 //
357 //if( *StartRunTime<fFirstDriveTime || *StartRunTime>fLastDriveTime)
358 if( StartRunTime>fLastDriveTime)
359 {
360 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
361 << " outside range of drive reports (" << fFirstDriveTime
362 << ", " << fLastDriveTime << ")" << endl;
363 return kERROR;
364 }
365
366 const Double_t zd = fSplineZd->Eval( time );
367 const Double_t az = fSplineAz->Eval( time );
368
369 if(TMath::Abs(zd)>90 || TMath::Abs(az)>360)
370 {
371 *fLog << warn << GetDescriptor() << ": Wrong Interpolated Pointing Position." << endl;
372 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = ("
373 << zd << ", " << az << "," << fRa << "," << fDec << ")" << endl;
374 }
375
376 fPointingPos->SetLocalPosition( zd, az );
377 //fPointingPos->SetSkyPosition( fRa*TMath::DegToRad()/15, fDec*TMath::DegToRad());
378 fPointingPos->SetSkyPosition( fRa, fDec);
379
380 return kTRUE;
381}
382
383
384
385Int_t MPointingPosInterpolate::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
386{
387 Bool_t rc = kFALSE;
388
389 if (IsEnvDefined(env, prefix, "NumStartEvents", print))
390 {
391 SetNumStartEvents(GetEnvValue(env, prefix, "NumStartEvents", fNumStartEvents));
392 rc = kTRUE;
393 }
394
395 if (IsEnvDefined(env, prefix, "TimeMode", print))
396 {
397
398 TString dat(GetEnvValue(env, prefix, "TimeMode", ""));
399 dat.ToLower();
400
401 rc = kTRUE;
402
403 if (dat.Contains("eventtime"))
404 SetTimeMode(kEventTime);
405 else if (dat.Contains("runtime"))
406 SetTimeMode(kRunTime);
407 else
408 rc = kFALSE;
409 }
410 return rc;
411}
412
Note: See TracBrowser for help on using the repository browser.