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

Last change on this file since 6289 was 6289, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 10.1 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 (reporttime->GetTime() == reportTime[n-1])
180 {
181 *fLog << warn <<"["<< GetName()
182 << "]: Warning: this report has the same time that the previous one...skipping it " << endl;
183 continue;
184 }
185
186 reportTime[n-1] = reporttime->GetTime();
187 currentZd [n-1] = report->GetCurrentZd();
188 currentAz [n-1] = report->GetCurrentAz();
189 nominalZd [n-1] = report->GetNominalZd();
190 nominalAz [n-1] = report->GetNominalAz();
191 fRa = report->GetRa();
192 fDec = report->GetDec();
193
194 if (fDebug)
195 {
196 *fLog << " GetTime(): " << reporttime->GetTime() << endl;
197 *fLog << " GetCurrentZd(): " << report->GetCurrentZd() << endl;
198 }
199 n++;
200 }
201
202 tlist.PrintStatistics();
203
204 gLog.SetNullOutput(kFALSE);
205
206 *fLog << inf << GetDescriptor() << ": loaded " << n-1 << " ReportDrive from "
207 << fFirstDriveTime << " to " << fLastDriveTime << endl;
208
209 if (fDebug)
210 {
211 for (int i=0;i<reportTime.GetSize();i++)
212 *fLog << i << " " << reportTime[i] << " "
213 << currentZd[i] << " " << currentAz[i] << " "
214 << nominalZd[i] << " " << nominalAz[i] << endl;
215 }
216
217 fSplineZd = new TSpline3("zenith",
218 reportTime.GetArray(), nominalZd.GetArray(), n-1);
219 fSplineAz = new TSpline3("azimuth",
220 reportTime.GetArray(), nominalAz.GetArray(), n-1);
221
222 if (fDebug)
223 {
224 TCanvas* c = new TCanvas();
225 c->Divide(2,1);
226 c->cd(1);
227 fSplineZd->Draw();
228 c->cd(2);
229 fSplineAz->Draw();
230 c->Modified();
231 c->Update();
232 c->SaveAs("pointing.root");
233 }
234
235 return kTRUE;
236}
237
238
239// --------------------------------------------------------------------------
240//
241// Input:
242// - MTime
243//
244// Output:
245// - MPointingPos
246//
247Int_t MPointingPosInterpolate::PreProcess( MParList *pList )
248{
249
250 Clear();
251
252 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
253 if (!fRunHeader)
254 {
255 *fLog << err << "MRunHeader not found... aborting." << endl;
256 return kFALSE;
257 }
258
259 fEvtTime = (MTime*)pList->FindObject("MTime");
260 if (!fEvtTime)
261 {
262 *fLog << err << "MTime not found... aborting." << endl;
263 return kFALSE;
264 }
265
266 fPointingPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
267 if (!fPointingPos)
268 return kFALSE;
269
270 if (fFilename.IsNull() && !fDirIter)
271 {
272 *fLog << err << "File name is empty or no MDirIter has been handed over... aborting" << endl;
273 return kFALSE;
274 }
275
276 if( !ReadDriveReport() )
277 return kFALSE;
278
279 return kTRUE;
280}
281
282
283// --------------------------------------------------------------------------
284//
285// Get the run start time, and get the pointing position for that time
286//
287Int_t MPointingPosInterpolate::Process()
288{
289
290 //const Int_t run = fRunHeader->GetRunNumber();
291 const MTime &StartRunTime = fRunHeader->GetRunStart();
292 Int_t time = StartRunTime.GetTime();
293
294 switch(fTimeMode)
295 {
296 case kRunTime:
297 //
298 // Check that we have drive report for this time
299 //
300 if( StartRunTime<fFirstDriveTime || StartRunTime>fLastDriveTime)
301 {
302 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
303 << " outside range of drive reports (" << fFirstDriveTime
304 << ", " << fLastDriveTime << ")" << endl;
305
306 if ( *fEvtTime<fFirstDriveTime ) time = fFirstDriveTime.GetTime();
307 if ( *fEvtTime>fLastDriveTime ) time = fLastDriveTime.GetTime();
308
309 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
310 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
311 }
312 break;
313
314 case kEventTime:
315 if (fDebug)
316 {
317 *fLog << " real time : " << time
318 << " first time: " << fFirstDriveTime.GetTime()
319 << " last time: " << fLastDriveTime.GetTime() << endl;
320 }
321 //
322 // Check that we have drive report for this time
323 //
324 if( *fEvtTime<fFirstDriveTime || *fEvtTime>fLastDriveTime)
325 {
326 if (fDebug)
327 {
328 *fLog << err << GetDescriptor() << ": Run time = "
329 << *fEvtTime << " outside range of drive reports ("
330 << fFirstDriveTime << ", "<< fLastDriveTime << ")" << endl;
331 }
332
333 if ( *fEvtTime < (fFirstDriveTime) ) time = fFirstDriveTime.GetTime();
334 else time = fLastDriveTime.GetTime();
335
336 if (fDebug)
337 {
338 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
339 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
340 }
341 }
342 break;
343 }
344
345 //
346 // Check that we have drive report for this time
347 //
348 //if( *StartRunTime<fFirstDriveTime || *StartRunTime>fLastDriveTime)
349 if( StartRunTime>fLastDriveTime)
350 {
351 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
352 << " outside range of drive reports (" << fFirstDriveTime
353 << ", " << fLastDriveTime << ")" << endl;
354 return kERROR;
355 }
356
357 const Double_t zd = fSplineZd->Eval( time );
358 const Double_t az = fSplineAz->Eval( time );
359
360 if(TMath::Abs(zd)>90 || TMath::Abs(az)>360)
361 {
362 *fLog << warn << GetDescriptor() << ": Wrong Interpolated Pointing Position." << endl;
363 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = ("
364 << zd << ", " << az << "," << fRa << "," << fDec << ")" << endl;
365 }
366
367 fPointingPos->SetLocalPosition( zd, az );
368 //fPointingPos->SetSkyPosition( fRa*TMath::DegToRad()/15, fDec*TMath::DegToRad());
369 fPointingPos->SetSkyPosition( fRa, fDec);
370
371 return kTRUE;
372}
373
374
375
376Int_t MPointingPosInterpolate::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
377{
378 Bool_t rc = kFALSE;
379
380 if (IsEnvDefined(env, prefix, "NumStartEvents", print))
381 {
382 SetNumStartEvents(GetEnvValue(env, prefix, "NumStartEvents", fNumStartEvents));
383 rc = kTRUE;
384 }
385 return rc;
386}
Note: See TracBrowser for help on using the repository browser.