---

AcdReconDigiAlg.cxx

Go to the documentation of this file.

#include "AcdReconDigiAlg.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 "TkrRecon/SiRecObjs.h"
#include "TkrRecon/GFparticle.h"
#include "TkrRecon/GFdata.h"
#include "TkrRecon/GFgamma.h"

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

#include "idents/AcdId.h"
#include "xml/IFile.h"

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

static float maxDOCA = 200.0;

static const AlgFactory<AcdReconDigiAlg>  Factory;
const IAlgFactory& AcdReconDigiAlgFactory = Factory;

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

    declareProperty("tupleName",  m_tupleName="");
    declareProperty("xmlFile", m_xmlFileStr="");
    declareProperty("vetoThreshFile", m_vetoThreshFileStr="");

}

//------------------------------------------------------------------------------
StatusCode AcdReconDigiAlg::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 << "AcdReconDigiAlg failed to get the GlastDetSvc" << endreq;
        return sc;
    }

    getParameters();

    readThresholds();

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

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

    return sc;
}



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

    DataObject *tempObj;
    sc = eventSvc()->retrieveObject("/Event/Digi/AcdDigis", tempObj);
    if (sc.isFailure()) {
        log << MSG::DEBUG << "No ACD Data available" << endreq;
        return StatusCode::SUCCESS;
    }

    m_AcdData = dynamic_cast<AcdDigiVector*>(tempObj);
    
    if (m_AcdData == 0) {
        log << MSG::DEBUG << "No ACD Data available " << endreq;
        return StatusCode::SUCCESS;
    }

    // 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 AcdReconDigiAlg::finalize() {
    
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "finalize" << endreq;
    
    return StatusCode::SUCCESS;
}


void AcdReconDigiAlg::getParameters ()
{

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

    xml::IFile *m_ifile = new xml::IFile(m_xmlFileStr.c_str()); 

    std::vector<int> topIds = m_ifile->getIntVector("acdTop", "ids");
    std::vector<double> topX = m_ifile->getDoubleVector("acdTop", "x_center");
    std::vector<double> topY = m_ifile->getDoubleVector("acdTop", "y_center");
    std::vector<double> topZ = m_ifile->getDoubleVector("acdTop", "z_center");
    std::vector<double> topXlen = m_ifile->getDoubleVector("acdTop", "xlen");
    std::vector<double> topYlen = m_ifile->getDoubleVector("acdTop", "ylen");
    std::vector<double> topZlen = m_ifile->getDoubleVector("acdTop", "zlen");

    // Load the ids, centers and dimensions for the top tiles
    for (int i=0; i < topIds.size(); i++)  { 
        m_tileCenters[topIds[i]] = Point(topX[i], topY[i], topZ[i]);
        m_tileDims[topIds[i]] = Vector(topXlen[i], topYlen[i], topZlen[i]);
    }

    // Load the id and center of the bigtile
    int id = m_ifile->getInt("bigtile", "id");
    double bigX = m_ifile->getDouble("bigtile", "x_center");
    double bigY = m_ifile->getDouble("bigtile", "y_center");
    double bigZ = m_ifile->getDouble("bigtile", "z_center");
    double bigXlen = m_ifile->getDouble("bigtile", "xlen");
    double bigYlen = m_ifile->getDouble("bigtile", "ylen");
    double bigZlen = m_ifile->getDouble("bigtile", "zlen");

    m_tileCenters[id] = Point(bigX, bigY, bigZ);
    m_tileDims[id] = Vector(bigXlen, bigYlen, bigZlen);

    // Load the ids and centers for the side tiles
    std::vector<int> sideIds = m_ifile->getIntVector("acdSide", "ids");
    std::vector<double> sideX = m_ifile->getDoubleVector("acdSide", "x_center");
    std::vector<double> sideY = m_ifile->getDoubleVector("acdSide", "y_center");
    std::vector<double> sideZ = m_ifile->getDoubleVector("acdSide", "z_center");
    std::vector<double> sideXlen = m_ifile->getDoubleVector("acdSide", "xlen");
    std::vector<double> sideYlen = m_ifile->getDoubleVector("acdSide", "ylen");
    std::vector<double> sideZlen = m_ifile->getDoubleVector("acdSide", "zlen");

    for (int j=0; j < sideIds.size(); j++) {
        m_tileCenters[sideIds[j]] = Point(sideX[j], sideY[j], sideZ[j]);
        m_tileDims[sideIds[j]] = Vector(sideXlen[j], sideYlen[j], sideZlen[j]);
    }

    std::vector<int> xgtIds = m_ifile->getIntVector("xgt", "ids");
    std::vector<double> xgtX = m_ifile->getDoubleVector("xgt", "x_center");
    std::vector<double> xgtY = m_ifile->getDoubleVector("xgt", "y_center");
    std::vector<double> xgtZ = m_ifile->getDoubleVector("xgt", "z_center");

    for (int k = 0; k < xgtIds.size(); k++) 
        m_tileCenters[xgtIds[k]] = Point(xgtX[k], xgtY[k], xgtZ[k]);

    for( std::map<int,Point>::const_iterator it= m_tileCenters.begin(); it != m_tileCenters.end(); ++it) {
        Point cen = (*it).second;
        log << MSG::DEBUG << "id = " << (*it).first << " center = (" << cen.x() << "," << cen.y() << "," << cen.z() << ")" << endreq; 
    }

    m_numSideRows = m_ifile->getInt("acd", "num_side_tiles");

    log << MSG::DEBUG << "number of ACD rows = " << m_numSideRows << endreq;


}

