Message: Re: Finding (grand^n)mother particle of a track [solved] Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

None Re: Finding (grand^n)mother particle of a track [solved] 

Forum: Event and Track Management
Re: Question Finding "mother particle" (Kyrre Ness Sjøbæk)
Re: Feedback Re: Finding (Gumplinger Peter)
Re: Feedback Re: Finding (Gumplinger Peter)
Date: 01 Aug, 2009
From: Kyrre Ness Sjøbæk <Kyrre Ness Sjøbæk>

Thank you!

This seems kindof similar to my current solution, which is using UserTrackingAction to store the data internaly in the UTA class. However you avoids having making the filter having to access the UTA-class.

The solution I eventually came up with was:

//Main filters:

#ifndef FilterMotherIon_hh
#define FilterMotherIon_hh

#include "G4VSDFilter.hh"
#include "globals.hh"

/*
 * Filter for a simple detector, in order to filter
 * hits depending on the originally decaying ion.
 *
 * True if hit comes from MotherIon or a child of MotherIon, else false.
 */

class EdepSpectraTrackingAction;

class FilterMotherIon : public G4VSDFilter {
public:
  FilterMotherIon(G4String name, G4int z, G4int a) :
    G4VSDFilter(name), Z(z), A(a) {};

  G4bool Accept(const G4Step* origStep) const;

private:
  G4int Z;
  G4int A;

};

#endif

#include "FilterMotherIon.hh"
#include "G4Track.hh"
#include "G4ParticleDefinition.hh"

#include "EdepSpectraTrackingAction.hh"
#include "G4RunManager.hh"

G4bool FilterMotherIon::Accept(const G4Step* origStep) const {
  //We need the trackingAction to do the digging
  EdepSpectraTrackingAction* trkAct = (EdepSpectraTrackingAction*) G4RunManager::GetRunManager()->GetUserTrackingAction();

  /*
   * Run down the particle tree until either condition is met:
   * - Correct particle is found (accept)
   * - Pre-primary particle is hit (not found - reject)
   */

  //Initialize on current particle
  G4int thisID   = origStep->GetTrack()->GetTrackID();
  G4int parentID = origStep->GetTrack()->GetParentID();
  G4ParticleDefinition* thisDef = origStep->GetTrack()->GetDefinition();
  //G4cout << "StartSearch, Z = " << Z << ", A = " << A << G4endl;
  //Search!
  while (true) {
    if (thisDef->GetAtomicNumber() == Z && thisDef->GetAtomicMass() == A) {
      //Gold!
      //G4cout << "Gold! ID = " << thisID << G4endl;
      return true;
    }
    else {
      //No luck. Try to dig deeper    
      if (parentID == 0) {
	//G4cout << "Rock" << G4endl;
	return false; //Rock bottom, no gold.
      }
      //Move the definition one step deeper
      thisID = parentID;
      parentID = trkAct->GetParentID(thisID);
      thisDef = trkAct->GetDefByID(thisID);
      //G4cout << "diggin'..." << G4endl;
    }
  }
}

//Some simple boolean filters (I have only tested and used the Not filter:

#ifndef FilterBools_hh
#define FilterBools_hh

#include "G4VSDFilter.hh"
#include "globals.hh"

/*
 * Set of boolean operations between two filters
 */

class FilterAnd : public G4VSDFilter {
public:
  FilterAnd(G4String name, G4VSDFilter* filter1, G4VSDFilter* filter2) :
    G4VSDFilter(name), Filter1(filter1), Filter2(filter2) {};

  inline G4bool Accept(const G4Step* origStep) const {
    return (Filter1->Accept(origStep) && Filter2->Accept(origStep));
  }

private:
  G4VSDFilter* Filter1;
  G4VSDFilter* Filter2;

};
class FilterOr : public G4VSDFilter {
public:
  FilterOr(G4String name, G4VSDFilter* filter1, G4VSDFilter* filter2) :
    G4VSDFilter(name), Filter1(filter1), Filter2(filter2) {};

  inline G4bool Accept(const G4Step* origStep) const {
    return (Filter1->Accept(origStep) || Filter2->Accept(origStep));
  }

private:
  G4VSDFilter* Filter1;
  G4VSDFilter* Filter2;
};

class FilterNot : public G4VSDFilter {
public:
  FilterNot(G4String name, G4VSDFilter* filter) :
    G4VSDFilter(name), Filter(filter) {};

  inline G4bool Accept(const G4Step* origStep) const {
    return (not Filter->Accept(origStep));
  }

private:
  G4VSDFilter* Filter;

};

#endif

//UserTrackingAction class: ifndef EdepSpectraTrackingAction_hh #define EdepSpectraTrackingAction_hh

#include "G4UserTrackingAction.hh"
#include "G4ParticleDefinition.hh"

#include <map>

class EdepSpectraTrackingAction : public G4UserTrackingAction{
public:
  EdepSpectraTrackingAction();

  void PreUserTrackingAction(const G4Track*);

  //Clear the maps for a new event
  void ResetMaps();

  inline G4int GetParentID(G4int child) {
    return parentMap[child];
  }
  inline G4ParticleDefinition* GetDefByID(G4int ID) {
    return defMap[ID];
  }

private:

  //What trackID <value> is parent to trackID <key> ?
  std::map<G4int, G4int> parentMap;
  // trackID <key>, particleDefinition<value>
  std::map<G4int, G4ParticleDefinition*> defMap;

};

#endif include "EdepSpectraTrackingAction.hh"

#include "G4Track.hh"

EdepSpectraTrackingAction::EdepSpectraTrackingAction() {

}

void EdepSpectraTrackingAction::PreUserTrackingAction(const G4Track* track) {
  //Method only called once pr. track - no need to check if already added
  G4int trackID = track->GetTrackID();

  parentMap[trackID] = track->GetParentID();
  defMap[trackID] = track->GetDefinition();

}

void EdepSpectraTrackingAction::ResetMaps() {
  //DEBUG: print the maps
  /*
  std::map<G4int,G4int>::iterator it;
  for ( it=parentMap.begin() ; it != parentMap.end(); it++ ) {
    G4cout << (*it).first << " => " << (*it).second << G4endl;
  }
  */

  parentMap.clear();
  defMap.clear();  
}

This is a somewhat simpler solution, and it does what I need it to do. Note that I need to calle EdepSpectraTrackingAction::ResetMaps() before each event -- i do this in UserPrimaryGeneratorAction::GeneratePrimaries() just before actually generating any primaries.

--- Kyrre

 Add Message Add Message
to: "Re: Finding (grand^n)mother particle of a track [solved]"

 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 ]