EDepSim::PersistencyManager Class Reference

#include <EDepSimPersistencyManager.hh>

Inheritance diagram for EDepSim::PersistencyManager:

Public Member Functions
	PersistencyManager ()
virtual	~PersistencyManager ()
virtual G4bool	Store (const G4Event *anEvent)
	stores anEvent and the associated objects into database. More...
virtual G4bool	Store (const G4Run *aRun)
virtual G4bool	Store (const G4VPhysicalVolume *aWorld)
virtual G4bool	Retrieve (G4Event *&e)
	Retrieve information from a file. These are not implemented. More...
virtual G4bool	Retrieve (G4Run *&r)
virtual G4bool	Retrieve (G4VPhysicalVolume *&w)
const TG4Event &	GetEventSummary ()
const std::vector< TG4PrimaryVertex > &	GetPrimaries () const
const std::vector< TG4Trajectory > &	GetTrajectories () const
const TG4HitSegmentDetectors &	GetSegmentDetectors () const
virtual G4bool	Open (G4String filename)
virtual G4bool	Close (void)
	Make sure the output file is closed. More...
virtual G4String	GetFilename (void) const
	Return the output file name. More...
virtual void	SetLengthThreshold (G4double thresh)
virtual G4double	GetLengthThreshold () const
	Get the threshold for length in a sensitive detector. More...
virtual void	SetGammaThreshold (G4double thresh)
virtual G4double	GetGammaThreshold () const
virtual void	SetNeutronThreshold (G4double thresh)
virtual G4double	GetNeutronThreshold () const
virtual void	SetTrajectoryPointAccuracy (double acc)
virtual double	GetTrajectoryPointAccuracy (void) const
virtual void	SetTrajectoryPointDeposit (double dep)
virtual double	GetTrajectoryPointDeposit (void) const
virtual void	SetSaveAllPrimaryTrajectories (bool val)
virtual bool	GetSaveAllPrimaryTrajectories (void) const
virtual void	AddTrajectoryBoundary (const G4String &boundary)
virtual void	ClearTrajectoryBoundaries ()
void	SetDetectorPartition (int partition)
	Set the detector mask. More...
int	GetDetectorPartition () const
	Get the detector partition. More...

Protected Member Functions
void	SetFilename (G4String file)
void	UpdateSummaries (const G4Event *event)
	Update the event summary fields. More...

Protected Attributes
TG4Event	fEventSummary
	A summary of the primary vertices in the event. More...

Private Types
typedef std::map< int, int >	TrackIdMap

Private Member Functions
void	SummarizePrimaries (std::vector< TG4PrimaryVertex > &primaries, const G4PrimaryVertex *event)
void	SummarizeTrajectories (std::vector< TG4Trajectory > &trajectories, const G4Event *event)
	Fill the trajectory container. More...
void	MarkTrajectories (const G4Event *event)
	Mark the G4 Trajectories that should be saved. More...
void	SummarizeSegmentDetectors (TG4HitSegmentDetectors &segmentDetectors, const G4Event *event)
	sensitive detector. More...
void	SummarizeHitSegments (std::vector< TG4HitSegment > &segments, G4VHitsCollection *hits)
	Fill a container of hit segments. More...
void	CopyHitContributors (std::vector< int > &dest, const std::vector< int > &src)
void	CopyTrajectoryPoints (TG4Trajectory &traj, G4VTrajectory *g4Traj)
double	FindTrajectoryAccuracy (G4VTrajectory *traj, int point1, int point2)
int	SplitTrajectory (G4VTrajectory *traj, int point1, int point2)
void	SelectTrajectoryPoints (std::vector< int > &selected, G4VTrajectory *g4Traj)
bool	SaveTrajectoryBoundary (G4VTrajectory *g4Traj, G4StepStatus status, G4String currentVolume, G4String prevVolume)

Private Attributes
TrackIdMap	fTrackIdMap
G4String	fFilename
	The filename of the output file. More...
G4double	fLengthThreshold
G4double	fGammaThreshold
G4double	fNeutronThreshold
	The threshold momentum below which a neutron will be rejected. More...
double	fTrajectoryPointAccuracy
double	fTrajectoryPointDeposit
bool	fSaveAllPrimaryTrajectories
std::vector< TPRegexp * >	fTrajectoryBoundaries
int	fDetectorPartition
	The detector partition. More...
EDepSim::PersistencyMessenger *	fPersistencyMessenger

Detailed Description

Definition at line 44 of file EDepSimPersistencyManager.hh.

Member Typedef Documentation

typedef std::map<int,int> EDepSim::PersistencyManager::TrackIdMap

private

The mapping between the internal G4 TrackID number and the external TrackId number. The G4 TrackID value is based on the order that the particle is placed onto the stack, not all TrackID value exist and the ordering in the vector is not monotonic. The external TrackId is the same as the location of the trajectory in the output array and can be used as an index so parent trajectories are "easy" to find. This is filled in SummarizeTrajectories, and used to remap TrackId in other methods.

Definition at line 274 of file EDepSimPersistencyManager.hh.

Constructor & Destructor Documentation

EDepSim::PersistencyManager::PersistencyManager

(

)

Creates a root persistency manager. Through the "magic" of G4VPersistencyManager the ultimate base class, this declared to the G4 persistency management system. You can only have one active persistency class at any give moment.

Definition at line 45 of file EDepSimPersistencyManager.cc.

    : G4VPersistencyManager(), fFilename("/dev/null"),
      fLengthThreshold(10*mm),
      fGammaThreshold(5*MeV), fNeutronThreshold(50*MeV),
      fTrajectoryPointAccuracy(1.*mm), fTrajectoryPointDeposit(0*MeV),
      fSaveAllPrimaryTrajectories(true) {
    fPersistencyMessenger = new EDepSim::PersistencyMessenger(this);
}

