---

acdReconAlg.cxx

Go to the documentation of this file.

#include "acdReconAlg.h"

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

#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"


#include "GlastEvent/TopLevel/EventModel.h"
#include "GlastEvent/TopLevel/Event.h"
#include "GaudiKernel/ObjectVector.h"
#include "GlastEvent/data/TdGlastData.h"

#include "TkrRecon/SiRecObjs.h"
#include "TkrRecon/GFparticle.h"
#include "TkrRecon/GFdata.h"
#include "TkrRecon/GFgamma.h"

#include "geometry/Point.h"
#include "geometry/Vector.h"
#include "geometry/Box.h"

#include "xml/IFile.h"

#include "idents/AcdId.h"


#include <algorithm>
#include <cstdio>
#include <stdlib.h>

double acdReconAlg::threshold_energy;

int acdReconAlg::xNumTowers;
int acdReconAlg::yNumTowers;
int acdReconAlg::xNumTopTiles;
int acdReconAlg::yNumTopTiles;
int acdReconAlg::numSideRows;


static float maxDOCA = 200.0;

static const AlgFactory<acdReconAlg>  Factory;
const IAlgFactory& acdReconAlgFactory = Factory;

acdReconAlg::acdReconAlg(const std::string& name, ISvcLocator* pSvcLocator) :
Algorithm(name, pSvcLocator) {

    declareProperty("tupleName",  m_tupleName="");

}

//------------------------------------------------------------------------------
StatusCode acdReconAlg::initialize ( )
{

    StatusCode sc = StatusCode::SUCCESS;

    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "initialize" << endreq;
    
    // Use the Job options service to set the Algorithm's parameters
    setProperties();

    sc = serviceLocator()->getService("GlastDetSvc",
        IID_IGlastDetSvc, reinterpret_cast<IInterface*&>( m_glastDetSvc));

    if( sc.isFailure() ) {
        log << MSG::ERROR << "acdReconAlg failed to get the GlastDetSvc" << endreq;
        return sc;
    }

    m_tileParams = new TileParams;

    getParameters();

    // get a pointer to our ntupleWriterSvc
    sc = serviceLocator()->getService("ntupleWriterSvc", 
        IID_INTupleWriterSvc, reinterpret_cast<IInterface*&>( ntupleWriteSvc));

    if( sc.isFailure() ) {
        log << MSG::ERROR << "acdReconAlg failed to get the ntupleWriterSvc" << endreq;
        return sc;
    }

    /*
    if( service("GismoSvc", m_gismoSvc).isFailure()) {
        log << MSG::ERROR << "failed to find the GismoSvc" << endreq;
        return StatusCode::FAILURE;
    }*/


    return sc;
}



//------------------------------------------------------------------------------
StatusCode acdReconAlg::execute() {
    
    StatusCode  sc = StatusCode::SUCCESS;
    MsgStream   log( msgSvc(), name() );
    log << MSG::DEBUG << "execute" << endreq;

    SmartDataPtr<TdGlastData> glastData(eventSvc(),"/Event/TdGlastData");

    if (glastData == 0) {
        log << MSG::DEBUG  << "No simulation detector data available " << endreq;
        return sc;
    }

    // retrieve the ACD data pointer
    m_AcdData = glastData->getVetoData();
    
    if (m_AcdData == 0) {
        log << MSG::DEBUG << "No ACD Data available " << endreq;
        return sc;
    }

    // run the reconstruction
    reconstruct(m_AcdData);


    log << MSG::DEBUG << "exiting ACD Recon " << endreq;

    return sc;
}


//------------------------------------------------------------------------------
/*
Finalize is called once at the end of event processing.
Any cleanup to occur - will occur here.
*/
StatusCode acdReconAlg::finalize() {
    
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "finalize" << endreq;
    
    return StatusCode::SUCCESS;
}


void acdReconAlg::getParameters ()
{

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

    xml::IFile *m_ifile = const_cast<xml::IFile*>(m_glastDetSvc->iniFile()); //OOPS!

    threshold_energy = m_ifile->getDouble("veto", "threshold");

    yNumTowers = m_ifile->getInt("glast", "yNum");
    xNumTowers = m_ifile->getInt("glast", "xNum");

    xNumTopTiles = m_ifile->getInt("veto", "numXtiles");
    yNumTopTiles = m_ifile->getInt("veto", "numYtiles");

    numSideRows = m_ifile->getInt("veto", "num_side_tiles");
        
}

void acdReconAlg::clear() {
    m_totEnergy = 0.0;
    m_tileCount = 0.0;
    m_gammaDOCA = maxDOCA;
    m_DOCA = maxDOCA;
    m_rowDOCA_vec.clear();
    m_rowDOCA_vec.resize(numSideRows+1, maxDOCA);  // one for each side, plus one for the top
    m_act_dist = -200.0;
}

