Message: Re: How can I Get the Energy the incident Particle on the Event Action? Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Feedback Re: How can I Get the Energy the incident Particle on the Event Action? 

Keywords: Primary Particle Energy in the Event Action
Forum: Event and Track Management
Re: None How can I Get the Energy the incident Particle on the Event Action? (Francisco Garcia)
Re: Question Re: How can I Get the Energy the incident Particle on the Event Action? (Sean Kirkwood)
Date: 17 Mar, 2008
From: Gumplinger Peter <Gumplinger Peter>

Hi Sean,

I am not sure what I am providing you with is 'simple'. ---

I have a hit class that is analogous to what we were doing in Geant3; e.g. we only made one hit per track segment through a sensitive detector volume. In this case, I keep the incident position, time, energy, exit positing, time, and particle PDG code:

#ifndef ScintHitSegment_h
#define ScintHitSegment_h 1

#include "globals.hh"
#include "G4VHit.hh"

#include "G4ThreeVector.hh"

#include "G4THitsCollection.hh"
#include "G4Allocator.hh"

#include "G4TouchableHandle.hh"

class G4Step;

class ScintHitSegment : public G4VHit {

public:

  ScintHitSegment();
  ~ScintHitSegment();

  ScintHitSegment(const ScintHitSegment &right);

  inline void *operator new(size_t);
  inline void  operator delete(void*);

  ScintHitSegment operator=(const ScintHitSegment&);
  ScintHitSegment operator+=(const ScintHitSegment&);

  /// Return the TrackID of the trajectory that made this hit.
  G4int GetTrackID(void) const {return fTrackID;}

  /// Return the ParentID of the trajectory that made this hit.
  G4int GetParentID(void) const {return fParentID;}

  /// Get the total energy deposited in this hit.
  G4double GetEnergyDeposit(void) const {return fEnergyDeposit;}

  /// Get the energy weighted hit position.
  G4ThreeVector GetEnergyWeightedPosition(void) const {return fPosition;}

  /// Get the total charged track length in this hit.
  G4double GetTrackLength(void) const {return fTrackLength;}

  /// The X position of the hit starting point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStartX(void) const {return fStartX;}

  /// The Y position of the hit starting point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStartY(void) const {return fStartY;}

  /// The Z position of the hit starting point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStartZ(void) const {return fStartZ;}

  /// The time of the hit starting point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStartT(void) const {return fStartT;}

  /// The Total energy at the hit starting point. Note that a hit by
  /// definition is in a single volume.
  G4double GetStartTotalEnergy(void) const {return fStartTotalEnergy;}

  /// The X position of the hit stopping point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStopX(void) const {return fStopX;}

  /// The Y position of the hit stopping point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStopY(void) const {return fStopY;}

  /// The Z position of the hit stopping point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStopZ(void) const {return fStopZ;}

  /// The time of the hit stopping point.  Note that a hit by
  /// definition is in a single volume.
  G4double GetStopT(void) const {return fStopT;}

  /// The PDG Encoding of the particle producing the hit.
  G4double GetPDG(void) const {return fPDG;}

  /// Add the effects of a step to this hit.
  void AddStep(G4Step* theStep);

  /// Hits for the same track in the same phys. volume belong to the same hit
  G4bool SameHit(G4Step* theStep);

  void Draw();
  void Print();

private:

  /// The G4 touchable handle to the volume that contains the hit.
  G4TouchableHandle fHandle;

  /// The TrackID of the trajectory that made this hit.
  G4int fTrackID;

  /// The ParentID of the trajectory that made this hit.
  G4int fParentID;

  /// The total energy deposit in this hit.
  G4double fEnergyDeposit;

  /// The energy weighted hit position.
  G4ThreeVector fPosition;

  /// The total charged track length in this hit.
  G4double fTrackLength;

  /// The X position of the hit starting point.
  G4double fStartX;

  /// The Y position of the hit starting point.
  G4double fStartY;

