In This Package:

SimHeaderCnvHistory.cc

Go to the documentation of this file.

#include "SimHeaderCnv.h"

#include "Event/SimParticleHistory.h"
#include "Event/SimTrack.h"
#include "Event/SimVertex.h"

#include "Event/GenHeader.h"

#include "PerSimEvent/PerSimParticleHistory.h"
#include "PerSimEvent/PerSimTrack.h"
#include "PerSimEvent/PerSimVertex.h"

using namespace DayaBay;
using namespace std;

DayaBay::SimParticleHistory* SimHeaderCnv::convert(vector<SimTrack*>& timap,
                                                   DayaBay::SimHeader& /*sh*/, 
                                                   const PerSimParticleHistory& pph)
{
    MsgStream log(msgSvc(), "SimHeaderCnvHistory::convert(p2t)");

    SimParticleHistory* ph = new SimParticleHistory;
    
    typedef vector<SimVertex*> VertexIndexMap;
    VertexIndexMap vimap;
    { // Vertices
        for (size_t ind=0; ind < pph.vertices.size(); ++ind) {
            PerSimVertex *pvertex = pph.vertices[ind];
            SimVertex* vertex = new SimVertex;
            SimProcess::Type type = (SimProcess::Type) pvertex->processType;
            vertex->setProcess(SimProcess(type,pvertex->processName));
            vertex->setTime(pvertex->time);
            vertex->setTotalEnergy(pvertex->energy);
            vertex->setPosition(pvertex->position);
            vertex->setMomentum(pvertex->momentum);
            ph->addVertex(vertex);
            vimap.push_back(vertex);
        }
    }

    { // Tracks
        for (size_t ind=0; ind < pph.tracks.size(); ++ind) {
            PerSimTrack* ptrack = pph.tracks[ind];
            SimTrack* track = new SimTrack;
            track->setTrackId(ptrack->geantID);
            track->setParentParticle(ptrack->parentPDG);
            track->setParticle(ptrack->particlePDG);
            SimTrack::vertex_list vertices;
            for (size_t vind=0; vind<ptrack->vertexIndex.size(); ++vind) {
                int index = ptrack->vertexIndex[vind];
                vertices.push_back(vimap[index]);
            }
            track->setVertices(vertices);

            SimVertex* vertex = 0;
            if (ptrack->ancestorVertex.index >= 0) 
                vertex = vimap[ptrack->ancestorVertex.index];
            track->setAncestorVertex(SimVertexReference(vertex,
                                                        ptrack->ancestorVertex.count));

            track->setUnrecordedDescendants(ptrack->unrecordedDescendants);

            timap.push_back(track);
            ph->addTrack(track);
        }
    }

    { // Fix up track references in vertices 
        
        for (size_t ind=0; ind < pph.vertices.size(); ++ind) {
            PerSimVertex *pvertex = pph.vertices[ind];
            SimVertex* vertex = vimap[ind];
            vertex->setTrack(SimTrackReference(timap[pvertex->track.index],
                                               pvertex->track.count));
            for (size_t stind=0; stind < pvertex->secondaryTracks.size(); ++stind) {
                PerSimIndirection& psi = pvertex->secondaryTracks[stind];
                vertex->addSecondary(SimTrackReference(timap[psi.index],psi.count));
            }
        }

    }

    { // Fix primary tracks
        const vector<int>& pti = pph.primaryTrackIndex;
        for (size_t ind=0; ind<pti.size(); ++ind) {
            const SimTrack* track = timap[ind];
            if (!track) {
                log << MSG::ERROR
                    << "Persistent primary track #" << ind << " not found in map" << endreq;
            }
            else {
                ph->addPrimaryTrack(track);
            }
        }
    }


    { // Fill track to track references
        for (size_t ind=0; ind < pph.tracks.size(); ++ind) {
            PerSimTrack* ptrack = pph.tracks[ind];
            SimTrack* track = timap[ind];

            int count = ptrack->ancestorTrack.count;
            SimTrack* ancestorTrack = 0;
            if (ptrack->ancestorTrack.index >= 0)
                ancestorTrack = timap[ptrack->ancestorTrack.index];

            track->setAncestorTrack(SimTrackReference(ancestorTrack,count));
        }
    }

    return ph;
}

