Index: trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.cc (revision 5943) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.cc (revision 5943) @@ -0,0 +1,394 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Marcos Lopez 03/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MExtrapolatePointingPos +// +// In the PreProcess, read the drive report and store it in an TSpline. +// In the Process, use the TSpline to calculate the PointingPos for the +// time of the current event. +// +// Input Containers: +// MRawEvtData +// +// Output Containers: +// MPointingPos +// +// +///////////////////////////////////////////////////////////////////////////// + +#include "MInterpolatePointingPos.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include + +#include "MTaskList.h" +#include "MParList.h" +#include "MEvtLoop.h" +#include "MReadReports.h" +#include "MReportDrive.h" +#include "MPointingPos.h" +#include "MTime.h" +#include "MRawRunHeader.h" +#include "MDirIter.h" + +#include + +ClassImp(MInterpolatePointingPos); + +using namespace std; + +const Int_t MInterpolatePointingPos::fgNumStartEvents = 10000; +// -------------------------------------------------------------------------- +// +// Constructor +// +MInterpolatePointingPos::MInterpolatePointingPos(const char *name, const char *title) + : fEvtTime(NULL), fPointingPos(NULL), fRunHeader(NULL), fDirIter(NULL), + fSplineZd(NULL), fSplineAz(NULL), fSplineRa(NULL), fSplineDec(NULL), + fTimeMode(MInterpolatePointingPos::kEventTime) + +{ + fName = name ? name : "MInterpolatePointingPos"; + fTitle = title ? title : "Task to interpolate the pointing positons from drive reports"; + + fFilename = ""; + fDebug = kFALSE; + + SetNumStartEvents(); +} + + +MInterpolatePointingPos::~MInterpolatePointingPos() +{ + Clear(); +} + +void MInterpolatePointingPos::Clear(Option_t *o) +{ + if(fSplineZd) + delete fSplineZd; + if(fSplineAz) + delete fSplineAz; + if(fSplineRa) + delete fSplineRa; + if(fSplineDec) + delete fSplineDec; +} + +// --------------------------------------------------------------------------- +// +// Read the drive report file for the whole night, a build from it the splines +// +Bool_t MInterpolatePointingPos::ReadDriveReport() +{ + + *fLog << endl << "["<< GetName() + << "]: Loading report file \"" << fFilename << "\" into TSpline..." << endl; + + // + // ParList + // + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + + // + // TaskList + // + MReadReports read; + read.AddTree("Drive"); + if (fDirIter) + read.AddFiles(*fDirIter); + else + read.AddFile(fFilename); // after the reading of the trees!!! + + read.AddToBranchList("MReportDrive.*"); + + tlist.AddToList(&read); + + // + // EventLoop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + if (!evtloop.PreProcess()) + return kFALSE; + + TArrayD reportTime(fNumStartEvents); + TArrayD currentZd(fNumStartEvents); + TArrayD currentAz(fNumStartEvents); + TArrayD nominalZd(fNumStartEvents); + TArrayD nominalAz(fNumStartEvents); + TArrayD ra(fNumStartEvents); + TArrayD dec(fNumStartEvents); + + Int_t n=0; + while (tlist.Process()) + { + MReportDrive* report = (MReportDrive*)plist.FindObject("MReportDrive"); + MTime* reporttime = (MTime*)plist.FindObject("MTimeDrive"); + + if(n==0) + fFirstDriveTime = *reporttime; + else + fLastDriveTime = *reporttime; + + // + // Sometimes there are two reports with the same time + // + if (reporttime->GetTime() == reportTime[n-1]) + { + cout << warn <<"["<< GetName() + << "]: Warning: this report has the same time that the previous one...skipping it " << endl; + continue; + } + + reportTime[n] = reporttime->GetTime(); + currentZd[n] = report->GetCurrentZd(); + currentAz[n] = report->GetCurrentAz(); + nominalZd[n] = report->GetNominalZd(); + nominalAz[n] = report->GetNominalAz(); + ra[n] = report->GetRa(); + dec[n] = report->GetDec(); + +//cout << " GetTime(): " << reporttime->GetTime() << endl; +//cout << " GetCurrentZd(): " << report->GetCurrentZd() << endl; + n++; + } + + tlist.PrintStatistics(); + + *fLog << "["<< GetName() << "]: loaded " << n << " ReportDrive from " + << fFirstDriveTime << " to " << fLastDriveTime << endl << endl; + + // + // Update the number of entries + // + reportTime.Set(n); + currentZd.Set(n); + currentAz.Set(n); + nominalZd.Set(n); + nominalAz.Set(n); + ra.Set(n); + dec.Set(n); + + if (fDebug) + { + for (int i=0;iDivide(2,2); + c->cd(1); + fSplineZd->Draw(); + c->cd(2); + fSplineAz->Draw(); + c->cd(3); + fSplineRa->Draw(); + c->cd(4); + fSplineDec->Draw(); + c->Modified(); + c->Update(); + } + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Input: +// - MTime +// +// Output: +// - MPointingPos +// +Int_t MInterpolatePointingPos::PreProcess( MParList *pList ) +{ + + Clear(); + + fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRunHeader) + { + *fLog << err << "MRunHeader not found... aborting." << endl; + return kFALSE; + } + + fEvtTime = (MTime*)pList->FindObject("MTime"); + if (!fEvtTime) + { + *fLog << err << "MTime not found... aborting." << endl; + return kFALSE; + } + + fPointingPos = (MPointingPos*)pList->FindCreateObj("MPointingPos"); + if (!fPointingPos) + return kFALSE; + + if (fFilename.IsNull() && !fDirIter) + { + *fLog << err << "File name is empty or no MDirIter has been handed over... aborting" << endl; + return kFALSE; + } + + if( !ReadDriveReport() ) + return kFALSE; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Get the run start time, and get the pointing position for that time +// +Int_t MInterpolatePointingPos::Process() +{ + + //const Int_t run = fRunHeader->GetRunNumber(); + const MTime* StartRunTime = &fRunHeader->GetRunStart(); + Double_t time = StartRunTime->GetTime(); + + switch(fTimeMode) + { + case kRunTime: + time = StartRunTime->GetTime(); + + // + // Check that we have drive report for this time + // + if( *StartRunTimefLastDriveTime) + { + *fLog << err << GetDescriptor() << ": Run time " << *StartRunTime + << " outside range of drive reports (" << fFirstDriveTime + << ", " << fLastDriveTime << ")" << endl; + + *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = (" << fSplineZd->Eval( time )<< ", " <Eval( time )<< ", " <Eval( time ) << ", " <Eval( time )<< ")" << endl; + + //fError = kTRUE; + //return kFALSE; + } + break; + + case kEventTime: + time = fEvtTime->GetTime(); + + if (fDebug) + { + cout << " real time : " << time + << " first time: " << fFirstDriveTime.GetTime() + << " last time: " << fLastDriveTime.GetTime() << endl; + } + // + // Check that we have drive report for this time + // + if( *fEvtTimefLastDriveTime) + { + *fLog << err << GetDescriptor() << ": Run time = " + << *fEvtTime << " outside range of drive reports (" + << fFirstDriveTime << ", "<< fLastDriveTime << ")" << endl; + + *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = (" << fSplineZd->Eval( time )<< ", " <Eval( time )<< ", " <Eval( time ) << ", " <Eval( time )<< ")" << endl; + + if ( *fEvtTimefLastDriveTime ) time = fLastDriveTime.GetTime(); + //fError = kTRUE; + //return kFALSE; + //return kCONTINUE; + } + break; + } + + + + + // + // Check that we have drive report for this time + // + if( *StartRunTimefLastDriveTime) + { + *fLog << err << GetDescriptor() << ": Run time " << *StartRunTime + << " outside range of drive reports (" << fFirstDriveTime + << ", " << fLastDriveTime << ")" << endl; + return kERROR; + } + + //if(run < 20000) + // time = fEvtTime->GetTime(); + //else + // time = StartRunTime->GetTime(); + + + const Double_t zd = fSplineZd->Eval( time ); + const Double_t az = fSplineAz->Eval( time ); + const Double_t ra = fSplineRa->Eval( time ); + const Double_t dec = fSplineDec->Eval( time ); + + + if(TMath::Abs(zd)>90 || TMath::Abs(az)>360 || ra>24 || ra<0 || TMath::Abs(dec)>90){ + *fLog << err << GetDescriptor() << ": Wrong Interpolated Pointing Position." << endl; + *fLog << " PointingPos: time = " << time << " (" << *fEvtTime << ") (zd, az, ra, dec) = (" << zd << ", " << az << ", " << ra << ", " << dec << ")" << endl; + } + + fPointingPos->SetLocalPosition( zd, az ); + fPointingPos->SetSkyPosition( ra*TMath::DegToRad()/15, dec*TMath::DegToRad()); + + return kTRUE; +} + + + +Int_t MInterpolatePointingPos::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + + if (IsEnvDefined(env, prefix, "NumStartEvents", print)) + { + SetNumStartEvents(GetEnvValue(env, prefix, "NumStartEvents", fNumStartEvents)); + rc = kTRUE; + } + return rc; +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.h (revision 5943) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MInterpolatePointingPos.h (revision 5943) @@ -0,0 +1,85 @@ +#ifndef MARS_MInterpolatePointingPos +#define MARS_MInterpolatePointingPos + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TSpline +#include +#endif + +#ifndef MARS_MTime +#include "MTime.h" +#endif + +class MTime; +class MPointingPos; +class MRawRunHeader; +class MDirIter; + +class MInterpolatePointingPos : public MTask +{ +public: + + enum TimeMode_t { + kRunTime, + kEventTime + }; + +private: + + static const Int_t fgNumStartEvents; //! Default for fNumStartEvents (now set to: 10000) + + Int_t fNumStartEvents; // Start number of allowed events + + Bool_t fDebug; + + TString fFilename; + MTime fFirstDriveTime; + MTime fLastDriveTime; + + MTime *fEvtTime; //! Raw event time + MPointingPos *fPointingPos; //! Telescope pointing postion + MRawRunHeader *fRunHeader; //! Run Header + MDirIter *fDirIter; //! Dir Iter + + TSpline3* fSplineZd; //! Zd vs. time + TSpline3* fSplineAz; //! Az vs. time + TSpline3* fSplineRa; //! Ra vs. time + TSpline3* fSplineDec; //! Dec vs. time + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Bool_t ReadDriveReport(); + + TimeMode_t fTimeMode; + + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +public: + + MInterpolatePointingPos(const char *name=NULL, const char *title=NULL); + + ~MInterpolatePointingPos(); + + void AddFiles(MDirIter *dir) { fDirIter = dir; } + void AddFile(const char *name) { fFilename = name; } + + void SetTimeMode( TimeMode_t mode) { fTimeMode = mode; } + void SetDebug( const Bool_t b=kTRUE) { fDebug = b; } + + void Clear(Option_t *o=""); + + Int_t GetNumStartEvents() const { return fNumStartEvents; } + void SetNumStartEvents ( const Int_t i=fgNumStartEvents ) { fNumStartEvents = i; } + + ClassDef(MInterpolatePointingPos, 1) // Interpolate the drive pointing positions +}; + +#endif + + + +