  /// The Z position of the hit starting point.
  G4double fStartZ;

  /// The time of the hit starting point.
  G4double fStartT;

  /// TotalEnergy of the hit starting point.
  G4double fStartTotalEnergy;

  /// The X position of the hit stopping point.
  G4double fStopX;

  /// The Y position of the hit stopping point.
  G4double fStopY;

  /// The Z position of the hit stopping point.
  G4double fStopZ;

  /// The time of the hit stopping point.
  G4double fStopT;

  /// The PDG Encoding of the particle
  G4double fPDG;

};

typedef G4THitsCollection<ScintHitSegment> ScintHitSegmentCollection;

extern G4Allocator<ScintHitSegment> ScintHitSegmentAllocator;

inline void* ScintHitSegment::operator new(size_t)
{
  void *aHit;
  aHit = (void *) ScintHitSegmentAllocator.MallocSingle();
  return aHit;
}

inline void ScintHitSegment::operator delete(void *aHit)
{
  ScintHitSegmentAllocator.FreeSingle((ScintHitSegment*) aHit);
}

#endif

#include "G4Step.hh"
#include "G4TouchableHandle.hh"
#include "G4VProcess.hh"

#include "G4UnitsTable.hh"
#include "G4VVisManager.hh"
#include "G4VisAttributes.hh"
#include "G4Polyline.hh"
#include "G4Color.hh"

#include "ScintHitSegment.hh"

G4Allocator<ScintHitSegment> ScintHitSegmentAllocator;

ScintHitSegment::ScintHitSegment()
  : fHandle(NULL), fTrackID(-1), fParentID(-1),
    fEnergyDeposit(0), fPosition(0,0,0), fTrackLength(0),
    fStartX(0), fStartY(0), fStartZ(0), fStartT(0), fStartTotalEnergy(0),
    fStopX(0),  fStopY(0),  fStopZ(0),  fStopT(0) {}

ScintHitSegment::ScintHitSegment(const ScintHitSegment &right)
    : G4VHit(),
    fHandle(right.fHandle),
    fTrackID(right.fTrackID),
    fParentID(right.fParentID),
    fEnergyDeposit(right.fEnergyDeposit),
    fPosition(right.fPosition),
    fTrackLength(right.fTrackLength),
    fStartX(right.fStartX), fStartY(right.fStartY),
    fStartZ(right.fStartZ), fStartT(right.fStartT),
    fStartTotalEnergy(right.fStartTotalEnergy),
    fStopX(right.fStopX), fStopY(right.fStopY),
    fStopZ(right.fStopZ), fStopT(right.fStopT) {}

ScintHitSegment::~ScintHitSegment() {}

ScintHitSegment ScintHitSegment::operator=(const ScintHitSegment& op2)
{
    fHandle = op2.fHandle;
    fTrackID = op2.fTrackID;
    fParentID = op2.fParentID;
    fEnergyDeposit = op2.fEnergyDeposit;
    fPosition = op2.fPosition;
    fTrackLength = op2.fTrackLength;
    fStartX = op2.fStartX;
    fStartY = op2.fStartY;
    fStartZ = op2.fStartZ;
    fStartT = op2.fStartT;
    fStartTotalEnergy = op2.fStartTotalEnergy;
    fStopX  = op2.fStopX;
    fStopY  = op2.fStopY;
    fStopZ  = op2.fStopZ;
    fStopT  = op2.fStopT;
    fPDG    = op2.fPDG;
    return *this;
}