PerSimParticleHistory* SimHeaderCnv::convert(map<const SimTrack*,int>& trackMap,
                                             const DayaBay::SimParticleHistory& ph)
{
    MsgStream log(msgSvc(), "SimHeaderCnvHistory::convert(t2p)");

    PerSimParticleHistory* pph = new PerSimParticleHistory();

    typedef map<const SimVertex*,int> VertexMap_t;
    VertexMap_t vertexMap;
    vertexMap[0] = -1;//Hopefully trigger a SegV to catch bogus values
    { // Vertices
        const list<SimVertex*>& vertices = ph.vertices();
        log << MSG::VERBOSE 
            << "Adding " << vertices.size() << " vertices" << endreq;
        list<SimVertex*>::const_iterator it, done = vertices.end();
        for (it=vertices.begin(); it != done; ++it) {
            SimVertex* vertex = *it;
            PerSimVertex* pvertex = new PerSimVertex;
            pvertex->processType = vertex->process().type();
            pvertex->processName = vertex->process().name();
            pvertex->time = vertex->time();
            pvertex->energy = vertex->totalEnergy();
            pvertex->position = vertex->position();
            pvertex->momentum = vertex->momentum();

            vertexMap[vertex] = pph->vertices.size();
            log << MSG::VERBOSE 
                << "Adding vertex from track " << vertex->track().track()->trackId()
                << " @ " << (void*) vertex << endreq;
            pph->vertices.push_back(pvertex);
        }
    }

    trackMap[0] = -1;//Hopefully trigger a SegV to catch any bogus values
    { // Tracks
        const list<SimTrack*>& tracks = ph.tracks();
        list<SimTrack*>::const_iterator it, done = tracks.end();
        for (it=tracks.begin(); it != done; ++it) {
            SimTrack* track = *it;
            PerSimTrack* ptrack = new PerSimTrack;
            ptrack->geantID = track->trackId();
            ptrack->parentPDG = track->parentParticle();
            ptrack->particlePDG = track->particle();
            for (size_t ind=0; ind<track->vertices().size(); ++ind) {
                VertexMap_t::iterator vit = vertexMap.find(track->vertices()[ind]);
                if (vit == vertexMap.end()) {
                    log << MSG::ERROR
                        << "Track references unknown vertex #" << ind << endreq;
                    continue;
                }
                ptrack->vertexIndex.push_back(vit->second);
            }

            const SimVertexReference& ancestorVertex = track->ancestorVertex();
            if (ancestorVertex.vertex()) {
                VertexMap_t::iterator vit = vertexMap.find(ancestorVertex.vertex());
                if (vit == vertexMap.end()) {
                    log << MSG::ERROR
                        << "Track's ancestor vertex is not found, is from track ID "
                        << ancestorVertex.vertex()->track().track()->trackId() 
                        << " and @ " << (void*) ancestorVertex.vertex()
                        << endreq;
                    // Historian uses non-default values for NULL
                    // SimVertexReferences so pass them along.
                    ptrack->ancestorVertex.count = ancestorVertex.indirection();
                }
                else {
                    ptrack->ancestorVertex = PerSimIndirection(vit->second,
                                                               ancestorVertex.indirection());
                }
            }
            else ptrack->ancestorVertex.count = ancestorVertex.indirection();

            // Rely on map<>'s default assignment
            ptrack->unrecordedDescendants = track->unrecordedDescendants();

            ptrack->primaryParticle = -1;
            if (track->primaryParticle()) {
                ptrack->primaryParticle = track->primaryParticle()->barcode();
            }

            trackMap[track] = pph->tracks.size();
            pph->tracks.push_back(ptrack);
        }
    }

    { // Fix up Vertices
        const list<SimVertex*>& vertices = ph.vertices();
        list<SimVertex*>::const_iterator it, done = vertices.end();
        for (it=vertices.begin(); it != done; ++it) {
            SimVertex* vertex = *it;
            PerSimVertex* pvertex = pph->vertices[vertexMap[vertex]];

            // Fill from SimTrackReference "track"
            if (vertex->track().track()) {
                int index = trackMap[vertex->track().track()];
                int count = vertex->track().indirection();
                pvertex->track = PerSimIndirection(index,count);
            }
            
            // Fill from SimTrackReferences "secondaries"
            const vector<SimTrackReference>& secondaries = vertex->secondaries();
            for (size_t ind=0; ind < secondaries.size(); ++ind) {
                if (!secondaries[ind].track()) continue;
                int index = trackMap[secondaries[ind].track()];
                int count = secondaries[ind].indirection();
                pvertex->secondaryTracks.push_back(PerSimIndirection(index,count));
            }
        }
    }

    { // Fill primary track indices
        const std::list<const SimTrack*>& primaries = ph.primaryTracks();
        std::list<const SimTrack*>::const_iterator it, done = primaries.end();
        for (it = primaries.begin(); it != done; ++it) {
            const SimTrack* track = *it;
            if (!track) continue;
            pph->primaryTrackIndex.push_back(trackMap[track]);
        }
    }

    { // Fix up missing track->ancestorTrack references
        const list<SimTrack*>& tracks = ph.tracks();
        list<SimTrack*>::const_iterator it, done = tracks.end();
        for (it=tracks.begin(); it != done; ++it) {
            SimTrack* track = *it;
            PerSimTrack* ptrack = pph->tracks[trackMap[track]];

            int count = track->ancestorTrack().indirection();
            const SimTrack* ancestorTrack = track->ancestorTrack().track();
            if (!ancestorTrack) {
                ptrack->ancestorTrack.count = count;
            }
            else {
                ptrack->ancestorTrack = PerSimIndirection(trackMap[ancestorTrack],count);
            }
        }
    }

    // Note: DO NOT fill idToTrack lookup.  It is an internal lookup
    // that will be refilled based on the the geantID stored in each
    // track....

    return pph;
}

