/* ======================================================================== *\ ! ! * ! * 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 07/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MTelAxisFromStars // // This task // - determines the transformation from expected positions of stars // in the camera to measured positions of these stars in the camera // - applies this transformation to expected positions of other objects // to obtain the estimated positions of these objects in the camera // - puts the estimated positions into the relevant containers // // Input Containers : // MStarLocalCam[MStarLocalCam], MStarLocalCamSource[MStarLocalCam] // // Output Containers : // MSkyCamTrans, MSrcPosCam // ///////////////////////////////////////////////////////////////////////////// #include #include #include #include "MTelAxisFromStars.h" #include "MParList.h" #include "MSrcPosCam.h" #include "MLog.h" #include "MLogManip.h" #include "MStarLocalCam.h" #include "MStarLocalPos.h" #include "MSkyCamTrans.h" ClassImp(MTelAxisFromStars); using namespace std; // -------------------------------------------------------------------------- // // Constructor // MTelAxisFromStars::MTelAxisFromStars(const char *name, const char *title) { fName = name ? name : "MTelAxisFromStars"; fTitle = title ? title : "Calculate source position from star positions"; // if scale factor fLambda should NOT be fixed set fFixdScaleFactor to // -1.0; otherwise set it to the value requested fFixedScaleFactor = 1.0; // if rotation angle fAlfa should NOT be fixed set fFixdRotationAngle to // -1.0; otherwise set it to the requested value fFixedRotationAngle = 0.0; // default type of input is : the result of the correlated Gauss fit // type 0 : result from the weighted average // type 1 : result from the uncorrelated Gauss fit fInputType = 2; } // -------------------------------------------------------------------------- // // Destructor // MTelAxisFromStars::~MTelAxisFromStars() { } // -------------------------------------------------------------------------- // // Set links to containers // Int_t MTelAxisFromStars::PreProcess(MParList *pList) { fStarLocalCam = (MStarLocalCam*)pList->FindObject("MStarLocalCam", "MStarLocalCam"); if (!fStarLocalCam) { *fLog << err << "MTelAxisFromStars::PreProcess; container 'MStarLocalCam' not found... aborting." << endl; return kFALSE; } fSourceLocalCam = (MStarLocalCam*)pList->FindObject("MSourceLocalCam", "MStarLocalCam"); if (!fSourceLocalCam) { *fLog << err << "MTelAxisFromStars::PreProcess; container 'MSourceLocalCam' not found... continue " << endl; } fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam")); if (!fSrcPos) { *fLog << err << "MTelAxisFromStars::PreProcess; MSrcPosCam not found... aborting" << endl; return kFALSE; } fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans")); if (!fSkyCamTrans) { *fLog << err << "MTelAxisFromStars::PreProcess; MSkyCamTrans not found... aborting" << endl; return kFALSE; } return kTRUE; } // -------------------------------------------------------------------------- // // Set the type of the input // // type = 0 calculated star positions (by averaging) // type = 1 fitted star positions (by uncorrelated Gauss fit) // type = 2 fitted star positions (by correlated Gauss fit) // void MTelAxisFromStars::SetInputType(Int_t type) { *fLog << all << "MTelAxisFromStars::SetInputType; input type is set equal to : " << type << endl; fInputType = type; } // -------------------------------------------------------------------------- // // Fix the scale factor fLambda // // void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda) { *fLog << all << "MTelAxisFromStars::FixScaleFactorAt; scale factor will be fixed at : " << lambda << endl; fFixedScaleFactor = lambda; } // -------------------------------------------------------------------------- // // Fix rotation angle fAlfa // // void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa) { *fLog << all << "MTelAxisFromStars::FixRotationAngleAt; rotation angle will be fixed at : " << alfa << endl; fFixedRotationAngle = alfa; // [degrees] } // -------------------------------------------------------------------------- // // Process // // call FindSkyCamTrans to find the Sky-Camera transformation // call TransSkyCam to transform some sky directions // into the camera system // put the estimated source position into MSrcPosCam // Int_t MTelAxisFromStars::Process() { //Int_t run = fRun->GetRunNumber(); //*fLog << "MTelAxisFromStars::Process; run = " << run << endl; //-------------------------------------- // Define the input for FindSkyCamTrans // // get the expected (axy[0], axy[1]) and the measured positions // (bxy[0], bxy[1]) of stars in the camera from MStarLocalCam Int_t fNumStars = fStarLocalCam->GetNumStars(); if (fNumStars <= 0) return kTRUE; TArrayD axy[2]; axy[0].Set(fNumStars); axy[1].Set(fNumStars); TArrayD bxy[2]; bxy[0].Set(fNumStars); bxy[1].Set(fNumStars); // error matrix of bxy TArrayD exy[2][2]; exy[0][0].Set(fNumStars); exy[0][1].Set(fNumStars); exy[1][0].Set(fNumStars); exy[1][1].Set(fNumStars); // transformation parameters Double_t fLambda; Double_t fAlfa; Double_t fA[2][2]; Double_t fD[2]; Double_t fErrD[2][2]; Int_t fNumIter; Int_t fNdof; Double_t fChi2; Double_t fChi2Prob; MStarLocalPos *star = 0; TIter next(fStarLocalCam->GetList()); Int_t ix = 0; // loop over all stars while ( (star = (MStarLocalPos*)next()) ) { axy[0][ix] = star->GetXExp(); axy[1][ix] = star->GetYExp(); if (fInputType == 0) { // values from averaging bxy[0][ix] = star->GetMeanXCalc(); bxy[1][ix] = star->GetMeanYCalc(); exy[0][0][ix]= star->GetSigmaMinorAxisCalc()*star->GetSigmaMinorAxisCalc(); exy[0][1][ix] = 0.0; exy[1][0][ix] = 0.0; exy[1][1][ix]= star->GetSigmaMajorAxisCalc()*star->GetSigmaMajorAxisCalc(); } else if (fInputType == 1) { // values from uncorrelatd Gauss fit bxy[0][ix] = star->GetMeanXFit(); bxy[1][ix] = star->GetMeanYFit(); exy[0][0][ix]= star->GetSigmaMinorAxisFit()*star->GetSigmaMinorAxisFit(); exy[0][1][ix] = 0.0; exy[1][0][ix] = 0.0; exy[1][1][ix]= star->GetSigmaMajorAxisFit()*star->GetSigmaMajorAxisFit(); } else if (fInputType == 2) { // values from correlatd Gauss fit bxy[0][ix] = star->GetMeanXCGFit(); bxy[1][ix] = star->GetMeanYCGFit(); // this is the error matrix for (MeanXCGFit, MeanYCGFit); // this is the error matrix which should be used exy[0][0][ix] = star->GetXXErrCGFit(); exy[0][1][ix] = star->GetXYErrCGFit(); exy[1][0][ix] = star->GetXYErrCGFit(); exy[1][1][ix] = star->GetYYErrCGFit(); // this is the error matrix constructed from SigmaXCGFit and SigmaYCGFit; // it is used because the errors above are too small, at present //exy[0][0][ix] = star->GetSigmaXCGFit() * star->GetSigmaXCGFit(); //exy[0][1][ix] = star->GetCorrXYCGFit() * // star->GetSigmaXCGFit() * star->GetSigmaYCGFit(); //exy[1][0][ix] = exy[0][1][ix]; //exy[1][1][ix] = star->GetSigmaYCGFit() * star->GetSigmaYCGFit(); } else { *fLog << err << "MTelAxisFromStars::Process; type of input is not defined" << endl; return kFALSE; } // don't include stars with undefined error Double_t deter = exy[0][0][ix]*exy[1][1][ix] - exy[0][1][ix]*exy[1][0][ix]; //*fLog << "ix ,deter, xx, xy, yy = " << ix << ": " // << deter << ", " << exy[0][0][ix] << ", " // << exy[0][1][ix] << ", " << exy[1][1][ix] << endl; if (deter <= 0.0) continue; //*fLog << "MTelAxisFromStars : " << endl; //*fLog << " ix, XExp, YExp, XFit, YFit, SigmaX2, SigmaXY, SigmaY2 = " // << ix << " : " // << axy[0][ix] << ", " << axy[1][ix] << ", " // << bxy[0][ix] << ", " << bxy[1][ix] << ", " // << exy[0][0][ix] << ", " << exy[0][1][ix] << ", " // << exy[1][1][ix] << endl; ix++; } //-------------------------------------- // Find the transformation from expected positions (axy[1], axy[2]) // to measured positions (bxy[1], bxy[2]) in the camera Int_t fNStars = ix; if (ix < fNumStars) { // reset the sizes of the arrays Int_t fNStars = ix; axy[0].Set(fNStars); axy[1].Set(fNStars); bxy[0].Set(fNStars); bxy[1].Set(fNStars); exy[0][0].Set(fNStars); exy[0][1].Set(fNStars); exy[1][0].Set(fNStars); exy[1][1].Set(fNStars); } Bool_t fitOK; if (fNStars < 1) { *fLog << "MTelAxisFromStars::Process; no star for MTelAxisFromStars" << endl; fitOK = kFALSE; } else { fitOK = FindSkyCamTrans(axy, bxy, exy, fFixedRotationAngle, fFixedScaleFactor, fLambda, fAlfa , fA, fD, fErrD, fNumIter, fNdof, fChi2, fChi2Prob); } if (!fitOK) { *fLog << err << "MTelAxisFromStars::Process; Fit to find transformation from star to camera system failed" << endl; if (fNStars > 0) { *fLog << err << " fNumIter, fNdof, fChi2, fChi2Prob = " << fNumIter << ", " << fNdof << ", " << fChi2 << ", " << fChi2Prob << endl; } return kTRUE; } //-------------------------------------- // Put the transformation parameters into the MSkyCamTrans container fSkyCamTrans->SetParameters(fLambda, fAlfa, fA, fD, fErrD, fNumStars, fNumIter, fNdof, fChi2, fChi2Prob); fSkyCamTrans->SetReadyToSave(); //-------------------------------------- // Put the estimated position, obtained by transforming the expected // position (0, 0), into SrcPosCam fSrcPos->SetXY(fD[0], fD[1]); fSrcPos->SetReadyToSave(); //-------------------------------------- // Apply the transformation to some expected positions (asxy[1], asxy[2]) // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera // and their error matrices esxy[2][2] // get the expected positions (asxy[1], asxy[2]) from another MStarLocalCam // container (with the name "MSourceLocalCam") Int_t fNumStarsSource = fSourceLocalCam->GetNumStars(); //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = " // << fNumStarsSource << endl; if (fNumStarsSource > 0) { TArrayD asxy[2]; asxy[0].Set(fNumStarsSource); asxy[1].Set(fNumStarsSource); TArrayD bsxy[2]; bsxy[0].Set(fNumStarsSource); bsxy[1].Set(fNumStarsSource); TArrayD esxy[2][2]; esxy[0][0].Set(fNumStarsSource); esxy[0][1].Set(fNumStarsSource); esxy[1][0].Set(fNumStarsSource); esxy[1][1].Set(fNumStarsSource); MStarLocalPos *starSource = 0; TIter nextSource(fSourceLocalCam->GetList()); ix = 0; while ( (starSource = (MStarLocalPos*)nextSource()) ) { asxy[0][ix] = starSource->GetXExp(); asxy[1][ix] = starSource->GetYExp(); ix++; } TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy); // put the estimated positions into the MStarLocalPos container TIter setnextSource(fSourceLocalCam->GetList()); ix = 0; while ( (starSource = (MStarLocalPos*)setnextSource()) ) { if (fInputType == 1) { starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), fChi2, fNdof); } else if (fInputType == 2) { Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] ); starSource->SetCGFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr, esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix], fChi2, fNdof); } ix++; } } //-------------------------------------- return kTRUE; } //--------------------------------------------------------------------------- // // FindSkyCamTrans // // This routine determines the transformation // // ( cos(alfa) -sin(alfa) ) // b = lambda * A * a + d A = ( ) // ^ ^ ^ ( sin(alfa) cos(alfa) ) // | | | // scale rotation shift // factor matrix // // from sky coordinates 'a' (projected onto the camera) to camera // coordinates 'b', using the positions of known stars in the camera. // The latter positions may have been determined by analysing the // DC currents in the different pixels. // // Input : a[2] x and y coordinates of stars projected onto the camera; // they were obtained from (RA, dec) of the stars and // (ThetaTel, PhiTel) and the time of observation; // these are the 'expected positions' of stars in the camera // b[2] 'measured positions' of these stars in the camera; // they may have been obtained from the DC currents // e[2][2] error matrix of b[2] // fixedrotationangle value [in degrees] at which rotation angle // alfa should be fixed; -1 means don't fix // fixedscalefactor value at which scale factor lambda // should be fixed; -1 means don't fix // // Output : lambda, alfadeg, A[2][2], d[2] fit results; // parameters describing the transformation // from 'expected positions' to the 'measured // positions' in the camera // errd[2][2] error matrix of d[2] // fNumIter number of iterations // fNdoF number of degrees of freedom // fChi2 chi-square value // fChi2Prob chi-square probability // // The units are assumed to be // [degrees] for alfadeg // [mm] for a, b, d // [1] for lambda Bool_t MTelAxisFromStars::FindSkyCamTrans( TArrayD a[2], TArrayD b[2], TArrayD e[2][2], Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda, Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2], Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob) { Int_t fNumStars = a[0].GetSize(); //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :" // << endl; for (Int_t ix=0; ix 1.1 ) lambda = 1.1; dlambda = lambda - lambdaold; lambdaold = lambda; //*fLog << "num, denom, lambda, dlambda = " << num // << ", " << denom << ", " << lambda << ", " // << dlambda << endl; } } //------- end of iteration ***************************************** alfadeg = alfa * kRad2Deg; // calculate error matrix of d[2] errd[0][0] = sumcinv[0][0]; errd[0][1] = sumcinv[0][1]; errd[1][0] = sumcinv[1][0]; errd[1][1] = sumcinv[1][1]; // evaluate quality of fit // calculate chi2 for (Int_t ix=0; ix