AFP_TDSensitiveDetector Class Reference

#include <AFP_TDSensitiveDetector.h>

Inheritance diagram for AFP_TDSensitiveDetector:

Collaboration diagram for AFP_TDSensitiveDetector:

Public Member Functions
	AFP_TDSensitiveDetector (const std::string &name, const std::string &hitCollectionName)
	~AFP_TDSensitiveDetector ()
void	Initialize (G4HCofThisEvent *) override final
G4bool	ProcessHits (G4Step *, G4TouchableHistory *) override final
template<class... Args>
void	AddHit (Args &&... args)
	Templated method to stuff a single hit into the sensitive detector class.

Static Public Attributes
static constexpr double	TDMaxQEff = 0.15

Private Member Functions
	FRIEND_TEST (AFP_TDSensitiveDetectortest, ProcessHits)
	FRIEND_TEST (AFP_TDSensitiveDetectortest, AddHit)
AFP_TDSimHitCollectionBuilder *	getHitCollection () const

Private Attributes
int	m_nHitID
std::string	m_hitCollectionName
AFP_TDSimHitCollectionBuilder *	m_HitColl {}

Detailed Description

Definition at line 22 of file AFP_TDSensitiveDetector.h.

Constructor & Destructor Documentation

AFP_TDSensitiveDetector()

AFP_TDSensitiveDetector::AFP_TDSensitiveDetector	(	const std::string &	name,
		const std::string &	hitCollectionName )

Definition at line 32 of file AFP_TDSensitiveDetector.cxx.

  : G4VSensitiveDetector( name )
  , m_nHitID(-1)
  , m_hitCollectionName(hitCollectionName)
{
}

~AFP_TDSensitiveDetector()

AFP_TDSensitiveDetector::~AFP_TDSensitiveDetector ( )

inline

Definition at line 32 of file AFP_TDSensitiveDetector.h.

{ /* I don't own myHitColl if all has gone well */ }

Member Function Documentation

AddHit()

template<class... Args>

void AFP_TDSensitiveDetector::AddHit ( Args &&... args )

inline

Templated method to stuff a single hit into the sensitive detector class.

This could get rather tricky, but the idea is to allow fast simulations to use the very same SD classes as the standard simulation.

Definition at line 41 of file AFP_TDSensitiveDetector.h.

  {
    if (m_HitColl) {
      m_HitColl->Emplace(std::forward<Args>(args)...);
    }
  }

FRIEND_TEST() [1/2]

AFP_TDSensitiveDetector::FRIEND_TEST ( AFP_TDSensitiveDetectortest , AddHit  )

private

FRIEND_TEST() [2/2]

AFP_TDSensitiveDetector::FRIEND_TEST ( AFP_TDSensitiveDetectortest , ProcessHits  )

private

getHitCollection()

AFP_TDSimHitCollectionBuilder * AFP_TDSensitiveDetector::getHitCollection ( ) const

private

Definition at line 413 of file AFP_TDSensitiveDetector.cxx.

{
  auto* eventInfo = AtlasG4EventUserInfo::GetEventUserInfo();
  if (!eventInfo) {
    return nullptr;
  }
  auto hitCollections = eventInfo->GetHitCollectionMap();
  return hitCollections ? hitCollections->Find<AFP_TDSimHitCollectionBuilder>(m_hitCollectionName) : nullptr;
}

Initialize()

void AFP_TDSensitiveDetector::Initialize ( G4HCofThisEvent * )

finaloverride

Definition at line 41 of file AFP_TDSensitiveDetector.cxx.

{
  m_HitColl = getHitCollection();
}

ProcessHits()

bool AFP_TDSensitiveDetector::ProcessHits ( G4Step * pStep, G4TouchableHistory *  )

finaloverride

Definition at line 48 of file AFP_TDSensitiveDetector.cxx.

