Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4740) +++ /trunk/MagicSoft/Mars/Changelog (revision 4741) @@ -19,4 +19,12 @@ -*-*- END OF LINE -*-*- + 2004/08/26 : Wolfgang Wittek + + * mstarcam/MStarCamTrans.[h,cc] + - include optical aberration when plotting the grid + + * mtemp/MTelAxisFromStars.[h,cc] + - add member functions SetPointingPosition + and SetSourcePosition Index: /trunk/MagicSoft/Mars/mstarcam/MStarCamTrans.cc =================================================================== --- /trunk/MagicSoft/Mars/mstarcam/MStarCamTrans.cc (revision 4740) +++ /trunk/MagicSoft/Mars/mstarcam/MStarCamTrans.cc (revision 4741) @@ -829,4 +829,12 @@ //gLog << "mmtodeg = " << mmtodeg << endl; + // aberr is the ratio r_optaberr/r_ideal between + // the distance from the camera center with optical aberration and + // the distance from the camera center with an ideal imaging + // the value 1.07 is valid if the expected position (with aberration) + // in the camera is calculated as the average of the positions obtained + // from the reflections at the individual mirrors + Double_t aberr = 1.07; + //-------------------------------------------------------------------- @@ -937,6 +945,6 @@ Loc0LocToCam(theta0, phi0, theta[j], phi[k], xx, yy); - x[k] = xx * mmtodeg; - y[k] = yy * mmtodeg; + x[k] = xx * mmtodeg * aberr; + y[k] = yy * mmtodeg * aberr; //gLog << "theta0, phi0 = " << theta0 << ", " << phi0 @@ -1003,6 +1011,6 @@ Loc0LocToCam(theta0, phi0, theta1[j], phi1[k], xx, yy); - x1[count] = xx * mmtodeg; - y1[count] = yy * mmtodeg; + x1[count] = xx * mmtodeg * aberr; + y1[count] = yy * mmtodeg * aberr; //gLog << "theta0, phi0 = " << theta0 << ", " << phi0 @@ -1103,6 +1111,6 @@ Cel0CelToCam(dec0, hh0, dec[j], hx, xx, yy); - xh[k] = xx * mmtodeg; - yh[k] = yy * mmtodeg; + xh[k] = xx * mmtodeg * aberr; + yh[k] = yy * mmtodeg * aberr; //gLog << "dec0, h0 = " << dec0 << ", " << h0 @@ -1172,6 +1180,6 @@ Cel0CelToCam(dec0, hh0, dec1[j], hhx, xx, yy); - xd[count] = xx * mmtodeg; - yd[count] = yy * mmtodeg; + xd[count] = xx * mmtodeg * aberr; + yd[count] = yy * mmtodeg * aberr; //gLog << "dec0, h0 = " << dec0 << ", " << h0 @@ -1224,2 +1232,3 @@ return kTRUE; } + Index: /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc (revision 4740) +++ /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc (revision 4741) @@ -48,12 +48,17 @@ #include "MParList.h" -#include "MSrcPosCam.h" + #include "MLog.h" #include "MLogManip.h" + +#include "MReportDrive.h" +#include "MPointingPos.h" +#include "MSrcPosCam.h" #include "MStarCam.h" #include "MStarPos.h" #include "MSkyCamTrans.h" +#include "MStarCamTrans.h" ClassImp(MTelAxisFromStars); @@ -82,4 +87,11 @@ // type 1 : result from the Gauss fit fInputType = 1; + + // default value of fAberr + // the value 1.07 is valid if the expected position (with aberration) + // in the camera is calculated as the average of the positions obtained + // from the reflections at the individual mirrors + fAberr = 1.07; + } @@ -90,4 +102,5 @@ MTelAxisFromStars::~MTelAxisFromStars() { + delete fStarCamTrans; } @@ -99,4 +112,12 @@ Int_t MTelAxisFromStars::PreProcess(MParList *pList) { + fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive")); + if (!fDrive) + { + *fLog << err << AddSerialNumber("MReportDrive") + << " not found... aborting." << endl; + return kFALSE; + } + fStarCam = (MStarCam*)pList->FindObject("MStarCam", "MStarCam"); @@ -115,4 +136,5 @@ + //----------------------------------------------------------------- fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam")); if (!fSrcPos) @@ -122,4 +144,18 @@ } + fPointPos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MPointingPos"); + if (!fPointPos) + { + *fLog << err << "MTelAxisFromStars::PreProcess; MPointingPos not found... aborting" << endl; + return kFALSE; + } + + fSourcePos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MSourcePos"); + if (!fSourcePos) + { + *fLog << err << "MTelAxisFromStars::PreProcess; MSourcePos[MPointingPos] not found... aborting" << endl; + return kFALSE; + } + fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans")); if (!fSkyCamTrans) @@ -129,6 +165,45 @@ } + //----------------------------------------------------------------- + // book an MStarCamTrans object + // this is needed when calling one of the member functions of MStarCamTrans + + MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!geom) + { + *fLog << err << AddSerialNumber("MGeomCam") + << " not found... aborting." << endl; + return kFALSE; + } + + MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory")); + if (!obs) + { + *fLog << err << AddSerialNumber("MObservatory") + << " not found... aborting." << endl; + return kFALSE; + } + + fStarCamTrans = new MStarCamTrans(*geom, *obs); return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set optical aberration factor +// +// fAberr is the ratio between +// the distance from the camera center with optical aberration and +// the distance from the camera center with an ideal imaging +// +// fAberr = r_real/r_ideal +// +void MTelAxisFromStars::SetOpticalAberr(Double_t aberr) +{ + *fLog << all << "MTelAxisFromStars::SetOpticalAberr; the optical aberration factor is set equal to : " + << aberr ; + + fAberr = aberr; } @@ -374,12 +449,29 @@ //-------------------------------------- - // Put the estimated position, obtained by transforming the expected - // position (0,0), into SrcPosCam + // Put the corrected pointing position into MPointingPos + // + SetPointingPosition(fStarCamTrans, fDrive, fSkyCamTrans, fPointPos); + + + //-------------------------------------- + // Put the estimated position of the source into SrcPosCam + // + // get the source direction from MReportDrive // Note : this has to be changed for the wobble mode, where the source - // won't be in the center of the camera - - fSrcPos->SetXY(fD[0], fD[1]); - fSrcPos->SetReadyToSave(); - + // isn't in the center of the camera + Double_t decsource = fDrive->GetDec(); + Double_t rasource = fDrive->GetRa(); + // + Double_t radrive = fDrive->GetRa(); + Double_t hdrive = fDrive->GetHa(); + Double_t hsource = hdrive + radrive - rasource; + fSourcePos->SetSkyPosition(rasource, decsource, hsource); + + SetSourcePosition(fStarCamTrans, fPointPos, fSourcePos, fSrcPos); + //fSrcPos->SetXY(fD[0], fD[1]); + //fSrcPos->SetReadyToSave(); + *fLog << "after calling SetSourcePosition : , X, Y ,fD = " + << fSrcPos->GetX() << ", " << fSrcPos->GetY() << ", " + << fD[0] << ", " << fD[1] << endl; //-------------------------------------- @@ -448,4 +540,6 @@ return kTRUE; } + + //--------------------------------------------------------------------------- @@ -1045,4 +1139,105 @@ // -------------------------------------------------------------------------- // +// SetPointingPosition +// +// put the corrected pointing direction into MPointingPos[MPointingPos]; +// this direction corresponds to the position (0,0) in the camera +// + +Bool_t MTelAxisFromStars::SetPointingPosition(MStarCamTrans *fstarcamtrans, + MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos) +{ + Double_t decdrive = fdrive->GetDec(); + Double_t hdrive = fdrive->GetHa(); + Double_t radrive = fdrive->GetRa(); + + // this is the estimated position (with optical aberration) in the camera + // corresponding to the direction in MReportDrive + Double_t Xpoint = (ftrans->GetShiftD())[0]; + Double_t Ypoint = (ftrans->GetShiftD())[1]; + + // get the sky direction corresponding to the position (0,0) in the camera + Double_t decpoint = 0.0; + Double_t hpoint = 0.0; + fstarcamtrans->CelCamToCel0(decdrive, hdrive, + Xpoint/fAberr, Ypoint/fAberr, decpoint, hpoint); + Double_t rapoint = radrive - hpoint + hdrive; + fpointpos->SetSkyPosition(rapoint, decpoint, hpoint); + + // get the local direction corresponding to the position (0,0) in the camera + Double_t thetadrive = fdrive->GetNominalZd(); + Double_t phidrive = fdrive->GetNominalAz(); + Double_t thetapoint = 0.0; + Double_t phipoint = 0.0; + fstarcamtrans->LocCamToLoc0(thetadrive, phidrive, + Xpoint/fAberr, Ypoint/fAberr, thetapoint, phipoint); + fpointpos->SetLocalPosition(thetapoint, phipoint); + fpointpos->SetReadyToSave(); + + *fLog << "SetPointingPosition : decdrive, hdrive, radrive Xpoint, Ypoint = " + << decdrive << ", " << hdrive << ", " << radrive << ", " + << Xpoint << ", " << Ypoint << endl; + + *fLog << "SetPointingPosition : thetadrive, phidrive, thetapoint, phipoint = " + << thetadrive << ", " << phidrive << ", " << thetapoint << ", " + << phipoint << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// SetSourcePosition +// +// put the estimated position of the source in the camera into +// MSrcPosCam[MSrcPosCam] +// +// and the estimated local direction of the source into +// MSourcePos[MPointingPos] +// + +Bool_t MTelAxisFromStars::SetSourcePosition(MStarCamTrans *fstarcamtrans, + MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos) +{ + // get the corrected pointing direction + // corresponding to the position (0,0) in the camera + Double_t decpoint = fpointpos->GetDec(); + Double_t hpoint = fpointpos->GetHa(); + + // get the sky direction of the source + Double_t decsource = fsourcepos->GetDec(); + Double_t hsource = fsourcepos->GetHa(); + + // get the estimated position (Xsource, Ysource) of the source in the camera; + // this is a position for an ideal imaging, without optical aberration + Double_t Xsource = 0.0; + Double_t Ysource = 0.0; + fstarcamtrans->Cel0CelToCam(decpoint, hpoint, + decsource, hsource, Xsource, Ysource); + fsrcpos->SetXY(Xsource*fAberr, Ysource*fAberr); + fsrcpos->SetReadyToSave(); + + // get the estimated local direction of the source + Double_t thetapoint = fpointpos->GetZd(); + Double_t phipoint = fpointpos->GetAz(); + Double_t thetasource = 0.0; + Double_t phisource = 0.0; + fstarcamtrans->Loc0CamToLoc(thetapoint, phipoint, + Xsource, Ysource, thetasource, phisource); + fsourcepos->SetLocalPosition(thetasource, phisource); + fsourcepos->SetReadyToSave(); + + *fLog << "SetSourcePosition : decpoint, hpoint, decsource, hsource, Xsource, Ysource = " + << decpoint << ", " << hpoint << ", " << decsource << ", " + << hsource << ", " << Xsource << ", " << Ysource << endl; + *fLog << "SetSourcePosition : thetapoint, phipoint, thetasource, phisource = " + << thetapoint << ", " << phipoint << ", " << thetasource << ", " + << phisource << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// // Int_t MTelAxisFromStars::PostProcess() Index: /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.h (revision 4740) +++ /trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.h (revision 4741) @@ -20,7 +20,10 @@ #endif +class MReportDrive; +class MPointingPos; class MSrcPosCam; class MStarCam; class MSkyCamTrans; +class MStarCamTrans; class MTelAxisFromStars : public MTask @@ -31,9 +34,14 @@ const MStarCam *fSourceCam; //! + MReportDrive *fDrive; //! + MPointingPos *fPointPos; //! + MPointingPos *fSourcePos; //! MSrcPosCam *fSrcPos; //! MSkyCamTrans *fSkyCamTrans; //! + MStarCamTrans *fStarCamTrans; //! Double_t fFixedScaleFactor; //! Double_t fFixedRotationAngle; //! [degrees] + Double_t fAberr; //! optical aberration factor Int_t fInputType; //! type of input @@ -50,6 +58,7 @@ void FixScaleFactorAt(Double_t lambda = 1.0); void FixRotationAngleAt(Double_t alfa = 0.0); // alfa in [degrees] - - void SetInputType(Int_t type = 2); + + void SetOpticalAberr(Double_t aberr); + void SetInputType(Int_t type = 1); Bool_t FindSkyCamTrans(TArrayD[2], TArrayD[2], TArrayD[2][2], @@ -60,4 +69,9 @@ void TransSkyCam(Double_t &, Double_t[2][2], Double_t[2], Double_t[2][2], TArrayD[2], TArrayD[2], TArrayD[2][2]); + + Bool_t SetPointingPosition(MStarCamTrans *fstarcamtrans, + MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos); + Bool_t SetSourcePosition(MStarCamTrans *fstarcamtrans, + MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos); ClassDef(MTelAxisFromStars, 0) // Task to calculate the source position from star positions Index: /trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C (revision 4740) +++ /trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C (revision 4741) @@ -15,8 +15,34 @@ } +//-------------------------------------------------------------------------- +// +// findTelAxisFromStars +// +// macro to +// - determine the positions of stars in the camera from the DC currents +// - and to find the expected position of the source in the camera +// + void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0) { gLog.SetNoColors(); + //------------------------------------------------------------------- + // book some histograms which are to be filled in the event loop + // + + TH2D *aberr = new TH2D("OpticalAberr", "Opt_Aberr", + 100, 0.0, 650.0, 100, 0.0, 650.0); + aberr->SetXTitle("Ideal distance from center [mm]"); + aberr->SetYTitle("Expected distance from center [mm]"); + + TProfile *reldiff = new TProfile("RelativeDiff", "Rel_Diff", + 100, 0.0, 650.0, -0.2, 0.2, "S"); + reldiff->SetXTitle("Ideal distance from center [mm]"); + reldiff->SetYTitle("(R - R0)/R0"); + + //------------------------------------------------------------------- + + MParList plist; MTaskList tlist; @@ -28,8 +54,11 @@ //$$$$$$$$$$$$$$$$ ww + MObservatory obs; MStarCam sourcecam; sourcecam.SetName("MSourceCam"); MHTelAxisFromStars htelaxis; + + plist.AddToList(&obs); plist.AddToList(&htelaxis); //$$$$$$$$$$$$$$$$ ww @@ -99,4 +128,61 @@ const Double_t dec = MAstro::Dms2Rad(38, 14, 29); sdirs.AddDirection(ra,dec,1,"My_UMa 51"); + + // dummy source dirtections + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 20, 0); + sdirs.AddDirection(ra,dec,1,"Dummy1"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 20, 0); + sdirs.AddDirection(ra,dec,1,"Dummy2"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 30, 0); + sdirs.AddDirection(ra,dec,1,"Dummy3"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 40, 0); + sdirs.AddDirection(ra,dec,1,"Dummy4"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 50, 0); + sdirs.AddDirection(ra,dec,1,"Dummy5"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 0, 0); + sdirs.AddDirection(ra,dec,1,"Dummy6"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 10, 0); + sdirs.AddDirection(ra,dec,1,"Dummy7"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 20, 0); + sdirs.AddDirection(ra,dec,1,"Dummy8"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 30, 0); + sdirs.AddDirection(ra,dec,1,"Dummy9"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 40, 0); + sdirs.AddDirection(ra,dec,1,"Dummy10"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(39, 50, 0); + sdirs.AddDirection(ra,dec,1,"Dummy11"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(40, 0, 0); + sdirs.AddDirection(ra,dec,1,"Dummy12"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(40, 10, 0); + sdirs.AddDirection(ra,dec,1,"Dummy13"); + + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(40, 20, 0); + sdirs.AddDirection(ra,dec,1,"Dummy14"); @@ -195,24 +281,44 @@ // This line prints the results: // starcam.Print(); + // This is how to access the TList of stars: - // TList* starlist = starcam.GetList(); + // TList* starlist = starcam.GetList(); // This is how to iterate over stars found: - // TIter Next(starlist); + // TIter Next(starlist); // MStarPos* star; // UInt_t starnum = 0; - // cout << filename << " "; - // cout << "Iterating over list" << endl; - // while ((star=(MStarPos*)Next())) - // { - // cout << "star[" << starnum << "] "; - // cout << star->GetMeanX() << " " - // << star->GetMeanY() << " "; - // starnum++; - // } - // cout << endl; + // cout << filename << " "; + // cout << "Iterating over list" << endl; + // while ( (star=(MStarPos*)Next()) ) + // { + // cout << "star[" << starnum << "] "; + // cout << star->GetMeanX() << " " + // << star->GetMeanY() << " "; + // starnum++; + // } + // cout << endl; + } + + // fill the histograms + // loop over stars in "MSourceCam" + TList* sourcelist = sourcecam.GetList(); + TIter Nextsource(sourcelist); + MStarPos* source; + while ( (source=(MStarPos*)Nextsource()) ) + { + Double_t R = sqrt( source->GetXExp()*source->GetXExp() + + source->GetYExp()*source->GetYExp() ); + + Double_t R0= sqrt( source->GetXIdeal()*source->GetXIdeal() + + source->GetYIdeal()*source->GetYIdeal() ); + + cout << "R0, R = " << R0 << ", " << R << endl; + + aberr->Fill(R0, R, 1.0); + reldiff->Fill(R0, (R-R0)/R0), 1.0); + } - }//integratedevents - + //-------------------------------------- MTime time; time.Set(2004, 4, 22, 21, 51, 15); @@ -225,6 +331,6 @@ // wait after each event - if (!HandleInput()) - break; + //if (!HandleInput()) + // break; } @@ -256,8 +362,21 @@ //$$$$$$$$$$$$$$$$ ww - + // plot the histograms + TCanvas *canv = new TCanvas("Optical_Aberration", "Optical aberration", + 0, 0, 600, 300); + canv->Divide(2,1); + canv->SetBorderMode(0); + + canv->cd(1); + aberr->Draw(); + + canv->cd(2); + reldiff->Draw(); } +//---------------------------------------------------------------------- +// +// Bool_t HandleInput() @@ -283,2 +402,12 @@ return kFALSE; } +//========================================================================= + + + + + + + + +