Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2566) +++ trunk/MagicSoft/Mars/Changelog (revision 2567) @@ -4,4 +4,14 @@ -*-*- END OF LINE -*-*- + 2003/11/25: Abelardo Moralejo + + * mimage/MStereoPar.[h,cc] + * mimage/MStereoCalc.[h,cc] + - Added. First version of some calculations for the data analysis + of two-telescope configurations (estimation of shower direction + and core position). + * mimage/Makefile, ImageLinkDef.h + - Updated accordingly. + 2003/11/23: Thomas Bretz Index: trunk/MagicSoft/Mars/mimage/ImageLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mimage/ImageLinkDef.h (revision 2566) +++ trunk/MagicSoft/Mars/mimage/ImageLinkDef.h (revision 2567) @@ -21,4 +21,6 @@ #pragma link C++ class MHHillasExt+; #pragma link C++ class MHNewImagePar+; +#pragma link C++ class MStereoPar+; +#pragma link C++ class MStereoCalc+; #endif Index: trunk/MagicSoft/Mars/mimage/MStereoCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoCalc.cc (revision 2567) +++ trunk/MagicSoft/Mars/mimage/MStereoCalc.cc (revision 2567) @@ -0,0 +1,167 @@ +/* ======================================================================== *\ +! +! * +! * 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): Abelardo Moralejo, 11/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MStereoCalc +// +// This is a task to calculate some shower parameters from the images of +// two telescopes in stereo mode. +// +// Input Containers: +// MGeomCam +// MHillas +// MMcEvt +// +// Output Containers: +// MStereoPar +// +///////////////////////////////////////////////////////////////////////////// +#include "MStereoCalc.h" + +#include "MParList.h" + +#include "MHillas.h" +#include "MMcEvt.hxx" +#include "MStereoPar.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MStereoCalc); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MStereoCalc::MStereoCalc(const char *name, const char *title) + : fStereoParName("MStereoPar") +{ + fName = name ? name : "MStereoCalc"; + fTitle = title ? title : "Calculate shower parameters in stereo mode"; + +} + +// -------------------------------------------------------------------------- +// +// Check for MMcEvt and MHillas containers. +// Try to find the Geometry conatiner. +// Try to find (and maybe create) the container MStereoPar. +// +Int_t MStereoCalc::PreProcess(MParList *pList) +{ + // necessary + TString mcname = "MMcEvt;"; + mcname += fCT1_id; + + fmcevt1 = (MMcEvt*)pList->FindObject(mcname); + if (!fmcevt1) + { + *fLog << err << mcname << " not found... aborting." << endl; + return kFALSE; + } + + // necessary + mcname = "MMcEvt;"; + mcname += fCT2_id; + + fmcevt2 = (MMcEvt*)pList->FindObject(mcname); + if (!fmcevt2) + { + *fLog << err << mcname << " not found... aborting." << endl; + return kFALSE; + } + + // necessary + TString geomname = "MGeomCam;"; + geomname += fCT1_id; + fGeomCam1 = (MGeomCam*)pList->FindObject(geomname); + if (!fGeomCam1) + { + *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + // necessary + geomname = "MGeomCam;"; + geomname += fCT2_id; + fGeomCam2 = (MGeomCam*)pList->FindObject(geomname); + if (!fGeomCam2) + { + *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + // necessary + TString hillasname = "MHillas;"; + hillasname += fCT1_id; + fHillas1 = (MHillas*)pList->FindObject(hillasname); + if (!fHillas1) + { + *fLog << err << hillasname << " missing in Parameter List... aborting." << endl; + return kFALSE; + } + + // necessary + hillasname = "MHillas;"; + hillasname += fCT2_id; + fHillas2 = (MHillas*)pList->FindObject(hillasname); + if (!fHillas2) + { + *fLog << err << hillasname << " missing in Parameter List... aborting." << endl; + return kFALSE; + } + + fStereoPar = (MStereoPar*)pList->FindCreateObj("MStereoPar"); + if (!fStereoPar) + { + *fLog << err << "Could not create MStereoPar... aborting" << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Call the Calc procedure of the MStereoPar object, where the +// calculations combining the data from the two telescopes are performed. +// +Int_t MStereoCalc::Process() +{ + fStereoPar->Calc(*fHillas1, *fmcevt1, *fGeomCam1, fCT1_x, fCT1_y, *fHillas2, *fmcevt2, *fGeomCam2, fCT2_x, fCT2_y); + + 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 2567) +++ trunk/MagicSoft/Mars/mimage/MStereoCalc.h (revision 2567) @@ -0,0 +1,65 @@ +#ifndef MARS_MStereoCalc +#define MARS_MStereoCalc + +///////////////////////////////////////////////////////////////////////////// +// // +// MStereoCalc // +// // +// Task to calculate some shower parameters in stereo mode // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif +#ifndef ROOT_TArrayL +#include +#endif + +class MGeomCam; +class MHillas; +class MMcEvt; +class MStereoPar; + +class MStereoCalc : public MTask +{ + const MGeomCam *fGeomCam1; //! Camera Geometry CT1 + const MHillas *fHillas1; //! input + const MMcEvt *fmcevt1; //! input + + const MGeomCam *fGeomCam2; //! Camera Geometry CT2 + const MHillas *fHillas2; //! input + const MMcEvt *fmcevt2; //! input + + Int_t fCT1_id; //! + Int_t fCT2_id; //! Identifiers of the two analyzed telescopes. + + Float_t fCT1_x; //! + Float_t fCT1_y; //! Position of first telescope + Float_t fCT2_x; //! + Float_t fCT2_y; //! Position of second telescope + + MStereoPar *fStereoPar; //! output container to store result + TString fStereoParName; // name of the 'MStereoPar' container + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + +public: + + MStereoCalc(const char *name=NULL, const char *title=NULL); + + void SetNameStereoPar(const char *name) { fStereoParName = name; } + + void SetCTids(Int_t i, Int_t j) { fCT1_id = i; fCT2_id = j; } + void SetCT1coor(Float_t x, Float_t y) { fCT1_x = x; fCT1_y = y; } // in m + void SetCT2coor(Float_t x, Float_t y) { fCT2_x = x; fCT2_y = y; } // in m + + ClassDef(MStereoCalc, 0) // Task to calculate some shower parameters in stereo mode +}; + +#endif + + Index: trunk/MagicSoft/Mars/mimage/MStereoPar.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MStereoPar.cc (revision 2567) +++ trunk/MagicSoft/Mars/mimage/MStereoPar.cc (revision 2567) @@ -0,0 +1,300 @@ +/* ======================================================================== *\ +! +! * +! * 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): Abelardo Moralejo 11/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MStereoPar +// +// Storage Container for shower parameters estimated using the information +// from two telescopes (presently for MC studies) +// +// +///////////////////////////////////////////////////////////////////////////// +#include "MStereoPar.h" +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHillas.h" +#include "MMcEvt.hxx" +#include "MGeomCam.h" + + +ClassImp(MStereoPar); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MStereoPar::MStereoPar(const char *name, const char *title) +{ + fName = name ? name : "MStereoPar"; + fTitle = title ? title : "Stereo image parameters"; + + +} + + +// -------------------------------------------------------------------------- +// +void MStereoPar::Reset() +{ + fCoreX = 0.; + fCoreY = 0.; + fSourceX = 0.; + fSourceY = 0.; +} + + +// -------------------------------------------------------------------------- +// +// 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) +{ + // + // 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); + + // + // Now we get another (arbitrary) point along the image long axis, + // fMeanX + cosdelta, fMeanY + sindelta, and calculate the direction + // 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); + + // + // The vectorial product of the latter two vectors is a vector + // perpendicular to the plane which contains the shower axis and + // passes through the telescope center (center of reflector). + // The vectorial product of that vector and (0,0,1) is a vector on + // the horizontal plane pointing from the telescope center to the + // 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; + + // + // Now we calculate again the versor, now assuming that the source + // direction is paralel to the telescope axis (camera position 0,0) + // This increases the precision of the core determination if the showers + // 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; + + // + // 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; + + // + // Be careful, the coordinates in MMcEvt.fCoreX,fCoreY are actually + // those of the vector going *from the shower core to the telescope*. + // Ours are those of the vector which goes from telescope 1 to the + // core estimated core. + // + + ///////////////////////////////////////////////////////////////////// + // + // Now estimate the source location on the camera by intersecting + // major axis of the ellipses. This assumes both telescopes are + // pointing paralel! We introduce the camera scale to account for + // 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 = scale2*hillas1.GetMeanY() + t * hillas1.GetSinDelta(); + + // + // Squared angular distance from reconstructed source position to + // camera center. + // + fTheta2 = fSourceX*fSourceX+fSourceY*fSourceY; + + // + // 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 ); + + + + 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 2567) +++ trunk/MagicSoft/Mars/mimage/MStereoPar.h (revision 2567) @@ -0,0 +1,92 @@ +#ifndef MARS_MStereoPar +#define MARS_MStereoPar + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MHillas; +class MGeomCam; +class MMcEvt; + +class MStereoPar : public MParContainer +{ +private: + + Float_t fCoreX; + Float_t fCoreY; // Estimated core position on ground + + Float_t fCoreX2; // Estimated core position on ground assuming + Float_t fCoreY2; // that the source direction is paralel to the + // telescope 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 fCT1Impact; // Estimated shower impact parameter from CT1 + Float_t fCT2Impact; // Estimated shower impact parameter from CT2 + + Float_t fCT1Impact2; // Estimated shower impact parameter from CT1 + // assuming that the source direction is paralel + // to the telescope axis. + + Float_t fCT2Impact2; // Estimated shower impact parameter from CT2 + // assuming that the source direction is paralel + // to the telescope axis. + + + 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); + + +public: + MStereoPar(const char *name=NULL, const char *title=NULL); + + void Reset(); + + Float_t GetCoreX() const { return fCoreX; } + Float_t GetCoreY() const { return fCoreY; } + Float_t GetSourceX() const { return fSourceX; } + Float_t GetSourceY() const { return fSourceY; } + Float_t GetTheta2() const { return fTheta2; } + Float_t GetCT1Impact() const { return fCT1Impact; } + Float_t GetCT2Impact() const { return fCT2Impact; } + Float_t GetCT1Impact2() const { return fCT1Impact2; } + 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); + + ClassDef(MStereoPar, 1) // Container to hold new image parameters +}; + +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mimage/Makefile =================================================================== --- trunk/MagicSoft/Mars/mimage/Makefile (revision 2566) +++ trunk/MagicSoft/Mars/mimage/Makefile (revision 2567) @@ -22,5 +22,5 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase -I../mhist -I../mgeom -I../manalysis -I../mgui +INCLUDES = -I. -I../mbase -I../mhist -I../mgeom -I../manalysis -I../mgui -I../mmc #------------------------------------------------------------------------------ @@ -40,5 +40,7 @@ MHHillasSrc.cc \ MHHillasExt.cc \ - MHNewImagePar.cc + MHNewImagePar.cc \ + MStereoPar.cc \ + MStereoCalc.cc SRCS = $(SRCFILES)