EDepSim::PersistencyManager::~PersistencyManager ( )

virtual

Definition at line 55 of file EDepSimPersistencyManager.cc.

                                               {
    ClearTrajectoryBoundaries();
    delete fPersistencyMessenger;
}

Member Function Documentation

void EDepSim::PersistencyManager::AddTrajectoryBoundary ( const G4String &  boundary )

virtual

Add a regular expression to define a volume boundary at which a trajectory point should be saved. The volume names are defined by the constructors which built the volume. Most of the names can be found by looking at the geant command list. The names are not the same as the root volume names.

Definition at line 94 of file EDepSimPersistencyManager.cc.

                                                                       {
    TPRegexp* bound = new TPRegexp(b.c_str());
    fTrajectoryBoundaries.push_back(bound);
}

void EDepSim::PersistencyManager::ClearTrajectoryBoundaries ( )

virtual

Clear the list of volume boundaries which will cause a trajectory point to be saved.

Definition at line 99 of file EDepSimPersistencyManager.cc.

                                                          {
    for (std::vector<TPRegexp*>::iterator r = fTrajectoryBoundaries.begin();
         r != fTrajectoryBoundaries.end();
         ++r) {
        delete (*r);
    }
    fTrajectoryBoundaries.clear();
}

G4bool EDepSim::PersistencyManager::Close ( void  )

virtual

Make sure the output file is closed.

Make sure the output file is closed. This is used to make sure that any information being summarized has been saved.

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 67 of file EDepSimPersistencyManager.cc.

                                            {
    EDepSimSevere(" -- Close is not implimented.");
    return false;
}

void EDepSim::PersistencyManager::CopyHitContributors ( std::vector< int > &  dest, const std::vector< int > &  src  )

private

Fill the map of sensitive detectors which use hit segments as the Copy and map the contributing trajectories from the vector of contributors for the EDepSim::HitSegment into the vector of contributors for the TG4HitSegment.

Definition at line 537 of file EDepSimPersistencyManager.cc.

                                                                            {

    dest.clear();

    for (std::vector<int>::const_iterator c = src.begin();
         c != src.end(); ++c) {
        // Check each contributor to make sure that it is a valid
        // trajectory.  If it isn't in the trajectory map, then set it
        // to a parent that is.
        EDepSim::Trajectory* ndTraj
            = dynamic_cast<EDepSim::Trajectory*>(
                EDepSim::TrajectoryMap::Get(*c));
        while (ndTraj && !ndTraj->SaveTrajectory()) {
            ndTraj = dynamic_cast<EDepSim::Trajectory*>(
                EDepSim::TrajectoryMap::Get(ndTraj->GetParentID()));
        }
        if (!ndTraj) {
            dest.push_back(-1);
            continue;
        }
        if (fTrackIdMap.find(ndTraj->GetTrackID()) != fTrackIdMap.end()) {
            dest.push_back(fTrackIdMap[ndTraj->GetTrackID()]);
        }
        else {
            EDepSimError("Contributor with unknown trajectory: "
                      << ndTraj->GetTrackID());
            dest.push_back(-2);
        }
    }

    // Remove the duplicate entries.
    std::sort(dest.begin(),dest.end());
    dest.erase(std::unique(dest.begin(),dest.end()),dest.end());
}

void EDepSim::PersistencyManager::CopyTrajectoryPoints ( TG4Trajectory &  traj, G4VTrajectory *  g4Traj  )

private

Copy the trajectory points into the trajectory summary. This uses SelectTrajectoryPoints to pick out the points that should appear in the summary.

Definition at line 452 of file EDepSimPersistencyManager.cc.

                                                                         {
    std::vector<int> selected;

    // Choose the trajectory points that are going to be saved.
    SelectTrajectoryPoints(selected, g4Traj);

    // Make sure the selected trajectory points are in order and unique.
    std::sort(selected.begin(),selected.end());
    selected.erase(std::unique(selected.begin(), selected.end()),
                   selected.end());

    ////////////////////////////////////
    // Save the trajectories.
    ////////////////////////////////////
    for (std::vector<int>::iterator tp = selected.begin();
         tp != selected.end(); ++tp) {
        EDepSim::TrajectoryPoint* edepPoint
            = dynamic_cast<EDepSim::TrajectoryPoint*>(g4Traj->GetPoint(*tp));
        TG4TrajectoryPoint point;
        point.Position.SetXYZT(edepPoint->GetPosition().x(),
                               edepPoint->GetPosition().y(),
                               edepPoint->GetPosition().z(),
                               edepPoint->GetTime());
        point.Momentum.SetXYZ(edepPoint->GetMomentum().x(),
                              edepPoint->GetMomentum().y(),
                              edepPoint->GetMomentum().z());
        point.Process = edepPoint->GetProcessType();
        point.Subprocess = edepPoint->GetProcessSubType();
        traj.Points.push_back(point);
    }
}

double EDepSim::PersistencyManager::FindTrajectoryAccuracy ( G4VTrajectory *  traj, int  point1, int  point2  )

private

Find the maximum deviation of a trajectory point from the interpolated path between point 1 and point 2

Definition at line 573 of file EDepSimPersistencyManager.cc.

                                                   {
    if ((point2-point1) < 2) return 0;

    G4ThreeVector p1 = g4Traj->GetPoint(point1)->GetPosition();
    G4ThreeVector p2 = g4Traj->GetPoint(point2)->GetPosition();

    if ( (p2-p1).mag() < fTrajectoryPointAccuracy) return 0;

    G4ThreeVector dir = (p2-p1).unit();

    int step = (point2-point1)/10 + 1;

    double approach = 0.0;
    for (int p = point1+1; p<point2; p = p + step) {
        p2 = g4Traj->GetPoint(p)->GetPosition() - p1;
        approach = std::max((p2 - (dir*p2)*dir).mag(), approach);
    }

    return approach;
}

