Changeset 8916 for trunk


Ignore:
Timestamp:
06/02/08 18:22:03 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r8913 r8916  
    3737     mimage/MHillas.[h,cc], mmuon/MMuonSearchPar..[h,cc],
    3838     mpedestal/MPedestalPix.[h,cc], mpointing/MPointingDev.[h,cc],
    39      mpointing/MSrcPosCam.[h,cc]:
     39     mpointing/MSrcPosCam.[h,cc], mpointing/MPointingPos.[h,cc],
     40     mpointing/MPointing.[h,cc]:
    4041     - moved some code to source file to prevent TMath inclusion in
    4142       header (root 5.18)
    … …  
    9091       Otherwise it is ambiguous in root 5.18
    9192
     93   * mhbase/MH.cc:
     94     - adde missing includes of TColor, TMath and TClass (root 5.18)
     95     - implemented a workaround which always uses the correct
     96       CreateGradientColorTable (root 5.18)
     97
    9298   * Makefile:
    9399     - linking of the shared object is now done in /tmp
    94100     - replaced = by := where possible
     101
     102   * mjobs/MJCalibrateSignal.cc:
     103     - do not invert contcoscal, that's wrong
     104
     105   * mmovie/MMovieWrite.cc:
     106     - added a #if-directive to use either gStyle or TColor
     107       for CreateGradientColorTable depending on root-version
     108
     109   * mimage/MStereoPar.[h,cc], mimage/MStereoCal.[h,cc]:
     110     - replaced Monate Carlo container by MPointingPos
     111     - made every algorithm unique
    95112
    96113

  • trunk/MagicSoft/Mars/mhbase/MH.cc

    r8892 r8916  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.37 2008-05-19 14:04:12 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.38 2008-06-02 17:21:24 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    2020!   Author(s): Thomas Bretz  07/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
    2121!
    22 !   Copyright: MAGIC Software Development, 2000-2007
     22!   Copyright: MAGIC Software Development, 2000-2008
    2323!
    2424!
    … …  
    5757#include <TH2.h>
    5858#include <TH3.h>
     59#include <TColor.h>
     60#include <TMath.h>
     61#include <TClass.h>
    5962#include <TStyle.h>       // TStyle::GetScreenFactor
    6063#include <TGaxis.h>
    … …  
    6366#include <TPaveStats.h>
    6467#include <TBaseClass.h>
     68#include <THashList.h>
    6569#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
    6670#include <THLimitsFinder.h>
    … …  
    15961600}
    15971601
     1602#ifdef CreateGradientColorTable
     1603#error CreateGradientColorTable already defined
     1604#endif
     1605
     1606#if ROOT_VERSION_CODE < ROOT_VERSION(5,18,00)
     1607#define CreateGradientColorTable gStyle->CreateGradientColorTable
     1608#else
     1609#define CreateGradientColorTable TColor::CreateGradientColorTable
     1610#endif
     1611
    15981612// --------------------------------------------------------------------------
    15991613//
    … …  
    16391653        Double_t g[5] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
    16401654        Double_t b[5] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
    1641         gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
     1655        CreateGradientColorTable(5, s, r, g, b, ncol);
    16421656        found=kTRUE;
    16431657    }
    … …  
    16491663        Double_t g[5] = { 0.00, 0.02, 0.40, 0.70, 1.00 };
    16501664        Double_t b[5] = { 0.00, 0.27, 0.51, 0.81, 1.00 };
    1651         gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
     1665        CreateGradientColorTable(5, s, r, g, b, ncol);
    16521666        found=kTRUE;
    16531667    }
    … …  
    16591673        double g[2] = {0.00, 1.00};
    16601674        double b[2] = {0.00, 1.00};
    1661         gStyle->CreateGradientColorTable(2, s, r, g, b, ncol);
     1675        CreateGradientColorTable(2, s, r, g, b, ncol);
    16621676        found=kTRUE;
    16631677    }
    … …  
    16691683        double g[5] = {0., 0.10, 0.20, 0.73, 1.00};
    16701684        double b[5] = {0., 0.03, 0.06, 0.00, 1.00};
    1671         gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
     1685        CreateGradientColorTable(5, s, r, g, b, ncol);
    16721686        found=kTRUE;
    16731687    }
    … …  
    16791693        double g[4] = {0.03, 0.04, 0.80, 1.00};
    16801694        double b[4] = {0.03, 0.04, 0.00, 1.00};
    1681         gStyle->CreateGradientColorTable(4, s, r, g, b, ncol);
     1695        CreateGradientColorTable(4, s, r, g, b, ncol);
    16821696        found=kTRUE;
    16831697    }
    … …  
    16891703        double g[8] = {0.70, 0.40, 0.02, 0.00, 0.10, 0.20, 0.73, 1.00};
    16901704        double b[8] = {0.81, 0.51, 0.27, 0.00, 0.03, 0.06, 0.00, 1.00};
    1691         gStyle->CreateGradientColorTable(8, s, r, g, b, ncol);
     1705        CreateGradientColorTable(8, s, r, g, b, ncol);
    16921706        found=kTRUE;
    16931707    }
    … …  
    16991713        double g[3] = {0., 0.0, 1.};
    17001714        double b[3] = {0., 0.0, 1.};
    1701         gStyle->CreateGradientColorTable(3, s, r, g, b, ncol);
     1715        CreateGradientColorTable(3, s, r, g, b, ncol);
    17021716        found=kTRUE;
    17031717    }
    … …  
    17091723        double g[3] = {0., 0.0, 1.};
    17101724        double b[3] = {0., 1.0, 1.};
    1711         gStyle->CreateGradientColorTable(3, s, r, g, b, ncol);
     1725        CreateGradientColorTable(3, s, r, g, b, ncol);
    17121726        found=kTRUE;
    17131727    }
    … …  
    17191733        double g[4] = {0.28, 0.93, 0.03, 1.};
    17201734        double b[4] = {0.79, 0.11, 0.03, 1.};
    1721         gStyle->CreateGradientColorTable(4, s, r, g, b, ncol);
     1735        CreateGradientColorTable(4, s, r, g, b, ncol);
    17221736        found=kTRUE;
    17231737    }
    … …  
    17281742        double g[5] = {0.1, 0.1, 0.0, 1.0, 1.0 };
    17291743        double b[5] = {0.2, 0.7, 0.0, 0.3, 0.9 };
    1730         gStyle->CreateGradientColorTable(5, s, r, g, b, ncol);
     1744        CreateGradientColorTable(5, s, r, g, b, ncol);
    17311745        found=kTRUE;
    17321746    }

  • trunk/MagicSoft/Mars/mimage/MStereoCalc.cc

    r2596 r8916  
    1818!   Author(s): Abelardo Moralejo, 11/2003 <mailto:moralejo@pd.infn.it>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     20!   Copyright: MAGIC Software Development, 2000-2008
    2121!
    2222!
    … …  
    3333//   MGeomCam
    3434//   MHillas
    35 //   MMcEvt
     35//   MPointingPos
    3636//
    3737//  Output Containers:
    … …  
    4444
    4545#include "MHillas.h"
    46 #include "MMcEvt.hxx"
    4746#include "MStereoPar.h"
     47#include "MPointingPos.h"
    4848
    4949#include "MLog.h"
    … …  
    6363    fName  = name  ? name  : "MStereoCalc";
    6464    fTitle = title ? title : "Calculate shower parameters in stereo mode";
    65 
    6665}
    6766
    … …  
    7473Int_t MStereoCalc::PreProcess(MParList *pList)
    7574{
    76     // necessary
    77 
    78     fmcevt1 = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt",fCT1id));
    79     if (!fmcevt1)
     75    fPointingPos1 = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos",fCT1id));
     76    if (!fPointingPos1)
    8077    {
    81         *fLog << err << AddSerialNumber("MMcEvt",fCT1id) << " not found... aborting." << endl;
     78        *fLog << err << AddSerialNumber("MPointingPos",fCT1id) << " not found... aborting." << endl;
    8279        return kFALSE;
    8380    }
    8481
    85     // necessary
    86     fmcevt2 = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt",fCT2id));
    87     if (!fmcevt2)
     82    fPointingPos2 = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos",fCT2id));
     83    if (!fPointingPos2)
    8884    {
    89       *fLog << err << AddSerialNumber("MMcEvt",fCT2id) << " not found... aborting." << endl;
     85        *fLog << err << AddSerialNumber("MPointingPos",fCT2id) << " not found... aborting." << endl;
    9086        return kFALSE;
    9187    }
    9288
    93     // necessary
    94     TString geomname = "MGeomCam;";
    95     geomname += fCT1id;
    96     fGeomCam1 = (MGeomCam*)pList->FindObject(geomname);
     89    fGeomCam1 = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam", fCT1id));
    9790    if (!fGeomCam1)
    9891    {
    99         *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl;
     92        *fLog << err << AddSerialNumber("MGeomCam", fCT1id) << " not found... aborting." << endl;
    10093        return kFALSE;
    10194    }
    10295
    103     // necessary
    104     geomname = "MGeomCam;";
    105     geomname += fCT2id;
    106     fGeomCam2 = (MGeomCam*)pList->FindObject(geomname);
     96    fGeomCam2 = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam", fCT2id));
    10797    if (!fGeomCam2)
    10898    {
    109         *fLog << err << geomname << " (Camera Geometry) missing in Parameter List... aborting." << endl;
     99        *fLog << err << AddSerialNumber("MGeomCam", fCT2id) << " not found... aborting." << endl;
    110100        return kFALSE;
    111101    }
    112102
    113     // necessary
    114     TString hillasname = "MHillas;";
    115     hillasname += fCT1id;
    116     fHillas1 = (MHillas*)pList->FindObject(hillasname);
     103    fHillas1 = (MHillas*)pList->FindObject(AddSerialNumber("MHillas", fCT1id));
    117104    if (!fHillas1)
    118105    {
    119         *fLog << err << hillasname << " missing in Parameter List... aborting." << endl;
     106        *fLog << err << AddSerialNumber("MHillas", fCT1id) << " not found... aborting." << endl;
    120107        return kFALSE;
    121108    }
    122109
    123     // necessary
    124     hillasname = "MHillas;";
    125     hillasname += fCT2id;
    126     fHillas2 = (MHillas*)pList->FindObject(hillasname);
     110    fHillas2 = (MHillas*)pList->FindObject(AddSerialNumber("MHillas", fCT2id));
    127111    if (!fHillas2)
    128112    {
    129         *fLog << err << hillasname << " missing in Parameter List... aborting." << endl;
     113        *fLog << err << AddSerialNumber("MHillas", fCT2id) << " not found... aborting." << endl;
    130114        return kFALSE;
    131115    }
    132116
    133     fStereoPar = (MStereoPar*)pList->FindCreateObj("MStereoPar");
     117    fStereoPar = (MStereoPar*)pList->FindCreateObj("MStereoPar", fStereoParName);
    134118    if (!fStereoPar)
    135     {
    136         *fLog << err << "Could not create MStereoPar... aborting" << endl;
    137119        return kFALSE;
    138     }
    139120
    140121    return kTRUE;
    … …  
    148129Int_t MStereoCalc::Process()
    149130{
    150     fStereoPar->Calc(*fHillas1, *fmcevt1, *fGeomCam1, fCT1x, fCT1y, *fHillas2, *fmcevt2, *fGeomCam2, fCT2x, fCT2y);
     131    fStereoPar->Calc(*fHillas1, *fPointingPos1, *fGeomCam1, fCT1x, fCT1y,
     132                     *fHillas2, *fPointingPos2, *fGeomCam2, fCT2x, fCT2y);
    151133
    152134    return kTRUE;
    153135}
    154 
    155 // --------------------------------------------------------------------------
    156 //
    157 // Does nothing at the moment.
    158 //
    159 Int_t MStereoCalc::PostProcess()
    160 {
    161     return kTRUE;
    162 }

  • trunk/MagicSoft/Mars/mimage/MStereoCalc.h

    r2595 r8916  
    1313#include "MTask.h"
    1414#endif
    15 #ifndef ROOT_TArrayL
    16 #include <TArrayL.h>
    17 #endif
    1815
    1916class MGeomCam;
    2017class MHillas;
    21 class MMcEvt;
    2218class MStereoPar;
     19class MPointingPos;
    2320
    2421class MStereoCalc : public MTask
    2522{
    26     const MGeomCam    *fGeomCam1;    //! Camera Geometry CT1
    27     const MHillas     *fHillas1;     //! input
    28     const MMcEvt      *fmcevt1;      //! input
     23    const MGeomCam     *fGeomCam1;     //! CT1: Camera Geometry
     24    const MHillas      *fHillas1;      //! CT1: Hillas parameters
     25    const MPointingPos *fPointingPos1; //! CT1: Pointing Direction
    2926
    30     const MGeomCam    *fGeomCam2;    //! Camera Geometry CT2
    31     const MHillas     *fHillas2;     //! input
    32     const MMcEvt      *fmcevt2;      //! input
     27    const MGeomCam     *fGeomCam2;     //! CT2: Camera Geometry
     28    const MHillas      *fHillas2;      //! CT2: Hillas parameters
     29    const MPointingPos *fPointingPos2; //! CT2: pointing direction
    3330
    34     Int_t fCT1id;   //!
    35     Int_t fCT2id;   //! Identifiers of the two analyzed telescopes.
     31    Int_t fCT1id;   // CT1: Identifier number
     32    Int_t fCT2id;   // CT2: Identifier number
    3633
    37     Float_t fCT1x;   //!
     34    Float_t fCT1x;   //! FIXME -> Move to parameter list
    3835    Float_t fCT1y;   //! Position of first telescope
    39     Float_t fCT2x;   //!
     36    Float_t fCT2x;   //! FIXME -> Move to parameter list
    4037    Float_t fCT2y;   //! Position of second telescope
    4138
    … …  
    4542    Int_t PreProcess(MParList *pList);
    4643    Int_t Process();
    47     Int_t PostProcess();
    48 
    4944
    5045public:

  • trunk/MagicSoft/Mars/mimage/MStereoPar.cc

    r2770 r8916  
    1818!   Author(s): Abelardo Moralejo 11/2003 <mailto:moralejo@pd.infn.it>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     20!   Copyright: MAGIC Software Development, 2000-2008
    2121!
    2222!
    … …  
    3333/////////////////////////////////////////////////////////////////////////////
    3434#include "MStereoPar.h"
     35
    3536#include <fstream>
     37
     38#include <TMath.h>
    3639
    3740#include "MLog.h"
    … …  
    3942
    4043#include "MHillas.h"
    41 #include "MMcEvt.hxx"
    4244#include "MGeomCam.h"
    43 
     45#include "MPointingPos.h"
    4446
    4547ClassImp(MStereoPar);
    … …  
    5557    fName  = name  ? name  : "MStereoPar";
    5658    fTitle = title ? title : "Stereo image parameters";
    57 
    58 
    59 }
    60 
     59}
    6160
    6261// --------------------------------------------------------------------------
    … …  
    7069}
    7170
    72 
    73 // --------------------------------------------------------------------------
    74 //
    75 //  Calculation of shower parameters
    76 //
    77 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)
     71// --------------------------------------------------------------------------
     72//
     73// Transformation of coordinates, from a point on the camera x, y , to
     74// the directon cosines of the corresponding direction, in the system of
     75// coordinates in which X-axis is North, Y-axis is west, and Z-axis
     76// points to the zenith. The transformation has been taken from TDAS 01-05,
     77// although the final system of coordinates is not the same defined there,
     78// but the one defined in Corsika (for convenience).
     79//
     80// rc is the distance from the reflector center to the camera. CTphi and
     81// CTtheta indicate the telescope orientation. The angle CTphi is the
     82// azimuth of the vector going along the telescope axis from the camera
     83// towards the reflector, measured from the North direction anticlockwise
     84// ( being West: phi=pi/2, South: phi=pi, East: phi=3pi/2 )
     85//
     86// rc and x,y must be given in the same units!
     87//
     88TVector3 MStereoPar::CamToDir(const MGeomCam &geom, const MPointingPos &pos, Float_t x, Float_t y) const
     89{
     90    const Double_t rc      = geom.GetCameraDist()*1e3;
     91
     92    const Double_t CTphi   = pos.GetAzRad();
     93    const Double_t CTtheta = pos.GetZdRad();
     94
     95    //
     96    // We convert phi to the convention defined in TDAS 01-05
     97    //
     98    const Double_t sinphi   = sin(TMath::TwoPi()-CTphi);
     99    const Double_t cosphi   = cos(CTphi);
     100
     101    const Double_t costheta = cos(CTtheta);
     102    const Double_t sintheta = sin(CTtheta);
     103
     104    const Double_t xc = x/rc;
     105    const Double_t yc = y/rc;
     106
     107    const Double_t norm = 1/sqrt(1+xc*xc+yc*yc);
     108
     109    const Double_t xref = xc * norm;
     110    const Double_t yref = yc * norm;
     111    const Double_t zref = -1 * norm;
     112
     113    const Double_t cosx =  xref * sinphi + yref * costheta*cosphi - zref * sintheta*cosphi;
     114    const Double_t cosy = -xref * cosphi + yref * costheta*sinphi - zref * sintheta*sinphi;
     115    const Double_t cosz =                  yref * sintheta        + zref * costheta;
     116
     117    //  Now change from system A of TDAS 01-05 to Corsika system:
     118    return TVector3(cosx, -cosy, -cosz);
     119}
     120
     121TVector3 MStereoPar::CamToDir(const MGeomCam &geom, const MPointingPos &pos, const TVector2 &p) const
     122{
     123    return CamToDir(geom, pos, p.X(), p.Y());
     124}
     125
     126void MStereoPar::CalcCT(const MHillas &h, const MPointingPos &p, const MGeomCam &g, TVector2 &cv1, TVector2 &cv2) const
    78127{
    79128    //
    80129    // Get the direction corresponding to the c.o.g. of the image on
    81     // the camera
    82     //
    83 
    84     Float_t ct1_cosx_a;
    85     Float_t ct1_cosy_a;
    86     Float_t ct1_cosz_a; // Direction from ct1 to the shower c.o.g.
    87 
    88     Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), hillas1.GetMeanX(), hillas1.GetMeanY(), &ct1_cosx_a, &ct1_cosy_a, &ct1_cosz_a);
     130    // the camera.
     131    //
     132    // ct1_a, Direction from ct1 to the shower c.o.g.
     133    //
     134    const TVector3 a = CamToDir(g, p, h.GetMeanX(), h.GetMeanY());
    89135
    90136    //
    … …  
    93139    // to which it corresponds.
    94140    //
    95    
    96     Float_t ct1_cosx_b;
    97     Float_t ct1_cosy_b;
    98     Float_t ct1_cosz_b;
    99 
    100     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);
     141    const TVector3 c = CamToDir(g, p, h.GetMeanX()+h.GetCosDelta(), h.GetMeanY()+h.GetSinDelta());
    101142
    102143    //
    … …  
    108149    // shower core position on the z=0 plane (ground).
    109150    //
    110 
    111     Float_t ct1_coreVersorX = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b;
    112     Float_t ct1_coreVersorY = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b;
     151    const Double_t coreVersorX = a.Z()*c.X() - a.X()*c.Z();
     152    const Double_t coreVersorY = a.Z()*c.Y() - a.X()*c.Z();
    113153
    114154    //
    … …  
    118158    // actually come from that direction (like for gammas from a point source)
    119159
    120     Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), 0., 0., &ct1_cosx_b, &ct1_cosy_b, &ct1_cosz_b);
    121 
    122     Float_t ct1_coreVersorX_best = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b;
    123     Float_t ct1_coreVersorY_best = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b;
    124    
    125     //
    126     // Now the second telescope
    127     //
    128 
    129     Float_t ct2_cosx_a;
    130     Float_t ct2_cosy_a;
    131     Float_t ct2_cosz_a; // Direction from ct2 to the shower c.o.g.
    132 
    133 
    134     Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), hillas2.GetMeanX(), hillas2.GetMeanY(), &ct2_cosx_a, &ct2_cosy_a, &ct2_cosz_a);
    135 
    136     Float_t ct2_cosx_b;
    137     Float_t ct2_cosy_b;
    138     Float_t ct2_cosz_b;
    139 
    140     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);
    141 
    142 
    143     Float_t ct2_coreVersorX = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b;
    144     Float_t ct2_coreVersorY = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b;
    145 
    146 
    147     Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), 0., 0., &ct2_cosx_b, &ct2_cosy_b, &ct2_cosz_b);
    148 
    149     Float_t ct2_coreVersorX_best = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b;
    150     Float_t ct2_coreVersorY_best = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b;
    151    
     160    const TVector3 b = CamToDir(g, p, 0, 0);
     161
     162    const Double_t coreVersorX_best = a.Z()*b.X() - a.X()*b.Z();
     163    const Double_t coreVersorY_best = a.Z()*b.Y() - a.Y()*b.Z();
     164
     165    cv1.Set(coreVersorX,      coreVersorY);
     166    cv2.Set(coreVersorX_best, coreVersorY_best);
     167}
     168
     169TVector2 MStereoPar::VersorToCore(const TVector2 &v1, const TVector2 &v2, const TVector2 &p1, const TVector2 &p2) const
     170{
     171    const TVector2 dp = p1 - p2;
     172
    152173    //
    153174    // Estimate core position:
    154175    //
    155     Float_t t = ct1_x - ct2_x - ct2_coreVersorX/ct2_coreVersorY*(ct1_y-ct2_y);
    156     t /= (ct2_coreVersorX/ct2_coreVersorY*ct1_coreVersorY - ct1_coreVersorX);
    157 
    158     fCoreX = ct1_x + t * ct1_coreVersorX;
    159     fCoreY = ct1_y + t * ct1_coreVersorY;
    160 
    161     // fCoreX, fCoreY, fCoreX2, fCoreY2 will have the same units
    162     // as ct1_x, ct1_y, ct2_x, ct2_y
    163 
    164 
    165     //
    166     // Now the estimated core position assuming the source is located in
    167     // the center of the camera:
    168     //
    169     t = ct1_x - ct2_x - ct2_coreVersorX_best /
    170         ct2_coreVersorY_best*(ct1_y-ct2_y);
    171     t /= (ct2_coreVersorX_best/ct2_coreVersorY_best*ct1_coreVersorY_best -
    172           ct1_coreVersorX_best);
    173 
    174     fCoreX2 = ct1_x + t * ct1_coreVersorX_best;
    175     fCoreY2 = ct1_y + t * ct1_coreVersorY_best;
    176 
     176    const Double_t t =
     177        (dp.X() - v2.X()/v2.Y()*dp.Y())/(v2.X()/v2.Y()*v1.Y() - v1.X());
     178
     179    // Core will have the same units as p1/p2
     180    return p1 + v1*t;
     181}
     182
     183Double_t MStereoPar::CalcImpact(const TVector2 &w, const TVector2 &v, const TVector2 &p) const
     184{
     185    const TVector2 d = v-p;
     186
     187    const Double_t f = d*w;
     188
     189    return TMath::Sqrt( d.Mod2() - f*f );
     190}
     191
     192Double_t MStereoPar::CalcImpact(const TVector3 &v, const TVector2 &p) const
     193{
     194    const TVector2 w = v.XYvector();
     195    return CalcImpact(w, w, p);
     196}
     197
     198Double_t MStereoPar::CalcImpact(const TVector2 &core, const TVector2 &p, const MPointingPos &point) const
     199{
     200    const TVector2 pt(-sin(point.GetZdRad()) * cos(point.GetAzRad()),
     201                      -sin(point.GetZdRad()) * sin(point.GetAzRad()));
     202
     203    return CalcImpact(pt, core, p);
     204}
     205
     206// --------------------------------------------------------------------------
     207//
     208//  Calculation of shower parameters
     209//
     210void MStereoPar::Calc(const MHillas &hillas1, const MPointingPos &pos1, const MGeomCam &geom1, const Float_t ct1_x, const Float_t ct1_y,
     211                      const MHillas &hillas2, const MPointingPos &pos2, const MGeomCam &geom2, const Float_t ct2_x, const Float_t ct2_y)
     212{
     213    TVector2 coreVersor1, coreVersor1_best;
     214    CalcCT(hillas1, pos1, geom1, coreVersor1, coreVersor1_best);
     215
     216    TVector2 coreVersor2, coreVersor2_best;
     217    CalcCT(hillas2, pos2, geom2, coreVersor2, coreVersor2_best);
     218
     219    const TVector2 p1(ct1_x, ct1_y);
     220    const TVector2 p2(ct2_x, ct2_y);
     221
     222    // Estimate core position:
     223    const TVector2 core = VersorToCore(coreVersor1, coreVersor2, p1, p2);
     224
     225    // Now the estimated core position assuming the source is
     226    // located in the center of the camera
     227    const TVector2 core2 = VersorToCore(coreVersor1_best, coreVersor2_best, p1, p2);
     228
     229    // Copy results to data members
    177230    //
    178231    // Be careful, the coordinates in MMcEvt.fCoreX,fCoreY are actually
    … …  
    181234    // core estimated core.
    182235    //
    183 
    184     /////////////////////////////////////////////////////////////////////
    185     //
     236    fCoreX  = core.X();
     237    fCoreY  = core.Y();
     238
     239    fCoreX2 = core2.X();
     240    fCoreY2 = core2.Y();
     241
    186242    // Now estimate the source location on the camera by intersecting
    187243    // major axis of the ellipses. This assumes both telescopes are
    … …  
    189245    // the use of telescopes with different focal distances.
    190246    //
    191 
    192     Float_t scale1 = mgeom1.GetConvMm2Deg();
    193     Float_t scale2 = mgeom2.GetConvMm2Deg();
    194 
    195     t = scale2*hillas2.GetMeanY() - scale1*hillas1.GetMeanY() +
    196         (scale1*hillas1.GetMeanX() - scale2*hillas2.GetMeanX()) *
    197         hillas2.GetSinDelta() / hillas2.GetCosDelta();
    198 
    199     t /= (hillas1.GetSinDelta() -
    200           hillas2.GetSinDelta()/hillas2.GetCosDelta()*hillas1.GetCosDelta());
    201 
    202     fSourceX = scale1*hillas1.GetMeanX() + t * hillas1.GetCosDelta();
    203     fSourceY = scale1*hillas1.GetMeanY() + t * hillas1.GetSinDelta();
    204 
    205     //
    206     // Squared angular distance from reconstructed source position to
    207     // camera center.
    208     //
    209     fTheta2 = fSourceX*fSourceX+fSourceY*fSourceY;
    210 
    211     //
     247    const TVector2 v1(hillas1.GetSinDelta(), hillas1.GetCosDelta());
     248    const TVector2 v2(hillas2.GetSinDelta(), hillas2.GetCosDelta());
     249
     250    const TVector2 h1 = hillas1.GetMean()*geom1.GetConvMm2Deg();
     251    const TVector2 h2 = hillas2.GetMean()*geom2.GetConvMm2Deg();
     252
     253    const TVector2 src = VersorToCore(v1, v2, h1, h2);
     254
     255    // Reconstructed source positon
     256    fSourceX = src.X();
     257    fSourceY = src.Y();
     258
     259    // Squared angular distance from reconstr. src pos to camera center.
     260    fTheta2 = src.Mod2();
     261
    212262    // Get the direction corresponding to the intersection of axes
    213     //
    214 
    215     Float_t source_direction[3];
    216 
    217     Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), fSourceX/scale1, fSourceY/scale1, &(source_direction[0]), &(source_direction[1]), &(source_direction[2]));
    218 
    219 
    220     //
    221     // Calculate impact parameters
    222     //
    223 
    224     Float_t scalar = (fCoreX-ct1_x)*source_direction[0] +
    225         (fCoreY-ct1_y)*source_direction[1];
    226     fCT1Impact = sqrt( (fCoreX-ct1_x)*(fCoreX-ct1_x) +
    227                        (fCoreY-ct1_y)*(fCoreY-ct1_y) -
    228                        scalar * scalar );
    229 
    230     scalar = (fCoreX-ct2_x)*source_direction[0] +
    231         (fCoreY-ct2_y)*source_direction[1];
    232     fCT2Impact = sqrt( (fCoreX-ct2_x)*(fCoreX-ct2_x) +
    233                        (fCoreY-ct2_y)*(fCoreY-ct2_y) -
    234                        scalar * scalar );
    235 
    236     //
     263    const TVector3 srcdir = CamToDir(geom1, pos1, src/geom1.GetConvMm2Deg());
     264
     265    fCT1Impact = CalcImpact(srcdir, p1);
     266    fCT2Impact = CalcImpact(srcdir, p2);
     267
    237268    // Now calculate i.p. assuming source is point-like and placed in
    238269    // the center of the camera.
    239     //
    240     scalar = (fCoreX2-ct1_x)*(-sin(mcevt1.GetTelescopeTheta())*
    241                              cos(mcevt1.GetTelescopePhi()))  +
    242       (fCoreY2-ct1_y)*(-sin(mcevt1.GetTelescopeTheta())*
    243                       sin(mcevt1.GetTelescopePhi()));
    244 
    245     fCT1Impact2 = sqrt( (fCoreX2-ct1_x)*(fCoreX2-ct1_x) +
    246                        (fCoreY2-ct1_y)*(fCoreY2-ct1_y) -
    247                        scalar * scalar );
    248 
    249     scalar = (fCoreX2-ct2_x)*(-sin(mcevt2.GetTelescopeTheta())*
    250                              cos(mcevt2.GetTelescopePhi()))  +
    251       (fCoreY2-ct2_y)*(-sin(mcevt2.GetTelescopeTheta())*
    252                       sin(mcevt2.GetTelescopePhi()));
    253 
    254     fCT2Impact2 = sqrt( (fCoreX2-ct2_x)*(fCoreX2-ct2_x) +
    255                        (fCoreY2-ct2_y)*(fCoreY2-ct2_y) -
    256                        scalar * scalar );
    257 
    258  
     270    fCT1Impact2 = CalcImpact(core2, p1, pos1);
     271    fCT2Impact2 = CalcImpact(core2, p2, pos2);
     272
    259273    SetReadyToSave();
    260274}
    261 
    262 // --------------------------------------------------------------------------
    263 //
    264 // Transformation of coordinates, from a point on the camera x, y , to
    265 // the director cosines of the corresponding direction, in the system of
    266 // coordinates in which X-axis is North, Y-axis is west, and Z-axis
    267 // points to the zenith. The transformation has been taken from TDAS 01-05,
    268 // although the final system of coordinates is not the same defined there,
    269 // but the one defined in Corsika (for convenience).
    270 //
    271 // rc is the distance from the reflector center to the camera. CTphi and
    272 // CTtheta indicate the telescope orientation. The angle CTphi is the
    273 // azimuth of the vector going along the telescope axis from the camera
    274 // towards the reflector, measured from the North direction anticlockwise
    275 // ( being West: phi=pi/2, South: phi=pi, East: phi=3pi/2 )
    276 //
    277 // rc and x,y must be given in the same units!
    278 //
    279  
    280 
    281 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)
    282 {
    283     //
    284     // We convert phi to the convention defined in TDAS 01-05
    285     //
    286     Float_t sinphi = sin(2*TMath::Pi()-CTphi);
    287     Float_t cosphi = cos(CTphi);
    288     Float_t costheta = cos(CTtheta);
    289     Float_t sintheta = sin(CTtheta);
    290 
    291     Float_t xc = x/rc;
    292     Float_t yc = y/rc;
    293 
    294     Float_t norm = 1/sqrt(1+xc*xc+yc*yc);
    295 
    296     Float_t xref = xc * norm;
    297     Float_t yref = yc * norm;
    298     Float_t zref = -1 * norm;
    299 
    300     *cosx =  xref * sinphi + yref * costheta*cosphi - zref * sintheta*cosphi;
    301     *cosy = -xref * cosphi + yref * costheta*sinphi - zref * sintheta*sinphi;
    302     *cosz =                  yref * sintheta        + zref * costheta;
    303 
    304     //  Now change from system A of TDAS 01-05 to Corsika system:
    305 
    306     *cosy *= -1;
    307     *cosz *= -1;
    308 
    309 }

  • trunk/MagicSoft/Mars/mimage/MStereoPar.h

    r2567 r8916  
    66#endif
    77
     8#ifndef ROOT_TVector3
     9#include <TVector3.h>
     10#endif
     11
    812class MHillas;
    913class MGeomCam;
    10 class MMcEvt;
     14class MPointingPos;
    1115
    1216class MStereoPar : public MParContainer
    … …  
    1418private:
    1519
    16     Float_t fCoreX;
    17     Float_t fCoreY;   // Estimated core position on ground
     20    Float_t fCoreX;   // Estimated core position on ground x
     21    Float_t fCoreY;   // Estimated core position on ground y
    1822
    19     Float_t fCoreX2;  // Estimated core position on ground assuming
    20     Float_t fCoreY2;  // that the source direction is paralel to the
    21                       // telescope axis.
     23    Float_t fCoreX2;  // Estimated core position on ground assuming that
     24    Float_t fCoreY2;  // the source direction is paralel to the tel. axis.
    2225
    2326    Float_t fSourceX; // Estimated source position on the camera
    2427    Float_t fSourceY; // Units are degrees!
    2528
    26     Float_t fTheta2;  // deg^2; Squared angular distance of estimated
    27                       // source position to cameracenter.
     29    Float_t fTheta2;  // deg^2; Squared angular distance of estimated source position to cameracenter.
    2830
    2931    Float_t fCT1Impact; // Estimated shower impact parameter from CT1
    … …  
    3840                         // to the telescope axis.
    3941
     42    TVector3 CamToDir(const MGeomCam &geom, const MPointingPos &pos, Float_t x, Float_t y) const;
     43    TVector3 CamToDir(const MGeomCam &geom, const MPointingPos &pos, const TVector2 &p) const;
    4044
    41     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);
     45    void     CalcCT(const MHillas &h, const MPointingPos &p, const MGeomCam &g, TVector2 &cv1, TVector2 &cv2) const;
    4246
     47    TVector2 VersorToCore(const TVector2 &v1, const TVector2 &v2, const TVector2 &p1, const TVector2 &p2) const;
     48
     49    Double_t CalcImpact(const TVector2 &w, const TVector2 &v, const TVector2 &p) const;
     50    Double_t CalcImpact(const TVector3 &v, const TVector2 &p) const;
     51    Double_t CalcImpact(const TVector2 &core, const TVector2 &p, const MPointingPos &point) const;
    4352
    4453public:
    … …  
    5766    Float_t GetCT2Impact2() const { return fCT2Impact2; }
    5867
    59 
    60     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);
     68    void Calc(const MHillas &h1, const MPointingPos &p1, const MGeomCam &g1, const Float_t ct1_x, const Float_t ct1_y,
     69              const MHillas &h2, const MPointingPos &p2, const MGeomCam &g2, const Float_t ct2_x, const Float_t ct2_y);
    6170
    6271    ClassDef(MStereoPar, 1) // Container to hold new image parameters
    … …  
    6473
    6574#endif
    66 
    67 
    68 
    69 
    70 
    71 
    72 
    73 
    74 
    75 
    76 
    77 
    78 
    79 
    80 
    81 
    82 
    83 
    84 
    85 
    86 
    87 
    88 
    89 
    90 
    91 
    92 
Note: See TracChangeset for help on using the changeset viewer.