StatusCode acdReconAlg::reconstruct (const IVetoData* v)
{

    StatusCode sc = StatusCode::SUCCESS;

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

    // reset all member variables to their defaults
    clear();

    // iterate through all hit ACD tiles
    for( IVetoData::const_iterator it= v->begin(); it != v->end(); ++it) {
        
        if ((*it).energy() < threshold_energy) continue; // toss out hits below threshold

        m_tileCount++;
        double tileEnergy = (*it).energy();
        m_totEnergy += tileEnergy;
        idents::AcdId id = (*it).type();
    }

    log << MSG::DEBUG << "num Tiles = " << m_tileCount << endreq;
    log << MSG::DEBUG << "total energy = " << m_totEnergy << endreq;

    m_glastDetSvc->accept(*m_tileParams);

    acdTileDOCA();

    sc = writeNTuple();

    return sc;
}


StatusCode acdReconAlg::acdTileDOCA() {

    StatusCode sc = StatusCode::SUCCESS;

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

    // Retrieve the TKR Reconstruction data
    SmartDataPtr<SiRecObjs> tkrRecData(eventSvc(),"/Event/TkrRecon/SiRecObjs");
    if (tkrRecData == 0) {
        log << MSG::DEBUG << "No TKR Reconstruction available " << endreq;
        return sc;
    }

    std::vector<const GFtrack*> xtracks;
    std::vector<const GFtrack*> ytracks;

    // Now retrieve all of the tracks

    int nparticles = tkrRecData->numParticles();
    if (nparticles > 0) {

        log << MSG::DEBUG << "Number of particles " << nparticles << endreq;

        unsigned int iParticle;
        for (iParticle = 0; iParticle < nparticles; iParticle++) {
            GFparticle *particle = tkrRecData->Particle(iParticle);
            xtracks.push_back(particle->getXGFtrack());
            ytracks.push_back(particle->getYGFtrack());
        }
    } else {
        log << MSG::DEBUG << "No reconstructed particles " << endreq;
    }

    int ngammas = tkrRecData->numGammas();
    log << MSG::DEBUG << "number of gammas = " << ngammas << endreq;

    if (ngammas > 0) {
        unsigned int iGamma;
        // Create a temporary vector to store the DOCA's for the 
        // top and side tiles - we don't want to use the reconstructed
        // gamma direction & vector for these DOCAs...
        std::vector<double> temprowDOCA_vec;
        temprowDOCA_vec.resize(numSideRows+1, maxDOCA); 
        for (iGamma = 0; iGamma < ngammas; iGamma++) {
            GFgamma *gamma = tkrRecData->Gamma(iGamma);
            Point gammaVertex = gamma->vertex();
            Vector gammaDirection = gamma->direction();
            double tempDOCA = DOCA(gammaVertex, gammaDirection, m_rowDOCA_vec);
            if (tempDOCA < m_gammaDOCA) m_gammaDOCA = tempDOCA;
            
            // Store the tracks associated with the gamma
            if (!gamma->getPair(SiCluster::X)->empty()) 
                xtracks.push_back(gamma->getPair(SiCluster::X));
            if (!gamma->getBest(SiCluster::X)->empty()) 
                xtracks.push_back(gamma->getBest(SiCluster::X));
            if (!gamma->getPair(SiCluster::Y)->empty())
                ytracks.push_back(gamma->getPair(SiCluster::Y)); 
            if (!gamma->getBest(SiCluster::X)->empty())
                ytracks.push_back(gamma->getBest(SiCluster::Y));
        }
    } else {
        log << MSG::DEBUG << "No reconstructed gammas " << endreq;
    }

    // Now we should have a list of x and y tracks
    // iterate through the list of tracks, and combine to create all possible
    // pairings of X and Y tracks.  The min DOCA is then reported as ACD_DOCA in
    // the ntuple

    for (std::vector<const GFtrack*>::const_iterator xtrk = xtracks.begin(); xtrk != xtracks.end(); xtrk++) {
        for(std::vector<const GFtrack*>::const_iterator ytrk = ytracks.begin(); ytrk != ytracks.end(); ytrk++) {
            Vector dir = Vector((*xtrk)->slope(), (*ytrk)->slope(), 1.).unit();
            Point x0((*xtrk)->positionAtZ(0), (*ytrk)->positionAtZ(0), 0.);

            float testDOCA = DOCA(x0, dir, m_rowDOCA_vec);
            if(testDOCA < m_DOCA) m_DOCA = testDOCA;
            
            float test_dist= hitTileDist(x0, dir);
            if(test_dist > m_act_dist) m_act_dist = test_dist;

        }
    }
   
    log << MSG::DEBUG << "ACD_DOCA = " << m_DOCA << endreq;
    log << MSG::DEBUG << "ACD_Act_dist = " << m_act_dist << endreq;

    return sc;

}

