ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGPiMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGPiMinusInelastic.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 "G4RPGPiMinusInelastic.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "Randomize.hh"
31 
34  G4Nucleus& targetNucleus)
35 {
36  const G4HadProjectile* originalIncident = &aTrack;
37 
38  if (originalIncident->GetKineticEnergy()<= 0.1) {
42  return &theParticleChange;
43  }
44 
45  // create the target particle
46 
47  G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
48  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
49 
50  G4ReactionProduct currentParticle(originalIncident->GetDefinition() );
51  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
52  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
53 
54  // Fermi motion and evaporation
55  // As of Geant3, the Fermi energy calculation had not been Done
56 
57  G4double ek = originalIncident->GetKineticEnergy();
58  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
59 
60  G4double tkin = targetNucleus.Cinema( ek );
61  ek += tkin;
62  currentParticle.SetKineticEnergy( ek );
63  G4double et = ek + amas;
64  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
65  G4double pp = currentParticle.GetMomentum().mag();
66  if( pp > 0.0 ) {
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  G4ThreeVector momentum = currentParticle.GetMomentum();
81  currentParticle.SetMomentum( momentum * (p/pp) );
82  }
83 
84  G4ReactionProduct modifiedOriginal = currentParticle;
85 
86  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
87  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
88  G4bool incidentHasChanged = false;
89  G4bool targetHasChanged = false;
90  G4bool quasiElastic = false;
91  G4FastVector<G4ReactionProduct,256> vec; // vec will contain the secondary particles
92  G4int vecLen = 0;
93  vec.Initialize( 0 );
94 
95  const G4double cutOff = 0.1;
96  if( currentParticle.GetKineticEnergy() > cutOff )
97  InitialCollision(vec, vecLen, currentParticle, targetParticle,
98  incidentHasChanged, targetHasChanged);
99 
100  CalculateMomenta(vec, vecLen,
101  originalIncident, originalTarget, modifiedOriginal,
102  targetNucleus, currentParticle, targetParticle,
103  incidentHasChanged, targetHasChanged, quasiElastic);
104 
105  SetUpChange(vec, vecLen,
106  currentParticle, targetParticle,
107  incidentHasChanged);
108 
109  delete originalTarget;
110  return &theParticleChange;
111 }
112 
113 
114 // Initial Collision
115 // selects the particle types arising from the initial collision of
116 // the projectile and target nucleon. Secondaries are assigned to
117 // forward and backward reaction hemispheres, but final state energies
118 // and momenta are not calculated here.
119 
120 void
122  G4int& vecLen,
123  G4ReactionProduct& currentParticle,
124  G4ReactionProduct& targetParticle,
125  G4bool& incidentHasChanged,
126  G4bool& targetHasChanged)
127 {
128  G4double KE = currentParticle.GetKineticEnergy()/GeV;
129 
130  G4int mult;
131  G4int partType;
132  std::vector<G4int> fsTypes;
133 
134  G4double testCharge;
135  G4double testBaryon;
136  G4double testStrange;
137 
138  // Get particle types according to incident and target types
139 
140  if (targetParticle.GetDefinition() == particleDef[pro]) {
141  mult = GetMultiplicityT12(KE);
142  fsTypes = GetFSPartTypesForPimP(mult, KE);
143  partType = fsTypes[0];
144  if (partType != pro) {
145  targetHasChanged = true;
146  targetParticle.SetDefinition(particleDef[partType]);
147  }
148 
149  testCharge = 0.0;
150  testBaryon = 1.0;
151  testStrange = 0.0;
152 
153  } else { // target was a neutron
154  mult = GetMultiplicityT32(KE);
155  fsTypes = GetFSPartTypesForPimN(mult, KE);
156  partType = fsTypes[0];
157  if (partType != neu) {
158  targetHasChanged = true;
159  targetParticle.SetDefinition(particleDef[partType]);
160  }
161 
162  testCharge = -1.0;
163  testBaryon = 1.0;
164  testStrange = 0.0;
165  }
166 
167  // Remove target particle from list
168 
169  fsTypes.erase(fsTypes.begin());
170 
171  // See if the incident particle changed type
172 
173  G4int choose = -1;
174  for(G4int i=0; i < mult-1; ++i ) {
175  partType = fsTypes[i];
176  if (partType == pim) {
177  choose = i;
178  break;
179  }
180  }
181  if (choose == -1) {
182  incidentHasChanged = true;
183  choose = G4int(G4UniformRand()*(mult-1) );
184  partType = fsTypes[choose];
185  currentParticle.SetDefinition(particleDef[partType]);
186  }
187 
188  fsTypes.erase(fsTypes.begin()+choose);
189 
190  // Remaining particles are secondaries. Put them into vec.
191 
192  G4ReactionProduct* rp(0);
193  for(G4int i=0; i < mult-2; ++i ) {
194  partType = fsTypes[i];
195  rp = new G4ReactionProduct();
196  rp->SetDefinition(particleDef[partType]);
197  (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
198  if (partType > pim && partType < pro) rp->SetMayBeKilled(false); // kaons
199  vec.SetElement(vecLen++, rp);
200  }
201 
202  // if (mult == 2 && !incidentHasChanged && !targetHasChanged)
203  // quasiElastic = true;
204 
205  // Check conservation of charge, strangeness, baryon number
206 
207  CheckQnums(vec, vecLen, currentParticle, targetParticle,
208  testCharge, testBaryon, testStrange);
209 
210  return;
211 }