CustomPhysics/src/DummyChargeFlipProcess.cc

0001 #include "SimG4Core/CustomPhysics/interface/DummyChargeFlipProcess.h"
0002
0003 #include <iostream>
0004 #include "G4ParticleTable.hh"
0005 #include "Randomize.hh"
0006 #include "G4NeutronElasticXS.hh"
0007 #include "G4Step.hh"
0008 #include "G4TrackStatus.hh"
0009 #include "G4Element.hh"
0010
0011 using namespace CLHEP;
0012
0013 DummyChargeFlipProcess::DummyChargeFlipProcess(const G4String& pname) : G4HadronicProcess(pname, fHadronic) {
0014   AddDataSet(new G4NeutronElasticXS());
0015   fPartChange = new G4ParticleChange();
0016 }
0017
0018 DummyChargeFlipProcess::~DummyChargeFlipProcess() { delete fPartChange; }
0019
0020 G4bool DummyChargeFlipProcess::IsApplicable(const G4ParticleDefinition& aParticleType) {
0021   return (aParticleType.GetParticleType() == "rhadron");
0022 }
0023
0024 G4VParticleChange* DummyChargeFlipProcess::PostStepDoIt(const G4Track& aTrack, const G4Step&) {
0025   fPartChange->Initialize(aTrack);
0026   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0027
0028   G4double ParentEnergy = aParticle->GetTotalEnergy();
0029   const G4ThreeVector& ParentDirection(aParticle->GetMomentumDirection());
0030
0031   G4double energyDeposit = 0.0;
0032   G4double finalGlobalTime = aTrack.GetGlobalTime();
0033
0034   G4int numberOfSecondaries = 1;
0035   fPartChange->SetNumberOfSecondaries(numberOfSecondaries);
0036   const G4TouchableHandle& thand = aTrack.GetTouchableHandle();
0037
0038   // get current position of the track
0039   aTrack.GetPosition();
0040   // create a new track object
0041   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0042   float randomParticle = G4UniformRand();
0043   G4ParticleDefinition* newType;
0044   if (randomParticle < 0.333)
0045     newType = particleTable->FindParticle(1009213);
0046   else if (randomParticle > 0.667)
0047     newType = particleTable->FindParticle(-1009213);
0048   else
0049     newType = particleTable->FindParticle(1009113);
0050
0051   //G4cout << "RHADRON: New charge = " << newType->GetPDGCharge() << G4endl;
0052
0053   G4DynamicParticle* newP = new G4DynamicParticle(newType, ParentDirection, ParentEnergy);
0054   G4Track* secondary = new G4Track(newP, finalGlobalTime, aTrack.GetPosition());
0055   // switch on good for tracking flag
0056   secondary->SetGoodForTrackingFlag();
0057   secondary->SetTouchableHandle(thand);
0058   // add the secondary track in the List
0059   fPartChange->AddSecondary(secondary);
0060
0061   // Kill the parent particle
0062   fPartChange->ProposeTrackStatus(fStopAndKill);
0063   fPartChange->ProposeLocalEnergyDeposit(energyDeposit);
0064   fPartChange->ProposeGlobalTime(finalGlobalTime);
0065
0066   ClearNumberOfInteractionLengthLeft();
0067
0068   return fPartChange;
0069 }