// Distance of Closest Approach (DOCA)
double acdReconAlg::DOCA (const Point &x0, const Vector &t0, std::vector<double> &doca_values)
{
    // This section looks for close-by hits to the ACD tiles

    float minDOCA = maxDOCA;
    float dist;

    // iterate over all tiles
    for( IVetoData::const_iterator it= m_AcdData->begin(); it != m_AcdData->end(); ++it) {

        float tile_energy = (*it).energy();
        if(tile_energy < threshold_energy) continue;

        Vector dX = (*it).position() - x0;

        float prod = dX * t0;
        dist = sqrt(dX.mag2() - prod*prod);
        if(dist < minDOCA){
            minDOCA = dist;
            m_hit_tile = *it;
        }
        
        idents::AcdId type = (*it).type();

        // Pick up the min. distance from each type of tile
        // i.e. top, and each type of side row tile
        if (type.top() && dist < doca_values[0]) doca_values[0] = dist;
        if (type.side()) {
            if (dist < doca_values[type.row()+1]) doca_values[type.row()+1] = dist;
        }
    }

    return minDOCA;
}

// Bill Atwood's new edge DOCA algorithm
double acdReconAlg::hitTileDist(const Point &x0, const Vector &t0)
{
    double return_dist = -200.;
    
    if(!m_AcdData) return return_dist;
    
    // iterate over all tiles
    for( IVetoData::const_iterator it= m_AcdData->begin(); it != m_AcdData->end(); ++it) {
        
        float eT = (*it).energy();
        if(eT < threshold_energy) continue;
        
        Point xT = (*it).position();
        idents::AcdId tileType((*it).type());
        int iFace = tileType.face();

        float dX = m_tileParams->length((*it).type());
        float dY = m_tileParams->width((*it).type());
        float dZ = m_tileParams->height((*it).type());
        
        // Figure out where in the plane of this face the trajectory hits
        double arc_dist = -1.; 
        if(iFace == 0) {// Top Tile. 
            arc_dist = (xT.z()-x0.z())/t0.z();                  
        }
        else if(iFace == 1 || iFace == 3) {// X Side Tile 
            arc_dist = (xT.x()-x0.x())/t0.x();
        }
        else if(iFace == 2 || iFace == 4) {// Y Side Tile
            arc_dist = (xT.y()-x0.y())/t0.y();
        }
        // If arc_dist is negative... had to go backwards to hit plane... 
        if(arc_dist < 0.) continue;
        
        Point x_isec = x0 + arc_dist*t0;
        
        Vector local_x0 = xT - x_isec; 
        if(iFace == 0) {// Top Tile
            double dist_x = dX/2. - fabs(local_x0.x());
            double dist_y = dY/2. - fabs(local_x0.y());                 
            double test_dist = (dist_x < dist_y) ? dist_x : dist_y;
            if(test_dist > return_dist) return_dist = test_dist;
        }
        else if(iFace == 1 || iFace == 3) {// X Side Tile
            double dist_z = dZ/2. - fabs(local_x0.z());
            double dist_y = dY/2. - fabs(local_x0.y());                 
            double test_dist = (dist_z < dist_y) ? dist_z : dist_y;
            if(test_dist > return_dist) return_dist = test_dist;
        }
        else if(iFace == 2 || iFace == 4) {// Y Side Tile
            double dist_z = dZ/2. - fabs(local_x0.z());
            double dist_x = dY/2. - fabs(local_x0.x());                 
            double test_dist = (dist_z < dist_x) ? dist_z : dist_x;
            if(test_dist > return_dist) return_dist = test_dist;
        }
    }
    return return_dist;
}


StatusCode acdReconAlg::writeNTuple() {

    StatusCode sc = StatusCode::SUCCESS;
    MsgStream   log( msgSvc(), name() );

    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_TotEnergy", m_totEnergy);
    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_TileCount", m_tileCount);
    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_DOCA", m_DOCA);
    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_GammaDOCA", m_gammaDOCA);
    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_TopDOCA", m_rowDOCA_vec[0]);
    sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_Act_Dist", m_act_dist);


    // loop over all DOCA values for the sides
    unsigned int iDOCA;
    for (iDOCA = 1; iDOCA <= numSideRows; iDOCA++) {
        std::string label("ACD_S");
        char num[5];
        _itoa((iDOCA-1), num, 10);
        label += num;
        label += "DOCA";
        sc = ntupleWriteSvc->addItem(m_tupleName.c_str(), label.c_str(), m_rowDOCA_vec[iDOCA]);
    }

    return sc;
}

---