void AcdReconDigiAlg::readThresholds() {

    xml::IFile::extractEnvVar(&m_vetoThreshFileStr);
    std::ifstream threshFile( m_vetoThreshFileStr.c_str() );
    
    if (threshFile.eof() || !threshFile.good()) return;
    int id = 0;
    double thresh = 0.0;

    threshFile >> id >> thresh;

    while ( ( !threshFile.eof() ) && ( id >= 0 )) {
        m_vetoThresh[id] = thresh;
        threshFile >> id >> thresh;
    }

    threshFile.close();
}


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

}

StatusCode AcdReconDigiAlg::reconstruct (const AcdDigiVector* v)
{

    StatusCode sc = StatusCode::SUCCESS;

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

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

    // iterate through all hit ACD tiles
    for( AcdDigiVector::const_iterator it= v->begin(); it != v->end(); ++it) {
        

        idents::AcdId id = (*it)->ID();
        
        //if ((*it)->Veto() == false) continue; // toss out hits below veto threshold
        // Use the pulse height and the veto thresholds to determine whether to continue
        if ((*it)->PulseHeight() < m_vetoThresh[id.id()]) continue;
        
        m_tileCount++;

    }

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

    acdTileDOCA();

    sc = writeNTuple();

    return sc;
}


StatusCode AcdReconDigiAlg::acdTileDOCA() {

    StatusCode sc = StatusCode::SUCCESS;

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

    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;

    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(m_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
            xtracks.push_back(gamma->getPair(SiCluster::X));
            xtracks.push_back(gamma->getBest(SiCluster::X));
            ytracks.push_back(gamma->getPair(SiCluster::Y));            
            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;

    return sc;

}

// Distance of Closest Approach (DOCA)
double AcdReconDigiAlg::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 of the tiles
    for( AcdDigiVector::const_iterator it= m_AcdData->begin(); it != m_AcdData->end(); ++it) {

        idents::AcdId type = (*it)->ID();

        //      if((*it)->Veto() == false) continue;
        // Use the pulse height and the veto thresholds to determine whether to continue
        if ((*it)->PulseHeight() < m_vetoThresh[type.id()]) continue;

        // Need to grab the center position of this tile
        Point tileCenter = m_tileCenters[type.id()];
        Vector dX = tileCenter - x0;

        float prod = dX * t0;
        dist = sqrt(dX.mag2() - prod*prod);
        if(dist < minDOCA){
            minDOCA = dist;
        }
        

        // 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 AcdReconDigiAlg::hitTileDist(const Point &x0, const Vector &t0)
{
    double return_dist = -200.;
    
    if(!m_AcdData) return return_dist;
    
    // iterate over all of the tiles
    for( AcdDigiVector::const_iterator it= m_AcdData->begin(); it != m_AcdData->end(); ++it) {


        idents::AcdId type = (*it)->ID();
        // Use the pulse height and the veto thresholds to determine whether to continue
        if ((*it)->PulseHeight() < m_vetoThresh[type.id()]) continue;
        //if((*it)->Veto() == false) continue;

        // skip XGTs for now
        if (type.layer() > 1) continue;

        int iFace = type.face();

        // Need to grab the center position of this tile
        Point xT = m_tileCenters[type.id()];
        
        float dX = m_tileDims[type.id()].x();
        float dY = m_tileDims[type.id()].y();
        float dZ = m_tileDims[type.id()].z();
        
        // 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 AcdReconDigiAlg::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 <= m_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;
}

---