Using RAT::TrackNav to find primaries associated with PMT hit photons

rat_backtrack.cpp

#include <iostream>
#include <string>
#include <assert.h>

#include <RAT/DSReader.hh>
#include <RAT/DS/Root.hh>
#include <RAT/TrackNav.hh>
#include <RAT/TrackCursor.hh>
#include <RAT/TrackNode.hh>

/**
 *  RAT::TrackNav provides a user interface for navigating through stored
 *  GEANT4 tracks.
 *
 *  Here, we use this class to find each detected photon's parent primary.
 *
 *  To compile:
 *      g++ -o backtrack backtrack.cpp `root-config --libs` \
 *      -I$RATROOT/include -L$RATROOT/lib -I$ROOTSYS/include \
 *      -lRATEvent_Linux-g++
 *
 *  A. Mastbaum (mastbaum@hep.upenn.edu), 6 August 2011 
 */

int main(int argc, char* argv[])
{
    std::string filename = argv[1];
    RAT::DSReader dsr(filename.c_str());
    RAT::DS::Root* ds;

    while(ds = dsr.NextEvent()) {
        RAT::DS::MC* mc = ds->GetMC();
        RAT::TrackNav nav(mc, false);
        RAT::TrackCursor c = nav.Cursor(false);                
        std::cout << "ds.mc has " <<  mc->GetMCParticleCount() << " primaries\n";

        // loop over mcpmts
        for(int i=0; i<mc->GetMCPMTCount(); i++) {
            RAT::DS::MCPMT* pmt = mc->GetMCPMT(i);

            // loop over hit photons on this mcpmt
            for(int j=0; j<pmt->GetMCPhotonCount(); j++) {
                // geant4 track id of hit photon is stored
                long track_id = pmt->GetMCPhoton(j)->GetTrackID();
                std::cout << "pmt " << i << ": "
                          << "photon " << j << " "
                          << "has track id " << track_id << ", ";
                if(track_id > 1e8) {
                    std::cout << "invalid track id " << track_id << std::endl;
                    continue;
                }

                // use rat track navigator to find the start of the track
                RAT::TrackNode* node = nav.FindID(track_id);
                c.Go(node);

                std::cout << "parent: "  << c.Parent()->GetParticleName() << std::endl;
            }
        }
    }

    return 0;
}