---

MCRecon.cxx

Go to the documentation of this file.

//  $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/src/MCRecon.cxx,v 1.10 2001/08/31 07:50:19 lsrea Exp $ 

//THB #include "data/MCData.h"
#include "reconstruction/MCRecon.h"
#include "reconstruction/ReconVisitor.h"

#include "instrument/MCTruth.h"
#include "instrument/FluxGenerator.h"

MCRecon::MCRecon ()
: Recon()
{
    setTitle("MCRecon: ver 2.0");
    
    i_MCId         = data.addElement("MC_Id");

    i_MCEnergy     = data.addElement("MC_Energy");
    i_MC_logE      = data.addElement("MC_logE");
    
    i_MCxDirection = data.addElement("MC_xDir");
    i_MCyDirection = data.addElement("MC_yDir");
    i_MCzDirection = data.addElement("MC_zDir");

    i_MCxPoint     = data.addElement("MC_X0");
    i_MCyPoint     = data.addElement("MC_Y0");
    i_MCzPoint     = data.addElement("MC_Z0");
#if 0 //THB  not set yet 
    i_MCxVertex = data.addElement("MC_xVertex");
    i_MCyVertex = data.addElement("MC_yVertex");
    i_MCzVertex = data.addElement("MC_zVertex");
    
    i_MCEnergy1    = data.addElement("MC_Energy1");
    i_MCEnergy2    = data.addElement("MC_Energy2");
    i_MCAngle1     = data.addElement("MC_Angle1");
    i_MCAngle2     = data.addElement("MC_Angle2");
    i_MCAngle12    = data.addElement("MC_Angle12");
    
    i_MCEvType     = data.addElement("MC_Evt_Type");
    i_MCMaterial   = data.addElement("MC_Mat_Code");

    i_MCSrcCode    = data.addElement("MC_Src_Code");

    i_lat =          data.addElement("MC_lattitude");
    i_lon =          data.addElement("MC_longitude");
    i_time =         data.addElement("MC_time") ;

#endif
}


void MCRecon::accept (ReconVisitor& rv)
{
    rv.visit(this);
}
void MCRecon::reconstruct (const MCTruth* truth)
{

    //old style const MCTruth* truth = MCTruth::instance();
    if( truth==0 ) truth = MCTruth::instance();
    if( truth==0 ) return; // baad

    //const MCTruth::_Particle& particle = *truth->particle();
    const FluxGenerator& particle = dynamic_cast<const FluxGenerator&>(*truth->particle());

    double energy = particle.energy();
    Vector dir=particle.launchDir().unit();
    
    data[i_MCEnergy] = energy;
    data[i_MC_logE] = log10(energy); // silly, but useful for binning
    data[i_MCxDirection] = dir.x();
    data[i_MCyDirection] = dir.y();
    data[i_MCzDirection] = dir.z();

    FluxGenerator::ParticleData* primary = particle.m_particles;
    data[i_MCId] = primary->m_idcode;

    Point vertex = primary->m_kids[0]->m_pos;

    // particle end point or interaction vertex
    data[i_MCxPoint] = vertex.x();
    data[i_MCyPoint] = vertex.y();
    data[i_MCzPoint] = vertex.z();


#if 0 //Ti_MCzPoint ODO: finish this 
    Point pos = particle.launchPt ();
    data[i_MCxPoint ] = pos.x();
    data[i_MCyPoint ] = pos.y();
    data[i_MCzPoint ] = pos.z();
     
    if(state() == DONE) return;
    
    float energy = mcData->energy();
    
    data[i_MCId] = mcData->id();

    data[i_MCEnergy] = energy;
    data[i_MC_logE] = log10(energy); // silly, but useful for binning

    data.addData(i_MCEnergy1, mcData->energy1());
    data.addData(i_MCEnergy2, mcData->energy2());
    
    Vector v0 = mcData->direction();
    data.addData(i_MCxDirection, v0.x());
    data.addData(i_MCyDirection, v0.y());
    data.addData(i_MCzDirection, v0.z());
    
    data[i_MCxPoint] = mcData->entry().x();
    data[i_MCyPoint] = mcData->entry().y();
    data[i_MCzPoint] = mcData->entry().z();
    
    // these Vectors are used more than once, thus more optimal to declare them
    Vector v1(mcData->direction1());
    Vector v2(mcData->direction2());
    
    double dot1  = v0 * v1;
    double dot2  = v0 * v2;
    double dot12 = v1 * v2;
    
    data.addData(i_MCAngle1,  acos(dot1));
    data.addData(i_MCAngle2,  acos(dot2));
    data.addData(i_MCAngle12, acos(dot12));
    
    data.addData(i_MCEvType, mcData->eventType());
    data.addData(i_MCMaterial, mcData->materialType());


    data[i_lat]  = mcData->orbit_info().lat();
    data[i_lon]  = mcData->orbit_info().lon();
    data[i_time] = mcData->orbit_info().time();

    data[i_MCxVertex] = mcData->conversionPoint().x();
    data[i_MCyVertex] = mcData->conversionPoint().y();
    data[i_MCzVertex] = mcData->conversionPoint().z();

    data.addData(i_MCSrcCode, mcData->code());
    
    setState(DONE);
#endif
    return ;
}

void MCRecon::clear ()
{
    Recon::clear();
    data.clear();
}

---