Skip to content

Commit 82bdf5c

Browse files
ctdaxctdax
authored
code-format
1 parent 22f2e00 commit 82bdf5c

File tree

1 file changed

+80
-43
lines changed

1 file changed

+80
-43
lines changed

‎SimG4Core/CustomPhysics/src/FullModelHadronicProcess.cc‎

Lines changed: 80 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -54,8 +54,9 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
5454
// Initialize parameters
5555
aParticleChange.Initialize(aTrack);
5656
const G4DynamicParticle* incomingRhadron = aTrack.GetDynamicParticle();
57-
CustomParticle* customIncomingRhadron = static_cast<CustomParticle*>(incomingRhadron->GetDefinition()); // This is used to get the cloud particle
58-
const G4ThreeVector& aPosition = aTrack.GetPosition(); // Position of the track
57+
CustomParticle* customIncomingRhadron =
58+
static_cast<CustomParticle*>(incomingRhadron->GetDefinition()); // This is used to get the cloud particle
59+
const G4ThreeVector& aPosition = aTrack.GetPosition(); // Position of the track
5960
const G4int incomingRhadronPDG = incomingRhadron->GetDefinition()->GetPDGEncoding();
6061
G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
6162
G4bool incomingRhadronSurvives = false;
@@ -76,18 +77,23 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
7677

7778
// Define the gluino and quark cloud G4LorentzVector (momentum, total energy) based on the momentum of the R-hadron and the ratio of the masses
7879
double scale = cloudParticle->GetDefinition()->GetPDGMass() / incomingRhadron->GetDefinition()->GetPDGMass();
79-
G4LorentzVector cloudMomentum(incomingRhadron->GetMomentum() * scale, std::sqrt(incomingRhadron->GetMomentum().mag() * scale * incomingRhadron->GetMomentum().mag() * scale + cloudParticle->GetDefinition()->GetPDGMass() * cloudParticle->GetDefinition()->GetPDGMass()));
80+
G4LorentzVector cloudMomentum(
81+
incomingRhadron->GetMomentum() * scale,
82+
std::sqrt(incomingRhadron->GetMomentum().mag() * scale * incomingRhadron->GetMomentum().mag() * scale +
83+
cloudParticle->GetDefinition()->GetPDGMass() * cloudParticle->GetDefinition()->GetPDGMass()));
8084
cloudParticle->Set4Momentum(cloudMomentum);
81-
G4LorentzVector gluinoMomentum(incomingRhadron->GetMomentum() * (1. - scale), incomingRhadron->GetTotalEnergy() - cloudParticle->GetTotalEnergy());
85+
G4LorentzVector gluinoMomentum(incomingRhadron->GetMomentum() * (1. - scale),
86+
incomingRhadron->GetTotalEnergy() - cloudParticle->GetTotalEnergy());
8287

8388
// Update the cloud kinetic energy based on the target nucleus and evaporative effects
84-
G4double cloudKineticEnergy = cloudParticle->GetKineticEnergy();
89+
G4double cloudKineticEnergy = cloudParticle->GetKineticEnergy();
8590
G4double initialCloudEnergy = cloudParticle->GetTotalEnergy();
8691
cloudKineticEnergy += targetNucleus.Cinema(cloudKineticEnergy);
8792
cloudKineticEnergy -= targetNucleus.EvaporationEffects(cloudKineticEnergy);
8893

8994
G4ThreeVector cloud3MomentumDirection = cloudParticle->GetMomentum().unit();
90-
G4double cloud3MomentumMagnitudeAfterEvaporativeEffects = std::sqrt(cloudKineticEnergy * (cloudKineticEnergy + 2. * cloudParticle->GetDefinition()->GetPDGMass()));
95+
G4double cloud3MomentumMagnitudeAfterEvaporativeEffects =
96+
std::sqrt(cloudKineticEnergy * (cloudKineticEnergy + 2. * cloudParticle->GetDefinition()->GetPDGMass()));
9197

9298
// If the R-hadron kinetic energy is less than 0.1 MeV, or the cloud kinetic energy is less than or equal to 0, stop the track but keep it alive. This should be very rare.
9399
if (cloudKineticEnergy + gluinoMomentum.e() - gluinoMomentum.m() <= 0.1 * MeV || cloudKineticEnergy <= 0.) {
@@ -113,17 +119,20 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
113119
TargetSurvives = true;
114120
}
115121

