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

Last change on this file since 6207 was 6169, checked in by mazin, 20 years ago
*** empty log message ***
File size: 10.0 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
326 if ( *fEvtTime < (fFirstDriveTime) ) time = fFirstDriveTime.GetTime();
327 else time = fLastDriveTime.GetTime();
328
329 if (fDebug)
330 {
331 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
332 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
333 }
334 }
335 break;
336 }
337
338 //
339 // Check that we have drive report for this time
340 //
341 //if( *StartRunTime<fFirstDriveTime || *StartRunTime>fLastDriveTime)
342 if( *StartRunTime>fLastDriveTime)
343 {
344 *fLog << err << GetDescriptor() << ": Run time " << *StartRunTime
345 << " outside range of drive reports (" << fFirstDriveTime
346 << ", " << fLastDriveTime << ")" << endl;
347 return kERROR;
348 }
349
350 const Double_t zd = fSplineZd->Eval( time );
351 const Double_t az = fSplineAz->Eval( time );
352
353 if(TMath::Abs(zd)>90 || TMath::Abs(az)>360)
354 {
355 *fLog << warn << GetDescriptor() << ": Wrong Interpolated Pointing Position." << endl;
356 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = ("
357 << zd << ", " << az << "," << fRa << "," << fDec << ")" << endl;
358 }
359
360 fPointingPos->SetLocalPosition( zd, az );
361 //fPointingPos->SetSkyPosition( fRa*TMath::DegToRad()/15, fDec*TMath::DegToRad());
362 fPointingPos->SetSkyPosition( fRa, fDec);
363
364 return kTRUE;
365}
366
367
368
369Int_t MPointingPosInterpolate::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
370{
371 Bool_t rc = kFALSE;
372
373 if (IsEnvDefined(env, prefix, "NumStartEvents", print))
374 {
375 SetNumStartEvents(GetEnvValue(env, prefix, "NumStartEvents", fNumStartEvents));
376 rc = kTRUE;
377 }
378 return rc;
379}
Note: See TracBrowser for help on using the repository browser.