int EDepSim::PersistencyManager::GetDetectorPartition ( ) const

inline

Get the detector partition.

Definition at line 194 of file EDepSimPersistencyManager.hh.

{return fDetectorPartition;}

const TG4Event& EDepSim::PersistencyManager::GetEventSummary

(

)

A public accessor to the summarized event. The primaries are summarized during by a call to UpdateSummaries. The EDepSim::PersistencyManager::Store(event) method is called from G4RunManager::AnalyzeEvent and will call the UpdateSummaries() method. Alternatively, the UpdateSummarize method can be called by a method of a derived class (e.g. in its Store method). The fEventSummary field is protected so derived classes can directly access it.

virtual G4String EDepSim::PersistencyManager::GetFilename ( void  ) const

inlinevirtual

Return the output file name.

Definition at line 98 of file EDepSimPersistencyManager.hh.

{return fFilename;}

virtual G4double EDepSim::PersistencyManager::GetGammaThreshold ( ) const

inlinevirtual

Get the momentum threshold required to save a gamma-ray as a trajectory. Gamma rays are rejected if the gamma-ray momentum is below either GetGammaThreshold() or GetMomentumThreshold().

Definition at line 119 of file EDepSimPersistencyManager.hh.

{return fGammaThreshold;}

virtual G4double EDepSim::PersistencyManager::GetLengthThreshold ( ) const

inlinevirtual

Get the threshold for length in a sensitive detector.

Definition at line 108 of file EDepSimPersistencyManager.hh.

{return fLengthThreshold;}

virtual G4double EDepSim::PersistencyManager::GetNeutronThreshold ( ) const

inlinevirtual

Get the momentum threshold required to save a neutron as a trajectory.

Definition at line 129 of file EDepSimPersistencyManager.hh.

{return fNeutronThreshold;}

const std::vector<TG4PrimaryVertex>& EDepSim::PersistencyManager::GetPrimaries

(

)

const

A public accessor to the summarized primaries. The primaries are summarized during by a call to UpdateSummaries. The EDepSim::PersistencyManager::Store(event) method is called from G4RunManager::AnalyzeEvent and will call the UpdateSummaries() method. Alternatively, the UpdateSummarize method can be called by a method of a derived class (e.g. in its Store method).

virtual bool EDepSim::PersistencyManager::GetSaveAllPrimaryTrajectories ( void  ) const

inlinevirtual

Get the flag to save primary particle trajectories. If this flag is true, then the trajectories for primary particles are saved even if they don't ultimately cause an energy deposition in a sensitive detector. If this flag is false, the trajectories for primary particles are only saved if they, or one of their children, deposit energy in a sensitive detector.

Definition at line 175 of file EDepSimPersistencyManager.hh.

                                                           {
        return fSaveAllPrimaryTrajectories;
    }

const TG4HitSegmentDetectors& EDepSim::PersistencyManager::GetSegmentDetectors

(

)

const

A public accessor to the summarized hit segment detectors. See the documentation for the GetPrimaries() method.

const std::vector<TG4Trajectory>& EDepSim::PersistencyManager::GetTrajectories

(

)

const

A public accessor to the summarized trajectories. See the documentation for the GetPrimaries() method.

virtual double EDepSim::PersistencyManager::GetTrajectoryPointAccuracy ( void  ) const

inlinevirtual

Get the required trajectory accuracy. Trajectory points are saved so that the interpolated position of a trajectory is never more than this distance from the unsaved points.

Definition at line 141 of file EDepSimPersistencyManager.hh.

                                                          {
        return fTrajectoryPointAccuracy;
    }

virtual double EDepSim::PersistencyManager::GetTrajectoryPointDeposit ( void  ) const

inlinevirtual

Get the minimum energy deposition for which a trajectory point will be saved. This is the trajectory point process deposition (GetProcessDeposit()).

Definition at line 155 of file EDepSimPersistencyManager.hh.

                                                         {
        return fTrajectoryPointDeposit;
    }

void EDepSim::PersistencyManager::MarkTrajectories ( const G4Event *  event )

private

Mark the G4 Trajectories that should be saved.

