ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGPiPlusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGPiPlusInelastic.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include "G4RPGPiPlusInelastic.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "Randomize.hh"
31 
34  G4Nucleus& targetNucleus)
35 {
36  const G4HadProjectile *originalIncident = &aTrack;
37  if (originalIncident->GetKineticEnergy()<= 0.1) {
41  return &theParticleChange;
42  }
43 
44  // create the target particle
45 
46  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
47  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
48 
49  G4ReactionProduct currentParticle(originalIncident->GetDefinition() );
50  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
51  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
52 
53  // Fermi motion and evaporation
54  // As of Geant3, the Fermi energy calculation had not been Done
55 
56  G4double ek = originalIncident->GetKineticEnergy();
57  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
58 
59  G4double tkin = targetNucleus.Cinema( ek );
60  ek += tkin;
61  currentParticle.SetKineticEnergy( ek );
62  G4double et = ek + amas;
63  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
64  G4double pp = currentParticle.GetMomentum().mag();
65  if( pp > 0.0 )
66  {
67  G4ThreeVector momentum = currentParticle.GetMomentum();
68  currentParticle.SetMomentum( momentum * (p/pp) );
69  }
70 
71  // calculate black track energies
72 
73  tkin = targetNucleus.EvaporationEffects( ek );
74  ek -= tkin;
75  currentParticle.SetKineticEnergy( ek );
76  et = ek + amas;
77  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
78  pp = currentParticle.GetMomentum().mag();
79  if( pp > 0.0 )
80  {
81  G4ThreeVector momentum = currentParticle.GetMomentum();
82  currentParticle.SetMomentum( momentum * (p/pp) );
83  }
84 
85  G4ReactionProduct modifiedOriginal = currentParticle;
86 
87  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
88  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
89  G4bool incidentHasChanged = false;
90  G4bool targetHasChanged = false;
91  G4bool quasiElastic = false;
92  G4FastVector<G4ReactionProduct,256> vec; // vec will contain the secondary particles
93  G4int vecLen = 0;
94  vec.Initialize( 0 );
95 
96  const G4double cutOff = 0.1;
97  if( currentParticle.GetKineticEnergy() > cutOff )
98  InitialCollision(vec, vecLen, currentParticle, targetParticle,
99  incidentHasChanged, targetHasChanged);
100 
101  CalculateMomenta( vec, vecLen,
102  originalIncident, originalTarget, modifiedOriginal,
103  targetNucleus, currentParticle, targetParticle,
104  incidentHasChanged, targetHasChanged, quasiElastic );
105 
106  SetUpChange( vec, vecLen,
107  currentParticle, targetParticle,
108  incidentHasChanged );
109 
110  delete originalTarget;
111  return &theParticleChange;
112 }
113 
114 
115 // Initial Collision
116 // selects the particle types arising from the initial collision of
117 // the projectile and target nucleon. Secondaries are assigned to
118 // forward and backward reaction hemispheres, but final state energies
119 // and momenta are not calculated here.
120 
121 void
123  G4int& vecLen,
124  G4ReactionProduct& currentParticle,
125  G4ReactionProduct& targetParticle,
126  G4bool& incidentHasChanged,
127  G4bool& targetHasChanged)
128 {
129  G4double KE = currentParticle.GetKineticEnergy()/GeV;
130 
131  G4int mult;
132  G4int partType;
133  std::vector<G4int> fsTypes;
134 
135  G4double testCharge;
136  G4double testBaryon;
137  G4double testStrange;
138 
139  // Get particle types according to incident and target types
140 
141  if (targetParticle.GetDefinition() == particleDef[pro]) {
142  mult = GetMultiplicityT32(KE);
143  fsTypes = GetFSPartTypesForPipP(mult, KE);
144  partType = fsTypes[0];
145  if (partType != pro) {
146  targetHasChanged = true;
147  targetParticle.SetDefinition(particleDef[partType]);
148  }
149 
150  testCharge = 2.0;
151  testBaryon = 1.0;
152  testStrange = 0.0;
153 
154  } else { // target was a neutron
155  mult = GetMultiplicityT12(KE);
156  fsTypes = GetFSPartTypesForPipN(mult, KE);
157  partType = fsTypes[0];
158  if (partType != neu) {
159  targetHasChanged = true;
160  targetParticle.SetDefinition(particleDef[partType]);
161  }
162 
163  testCharge = 1.0;
164  testBaryon = 1.0;
165  testStrange = 0.0;
166  }
167 
168  // Remove target particle from list
169 
170  fsTypes.erase(fsTypes.begin());
171 
172  // See if the incident particle changed type
173 
174  G4int choose = -1;
175  for(G4int i=0; i < mult-1; ++i ) {
176  partType = fsTypes[i];
177  if (partType == pip) {
178  choose = i;
179  break;
180  }
181  }
182  if (choose == -1) {
183  incidentHasChanged = true;
184  choose = G4int(G4UniformRand()*(mult-1) );
185  partType = fsTypes[choose];
186  currentParticle.SetDefinition(particleDef[partType]);
187  }
188  fsTypes.erase(fsTypes.begin()+choose);
189 
190  // Remaining particles are secondaries. Put them into vec.
191  // Improve this by randomizing secondary order, then alternate
192  // which secondary is put into forward or backward hemisphere
193 
194  G4ReactionProduct* rp(0);
195  for(G4int i=0; i < mult-2; ++i ) {
196  partType = fsTypes[i];
197  rp = new G4ReactionProduct();
198  rp->SetDefinition(particleDef[partType]);
199  (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
200  if (partType > pim && partType < pro) rp->SetMayBeKilled(false); // kaons
201  vec.SetElement(vecLen++, rp);
202  }
203 
204  // if (mult == 2 && !incidentHasChanged && !targetHasChanged)
205  // quasiElastic = true;
206 
207  // Check conservation of charge, strangeness, baryon number
208 
209  CheckQnums(vec, vecLen, currentParticle, targetParticle,
210  testCharge, testBaryon, testStrange);
211 
212  return;
213 }