McCoordsAlg.cxx

Go to the documentation of this file.

// Include files

#include "GaudiKernel/Algorithm.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/SmartDataPtr.h"

#include "astro/GPS.h"
#include "geometry/Vector.h"

#include <cassert>
#include <map>

// forward declatation of the worker
class McCworker;

namespace { // anonymous namespace for file-global
    std::string treename("MeritTuple");
    astro::GPS* gps;
#include "Item.h"
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class McCoordsAlg : public Algorithm {

public:
    McCoordsAlg(const std::string& name, ISvcLocator* pSvcLocator); 

    StatusCode initialize();
    StatusCode execute();
    StatusCode finalize();
private:
    McCworker * m_worker;
    //counter
    int m_count;
};

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

static const AlgFactory<McCoordsAlg>  Factory;
const IAlgFactory& McCoordsAlgFactory = Factory;

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

class McCworker{ 
public:
    McCworker();

    void evaluate();

private:
    // tuple items expect to find
    //TypedItem<unsigned int, 'i'> EvtRun, EvtEventId;

    // these all float or double
    Item McXDir;
    Item McYDir;
    Item McZDir;

    //McCoords entries to create
    float m_mcRa,m_mcDec,m_mcL,m_mcB;
    float m_mcZen,m_mcAzim;
};

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
McCoordsAlg::McCoordsAlg(const std::string& name, ISvcLocator* pSvcLocator) :
Algorithm(name, pSvcLocator)
{
    declareProperty("TreeName",  treename="MeritTuple");
}

StatusCode McCoordsAlg::initialize()
{
    StatusCode  sc = StatusCode::SUCCESS;

    MsgStream log(msgSvc(), name());

    // Use the Job options service to get the Algorithm's parameters
    setProperties();

    // get a pointer to RootTupleSvc 
    if( (sc = service("RootTupleSvc", rootTupleSvc, true) ). isFailure() ) {
        log << MSG::ERROR << " failed to get the RootTupleSvc" << endreq;
        return sc;
    }
    m_worker = new McCworker();

    // get the GPS instance
    gps = astro::GPS::instance();

    m_count = 0;

    return sc;
}

StatusCode McCoordsAlg::execute()
{
    StatusCode sc = StatusCode::SUCCESS;

    MsgStream log(msgSvc(), name());

    m_count++;

    // now have the worker do it
    m_worker->evaluate();

    return sc;
}

StatusCode McCoordsAlg::finalize()
{
    return StatusCode::SUCCESS;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

McCworker::McCworker()
// initialize pointers to current items
: McXDir("McXDir")
, McYDir("McYDir")
, McZDir("McZDir")
{
    //now create new items 

    addItem( "McRa",            m_mcRa);
    addItem( "McDec",           m_mcDec);
    addItem( "McL",             m_mcL);
    addItem( "McB",             m_mcB);
    addItem( "McZenithTheta",   m_mcZen);
    addItem( "McEarthAzimuth",  m_mcAzim);
}

void McCworker::evaluate()
{
    // convert to (ra, dec)

    Vector Mc_t0(McXDir, McYDir, McZDir);
    // Mc particle direction points down... 
    // toSky converts a *particle* direction
    // into a direction on the sky, so the minus-sign is taken care of!

    astro::SkyDir mcdir = gps->toSky( Mc_t0 );
    m_mcRa   = mcdir.ra();
    m_mcDec  = mcdir.dec();
    m_mcL = mcdir.l();
    m_mcB = mcdir.b();

    // Local zenith coordinates
    astro::SkyDir zenith(gps->zenithDir());  // pointing direction
    double zenith_theta = mcdir.difference(zenith); 
    if( fabs(zenith_theta)<1e-8) zenith_theta=0;
    // all this to do the azimuth angle :-(
    Hep3Vector north_pole(0,0,1);
    Hep3Vector east_dir( north_pole.cross(zenith()).unit() ); // east is perp to north_pole and zenith
    Hep3Vector north_dir( zenith().cross(east_dir));
    double earth_azimuth=atan2( mcdir().dot(east_dir), mcdir().dot(north_dir) );
    if( earth_azimuth <0) earth_azimuth += 2*M_PI; // to 0-360 deg.
    if( fabs(earth_azimuth)<1e-8) earth_azimuth=0;
    m_mcZen  = zenith_theta*180/M_PI;;
    m_mcAzim = earth_azimuth*180/M_PI;

    return;
}
//------------------------------------------------------------------------------

---