Definition at line 324 of file EDepSimPersistencyManager.cc.

                                                                     {
    const G4TrajectoryContainer* trajectories = event->GetTrajectoryContainer();
    if (!trajectories) {
        EDepSimVerbose("No Trajectories ");
        return;
    }

    // Go through all of the trajectories and save:
    //
    //   ** Trajectories from primary particles if
    //       1) a daughter deposited energy in a sensitive detector
    //       2) or, SaveAllPrimaryTrajectories() is true
    //
    //   ** Trajectories created by a particle decay if
    //       1) a daughter deposited energy in a sensitve detector
    //       2) or, SaveAllPrimaryTrajectories() is true.
    //
    //   ** Charged particle trajectories that pass through a sensitive
    //         detector.
    //
    //   ** Gamma-rays and neutrons above a threshold which have daughters
    //         depositing energy in a sensitive detector.
    //
    for (TrajectoryVector::iterator t = trajectories->GetVector()->begin();
         t != trajectories->GetVector()->end();
         ++t) {
        EDepSim::Trajectory* ndTraj = dynamic_cast<EDepSim::Trajectory*>(*t);
        std::string particleName = ndTraj->GetParticleName();
        std::string processName = ndTraj->GetProcessName();
        double initialMomentum = ndTraj->GetInitialMomentum().mag();

        // Check if all primary particle trajectories should be saved.  The
        // primary particle should always be saved if it, or any of it's
        // children, deposited energy in a sensitive detector.  If the primary
        // didn't deposit energy in a sensitive detector, then it will only be
        // saved if SaveAllPrimaryTrajectories is true.
        if (ndTraj->GetParentID() == 0) {
            if (ndTraj->GetSDTotalEnergyDeposit()>1*eV
                || GetSaveAllPrimaryTrajectories()) {
                ndTraj->MarkTrajectory();
                continue;
            }
        }

        // Don't save the neutrinos
        if (particleName == "anti_nu_e") continue;
        if (particleName == "anti_nu_mu") continue;
        if (particleName == "anti_nu_tau") continue;
        if (particleName == "nu_e") continue;
        if (particleName == "nu_mu") continue;
        if (particleName == "nu_tau") continue;

        // Save any decay product if it caused any energy deposit.
        if (processName == "Decay") {
            if (ndTraj->GetSDTotalEnergyDeposit()>1*eV
                || GetSaveAllPrimaryTrajectories()) {
                ndTraj->MarkTrajectory(false);
                continue;
            }
        }

        // Save particles that produce charged track inside a sensitive
        // detector.  This doesn't automatically save, but the parents will be
        // automatically considered for saving by the next bit of code.
        if (ndTraj->GetSDLength() > GetLengthThreshold()) {
            ndTraj->MarkTrajectory(false);
            continue;
        }

        // For the next part, only consider particles where the children have
        // deposited energy in a sensitive detector.
        if (ndTraj->GetSDTotalEnergyDeposit()<1*eV) {
            continue;
        }

        // Save higher energy gamma rays that have descendents depositing
        // energy in a sensitive detector.  This only affects secondary
        // photons since primary photons are handled above.
        if (particleName == "gamma" && initialMomentum > GetGammaThreshold()) {
            ndTraj->MarkTrajectory(false);
            continue;
        }

        // Save higher energy neutrons that have descendents depositing energy
        // in a sensitive detector.  This only affects secondary neutrons
        // since primary neutrons are controlled above.
        if (particleName == "neutron"
            && initialMomentum > GetNeutronThreshold()) {
            ndTraj->MarkTrajectory(false);
            continue;
        }
    }

    // Go through all of the event hit collections and make sure that all
    // primary trajectories and trajectories contributing to a hit are saved.
    // These are mostly a sub-set of the trajectories marked in the previous
    // step, but there are a few corner cases where trajectories are not saved
    // because of theshold issues.
    G4HCofThisEvent* hitCollections = event->GetHCofThisEvent();
    if (!hitCollections) return;
    for (int i=0; i < hitCollections->GetNumberOfCollections(); ++i) {
        G4VHitsCollection* g4Hits = hitCollections->GetHC(i);
        if (g4Hits->GetSize()<1) continue;
        for (unsigned int h=0; h<g4Hits->GetSize(); ++h) {
            EDepSim::HitSegment* g4Hit
                = dynamic_cast<EDepSim::HitSegment*>(g4Hits->GetHit(h));
            if (!g4Hit) {
                EDepSimError("Not a hit segment");
                continue;
            }

            // Make sure that the primary trajectory associated with this hit
            // is saved.  The primary trajectories are defined by
            // EDepSim::TrajectoryMap::FindPrimaryId().
            int primaryId = g4Hit->GetPrimaryTrajectoryId();
            EDepSim::Trajectory* ndTraj
                = dynamic_cast<EDepSim::Trajectory*>(
                    EDepSim::TrajectoryMap::Get(primaryId));
            if (ndTraj) {
                ndTraj->MarkTrajectory(false);
            }
            else {
                EDepSimError("Primary trajectory not found");
            }
        }
    }
}

G4bool EDepSim::PersistencyManager::Open ( G4String  filename )

virtual

Open the output (ie database) file. This is used by the persistency messenger to open files using the G4 macro language. It can be an empty method.

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 60 of file EDepSimPersistencyManager.cc.

                                                        {
    EDepSimSevere(" -- Open is not implimented for " << filename);
    SetFilename(filename);
    return false;
}

virtual G4bool EDepSim::PersistencyManager::Retrieve ( G4Event *&  e )

inlinevirtual

Retrieve information from a file. These are not implemented.

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 59 of file EDepSimPersistencyManager.hh.

{e=NULL; return false;}

virtual G4bool EDepSim::PersistencyManager::Retrieve ( G4Run *&  r )

inlinevirtual

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 60 of file EDepSimPersistencyManager.hh.

{r=NULL; return false;}

virtual G4bool EDepSim::PersistencyManager::Retrieve ( G4VPhysicalVolume *&  w )

inlinevirtual

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 61 of file EDepSimPersistencyManager.hh.

{w=NULL; return false;}

bool EDepSim::PersistencyManager::SaveTrajectoryBoundary ( G4VTrajectory *  g4Traj, G4StepStatus  status, G4String  currentVolume, G4String  prevVolume  )

private

Return true if a trajectory point should be saved. The decision is based on the stepping status (must be fGeomBoundary), the current and the previous volume name (one must match a trajectory boundary regexp, the other must not).

Definition at line 108 of file EDepSimPersistencyManager.cc.

                                                                         {
    if (status != fGeomBoundary) return false;
    std::string particleInfo = ":" + g4Traj->GetParticleName();
    if (std::abs(g4Traj->GetCharge())<0.1) particleInfo += ":neutral";
    else particleInfo += ":charged";
    std::string current = particleInfo + ":" + currentVolume;
    std::string previous = particleInfo + ":" + prevVolume;
    for (std::vector<TPRegexp*>::iterator r = fTrajectoryBoundaries.begin();
         r != fTrajectoryBoundaries.end();
         ++r) {
        // Check if a watched volume is being entered.
        if ((*r)->Match(current)>0 && (*r)->Match(previous)<1) {
            EDepSimNamedDebug("boundary","Entering " << current);
            return true;
        }
        // Check if a watched volume is being exited.
        if ((*r)->Match(current)<1 && (*r)->Match(previous)>0) {
            EDepSimNamedDebug("boundary","Exiting " << current);
            return true;
        }
    }
    return false;
}