116-
if (finalStateParticle->GetParticleType()=="rhadron") {
122+
if (finalStateParticle->GetParticleType() == "rhadron") {
117123
outgoingRhadronDefinition = finalStateParticle;
118124
outgoingCloudDefinition = finalStateCustomParticle->GetCloud();
119125
if (outgoingCloudDefinition == nullptr) {
120-
G4cerr << "FullModelHadronicProcess::PostStepDoIt Definition of outgoing particle cloud not available!" << G4endl;
126+
G4cerr << "FullModelHadronicProcess::PostStepDoIt Definition of outgoing particle cloud not available!"
127+
<< G4endl;
121128
exit(EXIT_FAILURE);
122129
}
123130
}
124131

125-
if (finalStateParticle == G4Proton::Proton() || finalStateParticle == G4Neutron::Neutron()) outgoingTargetDefinition = finalStateParticle;
126-
if (finalStateCustomParticle == nullptr && reactionProduct.size() == 2) outgoingTargetDefinition = finalStateParticle;
132+
if (finalStateParticle == G4Proton::Proton() || finalStateParticle == G4Neutron::Neutron())
133+
outgoingTargetDefinition = finalStateParticle;
134+
if (finalStateCustomParticle == nullptr && reactionProduct.size() == 2)
135+
outgoingTargetDefinition = finalStateParticle;
127136
if (finalStateParticle->GetPDGEncoding() == incomingRhadronPDG) {
128137
incomingRhadronSurvives = true;
129138
} else {
@@ -146,8 +155,7 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
146155
if (TargetSurvives) {
147156
outgoingTargetG4Dynamic->SetDefinition(aTarget);
148157
outgoingTargetG4Reaction = G4ReactionProduct(aTarget);
149-
}
150-
else {
158+
} else {
151159
outgoingTargetG4Dynamic->SetDefinition(outgoingTargetDefinition);
152160
outgoingTargetG4Reaction = G4ReactionProduct(outgoingTargetDefinition);
153161
}
@@ -157,7 +165,8 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
157165
G4LorentzRotation cloudParticleToLabFrameRotation = incomingCloudG4HadProjectile->GetTrafoToLab();
158166

159167
//Create a G4ReactionProduct object for the outgoing cloud
160-
G4ReactionProduct outgoingCloudG4Reaction(const_cast<G4ParticleDefinition*>(incomingCloudG4HadProjectile->GetDefinition()));
168+
G4ReactionProduct outgoingCloudG4Reaction(
169+
const_cast<G4ParticleDefinition*>(incomingCloudG4HadProjectile->GetDefinition()));
161170
outgoingCloudG4Reaction.SetMomentum(incomingCloudG4HadProjectile->Get4Momentum().vect());
162171
outgoingCloudG4Reaction.SetTotalEnergy(incomingCloudG4HadProjectile->GetTotalEnergy());
163172
if (!incomingRhadronSurvives) {
@@ -167,7 +176,7 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
167176
G4ReactionProduct modifiedoutgoingCloudG4Reaction = outgoingCloudG4Reaction;
168177

169178
//Set the hemisphere of the current and target particles. Initialize an empty vector for the secondary particles
170-
outgoingCloudG4Reaction.SetSide(1); // incident always goes in forward hemisphere
179+
outgoingCloudG4Reaction.SetSide(1); // incident always goes in forward hemisphere
171180
outgoingTargetG4Reaction.SetSide(-1); // target always goes in backward hemisphere
172181
G4bool quasiElastic = false;
173182
if (reactionProduct.size() == 2)
@@ -178,11 +187,15 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
178187

179188
//Fill the vector with the secondary particles
180189
for (G4int i = 0; i != reactionProductSize; i++) {
181-
if (outgoingParticleDefinitions[i] != aTarget && outgoingParticleDefinitions[i] != incomingCloudG4HadProjectile->GetDefinition() &&
182-
outgoingParticleDefinitions[i] != outgoingRhadronDefinition && outgoingParticleDefinitions[i] != outgoingTargetDefinition) {
190+
if (outgoingParticleDefinitions[i] != aTarget &&
191+
outgoingParticleDefinitions[i] != incomingCloudG4HadProjectile->GetDefinition() &&
192+
outgoingParticleDefinitions[i] != outgoingRhadronDefinition &&
193+
outgoingParticleDefinitions[i] != outgoingTargetDefinition) {
183194
G4ReactionProduct* secondaryReactionProduct = new G4ReactionProduct;
184195
secondaryReactionProduct->SetDefinition(outgoingParticleDefinitions[i]);
185-
(G4UniformRand() < 0.5) ? secondaryReactionProduct->SetSide(-1) : secondaryReactionProduct->SetSide(1); //Here we randomly determine the hemisphere of the secondary particle
196+
(G4UniformRand() < 0.5)
197+
? secondaryReactionProduct->SetSide(-1)
198+
: secondaryReactionProduct->SetSide(1); //Here we randomly determine the hemisphere of the secondary particle
186199
secondaryParticleVector.SetElement(secondaryParticleVectorLen++, secondaryReactionProduct);
187200
}
188201
}
@@ -206,7 +219,7 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
206219
quasiElastic);
207220

208221
//Declare the Cloud 4-momentum after the interaction and propose an energy deposit of the difference between the incoming and outgoing quark cloud energies
209-
G4LorentzVector outgoingCloudp4Prime(outgoingCloudG4Reaction.GetMomentum(), outgoingCloudG4Reaction.GetTotalEnergy());
222+
G4LorentzVector outgoingCloudp4Prime(outgoingCloudG4Reaction.GetMomentum(), outgoingCloudG4Reaction.GetTotalEnergy());
210223
outgoingCloudp4Prime *= cloudParticleToLabFrameRotation;
211224
G4double proposedEnergyDeposit = initialCloudEnergy - outgoingCloudG4Reaction.GetTotalEnergy();
212225

@@ -229,14 +242,17 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
229242

230243
//Stop the old track
231244
aParticleChange.ProposeTrackStatus(fStopAndKill);
232-
}
245+
}
233246

234247
//If the incident particle survives update its momentum direction. Includes error handling for when the momentum is zero
235248
else {
236249
dynamicOutgoingRhadron->SetDefinition(incomingRhadron->GetDefinition());
237250
dynamicOutgoingRhadron->SetMomentum(gluinoMomentum.vect() + outgoingCloudp4Prime.vect());
238251
if (dynamicOutgoingRhadron->GetMomentum().mag() > DBL_MIN)
239-
aParticleChange.ProposeMomentumDirection(dynamicOutgoingRhadron->GetMomentum().x() / dynamicOutgoingRhadron->GetMomentum().mag(), dynamicOutgoingRhadron->GetMomentum().y() / dynamicOutgoingRhadron->GetMomentum().mag(), dynamicOutgoingRhadron->GetMomentum().z() / dynamicOutgoingRhadron->GetMomentum().mag());
252+
aParticleChange.ProposeMomentumDirection(
253+
dynamicOutgoingRhadron->GetMomentum().x() / dynamicOutgoingRhadron->GetMomentum().mag(),
254+
dynamicOutgoingRhadron->GetMomentum().y() / dynamicOutgoingRhadron->GetMomentum().mag(),
255+
dynamicOutgoingRhadron->GetMomentum().z() / dynamicOutgoingRhadron->GetMomentum().mag());
240256
else
241257
aParticleChange.ProposeMomentumDirection(1.0, 0.0, 0.0);
242258
aParticleChange.ProposeEnergy(dynamicOutgoingRhadron->GetKineticEnergy());
@@ -247,9 +263,13 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
247263
if (outgoingTargetG4Reaction.GetMass() > 0.0) // outgoingTargetG4Reaction can be eliminated in TwoBody
248264
{
249265
targetParticleG4DynamicAfterInteraction->SetDefinition(outgoingTargetG4Reaction.GetDefinition());
250-
targetParticleG4DynamicAfterInteraction->SetMomentum(outgoingTargetG4Reaction.GetMomentum().rotate(2. * pi * G4UniformRand(), incomingCloud3Momentum)); // rotate(const G4double angle, const ThreeVector &axis) const;
251-
targetParticleG4DynamicAfterInteraction->SetMomentum((cloudParticleToLabFrameRotation * targetParticleG4DynamicAfterInteraction->Get4Momentum()).vect());
252-
G4Track* targetTrackAfterInteraction = new G4Track(targetParticleG4DynamicAfterInteraction, aTrack.GetGlobalTime(), aPosition);
266+
targetParticleG4DynamicAfterInteraction->SetMomentum(outgoingTargetG4Reaction.GetMomentum().rotate(
267+
2. * pi * G4UniformRand(),
268+
incomingCloud3Momentum)); // rotate(const G4double angle, const ThreeVector &axis) const;
269+
targetParticleG4DynamicAfterInteraction->SetMomentum(
270+
(cloudParticleToLabFrameRotation * targetParticleG4DynamicAfterInteraction->Get4Momentum()).vect());
271+
G4Track* targetTrackAfterInteraction =
272+
new G4Track(targetParticleG4DynamicAfterInteraction, aTrack.GetGlobalTime(), aPosition);
253273
targetTrackAfterInteraction->SetTouchableHandle(thisTouchable);
254274
aParticleChange.AddSecondary(targetTrackAfterInteraction);
255275
}
@@ -259,8 +279,10 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
259279
G4DynamicParticle* secondaryParticleAfterInteraction = new G4DynamicParticle();
260280
secondaryParticleAfterInteraction->SetDefinition(secondaryParticleVector[i]->GetDefinition());
261281
secondaryParticleAfterInteraction->SetMomentum(secondaryParticleVector[i]->GetMomentum());
262-
secondaryParticleAfterInteraction->SetMomentum((cloudParticleToLabFrameRotation * secondaryParticleAfterInteraction->Get4Momentum()).vect());
263-
G4Track* secondaryTrackAfterInteraction = new G4Track(secondaryParticleAfterInteraction, aTrack.GetGlobalTime(), aPosition);
282+
secondaryParticleAfterInteraction->SetMomentum(
283+
(cloudParticleToLabFrameRotation * secondaryParticleAfterInteraction->Get4Momentum()).vect());
284+
G4Track* secondaryTrackAfterInteraction =
285+
new G4Track(secondaryParticleAfterInteraction, aTrack.GetGlobalTime(), aPosition);
264286
secondaryTrackAfterInteraction->SetTouchableHandle(thisTouchable);
265287
aParticleChange.AddSecondary(secondaryTrackAfterInteraction);
266288

@@ -276,38 +298,45 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
276298
}
277299

278300
void FullModelHadronicProcess::CalculateMomenta(
279-
G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& secondaryParticleVector, //Vector of secondary particles
280-
G4int& secondaryParticleVectorLen, //Length of the secondary particle vector
281-
const G4HadProjectile* incomingCloudG4HadProjectile, //The incoming cloud projectile
282-
const G4DynamicParticle* outgoingTargetG4Dynamic, //The original target particle
283-
G4ReactionProduct& modifiedoutgoingCloudG4Reaction, //Fermi motion and evap. effects included
284-
G4Nucleus& targetNucleus, //The target nucleus
285-
G4ReactionProduct& outgoingCloudG4Reaction, //The outgoing cloud G4 Reaction
286-
G4ReactionProduct& outgoingTargetG4Reaction, //The outgoing particle previously defined as original target
287-
G4bool& incomingRhadronHasChanged, //True if the R-Hadron type has changed
288-
G4bool& targetHasChanged, //True if the target particle has changed
289-
G4bool quasiElastic) //True if the reaction product size equals 2, false otherwise
301+
G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& secondaryParticleVector, //Vector of secondary particles
302+
G4int& secondaryParticleVectorLen, //Length of the secondary particle vector
303+
const G4HadProjectile* incomingCloudG4HadProjectile, //The incoming cloud projectile
304+
const G4DynamicParticle* outgoingTargetG4Dynamic, //The original target particle
305+
G4ReactionProduct& modifiedoutgoingCloudG4Reaction, //Fermi motion and evap. effects included
306+
G4Nucleus& targetNucleus, //The target nucleus
307+
G4ReactionProduct& outgoingCloudG4Reaction, //The outgoing cloud G4 Reaction
308+
G4ReactionProduct& outgoingTargetG4Reaction, //The outgoing particle previously defined as original target
309+
G4bool& incomingRhadronHasChanged, //True if the R-Hadron type has changed
310+
G4bool& targetHasChanged, //True if the target particle has changed
311+
G4bool quasiElastic) //True if the reaction product size equals 2, false otherwise
290312
{
291313
FullModelReactionDynamics theReactionDynamics;
292-
incomingCloud3Momentum = incomingCloudG4HadProjectile->Get4Momentum().v(); //Use this for rotations later
314+
incomingCloud3Momentum = incomingCloudG4HadProjectile->Get4Momentum().v(); //Use this for rotations later
293315

294316
//If the reaction is quasi-elastic, use the TwoBody method to calculate the momenta of the outgoing particles.
295317
if (quasiElastic) {
296-
theReactionDynamics.TwoBody(secondaryParticleVector, secondaryParticleVectorLen, modifiedoutgoingCloudG4Reaction, outgoingTargetG4Dynamic, outgoingCloudG4Reaction, outgoingTargetG4Reaction, targetNucleus, targetHasChanged);
318+
theReactionDynamics.TwoBody(secondaryParticleVector,
319+
secondaryParticleVectorLen,
320+
modifiedoutgoingCloudG4Reaction,
321+
outgoingTargetG4Dynamic,
322+
outgoingCloudG4Reaction,
323+
outgoingTargetG4Reaction,
324+
targetNucleus,
325+
targetHasChanged);
297326
return;
298327
}
299328

300329
//If the reaction is not quasi-elastic, update the outgoing particles momenta based on effects detailed in the functions below. Then call the TwoBody method afterwards
301330
G4ReactionProduct leadingStrangeParticle;
302-
G4bool leadFlag = MarkLeadingStrangeParticle(outgoingCloudG4Reaction, outgoingTargetG4Reaction, leadingStrangeParticle);
331+
G4bool leadFlag =
332+
MarkLeadingStrangeParticle(outgoingCloudG4Reaction, outgoingTargetG4Reaction, leadingStrangeParticle);
303333
G4bool finishedTwoClu = false;
304334
if (modifiedoutgoingCloudG4Reaction.GetTotalMomentum() / MeV < 1.0) {
305335
for (G4int i = 0; i < secondaryParticleVectorLen; i++) {
306336
delete secondaryParticleVector[i];
307337
}
308338
secondaryParticleVectorLen = 0;
309-
}
310-
else {
339+
} else {
311340
theReactionDynamics.SuppressChargedPions(secondaryParticleVector,
312341
secondaryParticleVectorLen,
313342
modifiedoutgoingCloudG4Reaction,
@@ -341,7 +370,14 @@ void FullModelHadronicProcess::CalculateMomenta(
341370
return;
342371
}
343372

344-
theReactionDynamics.TwoBody(secondaryParticleVector, secondaryParticleVectorLen, modifiedoutgoingCloudG4Reaction, outgoingTargetG4Dynamic, outgoingCloudG4Reaction, outgoingTargetG4Reaction, targetNucleus, targetHasChanged);
373+
theReactionDynamics.TwoBody(secondaryParticleVector,
374+
secondaryParticleVectorLen,
375+
modifiedoutgoingCloudG4Reaction,
376+
outgoingTargetG4Dynamic,
377+
outgoingCloudG4Reaction,
378+
outgoingTargetG4Reaction,
379+
targetNucleus,
380+
targetHasChanged);
345381
}
346382

347383
G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(const G4ReactionProduct& outgoingCloudG4Reaction,
@@ -363,7 +399,8 @@ G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(const G4ReactionProd
363399
return lead;
364400
}
365401

366-
void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& secondaryParticleVector, G4int& secondaryParticleVectorLen) {
402+
void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& secondaryParticleVector,
403+
G4int& secondaryParticleVectorLen) {
367404
G4int i;
368405
for (i = 0; i < secondaryParticleVectorLen; ++i) {
369406
G4ThreeVector momentum = secondaryParticleVector[i]->GetMomentum();

0 commit comments

Comments
 (0)