/* ======================================================================== *\ ! ! * ! * 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): Wolfgang Wittek 02/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MSourcePosfromStarPos // // This is a task which // - calculates the position of the source in the camera // from the position of a known star in the camera // - and puts the source position into the container MSrcPosCam // // Input : // ASCII file containing for each run // - the run number // - the direction (theta, phi) the telescope is pointing to in [deg] // - the position (xStar, yStar) of a known star in the camera in [mm] // - the error (dxStar, dyStar) of this position in [mm] // // Output Containers : // MSrcPosCam // ///////////////////////////////////////////////////////////////////////////// #include #include #include #include #include "MSourcePosfromStarPos.h" #include "MParList.h" #include "MRawRunHeader.h" #include "MGeomCam.h" #include "MSrcPosCam.h" #include "MMcEvt.hxx" #include "MLog.h" #include "MLogManip.h" ClassImp(MSourcePosfromStarPos); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // MSourcePosfromStarPos::MSourcePosfromStarPos( const char *name, const char *title) : fIn(NULL) { fName = name ? name : "MSourcePosfromStarPos"; fTitle = title ? title : "Calculate source position from star position"; fFileNames = new TList; fFileNames->SetOwner(); fRuns = 0; fSize = 100; fStars = 0; fRunNr.Set(fSize); fThetaTel.Set(fSize); fPhiTel.Set(fSize); fdThetaTel.Set(fSize); fdPhiTel.Set(fSize); fDecStar.Set(1); fRaStar.Set(1); fxStar.ResizeTo(1,fSize); fyStar.ResizeTo(1,fSize); fdxStar.ResizeTo(1,fSize); fdyStar.ResizeTo(1,fSize); } // -------------------------------------------------------------------------- // // Delete the filename list and the input stream if one exists. // MSourcePosfromStarPos::~MSourcePosfromStarPos() { delete fFileNames; if (fIn) delete fIn; } // -------------------------------------------------------------------------- // // Set the sky coordinates of the source and of the star // // Input : // declination in units of (Deg, Min, Sec) // right ascension in units of (Hour, Min, Sec) // void MSourcePosfromStarPos::SetSourceAndStarPosition( TString nameSource, Double_t decSourceDeg, Double_t decSourceMin, Double_t decSourceSec, Double_t raSourceHour, Double_t raSourceMin, Double_t raSourceSec, TString nameStar, Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec, Double_t raStarHour, Double_t raStarMin, Double_t raStarSec ) { *fLog << "MSourcePosfromStarPos::SetSourceAndStarPosition :" << endl; *fLog << " Source : " << nameSource << " " << decSourceDeg << ":" << decSourceMin << ":" << decSourceSec << endl; *fLog << " Star : " << nameStar << " " << decStarDeg << ":" << decStarMin << ":" << decStarSec << endl; // convert into radians fDecSource = (decSourceDeg + decSourceMin/60.0 + decSourceSec/3600.0) / kRad2Deg; fRaSource = (raSourceHour + raSourceMin/60.0 + raSourceSec/3600.0) * 360.0 / (24.0 * kRad2Deg); fStars += 1; fDecStar.Set(fStars); fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0) / kRad2Deg; fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0) * 360.0 / (24.0 * kRad2Deg); *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fDecSource, fRaSource, fDecStar, fRaStar were set to : [radians] " << fDecSource << ", " << fRaSource << ", " << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl; } // -------------------------------------------------------------------------- // // Set the sky coordinates of another star // // Input : // declination in units of (Deg, Min, Sec) // right ascension in units of (Hour, Min, Sec) // void MSourcePosfromStarPos::AddStar( TString nameStar, Double_t decStarDeg, Double_t decStarMin, Double_t decStarSec, Double_t raStarHour, Double_t raStarMin, Double_t raStarSec ) { *fLog << "MSourcePosfromStarPos::AddStar :" << endl; *fLog << " Star : " << nameStar << " " << decStarDeg << ":" << decStarMin << ":" << decStarSec << endl; // convert into radians fStars = fDecStar.GetSize() + 1; fDecStar.Set(fStars); fDecStar[fStars-1] = (decStarDeg + decStarMin/60.0 + decStarSec/3600.0) / kRad2Deg; fRaStar[fStars-1] = (raStarHour + raStarMin/60.0 + raStarSec/3600.0) * 360.0 / (24.0 * kRad2Deg); *fLog << all << "MSourcePosfromStarPos::AddStar; fDecStar, fRaStar were set to : [radians] " << fDecStar[fStars-1] << ", " << fRaStar[fStars-1] << endl; } // -------------------------------------------------------------------------- // // Int_t MSourcePosfromStarPos::PreProcess(MParList *pList) { MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); if (!geom) { *fLog << err << "MSourcePosfromStarPos : MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; return kFALSE; } fMm2Deg = geom->GetConvMm2Deg(); // fDistCameraReflector is the distance of the camera from the reflector center in [mm] fDistCameraReflector = kRad2Deg / fMm2Deg; *fLog << all << "MSourcePosfromStarPos::PreProcess; fMm2Deg, fDistCameraReflector = " << fMm2Deg << ", " << fDistCameraReflector << endl; fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!fRun) { *fLog << err << "MSourcePosfromStarPos::PreProcess; MRawRunHeader not found... aborting." << endl; return kFALSE; } fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt"); if (!fMcEvt) { *fLog << err << "MSourcePosfromStarPos::PreProcess; MMcEvt not found... aborting." << endl; return kFALSE; } fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber("MSrcPosCam")); if (!fSrcPos) { *fLog << err << "MSourcePosfromStarPos::PreProcess; MSrcPosCam not found... aborting" << endl; return kFALSE; } //--------------------------------------------------------------------- // read all files and call ReadData() to read and store the information // *fLog << all << "---------------------------------" << endl; while(1) { if (!OpenNextFile()) { *fLog << "there is no more file to open" << endl; break; } *fLog << "read data" << endl; while (1) { if (fIn->eof()) { *fLog << "eof encountered; open next file" << endl; if (!OpenNextFile()) break; } // FIXME! Set InputStreamID ReadData(); } } *fLog << "all data were read" << endl; FixSize(); if (fDecSource == 0.0 || fRaSource == 0.0 || fStars == 0) { *fLog << err << "MSourcePosfromStarPos::PreProcess; there are no sky coordinates defined for the source or from stars" << endl; } *fLog << all << "---------------------------------" << endl; //------------------------------------------------------------- return kTRUE; } //========================================================================= // // SourcefromStar // // this routine calculates the position of a source (for example Crab) in the camera // from the position of a star (for example ZetaTauri) in the camera. The latter // position may have been determined by analysing the DC currents in the different // pixels. // // Input : thetaTel, phiTel the direction the telescope is pointing to, // in local coordinates // f the distance between camera and reflector // decStar, raStar the position of the star in sky coordinates // decSource, raSource the position of the source in sky coordinates // xStar, yStar the position of the star in the camera // dxStar, dyStar error of the position of the star in the camera // // Output : xSource, ySource the calculated position of the source in the camera // dxSource, dySource error of the calculated position of the source in // the camera // // Useful formulas can be found in TDAS 00-11 and TDAS 01-05 // void MSourcePosfromStarPos::SourcefromStar(Double_t &f, TArrayD &decStar, TArrayD &raStar, Double_t &decSource, Double_t &raSource, Double_t &thetaTel, Double_t &phiTel, TArrayD &xStar, TArrayD &yStar, TArrayD &dxStar, TArrayD &dyStar, Double_t &xSource, Double_t &ySource, Double_t &dxSource, Double_t &dySource) { /* *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in degrees" << endl; *fLog << " decStar, raStar = " << decStar[0]*kRad2Deg << ", " << raStar[0]*kRad2Deg << endl; *fLog << " decSource, raSource = " << decSource*kRad2Deg << ", " << raSource*kRad2Deg << endl; *fLog << " thetaTel, phiTel = " << thetaTel*kRad2Deg << ", " << phiTel*kRad2Deg << endl; *fLog << " xStar, yStar = " << xStar[0]*fMm2Deg << ", " << yStar[0]*fMm2Deg << endl; *fLog << "MSourcePosfromStarPos::SourcefromStar : printout in radians and mm" << endl; *fLog << " decStar, raStar = " << decStar[0] << ", " << raStar[0] << endl; *fLog << " decSource, raSource = " << decSource << ", " << raSource << endl; *fLog << " thetaTel, phiTel = " << thetaTel << ", " << phiTel << endl; *fLog << " xStar, yStar = " << xStar[0] << ", " << yStar[0] << endl; */ // the units are assumed to be radians for theta, phi, dec and ra // and mm for f, x and y // calculate rotation angle alpha of sky image in camera // (see TDAS 00-11, eqs. (18) and (20)) // a1 = cos(Lat), a3 = -sin(Lat), where Lat is the geographical latitude of La Palma Double_t a1 = 0.876627; Double_t a3 = -0.481171; Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) + ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) * ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) ); Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom; Double_t sinal = a1 * sin(phiTel) * denom; //*fLog << "thetaTel, phiTel, cosal, sinal = " << thetaTel << ", " // << phiTel << ", " << cosal << ", " << sinal << endl; // calculate coordinates of source in system B (see TDAS 00-11, eqs. (2)) // note that right ascension and hour angle go into opposite directions Double_t xB0 = cos(decSource) * cos(-raSource); Double_t yB0 = cos(decSource) * sin(-raSource); Double_t zB0 = -sin(decSource); //*fLog << "xB0, yB0, zB0 = " << xB0 << ", " << yB0 << ", " // << zB0 << endl; //----------------------------------------------------- // loop over stars Double_t sumx = 0.0; Double_t sumy = 0.0; Double_t sumwx = 0.0; Double_t sumwy = 0.0; for (Int_t i=0; iGetRunNumber(); *fLog << all << "MSourcePosfromStarPos::SourcefromStar; run, xSource, ySource = " << run << " : " << xSource << " +- " << dxSource << ", " << ySource << " +- " << dySource << endl; */ } // -------------------------------------------------------------------------- // // Get the source position and put it into MSrcPosCam // // Bool_t MSourcePosfromStarPos::ReInit(MParList *pList) { Int_t run = fRun->GetRunNumber(); *fLog << all << "MSourcePosfromStarPos::ReInit; run = " << run << endl; // define default values for source position and (theta, phi) Double_t xSource_def = 51.0; Double_t ySource_def = -59.0; Double_t theta_def = 25.0; Double_t phi_def = 90.0; //------------------------------------------------------------------- // search this run in the list for (Int_t i=0; iSetTelescopeTheta(thetarad); fMcEvt->SetTelescopePhi(phirad); fMcEvt->SetReadyToSave(); *fLog << all << "theta, phi = " << thetarad*kRad2Deg << ", " << phirad*kRad2Deg << " deg" << endl; //----------------------------------------- // Get source position and put it into MSrcPosCam if (fStars > 0) { TArrayD xStar(fxStar.GetNrows()); TArrayD dxStar(fdxStar.GetNrows()); TArrayD yStar(fyStar.GetNrows()); TArrayD dyStar(fdyStar.GetNrows()); for (Int_t j=0; jSetXY(fxSource, fySource); *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = " << fRunNr[i] << ", " << fxSource << " +- " << fdxSource << ", " << fySource << " +- " << fdySource << endl; fSrcPos->SetReadyToSave(); } else { // set default values fxSource = xSource_def; fySource = ySource_def; fSrcPos->SetXY(fxSource, fySource); fSrcPos->SetReadyToSave(); *fLog << warn << "MSourcePosfromStarPos::ReInit; no information on source position for run number = " << run << endl; *fLog << warn << " set xSource, ySource = " << fxSource << ", " << fySource << " mm" << endl; } return kTRUE; } } //------------------------------------------------------------------- // set default values fxSource = xSource_def; fySource = ySource_def; fSrcPos->SetXY(fxSource, fySource); fSrcPos->SetReadyToSave(); Double_t thetadeg = theta_def; Double_t thetarad = thetadeg / kRad2Deg; fMcEvt->SetTelescopeTheta(thetarad); Double_t phideg = phi_def; Double_t phirad = phideg / kRad2Deg; fMcEvt->SetTelescopePhi(phirad); fMcEvt->SetReadyToSave(); *fLog << warn << "MSourcePosfromStarPos::ReInit; no information on theta, phi and source position for run number = " << run << endl; *fLog << warn << " set xSource, ySource = " << fxSource << ", " << fySource << " mm" << endl; *fLog << warn << " set theta, phi = " << thetadeg << ", " << phideg << " deg" << endl; return kTRUE; } // -------------------------------------------------------------------------- // // Int_t MSourcePosfromStarPos::Process() { //Int_t run = fRun->GetRunNumber(); //*fLog << "MSourcePosfromStarPos::Process; run = " << run << endl; return kTRUE; } // -------------------------------------------------------------------------- // // Int_t MSourcePosfromStarPos::PostProcess() { return kTRUE; } // -------------------------------------------------------------------------- // // read the data from the ASCII file and store them // void MSourcePosfromStarPos::FixSize() { *fLog << "MSourcePosfromStarPos::FixSize; fix size of arrays : fRuns = " << fRuns << ", fStars = " << fStars << endl; fSize = fRuns; fRunNr.Set(fSize); fThetaTel.Set(fSize); fPhiTel.Set(fSize); fdThetaTel.Set(fSize); fdPhiTel.Set(fSize); Int_t fRows = fxStar.GetNrows(); fxStar.ResizeTo(fRows, fSize); fyStar.ResizeTo(fRows, fSize); fdxStar.ResizeTo(fRows, fSize); fdyStar.ResizeTo(fRows, fSize); } // -------------------------------------------------------------------------- // // read the data from the ASCII file and store them // void MSourcePosfromStarPos::ReadData() { Float_t val; Int_t ival; // extend size of arrays if necessary if ( fRuns >= (fSize-1) ) { fSize += 100; fRunNr.Set(fSize); fThetaTel.Set(fSize); fPhiTel.Set(fSize); fdThetaTel.Set(fSize); fdPhiTel.Set(fSize); fStars = fxStar.GetNrows(); fxStar.ResizeTo(fStars, fSize); fyStar.ResizeTo(fStars, fSize); fdxStar.ResizeTo(fStars, fSize); fdyStar.ResizeTo(fStars, fSize); } //------------------- // read header line //*fIn >> val; //*fLog << "val =" << val << endl; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //*fIn >> val; //------------------- fRuns += 1; *fIn >> ival; //*fLog << fRuns <<"th run : " << ival << endl; fRunNr.AddAt(ival, fRuns-1); //*fLog << "check : fRuns, fRunNr[fRuns-1], fRunNr[fRuns] = " << fRuns << ", " // << fRunNr[fRuns-1] << ", " << fRunNr[fRuns] << endl; // read mjdS, hmsS, mjdE, hmsE *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; *fIn >> val; fThetaTel.AddAt(val/kRad2Deg, fRuns-1); //*fLog << "val, fThetaTel[fRuns-1] = " << val << ", " // << fThetaTel[fRuns-1] << endl; *fIn >> val; fPhiTel.AddAt(val/kRad2Deg, fRuns-1); //*fLog << "val, fPhiTel[fRuns-1] = " << val << ", " // << fPhiTel[fRuns-1] << endl; //*fIn >> val; //fdThetaTel.AddAt(val/kRad2Deg, fRuns-1); //*fIn >> val; //fdPhiTel.AddAt(val/kRad2Deg, fRuns-1); // input is in [deg], convert to [mm] fStars = fxStar.GetNrows(); for (Int_t i=0; i> val; fxStar(i, fRuns-1) = val / fMm2Deg;; //*fLog << "val, fxStar(i, fRuns-1) = " << val << ", " // << fxStar(i, fRuns-1) << endl; *fIn >> val; fyStar(i, fRuns-1) = val / fMm2Deg; //*fLog << "val, fyStar(i, fRuns-1) = " << val << ", " // << fyStar(i, fRuns-1) << endl; *fIn >> val; //*fLog << "y=dxStar = " << val << endl; fdxStar(i, fRuns-1) = val / fMm2Deg; *fIn >> val; //*fLog << "y=dyStar = " << val << endl; fdyStar(i, fRuns-1) = val / fMm2Deg; } } // -------------------------------------------------------------------------- // // Add this file as the last entry in the chain // Int_t MSourcePosfromStarPos::AddFile(const char *txt, Int_t) { TNamed *name = new TNamed(txt, ""); fFileNames->AddLast(name); *fLog << "MSourcePosfromStarPos::AddFile; add file '" << txt << "'" << endl; return 1; } // -------------------------------------------------------------------------- // // This opens the next file in the list and deletes its name from the list. // Bool_t MSourcePosfromStarPos::OpenNextFile() { // // open the input stream and check if it is really open (file exists?) // if (fIn) delete fIn; fIn = NULL; // // Check for the existence of a next file to read // TNamed *file = (TNamed*)fFileNames->First(); if (!file) return kFALSE; // // open the file which is the first one in the chain // const char *name = file->GetName(); const char *expname = gSystem->ExpandPathName(name); fIn = new ifstream(expname); delete [] expname; const Bool_t noexist = !(*fIn); if (noexist) *fLog << dbginf << "Cannot open file '" << name << "'" << endl; else *fLog << "Open file: '" << name << "'" << endl; // // Remove this file from the list of pending files // fFileNames->Remove(file); return !noexist; } // --------------------------------------------------------------------------