{
  if (!m_HitColl) {
    m_HitColl = getHitCollection();
    if (!m_HitColl) {
      return false;
    }
  }
 
  if(verboseLevel>5)
    {
      G4cout << "AFP_TDSensitiveDetector::ProcessHits" << G4endl;
    }
 
  bool bRes=false;
 
  int nTrackID=-1;
  int nParticleEncoding=-1;
  float fKineticEnergy=0.0;
  float fEnergyDeposit=0.0;
  float fWaveLength=0.0;
  float fPreStepX=0.0;
  float fPreStepY=0.0;
  float fPreStepZ=0.0;
  float fPostStepX=0.0;
  float fPostStepY=0.0;
  float fPostStepZ=0.0;
  float fGlobalTime=0.0;
  int nStationID=-1;
  int nDetectorID=-1;
  int nQuarticID=-1;
  // int nPixelRow=-1;
  // int nPixelCol=-1;
 
  // step, track and particle info
  G4Track* pTrack = pStep->GetTrack();
  G4ParticleDefinition* pParticleDefinition = pTrack->GetDefinition();
  G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
  G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
  G4ThreeVector PreStepPointPos = pPreStepPoint->GetPosition();
  G4ThreeVector PostStepPointPos = pPostStepPoint->GetPosition();
 
  nTrackID=pTrack->GetTrackID();
  fKineticEnergy = pPreStepPoint->GetKineticEnergy();
  fEnergyDeposit = pStep->GetTotalEnergyDeposit();
 
  fPreStepX = PreStepPointPos.x();
  fPreStepY = PreStepPointPos.y();
  fPreStepZ = PreStepPointPos.z();
  fPostStepX = PostStepPointPos.x();
  fPostStepY = PostStepPointPos.y();
  fPostStepZ = PostStepPointPos.z();
  nParticleEncoding = pParticleDefinition->GetPDGEncoding();
  fGlobalTime = pStep->GetPreStepPoint()->GetGlobalTime()/CLHEP::picosecond; // time w.r.t. prestep or poststep ??
 
  // name of physical volume
  G4TouchableHandle touch1 = pPreStepPoint->GetTouchableHandle();
  G4VPhysicalVolume* volume = touch1->GetVolume();
  G4String VolumeName = volume->GetName();
 
  //G4ThreeVector ph0 = pStep->GetDeltaPosition().unit();
  //if (fKineticEnergy<10000. && ph0.y()>.2) pTrack->SetTrackStatus(fKillTrackAndSecondaries);
  //if (VolumeName.contains("TDQuarticBar")) return 1;
 
  if(verboseLevel>5)
    {
      G4cout << "hit volume name is " << VolumeName << G4endl;
 
      G4cout << "global, x_pre:  " << fPreStepX  << ", y_pre:  " << fPreStepY  << ", z_pre:  " << fPreStepZ << G4endl;
      G4cout << "global, x_post: " << fPostStepX << ", y_post: " << fPostStepY << ", z_post: " << fPostStepZ << G4endl;
    }
  //scan station and detector id
  char* ppv1, *ppv2;
  char szbuff[32];
  memset(&szbuff[0],0,sizeof(szbuff));
  strncpy(szbuff,VolumeName.data(),sizeof(szbuff));
  szbuff[sizeof(szbuff)-1] = '\0'; // idiomatic use of strncpy...
  ppv1=strchr(szbuff,'[');
  ppv2=strchr(szbuff,']');
  if(!ppv2 || !ppv1){
    G4cout << "ERROR: Invalid format of volume name " << VolumeName << G4endl;
    return false;
  }
  else *ppv2='\0';
 
  nStationID=10*(szbuff[3]-0x30)+(szbuff[4]-0x30);
  nDetectorID=atoi(ppv1+1);
 
  m_nHitID++;
 
  /*
    if (VolumeName.contains("TDSensor") || (bRes=VolumeName.contains("TDQuarticBar["))){
    nQuarticID=szbuff[7]-0x30;
 
    if(VolumeName.contains("TDSensor") && pStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()!="OpAbsorption" )
    {
    //hit in TD sensor but with no OpAbsorption (transportation)
    }
    else{
 
    fWaveLength = 2.*M_PI*CLHEP::hbarc/(CLHEP::MeV*CLHEP::nm)/fKineticEnergy;
 
    if (fWaveLength > 800. || fWaveLength < 200.) return 1; // 200-800 nm cut
    AFP_TDSimHit* pHit = new AFP_TDSimHit();
    pHit->m_nHitID=m_nHitID;
    pHit->m_nTrackID=nTrackID;
    pHit->m_nParticleEncoding=nParticleEncoding;
    pHit->m_fKineticEnergy=fKineticEnergy;
    pHit->m_fEnergyDeposit=fEnergyDeposit;
    pHit->m_fWaveLength=fWaveLength;
    pHit->m_fPreStepX=fPreStepX;
    pHit->m_fPreStepY=fPreStepY;
    pHit->m_fPreStepZ=fPreStepZ;
    pHit->m_fPostStepX=fPostStepX;
    pHit->m_fPostStepY=fPostStepY;
    pHit->m_fPostStepZ=fPostStepZ;
    pHit->m_fGlobalTime=fGlobalTime;
 
    pHit->m_nStationID=nStationID;
    pHit->m_nDetectorID=nDetectorID;
    pHit->m_nSensitiveElementID=(bRes? 2:1)+2*nQuarticID;//Q1: 1-2, Q2: 3-4
    m_HitColl->Insert(*pHit);
    }
    }
  */
 
#if G4VERSION_NUMBER < 1100
  if ( (VolumeName.contains("TDQuarticBarVacBorder")) && pParticleDefinition->GetPDGCharge() !=0 )
#else
  if ( (G4StrUtil::contains(VolumeName, "TDQuarticBarVacBorder")) && pParticleDefinition->GetPDGCharge() !=0 )
#endif
    {
      nQuarticID=szbuff[7]-0x30;
      /*
        m_HitColl->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
        fWaveLength,fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,
        fPostStepZ,fGlobalTime,nStationID,nDetectorID,(2+2*nQuarticID));//Q1: 1-2, Q2: 3-4
        // m_HitColl->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
        //                                fWaveLength,fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,
        //                                fPostStepZ,fGlobalTime,nStationID,nDetectorID,((bRes? 2:1)+2*nQuarticID));//Q1: 1-2, Q2: 3-4
      */
    }

#if G4VERSION_NUMBER < 1100
  if ( (bRes=VolumeName.contains("TDQuarticBar[")) )
#else
  if ( (bRes=G4StrUtil::contains(VolumeName, "TDQuarticBar[")) )
#endif
    {
      nQuarticID=szbuff[7]-0x30;
 
      // Cut on maximum number of generated photons/bar
      if (m_HitColl->HasReachedLimit(nStationID, nQuarticID, nDetectorID)) return 1;
 
      // Get the Touchable History:
      const G4TouchableHistory* myTouch = static_cast<const G4TouchableHistory*>(pPreStepPoint->GetTouchable());
      // Calculate the local step position.
      // From a G4 FAQ:
      // http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/geometry/17/1.html
      const G4AffineTransform transformation = myTouch->GetHistory()->GetTopTransform();
      G4ThreeVector PreStepPointPos2  = transformation.TransformPoint(PreStepPointPos);
      G4ThreeVector PostStepPointPos2 = transformation.TransformPoint(PostStepPointPos);
 
      G4String shape( myTouch->GetSolid()->GetEntityType() );
 
      G4ThreeVector normpX( 1., 0., 0.);
      G4ThreeVector normnX(-1., 0., 0.);
      G4ThreeVector normpY( 0., 1., 0.);
      G4ThreeVector normnY( 0.,-1., 0.);
      G4ThreeVector normpZ( 0., 0., 1.);
      G4ThreeVector normnZ( 0., 0.,-1.);
 
      //G4double BarpX = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normpX);
      //G4double BarnX = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normnX);
      //G4double BarpY = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normpY);
      //G4double BarnY = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normnY);
      //G4double BarpZ = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normpZ);
      //G4double BarnZ = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PreStepPointPos2, normnZ);
 
      //G4double BarHalfX = .5 * (BarpX+BarnX);
      //G4double BarHalfY = .5 * (BarpY+BarnY);
      //G4double BarHalfZ = .5 * (BarpZ+BarnZ);
 
      G4double PreProtonX =  PreStepPointPos2.x();
      G4double PreProtonY =  PreStepPointPos2.y();
      G4double PreProtonZ =  PreStepPointPos2.z();
 
      G4double PostProtonX =  PostStepPointPos2.x();
      G4double PostProtonY =  PostStepPointPos2.y();
      G4double PostProtonZ =  PostStepPointPos2.z();
 
      G4ThreeVector p0 = pStep->GetDeltaPosition().unit();
 
      G4Material* mat = pStep->GetTrack()->GetMaterial();
      G4MaterialPropertiesTable *matPropTable = mat->GetMaterialPropertiesTable();
      // Refractive index
      G4MaterialPropertyVector* Rind = matPropTable->GetProperty("RINDEX");
      // Absorbtion length
      G4MaterialPropertyVector* Alen = matPropTable->GetProperty("ABSLENGTH");
 
      const G4double charge  = pParticleDefinition->GetPDGCharge();
      const G4double beta    = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
      G4double BetaInverse   = 1. / beta;
 
      //G4int Rsize   = Rind->Entries() - 1;
      //G4double Pmin = Rind->GetMinPhotonEnergy();     // 800 nm
#if G4VERSION_NUMBER < 1100
      G4double Pmin = Rind->GetMinLowEdgeEnergy();
      //G4double Pmax = Rind->GetMaxPhotonEnergy();  // 200 nm
      G4double Pmax = Rind->GetMaxLowEdgeEnergy();
#else
      G4double Pmin = Rind->GetMinEnergy();
      G4double Pmax = Rind->GetMaxEnergy();
#endif
      G4double dp   = Pmax - Pmin;
      //G4double maxCosTheta  = BetaInverse / Rind->GetMinProperty();
      G4double maxCosTheta  = BetaInverse / Rind->GetMinValue();
      G4double maxSin2Theta = (1.0 - maxCosTheta) * (1.0 + maxCosTheta);
 
      //G4double meanRI = .5*(Rind->GetMinProperty() + Rind->GetMaxProperty());
      G4double meanRI = .5*(Rind->GetMinValue() + Rind->GetMaxValue());
 
      // Formula taken from G4 docu (to be changed since the integral approximation)
      G4double MeanNumberOfPhotons = 370.*(charge/CLHEP::eplus)*(charge/CLHEP::eplus)* (1.0 - 1.0/(beta * meanRI * beta * meanRI)) / (CLHEP::cm*CLHEP::eV);
      if (MeanNumberOfPhotons <= 0.0) return 1;
 
      G4double step_length = pStep->GetStepLength();
 
      MeanNumberOfPhotons = MeanNumberOfPhotons * step_length * dp;
      G4int NumPhotons = (G4int) G4Poisson( MeanNumberOfPhotons );
      if (NumPhotons <= 0) return 1;
      //ATH_MSG_INFO("number of photons: " << NumPhotons);
 
      G4int NumPhotonsCuts=0;
      for (G4int I = 0; I < NumPhotons; I++) {
 
        G4double rand;
        G4double sampledEnergy, sampledRI;
        G4double cosTheta, sin2Theta;
 
        // Sample an energy for photon (using MC elimination method)
        do {
          rand = G4UniformRand();
          sampledEnergy = Pmin + rand * dp;
          //sampledRI = Rind->GetProperty(sampledEnergy);
          sampledRI = Rind->Value(sampledEnergy);
          cosTheta = BetaInverse / sampledRI;
 
          sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
          rand = G4UniformRand();
 
        } while (rand * maxSin2Theta > sin2Theta);
 
        // Generate random position of photon on the cone surface defined by Theta
        rand = G4UniformRand();
        G4double phi = 2.*M_PI*rand;
        G4double sinPhi = sin(phi);
        G4double cosPhi = cos(phi);
 
        // Calculate x,y,z coordinates of photon momentum
        // in coordinate system with primary particle direction aligned with the z-axis
        // + Rotate momentum direction back to the global coordinate system
        G4double sinTheta = sqrt(sin2Theta);
        G4double px = sinTheta*cosPhi;
        G4double py = sinTheta*sinPhi;
        G4double pz = cosTheta;
        G4ParticleMomentum photonMomentum(px, py, pz);
        photonMomentum.rotateUz(p0);
 
        G4double PX = photonMomentum.getX();
        G4double PY = photonMomentum.getY();
        G4double PZ = photonMomentum.getZ();
 
        // calculate projections coordinates
        //G4double PXp = PX/sqrt(PX*PX+PY*PY+PZ*PZ);
        G4double PYp = PY/sqrt(PX*PX+PY*PY+PZ*PZ);
        G4double PZp = PZ/sqrt(PX*PX+PY*PY+PZ*PZ);
 
        G4double PYt = PY/sqrt(PY*PY+PZ*PZ);
        G4double PZt = PZ/sqrt(PY*PY+PZ*PZ);
 
        // Cosines (alpha, delta)
        cosPhi = (PYp*PYt + PZp*PZt);
        if (nStationID == 0)    cosTheta = ( -PYt*sin(48.*CLHEP::deg) + PZt*cos(48.*CLHEP::deg) );
        else                    cosTheta = ( -PYt*sin(48.*CLHEP::deg) - PZt*cos(48.*CLHEP::deg) );
 
        // Total internal reflection conditions
        G4double cosThetaC = sqrt(1.-1./sampledRI/sampledRI);
        if (sqrt(1.-cosPhi*cosPhi)>cosThetaC) continue;
        if (sqrt(1.-cosTheta*cosTheta)*cosPhi>cosThetaC) continue;
 
        // Parametric equation of line where photons are generated
        rand = G4UniformRand();
        G4double PhotonX = PreProtonX + (PostProtonX-PreProtonX)*rand;
        G4double PhotonY = PreProtonY + (PostProtonY-PreProtonY)*rand;
        G4double PhotonZ = PreProtonZ + (PostProtonZ-PreProtonZ)*rand;
        G4ThreeVector PhotonPos(PhotonX,PhotonY,PhotonZ);
        G4double PhotonR = sqrt( (PreProtonX-PhotonX)*(PreProtonX-PhotonX) + (PreProtonY-PhotonY)*(PreProtonY-PhotonY) + (PreProtonZ-PhotonZ)*(PreProtonZ-PhotonZ) );
 
        G4double Y0;
        // including the scattering from the top edge of the bar (perfect absorber now)
        if (cosTheta>=0)
          {
            Y0 = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(PhotonPos, normnY);
            Y0 = Y0/cosTheta/cosPhi;
          }
        else
          {
            continue;
            //Y0 = 2.*((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PhotonPos, normpY) + ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(PhotonPos, normnY);
            //Y0 = -Y0/cosTheta/cosPhi;
          }
 
        // absorption of photons inside the crystal
        //float Pabs = 1. - exp( - Y0/Alen->GetProperty(sampledEnergy) );
        float Pabs = 1. - exp( - Y0/Alen->Value(sampledEnergy) );
        rand = G4UniformRand();
        if (Pabs>rand) continue;
 
        // maximum PMT efficiency cut (15%) to to avoid crashes due to the too large memory consumption
        rand = G4UniformRand();
        if (rand>TDMaxQEff) continue;
 
        NumPhotonsCuts++;
 
        float fGlobalTime2 = fGlobalTime;
        fGlobalTime2 += ( PhotonR * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
 
        // for group velocity of light: Edn/dE
        float EdndE;
        //if (sampledEnergy > (Pmin+.5*dp)) EdndE = (sampledRI - Rind->GetProperty(sampledEnergy-0.0001*CLHEP::eV))/0.0001*sampledEnergy/CLHEP::eV;
        if (sampledEnergy > (Pmin+.5*dp)) EdndE = (sampledRI - Rind->Value(sampledEnergy-0.0001*CLHEP::eV))/0.0001*sampledEnergy/CLHEP::eV;
        //else                            EdndE = (Rind->GetProperty(sampledEnergy+0.0001*CLHEP::eV) - sampledRI)/0.0001*sampledEnergy/CLHEP::eV;
        else                              EdndE = (Rind->Value(sampledEnergy+0.0001*CLHEP::eV) - sampledRI)/0.0001*sampledEnergy/CLHEP::eV;
        fGlobalTime2 += ( (sampledRI + EdndE)* Y0 * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
 
        if (verboseLevel>5)
          {
            G4cout << "FastCher EdndE: " << EdndE << G4endl;
          }
        fWaveLength = 2.*M_PI*CLHEP::hbarc/sampledEnergy/(CLHEP::MeV*CLHEP::nm);
 
        // Cut on maximum number of generated photons/bar
        if (m_HitColl->HasReachedLimit(nStationID, nQuarticID, nDetectorID)) return 1;
 
        int nSensitiveElementID=-1;
        if(nQuarticID==0) { nSensitiveElementID=1; }
        else if(nQuarticID==1) { nSensitiveElementID=3; }
 
        m_HitColl->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
                           fWaveLength,PhotonX,PhotonY,PhotonZ,(PhotonX+PX),(PhotonY+PY),(PhotonZ+PZ),
                           fGlobalTime2,nStationID,nDetectorID,nSensitiveElementID);
        m_HitColl->CountHit(nStationID, nQuarticID, nDetectorID);
      }
      if(verboseLevel>5)
        {
          G4cout << "FastCher number of photons: " << NumPhotonsCuts << G4endl;
        }
    }
  return true;
}

Member Data Documentation

m_HitColl

AFP_TDSimHitCollectionBuilder* AFP_TDSensitiveDetector::m_HitColl {}

private

Definition at line 57 of file AFP_TDSensitiveDetector.h.

{};

m_hitCollectionName

std::string AFP_TDSensitiveDetector::m_hitCollectionName

private

Definition at line 55 of file AFP_TDSensitiveDetector.h.

m_nHitID

int AFP_TDSensitiveDetector::m_nHitID

private

Definition at line 53 of file AFP_TDSensitiveDetector.h.

TDMaxQEff

double AFP_TDSensitiveDetector::TDMaxQEff = 0.15

staticconstexpr

Definition at line 48 of file AFP_TDSensitiveDetector.h.

---

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

• AFP_TDSensitiveDetector.h
• AFP_TDSensitiveDetector.cxx