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

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