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

Last change on this file since 6491 was 6488, checked in by mazin, 20 years ago
*** empty log message ***
File size: 10.4 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-2])
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
223 if (fDebug)
224 {
225 *fLog << n-1 << " " << reportTime.GetSize() << endl;
226 *fLog << n-1 << " " << nominalZd.GetSize() << endl;
227 *fLog << fFirstDriveTime.GetTime() << " " << fSplineZd->Eval((fFirstDriveTime.GetTime()+fLastDriveTime.GetTime())/2.) << endl;
228 TCanvas* c = new TCanvas();
229 c->Divide(2,1);
230 c->cd(1);
231 fSplineZd->Draw();
232 c->cd(2);
233 fSplineAz->Draw();
234 c->Modified();
235 c->Update();
236 c->SaveAs("pointing.root");
237 }
238
239 return kTRUE;
240}
241
242
243// --------------------------------------------------------------------------
244//
245// Input:
246// - MTime
247//
248// Output:
249// - MPointingPos
250//
251Int_t MPointingPosInterpolate::PreProcess( MParList *pList )
252{
253
254 Clear();
255
256 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
257 if (!fRunHeader)
258 {
259 *fLog << err << "MRunHeader not found... aborting." << endl;
260 return kFALSE;
261 }
262
263 fEvtTime = (MTime*)pList->FindObject("MTime");
264 if (!fEvtTime)
265 {
266 *fLog << err << "MTime not found... aborting." << endl;
267 return kFALSE;
268 }
269
270 fPointingPos = (MPointingPos*)pList->FindCreateObj("MPointingPos");
271 if (!fPointingPos)
272 return kFALSE;
273
274 if (fFilename.IsNull() && !fDirIter)
275 {
276 *fLog << err << "File name is empty or no MDirIter has been handed over... aborting" << endl;
277 return kFALSE;
278 }
279
280 if( !ReadDriveReport() )
281 return kFALSE;
282
283 return kTRUE;
284}
285
286
287// --------------------------------------------------------------------------
288//
289// Get the run start time, and get the pointing position for that time
290//
291Int_t MPointingPosInterpolate::Process()
292{
293
294 //const Int_t run = fRunHeader->GetRunNumber();
295 const MTime &StartRunTime = fRunHeader->GetRunStart();
296 Int_t time = StartRunTime.GetTime();
297
298 switch(fTimeMode)
299 {
300 case kRunTime:
301 //
302 // Check that we have drive report for this time
303 //
304 if( StartRunTime<fFirstDriveTime || StartRunTime>fLastDriveTime)
305 {
306 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
307 << " outside range of drive reports (" << fFirstDriveTime
308 << ", " << fLastDriveTime << ")" << endl;
309
310 if ( *fEvtTime<fFirstDriveTime ) time = fFirstDriveTime.GetTime();
311 if ( *fEvtTime>fLastDriveTime ) time = fLastDriveTime.GetTime();
312
313 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
314 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
315 }
316 break;
317
318 case kEventTime:
319 time = fEvtTime->GetTime();
320
321 if (fDebug)
322 {
323 *fLog << " real time : " << time
324 << " first time: " << fFirstDriveTime.GetTime()
325 << " last time: " << fLastDriveTime.GetTime() << endl;
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 if (fDebug)
343 {
344 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az) = ("
345 << fSplineZd->Eval( time )<< ", " <<fSplineAz->Eval( time )<< ")" << endl;
346 }
347 }
348 break;
349 }
350
351 //
352 // Check that we have drive report for this time
353 //
354 //if( *StartRunTime<fFirstDriveTime || *StartRunTime>fLastDriveTime)
355 if( StartRunTime>fLastDriveTime)
356 {
357 *fLog << err << GetDescriptor() << ": Run time " << StartRunTime
358 << " outside range of drive reports (" << fFirstDriveTime
359 << ", " << fLastDriveTime << ")" << endl;
360 return kERROR;
361 }
362
363 const Double_t zd = fSplineZd->Eval( time );
364 const Double_t az = fSplineAz->Eval( time );
365
366 if(TMath::Abs(zd)>90 || TMath::Abs(az)>360)
367 {
368 *fLog << warn << GetDescriptor() << ": Wrong Interpolated Pointing Position." << endl;
369 *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = ("
370 << zd << ", " << az << "," << fRa << "," << fDec << ")" << endl;
371 }
372
373 fPointingPos->SetLocalPosition( zd, az );
374 //fPointingPos->SetSkyPosition( fRa*TMath::DegToRad()/15, fDec*TMath::DegToRad());
375 fPointingPos->SetSkyPosition( fRa, fDec);
376
377 return kTRUE;
378}
379
380
381
382Int_t MPointingPosInterpolate::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
383{
384 Bool_t rc = kFALSE;
385
386 if (IsEnvDefined(env, prefix, "NumStartEvents", print))
387 {
388 SetNumStartEvents(GetEnvValue(env, prefix, "NumStartEvents", fNumStartEvents));
389 rc = kTRUE;
390 }
391 return rc;
392}
Note: See TracBrowser for help on using the repository browser.