ScintHitSegment ScintHitSegment::operator+=(const ScintHitSegment& op2)
{
   fPosition *= fEnergyDeposit;
   fPosition += op2.GetEnergyDeposit()*op2.GetEnergyWeightedPosition();

   fEnergyDeposit = fEnergyDeposit + op2.GetEnergyDeposit();
   fPosition /= fEnergyDeposit;

   fTrackLength += op2.GetTrackLength();

   if(op2.GetStopT() >fStopT )fStopT  = op2.GetStopT();
   if(op2.GetStartT()<fStartT)fStartT = op2.GetStartT();

   fStartTotalEnergy += op2.GetStartTotalEnergy();

   if(fParentID == op2.GetTrackID()){
     fParentID = op2.GetParentID();
     fTrackID  = op2.GetTrackID();
     fPDG      = op2.GetPDG();
     fStartX   = op2.GetStartX();
     fStartY   = op2.GetStartY();
     fStartZ   = op2.GetStartZ();
     fStopX    = op2.GetStopX();
     fStopY    = op2.GetStopY();
     fStopZ    = op2.GetStopZ();
   }

   return *this;

}

G4bool ScintHitSegment::SameHit(G4Step* theStep)
{
  // Check that the track IDs are the same

  G4int trackId = theStep->GetTrack()->GetTrackID();
  if (fTrackID != trackId) return false;

  // Check that the hit and new step are in the same volume

  const G4TouchableHandle& touchable
      = theStep->GetPreStepPoint()->GetTouchableHandle();

  G4int historyDepth1 = fHandle->GetHistoryDepth();
  G4int historyDepth2 = touchable ->GetHistoryDepth();

  if (historyDepth1 != historyDepth2) return false;

  for (G4int i = 0; i<historyDepth1; ++i) {
    if (fHandle->GetVolume(i) != touchable->GetVolume(i)) return false;
    if (fHandle->GetReplicaNumber(i) !=
        touchable->GetReplicaNumber(i)) return false;
  }

  return true;

}

void ScintHitSegment::AddStep(G4Step* theStep)
{
  G4TouchableHandle touchable
      = theStep->GetPreStepPoint()->GetTouchableHandle();

  G4ThreeVector prePos  = theStep->GetPreStepPoint()->GetPosition();
  G4ThreeVector postPos = theStep->GetPostStepPoint()->GetPosition();

  G4StepStatus stepStatus = theStep->GetPostStepPoint()->GetStepStatus();
  G4ParticleDefinition* particle =  theStep->GetTrack()->GetDefinition();

  G4double energyDeposit = theStep->GetTotalEnergyDeposit();
  G4double stepLength = theStep->GetStepLength();

  // Occasionally, a neutral particle will produce a particle below
  // threshold, and it will be recorded as generating the hit.  All of the
  // energy should be deposited at the stopping point of the track.

  if (stepStatus == fPostStepDoItProc
        && std::abs(particle->GetPDGCharge()) == 0.0) {
     prePos = postPos;
     stepLength = 0.0*cm;
  }

  /// First make sure the hit has been initialized.

  if (!fHandle) {

     fHandle = touchable;
     fStartX = prePos.x();
     fStartY = prePos.y();
     fStartZ = prePos.z();
     fStartT = theStep->GetPreStepPoint()->GetGlobalTime();
     fStartTotalEnergy = theStep->GetPreStepPoint()->GetTotalEnergy();
     fStopX = postPos.x();
     fStopY = postPos.y();
     fStopZ = postPos.z();
     fStopT = theStep->GetPostStepPoint()->GetGlobalTime();
     fTrackID = theStep->GetTrack()->GetTrackID();
     fParentID = theStep->GetTrack()->GetParentID();

  } else {

     // Record the track that contributes to this hit.
     G4int trackId = theStep->GetTrack()->GetTrackID();
     G4int parentId = theStep->GetTrack()->GetParentID();

     if (trackId == fTrackID) {
        fStopX = postPos.x();
        fStopY = postPos.y();
        fStopZ = postPos.z();
        fStopT = theStep->GetPostStepPoint()->GetGlobalTime();
     } else {
        fTrackID = trackId;
        fParentID = parentId;
     }
  }

  fPosition *= fEnergyDeposit;

  fEnergyDeposit += energyDeposit;
  fTrackLength   += stepLength;

  fPosition += energyDeposit*((postPos+prePos)/2.);
  fPosition /= fEnergyDeposit;

}

