Belle II Software  release-05-01-25
G4LongLivedNeutralDecay.cc
1 #include <simulation/longlivedneutral/G4LongLivedNeutral.h>
2 #include <simulation/longlivedneutral/G4LongLivedNeutralDecay.h>
3 
4 #include "G4PhysicalConstants.hh"
5 #include "G4SystemOfUnits.hh"
6 #include "G4DynamicParticle.hh"
7 #include "G4PrimaryParticle.hh"
8 #include "G4DecayProducts.hh"
9 #include "G4DecayTable.hh"
10 #include "G4VDecayChannel.hh"
11 #include "G4PhysicsLogVector.hh"
12 #include "G4ParticleChangeForDecay.hh"
13 
14 
15 using namespace Belle2;
16 
17 // constructor
19  : G4VRestDiscreteProcess(processName, fDecay),
20  verboseLevel(1),
21  HighestValue(20.0),
22  fRemainderLifeTime(-1.0)
23 
24 {
25  // set Process Sub Type
26  SetProcessSubType(static_cast<int>(DECAY));
27 
28 #ifdef G4VERBOSE
29  if (GetVerboseLevel() > 1) {
30  G4cout << "G4LongLivedNeutralDecay constructor " << " Name:" << processName << G4endl;
31  }
32 #endif
33 
34  pParticleChange = &fParticleChangeForDecay;
35 }
36 
38 {
39 
40 }
41 
42 G4bool G4LongLivedNeutralDecay::IsApplicable(const G4ParticleDefinition& aParticleType)
43 {
44  // original G4Decay checks for particle definition mass and lifetime, which are zero for generality
45  // the correct mass and lifetime are set for the G4PrimaryParticle
46  if (!aParticleType.GetPDGStable()) {
47  return true;
48  } else return false;
49 }
50 
51 
52 G4double G4LongLivedNeutralDecay::GetMeanLifeTime(const G4Track& aTrack ,
53  G4ForceCondition*)
54 {
55  const G4DynamicParticle* dynamicParticle = aTrack.GetDynamicParticle();
56  const G4PrimaryParticle* primaryParticle = dynamicParticle->GetPrimaryParticle();
57  G4double pTime = primaryParticle->GetProperTime();
58 
59  return pTime;
60 }
61 
62 G4double G4LongLivedNeutralDecay::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*)
63 {
64 
65  const G4DynamicParticle* dynamicParticle = aTrack.GetDynamicParticle();
66  const G4PrimaryParticle* primaryParticle = dynamicParticle->GetPrimaryParticle();
67  G4double pTime = primaryParticle->GetProperTime();
68  G4double aMass = primaryParticle->GetMass();
69  G4double aMomentum = primaryParticle->GetTotalMomentum();
70  return c_light * pTime * aMomentum / aMass;
71 }
72 
73 void G4LongLivedNeutralDecay::BuildPhysicsTable(const G4ParticleDefinition&)
74 {
75  return;
76 }
77 
78 G4VParticleChange* G4LongLivedNeutralDecay::DecayIt(const G4Track& aTrack, const G4Step&)
79 {
80  // The DecayIt() method returns by pointer a particle-change object.
81  // Units are expressed in GEANT4 internal units.
82  // Initialize ParticleChange
83  // all members of G4VParticleChange are set to equal to
84  // corresponding member in G4Track
85 //
86 
87  fParticleChangeForDecay.Initialize(aTrack);
88 
89  // LongLivedNeutral: work only with DynamicParticle, not ParticleDefinition
90  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
91  const G4PrimaryParticle* aPrimaryParticle = aParticle->GetPrimaryParticle();
92 
93  //check if thePreAssignedDecayProducts exists
94  const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
95  G4bool isPreAssigned = (o_products != nullptr);
96  G4DecayProducts* products = nullptr;
97 
98  if (isPreAssigned) {
99  // default: copy decay products
100  products = new G4DecayProducts(*o_products);
101  } else {
102  // Should not happen, since interaction and free path length are set to max
103  G4ExceptionDescription ed;
104  ed << "LongLivedNeutral particle:"
105  << " decay probability exist but no pre-assigned decay products defined "
106  << "- the particle will be killed;\n";
107  G4Exception("G4LongLivedNeutralDecay::DecayIt ",
108  "DECAY101", JustWarning, ed);
109 
110  fParticleChangeForDecay.SetNumberOfSecondaries(0);
111  // Kill the parent particle
112  fParticleChangeForDecay.ProposeTrackStatus(fStopAndKill) ;
113  fParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0);
114 
115  ClearNumberOfInteractionLengthLeft();
116  return &fParticleChangeForDecay ;
117  }
118 
119  // get parent particle information ...................................
120  G4double ParentEnergy = aPrimaryParticle->GetTotalEnergy();
121  G4double ParentMass = aPrimaryParticle->GetMass();
122  if (ParentEnergy < ParentMass) {
123  G4ExceptionDescription ed;
124  ed << "Total Energy is less than its mass - increased the energy"
125  << "\n Particle: " << "LongLivedNeutral"
126  << "\n Energy:" << ParentEnergy / MeV << "[MeV]"
127  << "\n Mass:" << ParentMass / MeV << "[MeV]";
128  G4Exception("G4LongLivedNeutralDecay::DecayIt ",
129  "DECAY102", JustWarning, ed);
130  ParentEnergy = ParentMass;
131  }
132 
133  G4ThreeVector ParentDirection(aPrimaryParticle->GetMomentumDirection());
134 
135  // assign a new parent to the products from info of the primary particle
136  const G4DynamicParticle dParticle(aPrimaryParticle->GetParticleDefinition(),
137  aPrimaryParticle->GetMomentumDirection(),
138  aPrimaryParticle->GetKineticEnergy(),
139  aPrimaryParticle->GetMass());
140 
141  products->SetParentParticle(dParticle);
142  //boost all decay products to laboratory frame
143  G4double energyDeposit = 0.0;
144  G4double finalGlobalTime = aTrack.GetGlobalTime();
145  G4double finalLocalTime = aTrack.GetLocalTime();
146  if (aTrack.GetTrackStatus() == fStopButAlive) {
147  // AtRest case
148  finalGlobalTime += fRemainderLifeTime;
149  finalLocalTime += fRemainderLifeTime;
150  energyDeposit += aPrimaryParticle->GetKineticEnergy();
151  if (isPreAssigned) products->Boost(ParentEnergy, ParentDirection);
152  } else {
153  // default for LongLivedNeutral: PostStep case (decay in flight)
154  products->Boost(ParentEnergy, ParentDirection);
155  }
156 
157  //add products in fParticleChangeForDecay
158  G4int numberOfSecondaries = products->entries();
159  fParticleChangeForDecay.SetNumberOfSecondaries(numberOfSecondaries);
160 #ifdef G4VERBOSE
161  if (GetVerboseLevel() > 1) {
162  G4cout << "G4LongLivedNeutralDecay::DoIt : Decay vertex :";
163  G4cout << " Time: " << finalGlobalTime / ns << "[ns]";
164  G4cout << " proper time: " << finalLocalTime / ns << "[ns]";
165  G4cout << " X:" << (aTrack.GetPosition()).x() / cm << "[cm]";
166  G4cout << " Y:" << (aTrack.GetPosition()).y() / cm << "[cm]";
167  G4cout << " Z:" << (aTrack.GetPosition()).z() / cm << "[cm]";
168  G4cout << G4endl;
169  G4cout << "G4LongLivedNeutralDecay::DoIt : decay products in Lab. Frame" << G4endl;
170  products->DumpInfo();
171  }
172 #endif
173 
174  G4int index;
175  G4ThreeVector currentPosition;
176  const G4TouchableHandle thand = aTrack.GetTouchableHandle();
177  for (index = 0; index < numberOfSecondaries; index++) {
178  // get current position of the track
179  currentPosition = aTrack.GetPosition();
180  // create a new track object
181  G4Track* secondary = new G4Track(products->PopProducts(),
182  finalGlobalTime ,
183  currentPosition);
184  // switch on good for tracking flag
185  secondary->SetGoodForTrackingFlag();
186  secondary->SetTouchableHandle(thand);
187  // add the secondary track in the List
188  fParticleChangeForDecay.AddSecondary(secondary);
189  }
190  delete products;
191 
192  // Kill the parent particle
193  fParticleChangeForDecay.ProposeTrackStatus(fStopAndKill) ;
194  fParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit);
195  fParticleChangeForDecay.ProposeLocalTime(finalLocalTime);
196 
197  // Clear NumberOfInteractionLengthLeft
198  ClearNumberOfInteractionLengthLeft();
199 
200  return &fParticleChangeForDecay ;
201 }
202 
203 
205 {
206  currentInteractionLength = -1.0;
207  ResetNumberOfInteractionLengthLeft();
208 
209 
210  fRemainderLifeTime = -1.0;
211 }
212 
214 {
215  // Clear NumberOfInteractionLengthLeft
216  ClearNumberOfInteractionLengthLeft();
217 
218  currentInteractionLength = -1.0;
219 }
220 
221 
223  const G4Track& track,
224  G4double previousStepSize,
225  G4ForceCondition* condition
226 )
227 {
228  // condition is set to "Not Forced"
229  *condition = NotForced;
230 
231  // pre-assigned Decay time
232  G4double pTime = track.GetDynamicParticle()->GetPrimaryParticle()->GetProperTime();
233  if (pTime >= 0.) {
234  // default for LongLived Particle: pre-assigned decay time case
235  fRemainderLifeTime = pTime - track.GetProperTime();
236  if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = 0.0;
237 
238  // use normalized kinetic energy (= Ekin/mass)
239  G4double rvalue = c_light * fRemainderLifeTime;
240  const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
241  const G4PrimaryParticle* primaryParticle = dynamicParticle->GetPrimaryParticle();
242  G4double aMass = primaryParticle->GetMass();
243  G4double aMomentum = primaryParticle->GetTotalMomentum();
244  rvalue *= aMomentum / aMass;
245  return rvalue;
246 
247  } else {
248  // no pre-assigned proper decay time --> stable
249  return previousStepSize * DBL_MAX;
250  }
251 }
252 
254  const G4Track& track,
255  G4ForceCondition* condition
256 )
257 {
258  // condition is set to "Not Forced"
259  *condition = NotForced;
260  G4double pTime = track.GetDynamicParticle()->GetPrimaryParticle()->GetProperTime();
261  if (pTime >= 0.) {
262  // general case for LongLivedNeutral
263  fRemainderLifeTime = pTime - track.GetProperTime();
264  if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
265  } else {
266  // LongLivedNeutral case: set stable if there is no pre-assigned decay time
267  fRemainderLifeTime = 1e24 * s;
268  }
269  return fRemainderLifeTime;
270 }
271 
272 
274  const G4Track& aTrack,
275  const G4Step& aStep
276 )
277 {
278  if ((aTrack.GetTrackStatus() == fStopButAlive) ||
279  (aTrack.GetTrackStatus() == fStopAndKill)) {
280  fParticleChangeForDecay.Initialize(aTrack);
281  return &fParticleChangeForDecay;
282  } else {
283  return DecayIt(aTrack, aStep);
284  }
285 }
286 
Belle2::G4LongLivedNeutralDecay::IsApplicable
virtual G4bool IsApplicable(const G4ParticleDefinition &) override
returns "true" if the decay process can be applied to the particle type.
Definition: G4LongLivedNeutralDecay.cc:42
Belle2::G4LongLivedNeutralDecay::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
G4VProcess::PostStepDoIt() implemention.
Definition: G4LongLivedNeutralDecay.cc:273
Belle2::G4LongLivedNeutralDecay::GetMeanLifeTime
virtual G4double GetMeanLifeTime(const G4Track &aTrack, G4ForceCondition *condition) override
GetMeanLifeTime returns ctau for decay at rest.
Definition: G4LongLivedNeutralDecay.cc:52
Belle2::G4LongLivedNeutralDecay::BuildPhysicsTable
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
In G4Decay, thePhysicsTable stores values of beta * std::sqrt( 1 - beta*beta) as a function of normal...
Definition: G4LongLivedNeutralDecay.cc:73
Belle2::G4LongLivedNeutralDecay::DecayIt
virtual G4VParticleChange * DecayIt(const G4Track &aTrack, const G4Step &aStep)
G4Decay::DecayIt() implemention.
Definition: G4LongLivedNeutralDecay.cc:78
Belle2::G4LongLivedNeutralDecay::GetMeanFreePath
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
GetMeanFreePath returns ctau*beta*gamma for decay in flight.
Definition: G4LongLivedNeutralDecay.cc:62
Belle2::G4LongLivedNeutralDecay::fParticleChangeForDecay
G4ParticleChangeForDecay fParticleChangeForDecay
ParticleChange for decay process.
Definition: G4LongLivedNeutralDecay.h:170
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::G4LongLivedNeutralDecay::PostStepGetPhysicalInteractionLength
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4VProcess::PostStepGetPhysicalInteractionLength() implemention.
Definition: G4LongLivedNeutralDecay.cc:222
Belle2::G4LongLivedNeutralDecay::EndTracking
virtual void EndTracking() override
inform End of tracking for each track to the physics process
Definition: G4LongLivedNeutralDecay.cc:213
Belle2::G4LongLivedNeutralDecay::G4LongLivedNeutralDecay
G4LongLivedNeutralDecay(const G4String &processName="LongLivedNeutralDecay")
Constructor.
Definition: G4LongLivedNeutralDecay.cc:18
Belle2::G4LongLivedNeutralDecay::~G4LongLivedNeutralDecay
virtual ~G4LongLivedNeutralDecay()
Destructor.
Definition: G4LongLivedNeutralDecay.cc:37
Belle2::G4LongLivedNeutralDecay::fRemainderLifeTime
G4double fRemainderLifeTime
ParticleChange for decay process.
Definition: G4LongLivedNeutralDecay.h:168
Belle2::G4LongLivedNeutralDecay::AtRestGetPhysicalInteractionLength
virtual G4double AtRestGetPhysicalInteractionLength(const G4Track &track, G4ForceCondition *condition) override
G4VProcess::AtRestGetPhysicalInteractionLength() implemention.
Definition: G4LongLivedNeutralDecay.cc:253
Belle2::G4LongLivedNeutralDecay::StartTracking
virtual void StartTracking(G4Track *) override
inform Start of tracking for each track to the physics process
Definition: G4LongLivedNeutralDecay.cc:204