void SimHeaderCnv::relate(DayaBay::SimHeader& sh)
{
    MsgStream log(msgSvc(), "SimHeaderCnvHistory::relate(p2t)");
    const HepMC::GenEvent* genEvent = 0;
    const std::vector<const DayaBay::IHeader*>& inputHeaders = sh.inputHeaders();
    log << MSG::DEBUG 
        << "Checking for GenHeader in " 
        << inputHeaders.size() << " input headers" 
        << endreq;
    for (size_t ind=0; ind<inputHeaders.size(); ++ind) {
        const GenHeader* genHeader = dynamic_cast<const GenHeader*>(inputHeaders[ind]);
        if (1) {
            const HeaderObject* ho = dynamic_cast<const HeaderObject*>(inputHeaders[ind]);
            log << MSG::DEBUG
                << "checking " << ind << "/" << inputHeaders.size() 
                << ": " << ho->name()
                << endreq;
        }
        if (!genHeader) continue;
        genEvent = genHeader->event();
    }
    if (!genEvent) {
        log << MSG::DEBUG
            << "Could not find the HepMC::GenEvent used to create me." 
            << endreq;
        return;
    }
    

    if (!m_perInObj->history) {
        log << MSG::DEBUG
            << "No persistent particle history"
            << endreq;
        return;
    }
    const PerSimParticleHistory& pph = *(m_perInObj->history);
    if (!sh.particleHistory()) {
        log << MSG::ERROR
            << "No transient particle history"
            << endreq;
        return;
    }

    // Sometimes the genEvent barcodes for particles and vertices has been pruned
    // in the production of the original file. Report this here and
    // suppress issuing an ERROR for every particle that does not have a barcodes.
    bool empty = genEvent->particles_empty() && genEvent->vertices_empty() ;
    if (empty) log << MSG::INFO << "No GenEvent particles or vertices exist"<< endreq;


    const list<SimTrack*>& tracks = sh.particleHistory()->tracks();
    list<SimTrack*>::const_iterator it, done = tracks.end();
    size_t ind=0;
    for (it=tracks.begin(); it !=done && ind < pph.tracks.size(); ++ind, ++it) {
        PerSimTrack* ptrack = pph.tracks[ind];
        SimTrack* track = *it;

        if (genEvent) {
            int barcode = ptrack->primaryParticle;

            if (barcode < 0) {
                track->setPrimaryParticle(0);
            }
            else {
                const HepMC::GenParticle* part = genEvent->barcode_to_particle(barcode);
                if (part) {
                    track->setPrimaryParticle(part);
                }
                else {
                  if (!empty){
                    log << MSG::ERROR
                        << "Failed to get a primary particle for barcode " << barcode << endreq;
                  }
                }
            }
        }
    }
}

---