void EDepSim::PersistencyManager::SelectTrajectoryPoints ( std::vector< int > &  selected, G4VTrajectory *  g4Traj  )

private

Fill a vector with the indices of trajectory points that should be copied to the output file.

Definition at line 619 of file EDepSimPersistencyManager.cc.

                                                                           {

    selected.clear();
    if (g4Traj->GetPointEntries() < 1) {
        EDepSimError("Trajectory with no points"
                     << " " << g4Traj->GetTrackID()
                     << " " << g4Traj->GetParentID()
                     << " " << g4Traj->GetParticleName());
        return;
    }

    ////////////////////////////////////
    // Save the first point of the trajectory.
    ////////////////////////////////////
    selected.push_back(0);

    /////////////////////////////////////
    // Save the last point of the trajectory.
    /////////////////////////////////////
    int lastIndex = g4Traj->GetPointEntries()-1;
    if (lastIndex < 1) {
        EDepSimError("Trajectory with one point"
                     << " " << g4Traj->GetTrackID()
                     << " " << g4Traj->GetParentID()
                     << " " << g4Traj->GetParticleName());
        return;
    }
    selected.push_back(lastIndex);

    //////////////////////////////////////////////
    // Short out trajectories that don't create any energy deposit in a
    // sensitive detector.  These are trajectories that disappear from the
    // detector, so we don't need to record the extra trajectory points.  The
    // starting and stopping point of the particle are recorded in the
    // trajectory.
    //////////////////////////////////////////////
    EDepSim::Trajectory* ndTraj = dynamic_cast<EDepSim::Trajectory*>(g4Traj);
    if (ndTraj->GetSDTotalEnergyDeposit() < 1*eV) return;

    // Find the trajectory points where particles are entering and leaving the
    // detectors.
    EDepSim::TrajectoryPoint* edepPoint
        = dynamic_cast<EDepSim::TrajectoryPoint*>(g4Traj->GetPoint(0));
    G4String prevVolumeName = edepPoint->GetPhysVolName();
    for (int tp = 1; tp < lastIndex; ++tp) {
        edepPoint
            = dynamic_cast<EDepSim::TrajectoryPoint*>(g4Traj->GetPoint(tp));
        G4String volumeName = edepPoint->GetPhysVolName();
        // Save the point on a boundary crossing for volumes where we are
        // saving the entry and exit points.
        if (SaveTrajectoryBoundary(g4Traj,edepPoint->GetStepStatus(),
                                   volumeName,prevVolumeName)) {
            selected.push_back(tp);
        }
        prevVolumeName = volumeName;
    }

    // Save trajectory points where there is a "big" interaction.
    for (int tp = 1; tp < lastIndex; ++tp) {
        edepPoint
            = dynamic_cast<EDepSim::TrajectoryPoint*>(g4Traj->GetPoint(tp));
        // Just navigation....
        if (edepPoint->GetProcessType() == fTransportation) continue;
        // Not much energy deposit...
        if (edepPoint->GetProcessDeposit() < GetTrajectoryPointDeposit())
            continue;
        // Don't save optical photons...
        if (edepPoint->GetProcessType() == fOptical) continue;
        // Not a physics step...
        if (edepPoint->GetProcessType() == fGeneral) continue;
        if (edepPoint->GetProcessType() == fUserDefined) continue;
        // Don't save continuous ionization steps.
        if (edepPoint->GetProcessType() == fElectromagnetic
            && edepPoint->GetProcessSubType() == fIonisation) continue;
        // Don't save multiple scattering.
        if (edepPoint->GetProcessType() == fElectromagnetic
            && edepPoint->GetProcessSubType() == fMultipleScattering) continue;
        selected.push_back(tp);
    }

    // Make sure there aren't any duplicates in the selected trajectory points.
    std::sort(selected.begin(), selected.end());
    selected.erase(std::unique(selected.begin(), selected.end()),
                   selected.end());

    if (ndTraj->GetSDEnergyDeposit() < 1*eV) return;

    double desiredAccuracy = GetTrajectoryPointAccuracy();
    // Make sure that the trajectory accuracy stays in tolerance.
    for (int throttle = 0; throttle < 1000; ++throttle) {
        bool addPoint = false;
        for (std::vector<int>::iterator p1 = selected.begin();
             p1 != selected.end();
             ++p1) {
            std::vector<int>::iterator p2 = p1+1;
            if (p2==selected.end()) break;
            double trajectoryAccuracy = FindTrajectoryAccuracy(g4Traj,*p1,*p2);
            if (trajectoryAccuracy <= desiredAccuracy) continue;
            int split = SplitTrajectory(g4Traj,*p1,*p2);
            if (split < 0) continue;
            selected.push_back(split);
            addPoint = true;
            break;
        }
        std::sort(selected.begin(), selected.end());
        selected.erase(std::unique(selected.begin(), selected.end()),
                       selected.end());
        if (!addPoint) break;
    }
}

void EDepSim::PersistencyManager::SetDetectorPartition ( int  partition )

inline

Set the detector mask.

Definition at line 191 of file EDepSimPersistencyManager.hh.

{fDetectorPartition = partition;}

void EDepSim::PersistencyManager::SetFilename ( G4String  file )

inlineprotected

Set the output filename. This can be used by the derived classes to inform the base class of the output file name.

Definition at line 199 of file EDepSimPersistencyManager.hh.

                                    {
        fFilename = file;
    }

virtual void EDepSim::PersistencyManager::SetGammaThreshold ( G4double  thresh )

inlinevirtual

Set the momentum threshold required to save a gamma-ray as a trajectory.