void ScintHitSegment::Draw(void)
{
    G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
    if(pVVisManager) {
        G4Polyline line;
        line.push_back(G4Point3D(fStartX, fStartY, fStartZ));
        line.push_back(G4Point3D(fStopX,  fStopY,  fStopZ ));
        G4Color color(1.,0.,1.);
        G4VisAttributes attribs(color);
        line.SetVisAttributes(attribs);
        pVVisManager->Draw(line);
    }
}

void ScintHitSegment::Print(void)
{
    G4cout << "Hit Energy Deposit: "
           << G4BestUnit(fEnergyDeposit,"Energy")
           << G4endl;
}

This hit class goes with this Sensitive Detector:

#ifndef ScintSD_h
#define ScintSD_h 1

#include "ScintHitSegment.hh"

#include "G4VSensitiveDetector.hh"

class G4Step; class G4HCofThisEvent; class G4TouchableHistory;

class ScintSD : public G4VSensitiveDetector {

public:

  ScintSD(G4String name);
  virtual ~ScintSD();

  void Initialize(G4HCofThisEvent*);
  G4bool ProcessHits(G4Step*, G4TouchableHistory*);
  void EndOfEvent(G4HCofThisEvent*);

private:

  /// The collection of hits that is being filled in the current event.  It
  /// is constructed in Initialize, filled in ProcessHits, and added to the
  /// event in EndOfEvent.

  ScintHitSegmentCollection* fHits;

  /// The hit collection id of fHits
  G4int fHCID;

  /// The last hit that was found.
  G4int fLastHit;

};

#endif

#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4VProcess.hh"
#include "G4TouchableHistory.hh"
#include "G4SDManager.hh"

#include "ScintSD.hh"
#include "ScintHitSegment.hh"

ScintSD::ScintSD(G4String name)
  : G4VSensitiveDetector(name), fHits(NULL), fHCID(-1), fLastHit(-1)
{
  collectionName.insert("scintCollection");
}

ScintSD::~ScintSD() {}

void ScintSD::Initialize(G4HCofThisEvent* HCE)
{
  fHits =
     new ScintHitSegmentCollection(GetName(), GetCollectionName(0));

  if (fHCID<0) {
      G4String hcName = GetName() + "/" + GetCollectionName(0);
      fHCID = G4SDManager::GetSDMpointer()->GetCollectionID(hcName);
  }

  HCE->AddHitsCollection(fHCID, fHits);

}

G4bool ScintSD::ProcessHits(G4Step* theStep, G4TouchableHistory*) {

  G4double energyDeposit = theStep->GetTotalEnergyDeposit();

  if (energyDeposit <= 0.) return true;

  ScintHitSegment* currentHit = NULL;

  G4StepPoint* preStepPoint = theStep->GetPreStepPoint();
  const G4VProcess* preProcess = preStepPoint->GetProcessDefinedStep();

  G4ProcessType preProcessType = fNotDefined;
  if (preProcess) preProcessType = preProcess->GetProcessType();

  // If the presStepPoint is of type fTransporation, then we need a new
  // hit since we are crossing a geometry boundary. On the other hand,
  // we look at the most recent hit stored in the collection to see if
  // it is the 'same hit' (same track ID and volume tree).

  if (preProcessType != fTransportation) {
     if (fHits) {
        fLastHit = fHits->entries()-1;
        if (fLastHit>-1) {
           ScintHitSegment *tmpHit = (*fHits)[fLastHit];
           if (tmpHit->SameHit(theStep)){
              currentHit = tmpHit;
           }
        }
     }
  }

  if (!currentHit) {
     currentHit = new ScintHitSegment();
     fHits->insert(currentHit);
     fLastHit = fHits->entries()-1;
  }

  currentHit->AddStep(theStep);

  return true;

}

void ScintSD::EndOfEvent(G4HCofThisEvent*){}

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

 Add Message Add Message
to: "Re: How can I Get the Energy the incident Particle on the Event Action?"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]