Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8915) +++ trunk/MagicSoft/Mars/Changelog (revision 8916) @@ -37,5 +37,6 @@ mimage/MHillas.[h,cc], mmuon/MMuonSearchPar..[h,cc], mpedestal/MPedestalPix.[h,cc], mpointing/MPointingDev.[h,cc], - mpointing/MSrcPosCam.[h,cc]: + mpointing/MSrcPosCam.[h,cc], mpointing/MPointingPos.[h,cc], + mpointing/MPointing.[h,cc]: - moved some code to source file to prevent TMath inclusion in header (root 5.18) @@ -90,7 +91,23 @@ Otherwise it is ambiguous in root 5.18 + * mhbase/MH.cc: + - adde missing includes of TColor, TMath and TClass (root 5.18) + - implemented a workaround which always uses the correct + CreateGradientColorTable (root 5.18) + * Makefile: - linking of the shared object is now done in /tmp - replaced = by := where possible + + * mjobs/MJCalibrateSignal.cc: + - do not invert contcoscal, that's wrong + + * mmovie/MMovieWrite.cc: + - added a #if-directive to use either gStyle or TColor + for CreateGradientColorTable depending on root-version + + * mimage/MStereoPar.[h,cc], mimage/MStereoCal.[h,cc]: + - replaced Monate Carlo container by MPointingPos + - made every algorithm unique Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8915) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8916) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.37 2008-05-19 14:04:12 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.38 2008-06-02 17:21:24 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -20,5 +20,5 @@ ! Author(s): Thomas Bretz 07/2001 ! -! Copyright: MAGIC Software Development, 2000-2007 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -57,4 +57,7 @@ #include #include +#include +#include +#include #include // TStyle::GetScreenFactor #include @@ -63,4 +66,5 @@ #include #include +#include #if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01) #include @@ -1596,4 +1600,14 @@ } +#ifdef CreateGradientColorTable +#error CreateGradientColorTable already defined +#endif + +#if ROOT_VERSION_CODE < ROOT_VERSION(5,18,00) +#define CreateGradientColorTable gStyle->CreateGradientColorTable +#else +#define CreateGradientColorTable TColor::CreateGradientColorTable +#endif + // -------------------------------------------------------------------------- // @@ -1639,5 +1653,5 @@ Double_t g[5] = { 0.01, 0.02, 0.39, 0.68, 0.97 }; Double_t b[5] = { 0.17, 0.39, 0.62, 0.79, 0.97 }; - gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + CreateGradientColorTable(5, s, r, g, b, ncol); found=kTRUE; } @@ -1649,5 +1663,5 @@ Double_t g[5] = { 0.00, 0.02, 0.40, 0.70, 1.00 }; Double_t b[5] = { 0.00, 0.27, 0.51, 0.81, 1.00 }; - gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + CreateGradientColorTable(5, s, r, g, b, ncol); found=kTRUE; } @@ -1659,5 +1673,5 @@ double g[2] = {0.00, 1.00}; double b[2] = {0.00, 1.00}; - gStyle->CreateGradientColorTable(2, s, r, g, b, ncol); + CreateGradientColorTable(2, s, r, g, b, ncol); found=kTRUE; } @@ -1669,5 +1683,5 @@ double g[5] = {0., 0.10, 0.20, 0.73, 1.00}; double b[5] = {0., 0.03, 0.06, 0.00, 1.00}; - gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + CreateGradientColorTable(5, s, r, g, b, ncol); found=kTRUE; } @@ -1679,5 +1693,5 @@ double g[4] = {0.03, 0.04, 0.80, 1.00}; double b[4] = {0.03, 0.04, 0.00, 1.00}; - gStyle->CreateGradientColorTable(4, s, r, g, b, ncol); + CreateGradientColorTable(4, s, r, g, b, ncol); found=kTRUE; } @@ -1689,5 +1703,5 @@ double g[8] = {0.70, 0.40, 0.02, 0.00, 0.10, 0.20, 0.73, 1.00}; double b[8] = {0.81, 0.51, 0.27, 0.00, 0.03, 0.06, 0.00, 1.00}; - gStyle->CreateGradientColorTable(8, s, r, g, b, ncol); + CreateGradientColorTable(8, s, r, g, b, ncol); found=kTRUE; } @@ -1699,5 +1713,5 @@ double g[3] = {0., 0.0, 1.}; double b[3] = {0., 0.0, 1.}; - gStyle->CreateGradientColorTable(3, s, r, g, b, ncol); + CreateGradientColorTable(3, s, r, g, b, ncol); found=kTRUE; } @@ -1709,5 +1723,5 @@ double g[3] = {0., 0.0, 1.}; double b[3] = {0., 1.0, 1.}; - gStyle->CreateGradientColorTable(3, s, r, g, b, ncol); + CreateGradientColorTable(3, s, r, g, b, ncol); found=kTRUE; } @@ -1719,5 +1733,5 @@ double g[4] = {0.28, 0.93, 0.03, 1.}; double b[4] = {0.79, 0.11, 0.03, 1.}; - gStyle->CreateGradientColorTable(4, s, r, g, b, ncol); + CreateGradientColorTable(4, s, r, g, b, ncol); found=kTRUE; } @@ -1728,5 +1742,5 @@ double g[5] = {0.1, 0.1, 0.0, 1.0, 1.0 }; double b[5] = {0.2, 0.7, 0.0, 0.3, 0.9 }; - gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + CreateGradientColorTable(5, s, r, g, b, ncol); found=kTRUE; } Index: trunk/MagicSoft/Mars/mimage/MStereoCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoCalc.cc (revision 8915) +++ trunk/MagicSoft/Mars/mimage/MStereoCalc.cc (revision 8916) @@ -18,5 +18,5 @@ ! Author(s): Abelardo Moralejo, 11/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -33,5 +33,5 @@ // MGeomCam // MHillas -// MMcEvt +// MPointingPos // // Output Containers: @@ -44,6 +44,6 @@ #include "MHillas.h" -#include "MMcEvt.hxx" #include "MStereoPar.h" +#include "MPointingPos.h" #include "MLog.h" @@ -63,5 +63,4 @@ fName = name ? name : "MStereoCalc"; fTitle = title ? title : "Calculate shower parameters in stereo mode"; - } @@ -74,67 +73,49 @@ Int_t MStereoCalc::PreProcess(MParList *pList) { - // necessary - - fmcevt1 = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt",fCT1id)); - if (!fmcevt1) + fPointingPos1 = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos",fCT1id)); + if (!fPointingPos1) { - *fLog << err << AddSerialNumber("MMcEvt",fCT1id) << " not found... aborting." << endl; + *fLog << err << AddSerialNumber("MPointingPos",fCT1id) << " not found... aborting." << endl; return kFALSE; } - // necessary - fmcevt2 = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt",fCT2id)); - if (!fmcevt2) + fPointingPos2 = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos",fCT2id)); + if (!fPointingPos2) { - *fLog << err << AddSerialNumber("MMcEvt",fCT2id) << " not found... aborting." << endl; + *fLog << err << AddSerialNumber("MPointingPos",fCT2id) << " not found... aborting." << endl; return kFALSE; } - // necessary - TString geomname = "MGeomCam;"; - geomname += fCT1id; - fGeomCam1 = (MGeomCam*)pList->FindObject(geomname); + fGeomCam1 = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam", fCT1id)); if (!fGeomCam1) { - *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl; + *fLog << err << AddSerialNumber("MGeomCam", fCT1id) << " not found... aborting." << endl; return kFALSE; } - // necessary - geomname = "MGeomCam;"; - geomname += fCT2id; - fGeomCam2 = (MGeomCam*)pList->FindObject(geomname); + fGeomCam2 = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam", fCT2id)); if (!fGeomCam2) { - *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl; + *fLog << err << AddSerialNumber("MGeomCam", fCT2id) << " not found... aborting." << endl; return kFALSE; } - // necessary - TString hillasname = "MHillas;"; - hillasname += fCT1id; - fHillas1 = (MHillas*)pList->FindObject(hillasname); + fHillas1 = (MHillas*)pList->FindObject(AddSerialNumber("MHillas", fCT1id)); if (!fHillas1) { - *fLog << err << hillasname << " missing in Parameter List... aborting." << endl; + *fLog << err << AddSerialNumber("MHillas", fCT1id) << " not found... aborting." << endl; return kFALSE; } - // necessary - hillasname = "MHillas;"; - hillasname += fCT2id; - fHillas2 = (MHillas*)pList->FindObject(hillasname); + fHillas2 = (MHillas*)pList->FindObject(AddSerialNumber("MHillas", fCT2id)); if (!fHillas2) { - *fLog << err << hillasname << " missing in Parameter List... aborting." << endl; + *fLog << err << AddSerialNumber("MHillas", fCT2id) << " not found... aborting." << endl; return kFALSE; } - fStereoPar = (MStereoPar*)pList->FindCreateObj("MStereoPar"); + fStereoPar = (MStereoPar*)pList->FindCreateObj("MStereoPar", fStereoParName); if (!fStereoPar) - { - *fLog << err << "Could not create MStereoPar... aborting" << endl; return kFALSE; - } return kTRUE; @@ -148,15 +129,7 @@ Int_t MStereoCalc::Process() { - fStereoPar->Calc(*fHillas1, *fmcevt1, *fGeomCam1, fCT1x, fCT1y, *fHillas2, *fmcevt2, *fGeomCam2, fCT2x, fCT2y); + fStereoPar->Calc(*fHillas1, *fPointingPos1, *fGeomCam1, fCT1x, fCT1y, + *fHillas2, *fPointingPos2, *fGeomCam2, fCT2x, fCT2y); return kTRUE; } - -// -------------------------------------------------------------------------- -// -// Does nothing at the moment. -// -Int_t MStereoCalc::PostProcess() -{ - return kTRUE; -} Index: trunk/MagicSoft/Mars/mimage/MStereoCalc.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoCalc.h (revision 8915) +++ trunk/MagicSoft/Mars/mimage/MStereoCalc.h (revision 8916) @@ -13,29 +13,26 @@ #include "MTask.h" #endif -#ifndef ROOT_TArrayL -#include -#endif class MGeomCam; class MHillas; -class MMcEvt; class MStereoPar; +class MPointingPos; class MStereoCalc : public MTask { - const MGeomCam *fGeomCam1; //! Camera Geometry CT1 - const MHillas *fHillas1; //! input - const MMcEvt *fmcevt1; //! input + const MGeomCam *fGeomCam1; //! CT1: Camera Geometry + const MHillas *fHillas1; //! CT1: Hillas parameters + const MPointingPos *fPointingPos1; //! CT1: Pointing Direction - const MGeomCam *fGeomCam2; //! Camera Geometry CT2 - const MHillas *fHillas2; //! input - const MMcEvt *fmcevt2; //! input + const MGeomCam *fGeomCam2; //! CT2: Camera Geometry + const MHillas *fHillas2; //! CT2: Hillas parameters + const MPointingPos *fPointingPos2; //! CT2: pointing direction - Int_t fCT1id; //! - Int_t fCT2id; //! Identifiers of the two analyzed telescopes. + Int_t fCT1id; // CT1: Identifier number + Int_t fCT2id; // CT2: Identifier number - Float_t fCT1x; //! + Float_t fCT1x; //! FIXME -> Move to parameter list Float_t fCT1y; //! Position of first telescope - Float_t fCT2x; //! + Float_t fCT2x; //! FIXME -> Move to parameter list Float_t fCT2y; //! Position of second telescope @@ -45,6 +42,4 @@ Int_t PreProcess(MParList *pList); Int_t Process(); - Int_t PostProcess(); - public: Index: trunk/MagicSoft/Mars/mimage/MStereoPar.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoPar.cc (revision 8915) +++ trunk/MagicSoft/Mars/mimage/MStereoPar.cc (revision 8916) @@ -18,5 +18,5 @@ ! Author(s): Abelardo Moralejo 11/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2008 ! ! @@ -33,5 +33,8 @@ ///////////////////////////////////////////////////////////////////////////// #include "MStereoPar.h" + #include + +#include #include "MLog.h" @@ -39,7 +42,6 @@ #include "MHillas.h" -#include "MMcEvt.hxx" #include "MGeomCam.h" - +#include "MPointingPos.h" ClassImp(MStereoPar); @@ -55,8 +57,5 @@ fName = name ? name : "MStereoPar"; fTitle = title ? title : "Stereo image parameters"; - - -} - +} // -------------------------------------------------------------------------- @@ -70,21 +69,68 @@ } - -// -------------------------------------------------------------------------- -// -// Calculation of shower parameters -// -void MStereoPar::Calc(const MHillas &hillas1, const MMcEvt &mcevt1, const MGeomCam &mgeom1, const Float_t ct1_x, const Float_t ct1_y, const MHillas &hillas2, const MMcEvt &mcevt2, const MGeomCam &mgeom2, const Float_t ct2_x, const Float_t ct2_y) +// -------------------------------------------------------------------------- +// +// Transformation of coordinates, from a point on the camera x, y , to +// the directon cosines of the corresponding direction, in the system of +// coordinates in which X-axis is North, Y-axis is west, and Z-axis +// points to the zenith. The transformation has been taken from TDAS 01-05, +// although the final system of coordinates is not the same defined there, +// but the one defined in Corsika (for convenience). +// +// rc is the distance from the reflector center to the camera. CTphi and +// CTtheta indicate the telescope orientation. The angle CTphi is the +// azimuth of the vector going along the telescope axis from the camera +// towards the reflector, measured from the North direction anticlockwise +// ( being West: phi=pi/2, South: phi=pi, East: phi=3pi/2 ) +// +// rc and x,y must be given in the same units! +// +TVector3 MStereoPar::CamToDir(const MGeomCam &geom, const MPointingPos &pos, Float_t x, Float_t y) const +{ + const Double_t rc = geom.GetCameraDist()*1e3; + + const Double_t CTphi = pos.GetAzRad(); + const Double_t CTtheta = pos.GetZdRad(); + + // + // We convert phi to the convention defined in TDAS 01-05 + // + const Double_t sinphi = sin(TMath::TwoPi()-CTphi); + const Double_t cosphi = cos(CTphi); + + const Double_t costheta = cos(CTtheta); + const Double_t sintheta = sin(CTtheta); + + const Double_t xc = x/rc; + const Double_t yc = y/rc; + + const Double_t norm = 1/sqrt(1+xc*xc+yc*yc); + + const Double_t xref = xc * norm; + const Double_t yref = yc * norm; + const Double_t zref = -1 * norm; + + const Double_t cosx = xref * sinphi + yref * costheta*cosphi - zref * sintheta*cosphi; + const Double_t cosy = -xref * cosphi + yref * costheta*sinphi - zref * sintheta*sinphi; + const Double_t cosz = yref * sintheta + zref * costheta; + + // Now change from system A of TDAS 01-05 to Corsika system: + return TVector3(cosx, -cosy, -cosz); +} + +TVector3 MStereoPar::CamToDir(const MGeomCam &geom, const MPointingPos &pos, const TVector2 &p) const +{ + return CamToDir(geom, pos, p.X(), p.Y()); +} + +void MStereoPar::CalcCT(const MHillas &h, const MPointingPos &p, const MGeomCam &g, TVector2 &cv1, TVector2 &cv2) const { // // Get the direction corresponding to the c.o.g. of the image on - // the camera - // - - Float_t ct1_cosx_a; - Float_t ct1_cosy_a; - Float_t ct1_cosz_a; // Direction from ct1 to the shower c.o.g. - - Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), hillas1.GetMeanX(), hillas1.GetMeanY(), &ct1_cosx_a, &ct1_cosy_a, &ct1_cosz_a); + // the camera. + // + // ct1_a, Direction from ct1 to the shower c.o.g. + // + const TVector3 a = CamToDir(g, p, h.GetMeanX(), h.GetMeanY()); // @@ -93,10 +139,5 @@ // to which it corresponds. // - - Float_t ct1_cosx_b; - Float_t ct1_cosy_b; - Float_t ct1_cosz_b; - - Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), hillas1.GetMeanX()+hillas1.GetCosDelta(), hillas1.GetMeanY()+hillas1.GetSinDelta(), &ct1_cosx_b, &ct1_cosy_b, &ct1_cosz_b); + const TVector3 c = CamToDir(g, p, h.GetMeanX()+h.GetCosDelta(), h.GetMeanY()+h.GetSinDelta()); // @@ -108,7 +149,6 @@ // shower core position on the z=0 plane (ground). // - - Float_t ct1_coreVersorX = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b; - Float_t ct1_coreVersorY = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b; + const Double_t coreVersorX = a.Z()*c.X() - a.X()*c.Z(); + const Double_t coreVersorY = a.Z()*c.Y() - a.X()*c.Z(); // @@ -118,61 +158,74 @@ // actually come from that direction (like for gammas from a point source) - Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), 0., 0., &ct1_cosx_b, &ct1_cosy_b, &ct1_cosz_b); - - Float_t ct1_coreVersorX_best = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b; - Float_t ct1_coreVersorY_best = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b; - - // - // Now the second telescope - // - - Float_t ct2_cosx_a; - Float_t ct2_cosy_a; - Float_t ct2_cosz_a; // Direction from ct2 to the shower c.o.g. - - - Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), hillas2.GetMeanX(), hillas2.GetMeanY(), &ct2_cosx_a, &ct2_cosy_a, &ct2_cosz_a); - - Float_t ct2_cosx_b; - Float_t ct2_cosy_b; - Float_t ct2_cosz_b; - - Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), hillas2.GetMeanX()+hillas2.GetCosDelta(), hillas2.GetMeanY()+hillas2.GetSinDelta(), &ct2_cosx_b, &ct2_cosy_b, &ct2_cosz_b); - - - Float_t ct2_coreVersorX = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b; - Float_t ct2_coreVersorY = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b; - - - Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), 0., 0., &ct2_cosx_b, &ct2_cosy_b, &ct2_cosz_b); - - Float_t ct2_coreVersorX_best = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b; - Float_t ct2_coreVersorY_best = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b; - + const TVector3 b = CamToDir(g, p, 0, 0); + + const Double_t coreVersorX_best = a.Z()*b.X() - a.X()*b.Z(); + const Double_t coreVersorY_best = a.Z()*b.Y() - a.Y()*b.Z(); + + cv1.Set(coreVersorX, coreVersorY); + cv2.Set(coreVersorX_best, coreVersorY_best); +} + +TVector2 MStereoPar::VersorToCore(const TVector2 &v1, const TVector2 &v2, const TVector2 &p1, const TVector2 &p2) const +{ + const TVector2 dp = p1 - p2; + // // Estimate core position: // - Float_t t = ct1_x - ct2_x - ct2_coreVersorX/ct2_coreVersorY*(ct1_y-ct2_y); - t /= (ct2_coreVersorX/ct2_coreVersorY*ct1_coreVersorY - ct1_coreVersorX); - - fCoreX = ct1_x + t * ct1_coreVersorX; - fCoreY = ct1_y + t * ct1_coreVersorY; - - // fCoreX, fCoreY, fCoreX2, fCoreY2 will have the same units - // as ct1_x, ct1_y, ct2_x, ct2_y - - - // - // Now the estimated core position assuming the source is located in - // the center of the camera: - // - t = ct1_x - ct2_x - ct2_coreVersorX_best / - ct2_coreVersorY_best*(ct1_y-ct2_y); - t /= (ct2_coreVersorX_best/ct2_coreVersorY_best*ct1_coreVersorY_best - - ct1_coreVersorX_best); - - fCoreX2 = ct1_x + t * ct1_coreVersorX_best; - fCoreY2 = ct1_y + t * ct1_coreVersorY_best; - + const Double_t t = + (dp.X() - v2.X()/v2.Y()*dp.Y())/(v2.X()/v2.Y()*v1.Y() - v1.X()); + + // Core will have the same units as p1/p2 + return p1 + v1*t; +} + +Double_t MStereoPar::CalcImpact(const TVector2 &w, const TVector2 &v, const TVector2 &p) const +{ + const TVector2 d = v-p; + + const Double_t f = d*w; + + return TMath::Sqrt( d.Mod2() - f*f ); +} + +Double_t MStereoPar::CalcImpact(const TVector3 &v, const TVector2 &p) const +{ + const TVector2 w = v.XYvector(); + return CalcImpact(w, w, p); +} + +Double_t MStereoPar::CalcImpact(const TVector2 &core, const TVector2 &p, const MPointingPos &point) const +{ + const TVector2 pt(-sin(point.GetZdRad()) * cos(point.GetAzRad()), + -sin(point.GetZdRad()) * sin(point.GetAzRad())); + + return CalcImpact(pt, core, p); +} + +// -------------------------------------------------------------------------- +// +// Calculation of shower parameters +// +void MStereoPar::Calc(const MHillas &hillas1, const MPointingPos &pos1, const MGeomCam &geom1, const Float_t ct1_x, const Float_t ct1_y, + const MHillas &hillas2, const MPointingPos &pos2, const MGeomCam &geom2, const Float_t ct2_x, const Float_t ct2_y) +{ + TVector2 coreVersor1, coreVersor1_best; + CalcCT(hillas1, pos1, geom1, coreVersor1, coreVersor1_best); + + TVector2 coreVersor2, coreVersor2_best; + CalcCT(hillas2, pos2, geom2, coreVersor2, coreVersor2_best); + + const TVector2 p1(ct1_x, ct1_y); + const TVector2 p2(ct2_x, ct2_y); + + // Estimate core position: + const TVector2 core = VersorToCore(coreVersor1, coreVersor2, p1, p2); + + // Now the estimated core position assuming the source is + // located in the center of the camera + const TVector2 core2 = VersorToCore(coreVersor1_best, coreVersor2_best, p1, p2); + + // Copy results to data members // // Be careful, the coordinates in MMcEvt.fCoreX,fCoreY are actually @@ -181,7 +234,10 @@ // core estimated core. // - - ///////////////////////////////////////////////////////////////////// - // + fCoreX = core.X(); + fCoreY = core.Y(); + + fCoreX2 = core2.X(); + fCoreY2 = core2.Y(); + // Now estimate the source location on the camera by intersecting // major axis of the ellipses. This assumes both telescopes are @@ -189,121 +245,30 @@ // the use of telescopes with different focal distances. // - - Float_t scale1 = mgeom1.GetConvMm2Deg(); - Float_t scale2 = mgeom2.GetConvMm2Deg(); - - t = scale2*hillas2.GetMeanY() - scale1*hillas1.GetMeanY() + - (scale1*hillas1.GetMeanX() - scale2*hillas2.GetMeanX()) * - hillas2.GetSinDelta() / hillas2.GetCosDelta(); - - t /= (hillas1.GetSinDelta() - - hillas2.GetSinDelta()/hillas2.GetCosDelta()*hillas1.GetCosDelta()); - - fSourceX = scale1*hillas1.GetMeanX() + t * hillas1.GetCosDelta(); - fSourceY = scale1*hillas1.GetMeanY() + t * hillas1.GetSinDelta(); - - // - // Squared angular distance from reconstructed source position to - // camera center. - // - fTheta2 = fSourceX*fSourceX+fSourceY*fSourceY; - - // + const TVector2 v1(hillas1.GetSinDelta(), hillas1.GetCosDelta()); + const TVector2 v2(hillas2.GetSinDelta(), hillas2.GetCosDelta()); + + const TVector2 h1 = hillas1.GetMean()*geom1.GetConvMm2Deg(); + const TVector2 h2 = hillas2.GetMean()*geom2.GetConvMm2Deg(); + + const TVector2 src = VersorToCore(v1, v2, h1, h2); + + // Reconstructed source positon + fSourceX = src.X(); + fSourceY = src.Y(); + + // Squared angular distance from reconstr. src pos to camera center. + fTheta2 = src.Mod2(); + // Get the direction corresponding to the intersection of axes - // - - Float_t source_direction[3]; - - Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), fSourceX/scale1, fSourceY/scale1, &(source_direction[0]), &(source_direction[1]), &(source_direction[2])); - - - // - // Calculate impact parameters - // - - Float_t scalar = (fCoreX-ct1_x)*source_direction[0] + - (fCoreY-ct1_y)*source_direction[1]; - fCT1Impact = sqrt( (fCoreX-ct1_x)*(fCoreX-ct1_x) + - (fCoreY-ct1_y)*(fCoreY-ct1_y) - - scalar * scalar ); - - scalar = (fCoreX-ct2_x)*source_direction[0] + - (fCoreY-ct2_y)*source_direction[1]; - fCT2Impact = sqrt( (fCoreX-ct2_x)*(fCoreX-ct2_x) + - (fCoreY-ct2_y)*(fCoreY-ct2_y) - - scalar * scalar ); - - // + const TVector3 srcdir = CamToDir(geom1, pos1, src/geom1.GetConvMm2Deg()); + + fCT1Impact = CalcImpact(srcdir, p1); + fCT2Impact = CalcImpact(srcdir, p2); + // Now calculate i.p. assuming source is point-like and placed in // the center of the camera. - // - scalar = (fCoreX2-ct1_x)*(-sin(mcevt1.GetTelescopeTheta())* - cos(mcevt1.GetTelescopePhi())) + - (fCoreY2-ct1_y)*(-sin(mcevt1.GetTelescopeTheta())* - sin(mcevt1.GetTelescopePhi())); - - fCT1Impact2 = sqrt( (fCoreX2-ct1_x)*(fCoreX2-ct1_x) + - (fCoreY2-ct1_y)*(fCoreY2-ct1_y) - - scalar * scalar ); - - scalar = (fCoreX2-ct2_x)*(-sin(mcevt2.GetTelescopeTheta())* - cos(mcevt2.GetTelescopePhi())) + - (fCoreY2-ct2_y)*(-sin(mcevt2.GetTelescopeTheta())* - sin(mcevt2.GetTelescopePhi())); - - fCT2Impact2 = sqrt( (fCoreX2-ct2_x)*(fCoreX2-ct2_x) + - (fCoreY2-ct2_y)*(fCoreY2-ct2_y) - - scalar * scalar ); - - + fCT1Impact2 = CalcImpact(core2, p1, pos1); + fCT2Impact2 = CalcImpact(core2, p2, pos2); + SetReadyToSave(); } - -// -------------------------------------------------------------------------- -// -// Transformation of coordinates, from a point on the camera x, y , to -// the director cosines of the corresponding direction, in the system of -// coordinates in which X-axis is North, Y-axis is west, and Z-axis -// points to the zenith. The transformation has been taken from TDAS 01-05, -// although the final system of coordinates is not the same defined there, -// but the one defined in Corsika (for convenience). -// -// rc is the distance from the reflector center to the camera. CTphi and -// CTtheta indicate the telescope orientation. The angle CTphi is the -// azimuth of the vector going along the telescope axis from the camera -// towards the reflector, measured from the North direction anticlockwise -// ( being West: phi=pi/2, South: phi=pi, East: phi=3pi/2 ) -// -// rc and x,y must be given in the same units! -// - - -void MStereoPar::Camera2direction(Float_t rc, Float_t CTphi, Float_t CTtheta, Float_t x, Float_t y, Float_t* cosx, Float_t* cosy, Float_t* cosz) -{ - // - // We convert phi to the convention defined in TDAS 01-05 - // - Float_t sinphi = sin(2*TMath::Pi()-CTphi); - Float_t cosphi = cos(CTphi); - Float_t costheta = cos(CTtheta); - Float_t sintheta = sin(CTtheta); - - Float_t xc = x/rc; - Float_t yc = y/rc; - - Float_t norm = 1/sqrt(1+xc*xc+yc*yc); - - Float_t xref = xc * norm; - Float_t yref = yc * norm; - Float_t zref = -1 * norm; - - *cosx = xref * sinphi + yref * costheta*cosphi - zref * sintheta*cosphi; - *cosy = -xref * cosphi + yref * costheta*sinphi - zref * sintheta*sinphi; - *cosz = yref * sintheta + zref * costheta; - - // Now change from system A of TDAS 01-05 to Corsika system: - - *cosy *= -1; - *cosz *= -1; - -} Index: trunk/MagicSoft/Mars/mimage/MStereoPar.h =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoPar.h (revision 8915) +++ trunk/MagicSoft/Mars/mimage/MStereoPar.h (revision 8916) @@ -6,7 +6,11 @@ #endif +#ifndef ROOT_TVector3 +#include +#endif + class MHillas; class MGeomCam; -class MMcEvt; +class MPointingPos; class MStereoPar : public MParContainer @@ -14,16 +18,14 @@ private: - Float_t fCoreX; - Float_t fCoreY; // Estimated core position on ground + Float_t fCoreX; // Estimated core position on ground x + Float_t fCoreY; // Estimated core position on ground y - Float_t fCoreX2; // Estimated core position on ground assuming - Float_t fCoreY2; // that the source direction is paralel to the - // telescope axis. + Float_t fCoreX2; // Estimated core position on ground assuming that + Float_t fCoreY2; // the source direction is paralel to the tel. axis. Float_t fSourceX; // Estimated source position on the camera Float_t fSourceY; // Units are degrees! - Float_t fTheta2; // deg^2; Squared angular distance of estimated - // source position to cameracenter. + Float_t fTheta2; // deg^2; Squared angular distance of estimated source position to cameracenter. Float_t fCT1Impact; // Estimated shower impact parameter from CT1 @@ -38,7 +40,14 @@ // to the telescope axis. + TVector3 CamToDir(const MGeomCam &geom, const MPointingPos &pos, Float_t x, Float_t y) const; + TVector3 CamToDir(const MGeomCam &geom, const MPointingPos &pos, const TVector2 &p) const; - void Camera2direction(Float_t rc, Float_t CTphi, Float_t CTtheta, Float_t x, Float_t y, Float_t* cosx, Float_t* cosy, Float_t* cosz); + void CalcCT(const MHillas &h, const MPointingPos &p, const MGeomCam &g, TVector2 &cv1, TVector2 &cv2) const; + TVector2 VersorToCore(const TVector2 &v1, const TVector2 &v2, const TVector2 &p1, const TVector2 &p2) const; + + Double_t CalcImpact(const TVector2 &w, const TVector2 &v, const TVector2 &p) const; + Double_t CalcImpact(const TVector3 &v, const TVector2 &p) const; + Double_t CalcImpact(const TVector2 &core, const TVector2 &p, const MPointingPos &point) const; public: @@ -57,6 +66,6 @@ Float_t GetCT2Impact2() const { return fCT2Impact2; } - - void Calc(const MHillas &hillas1, const MMcEvt &mcevt1, const MGeomCam &mgeom1, const Float_t ct1_x, const Float_t ct1_y, const MHillas &hillas2, const MMcEvt &mcevt2, const MGeomCam &mgeom2, const Float_t ct2_x, const Float_t ct2_y); + void Calc(const MHillas &h1, const MPointingPos &p1, const MGeomCam &g1, const Float_t ct1_x, const Float_t ct1_y, + const MHillas &h2, const MPointingPos &p2, const MGeomCam &g2, const Float_t ct2_x, const Float_t ct2_y); ClassDef(MStereoPar, 1) // Container to hold new image parameters @@ -64,29 +73,2 @@ #endif - - - - - - - - - - - - - - - - - - - - - - - - - - -