Definition at line 112 of file EDepSimPersistencyManager.hh.

                                                    {
        fGammaThreshold = thresh;
    }

virtual void EDepSim::PersistencyManager::SetLengthThreshold ( G4double  thresh )

inlinevirtual

Set the threshold for length in a sensitive detector above which a trajectory will be saved. If a trajectory created this much track inside a sensitive detector, then it is saved.

Definition at line 103 of file EDepSimPersistencyManager.hh.

                                                     {
        fLengthThreshold = thresh;
    }

virtual void EDepSim::PersistencyManager::SetNeutronThreshold ( G4double  thresh )

inlinevirtual

Set the momentum threshold required to save a neutron as a trajectory.

Definition at line 123 of file EDepSimPersistencyManager.hh.

                                                      {
        fNeutronThreshold = thresh;
    }

virtual void EDepSim::PersistencyManager::SetSaveAllPrimaryTrajectories ( bool  val )

inlinevirtual

Set the flag to save primary particle trajectories. If this flag is true, then the trajectories for primary particles are saved even if they don't ultimately cause an energy deposition in a sensitive detector. If this flag is false, the trajectories for primary particles are only saved if they, or one of their children, deposit energy in a sensitive detector.

Definition at line 165 of file EDepSimPersistencyManager.hh.

                                                         {
        fSaveAllPrimaryTrajectories = val;
    }

virtual void EDepSim::PersistencyManager::SetTrajectoryPointAccuracy ( double  acc )

inlinevirtual

Get the required trajectory accuracy. Trajectory points are saved so that the interpolated position of a trajectory is never more than this distance from the unsaved points.

Definition at line 134 of file EDepSimPersistencyManager.hh.

                                                        {
        fTrajectoryPointAccuracy = acc;
    }

virtual void EDepSim::PersistencyManager::SetTrajectoryPointDeposit ( double  dep )

inlinevirtual

Get the minimum energy deposition for which a trajectory point will be saved. This is the trajectory point process deposition (GetProcessDeposit()).

Definition at line 148 of file EDepSimPersistencyManager.hh.

                                                       {
        fTrajectoryPointDeposit = dep;
    }

int EDepSim::PersistencyManager::SplitTrajectory ( G4VTrajectory *  traj, int  point1, int  point2  )

private

Save another trajectory point between point 1 and point to optimize the trajectory accuracy.

Definition at line 595 of file EDepSimPersistencyManager.cc.

                                                                    {

    int point3 = 0.5*(point1 + point2);
    if (point3 <= point1) EDepSimThrow("Points too close to split");
    if (point2 <= point3) EDepSimThrow("Points too close to split");
    double bestAccuracy = FindTrajectoryAccuracy(g4Traj, point1, point3);

    bestAccuracy = std::max(FindTrajectoryAccuracy(g4Traj, point3, point2),
                            bestAccuracy);
    for (int p = point1+1; p<point2-1; ++p) {
        double a1 = FindTrajectoryAccuracy(g4Traj,point1, p);
        double a2 = FindTrajectoryAccuracy(g4Traj,p, point2);
        double accuracy = std::max(a1,a2);
        if (accuracy < bestAccuracy) {
            point3 = p;
            bestAccuracy  = accuracy;
        }
    }

    return point3;
}

G4bool EDepSim::PersistencyManager::Store ( const G4Event *  anEvent )

virtual

stores anEvent and the associated objects into database.

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 73 of file EDepSimPersistencyManager.cc.

                                                              {
    UpdateSummaries(anEvent);
    return false;
}

G4bool EDepSim::PersistencyManager::Store ( const G4Run *  aRun )

virtual

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 79 of file EDepSimPersistencyManager.cc.

                                                         {
    if (!aRun) return false;
    EDepSimSevere(" -- Run store called without a save method for "
               << "GEANT4 run " << aRun->GetRunID());
    return false;
}

G4bool EDepSim::PersistencyManager::Store ( const G4VPhysicalVolume *  aWorld )

virtual

Reimplemented in EDepSim::RootPersistencyManager.

Definition at line 87 of file EDepSimPersistencyManager.cc.

                                                                       {
    if (!aWorld) return false;
    EDepSimSevere(" -- Geometry store called without a save method for "
               << aWorld->GetName());
    return false;
}

void EDepSim::PersistencyManager::SummarizeHitSegments ( std::vector< TG4HitSegment > &  segments, G4VHitsCollection *  hits  )

private

Fill a container of hit segments.

Definition at line 510 of file EDepSimPersistencyManager.cc.

                                                                        {
    dest.clear();

    EDepSim::HitSegment* g4Hit = dynamic_cast<EDepSim::HitSegment*>(g4Hits->GetHit(0));
    if (!g4Hit) return;

    for (std::size_t h=0; h<g4Hits->GetSize(); ++h) {
        g4Hit = dynamic_cast<EDepSim::HitSegment*>(g4Hits->GetHit(h));
        TG4HitSegment hit;
        hit.PrimaryId = fTrackIdMap[g4Hit->GetPrimaryTrajectoryId()];
        hit.EnergyDeposit = g4Hit->GetEnergyDeposit();
        hit.SecondaryDeposit = g4Hit->GetSecondaryDeposit();
        hit.TrackLength = g4Hit->GetTrackLength();
        CopyHitContributors(hit.Contrib,g4Hit->GetContributors());
        hit.Start.SetXYZT(g4Hit->GetStart().x(),
                          g4Hit->GetStart().y(),
                          g4Hit->GetStart().z(),
                          g4Hit->GetStart().t());
        hit.Stop.SetXYZT(g4Hit->GetStop().x(),
                          g4Hit->GetStop().y(),
                          g4Hit->GetStop().z(),
                          g4Hit->GetStop().t());
        dest.push_back(hit);
    }
}

void EDepSim::PersistencyManager::SummarizePrimaries ( std::vector< TG4PrimaryVertex > &  primaries, const G4PrimaryVertex *  event  )

private

Definition at line 158 of file EDepSimPersistencyManager.cc.

                                {
    dest.clear();

    if (!src) return;

    while (src) {
        TG4PrimaryVertex vtx;

        vtx.Position.SetX(src->GetX0());
        vtx.Position.SetY(src->GetY0());
        vtx.Position.SetZ(src->GetZ0());
        vtx.Position.SetT(src->GetT0());

        // Add the particles associated with the vertex to the summary.
        for (int i=0; i< src->GetNumberOfParticle(); ++i) {
            TG4PrimaryParticle prim;
            G4PrimaryParticle *g4Prim = src->GetPrimary(i);
            double E = 0.0;
            if (g4Prim->GetG4code()) {
                prim.Name = g4Prim->GetG4code()->GetParticleName();
                E = pow(g4Prim->GetG4code()->GetPDGMass(),2);
            }
            else {
                prim.Name = "unknown";
            }
            prim.PDGCode = g4Prim->GetPDGcode();
            prim.TrackId = g4Prim->GetTrackID() - 1;
            prim.Momentum.SetX(g4Prim->GetPx());
            prim.Momentum.SetY(g4Prim->GetPy());
            prim.Momentum.SetZ(g4Prim->GetPz());
            E += pow(prim.Momentum.P(),2);
            if (E>0) E = std::sqrt(E);
            else E = 0;
            prim.Momentum.SetE(E);
            vtx.Particles.push_back(prim);
        }

        // Check to see if there is anyu user information associated with the
        // vertex.
        EDepSim::VertexInfo* srcInfo
            = dynamic_cast<EDepSim::VertexInfo*>(src->GetUserInformation());
        if (srcInfo) {
            vtx.GeneratorName = srcInfo->GetName();
            vtx.Reaction = srcInfo->GetReaction();
            vtx.Filename = srcInfo->GetFilename();
            vtx.InteractionNumber = srcInfo->GetInteractionNumber();
            vtx.CrossSection = srcInfo->GetCrossSection();
            vtx.DiffCrossSection = srcInfo->GetDiffCrossSection();
            vtx.Weight = srcInfo->GetWeight();
            vtx.Probability = srcInfo->GetProbability();

            const G4PrimaryVertex* infoVertex
                = srcInfo->GetInformationalVertex();
            if (infoVertex) {
                SummarizePrimaries(vtx.Informational,
                                   srcInfo->GetInformationalVertex());
            }
        }

        dest.push_back(vtx);
        src = src->GetNext();
    }
}

void EDepSim::PersistencyManager::SummarizeSegmentDetectors ( TG4HitSegmentDetectors &  segmentDetectors, const G4Event *  event  )

private

sensitive detector.

Definition at line 486 of file EDepSimPersistencyManager.cc.

                          {
    dest.clear();

    G4HCofThisEvent* HCofEvent = event->GetHCofThisEvent();
    if (!HCofEvent) return;
    G4SDManager *sdM = G4SDManager::GetSDMpointer();
    G4HCtable *hcT = sdM->GetHCtable();
    // Copy each of the hit categories into the output event.
    for (int i=0; i<hcT->entries(); ++i) {
        G4String SDname = hcT->GetSDname(i);
        G4String HCname = hcT->GetHCname(i);
        int HCId = sdM->GetCollectionID(SDname+"/"+HCname);
        G4VHitsCollection* g4Hits = HCofEvent->GetHC(HCId);
        if (g4Hits->GetSize()<1) continue;
        EDepSim::HitSegment* hitSeg
            = dynamic_cast<EDepSim::HitSegment*>(g4Hits->GetHit(0));
        if (!hitSeg) continue;
        SummarizeHitSegments(dest[SDname],g4Hits);
    }
}

void EDepSim::PersistencyManager::SummarizeTrajectories ( std::vector< TG4Trajectory > &  trajectories, const G4Event *  event  )

private

Fill the trajectory container.

Rewrite the track ids so that they are consecutive.

Definition at line 224 of file EDepSimPersistencyManager.cc.

                          {
    dest.clear();
    MarkTrajectories(event);

    // Build a map of the original G4 TrackID to the new relocated TrackId
    // (not capitalization).  This also uses the fact that maps are sorted as
    // a hack to prevent writing a predicate to sort TG4Trajectories (because
    // I'm lazy).
    fTrackIdMap.clear();

    const G4TrajectoryContainer* trajectories = event->GetTrajectoryContainer();
    if (!trajectories) {
        EDepSimError("No trajectory container");
        return;
    }
    if (!trajectories->GetVector()) {
        EDepSimError("No trajectory vector in trajectory container");
        return;
    }

    int index = 0;
    TG4TrajectoryContainer tempContainer;
    for (TrajectoryVector::iterator t = trajectories->GetVector()->begin();
         t != trajectories->GetVector()->end();
         ++t) {
        EDepSim::Trajectory* ndTraj = dynamic_cast<EDepSim::Trajectory*>(*t);

        // Check if the trajectory should be saved.
        if (!ndTraj->SaveTrajectory()) continue;

        // Set the particle type information.
        G4ParticleDefinition* part
            = G4ParticleTable::GetParticleTable()->FindParticle(
                ndTraj->GetParticleName());
        if (!part) {
            EDepSimError(std::string("EDepSim::RootPersistencyManager::")
                      + "No particle information for "
                      + ndTraj->GetParticleName());
        }

        fTrackIdMap[ndTraj->GetTrackID()] = index++;

        TG4Trajectory traj;
        traj.TrackId = ndTraj->GetTrackID();
        traj.Name = ndTraj->GetParticleName();
        traj.PDGCode = ndTraj->GetPDGEncoding();
        traj.ParentId = ndTraj->GetParentID();
        traj.InitialMomentum.SetXYZM(ndTraj->GetInitialMomentum().x(),
                                     ndTraj->GetInitialMomentum().y(),
                                     ndTraj->GetInitialMomentum().z(),
                                     part->GetPDGMass());
        CopyTrajectoryPoints(traj,ndTraj);
        int throttle = 999999;
        do {
            if (traj.ParentId == 0) break;
            EDepSim::Trajectory* pTraj
                = dynamic_cast<EDepSim::Trajectory*>(
                    EDepSim::TrajectoryMap::Get(traj.ParentId));
            if (!pTraj) {
                EDepSimError("Trajectory " << traj.ParentId << " does not exist");
                throw;
                break;
            }
            if (pTraj->SaveTrajectory()) break;
            traj.ParentId = pTraj->GetParentID();
        } while (--throttle > 0);
        tempContainer.push_back(traj);
    }

    // Reorder the trajectories using the sorted fTrackIdMap.  This uses the
    // fact that the map is going to have the right sort order.
    for (TrackIdMap::iterator i = fTrackIdMap.begin();
         i != fTrackIdMap.end(); ++i) {
        dest.push_back(tempContainer[i->second]);
    }

    // Now that the trajectories are reordered, reassign the mapping from the
    // original ndTraj->GetTrackID() value to the index in the Trajectories
    // vector.
    index = 0;
    for (TrackIdMap::iterator i = fTrackIdMap.begin();
         i != fTrackIdMap.end(); ++i) {
        i->second = index++;
    }

    // Make double sure that TrackID zero maps to the non-existent -1.
    fTrackIdMap[0] = -1;

    /// Rewrite the track ids so that they are consecutive.
    for (TG4TrajectoryContainer::iterator
             t = dest.begin();
         t != dest.end(); ++t) {
        t->TrackId = fTrackIdMap[t->TrackId];
        t->ParentId = fTrackIdMap[t->ParentId];
    }

}

void EDepSim::PersistencyManager::UpdateSummaries ( const G4Event *  event )

protected

Update the event summary fields.

Definition at line 135 of file EDepSimPersistencyManager.cc.

                                                                    {

    const G4Run* runInfo = G4RunManager::GetRunManager()->GetCurrentRun();

    fEventSummary.RunId = runInfo->GetRunID();
    fEventSummary.EventId = event->GetEventID();
    EDepSimLog("Event Summary for run " << fEventSummary.RunId
               << " event " << fEventSummary.EventId);

    // Summarize the trajectories first so that fTrackIdMap is filled.
    MarkTrajectories(event);

    SummarizePrimaries(fEventSummary.Primaries,event->GetPrimaryVertex());
    EDepSimLog("   Primaries " << fEventSummary.Primaries.size());

    SummarizeTrajectories(fEventSummary.Trajectories,event);
    EDepSimLog("   Trajectories " << fEventSummary.Trajectories.size());

    SummarizeSegmentDetectors(fEventSummary.SegmentDetectors, event);
    EDepSimLog("   Segment Detectors "
               << fEventSummary.SegmentDetectors.size());
}

Member Data Documentation

int EDepSim::PersistencyManager::fDetectorPartition

private

The detector partition.

Definition at line 311 of file EDepSimPersistencyManager.hh.

TG4Event EDepSim::PersistencyManager::fEventSummary

protected

A summary of the primary vertices in the event.

Definition at line 207 of file EDepSimPersistencyManager.hh.

G4String EDepSim::PersistencyManager::fFilename

private

The filename of the output file.

Definition at line 278 of file EDepSimPersistencyManager.hh.

G4double EDepSim::PersistencyManager::fGammaThreshold

private

Threshold momentum above which gamma-ray trajectories are written. Gamma-rays are rejected if the momentum is less than either fGammaThreshold or fMomentumThreshold.

Definition at line 287 of file EDepSimPersistencyManager.hh.

G4double EDepSim::PersistencyManager::fLengthThreshold

private

If a trajectory left this much track in a sensitive detector then it should be saved.

Definition at line 282 of file EDepSimPersistencyManager.hh.

G4double EDepSim::PersistencyManager::fNeutronThreshold

private

The threshold momentum below which a neutron will be rejected.

Definition at line 290 of file EDepSimPersistencyManager.hh.

EDepSim::PersistencyMessenger* EDepSim::PersistencyManager::fPersistencyMessenger

private

Definition at line 314 of file EDepSimPersistencyManager.hh.

bool EDepSim::PersistencyManager::fSaveAllPrimaryTrajectories

private

Flag to determine if all primary trajectories are saved, or only those that ultimately create energy in a sensitive detector. The primary particles are always saved.

Definition at line 304 of file EDepSimPersistencyManager.hh.

TrackIdMap EDepSim::PersistencyManager::fTrackIdMap

private

Definition at line 275 of file EDepSimPersistencyManager.hh.

std::vector<TPRegexp*> EDepSim::PersistencyManager::fTrajectoryBoundaries

private

A list of detectors for which trajectory points should be saved as particles enter and exit.

Definition at line 308 of file EDepSimPersistencyManager.hh.

double EDepSim::PersistencyManager::fTrajectoryPointAccuracy

private

The maximum distance between the interpolated trajectory position and an unsaved trajectory point.

Definition at line 294 of file EDepSimPersistencyManager.hh.

double EDepSim::PersistencyManager::fTrajectoryPointDeposit

private

The minimum energy for the trajectory point process deposit (GetProcessDeposit()) to be saved. Points with less than this energy are not selected.

Definition at line 299 of file EDepSimPersistencyManager.hh.

---

The documentation for this class was generated from the following files:

• edep-sim/src/EDepSimPersistencyManager.hh
• edep-sim/src/EDepSimPersistencyManager.cc