ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TheoFSGenerator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TheoFSGenerator.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 // G4TheoFSGenerator
28 //
29 // 20110307 M. Kelsey -- Add call to new theTransport->SetPrimaryProjectile()
30 // to provide access to full initial state (for Bertini)
31 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
32 
33 #include "G4DynamicParticle.hh"
34 #include "G4TheoFSGenerator.hh"
36 #include "G4ReactionProduct.hh"
37 #include "G4IonTable.hh"
38 
40  : G4HadronicInteraction(name)
41  , theTransport(0), theHighEnergyGenerator(0)
42  , theQuasielastic(0)
43  {
45 }
46 
48 {
49  delete theParticleChange;
50 }
51 
52 void G4TheoFSGenerator::ModelDescription(std::ostream& outFile) const
53 {
54  outFile << GetModelName() <<" consists of a " << theHighEnergyGenerator->GetModelName()
55  << " string model and a stage to de-excite the excited nuclear fragment."
56  << "\n<p>"
57  << "The string model simulates the interaction of\n"
58  << "an incident hadron with a nucleus, forming \n"
59  << "excited strings, decays these strings into hadrons,\n"
60  << "and leaves an excited nucleus. \n"
61  << "<p>The string model:\n";
63  outFile <<"\n<p>";
65 }
66 
67 
69 {
70  // init particle change
73  G4double timePrimary=thePrimary.GetGlobalTime();
74 
75  // check if models have been registered, and use default, in case this is not true @@
76 
77  const G4DynamicParticle aPart(thePrimary.GetDefinition(),thePrimary.Get4Momentum().vect());
78 
79  if ( theQuasielastic ) {
80 
81  if ( theQuasielastic->GetFraction(theNucleus, aPart) > G4UniformRand() )
82  {
83  //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl;
84  G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, aPart);
85  //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl;
86  if (result)
87  {
88  for(unsigned int i=0; i<result->size(); i++)
89  {
90  G4DynamicParticle * aNew =
91  new G4DynamicParticle(result->operator[](i)->GetDefinition(),
92  result->operator[](i)->Get4Momentum().e(),
93  result->operator[](i)->Get4Momentum().vect());
95  delete result->operator[](i);
96  }
97  delete result;
98 
99  } else
100  {
104 
105  }
106  return theParticleChange;
107  }
108  }
109 
110  // get result from high energy model
111 
112  G4KineticTrackVector * theInitialResult =
113  theHighEnergyGenerator->Scatter(theNucleus, aPart);
114 
115 //#define DEBUG_initial_result
116  #ifdef DEBUG_initial_result
117  G4double E_out(0);
119  std::vector<G4KineticTrack *>::iterator ir_iter;
120  for(ir_iter=theInitialResult->begin(); ir_iter!=theInitialResult->end(); ir_iter++)
121  {
122  //G4cout << "TheoFS secondary, mom " << (*ir_iter)->GetDefinition()->GetParticleName() << " " << (*ir_iter)->Get4Momentum() << G4endl;
123  E_out += (*ir_iter)->Get4Momentum().e();
124  }
125  G4double init_mass= ionTable->GetIonMass(theNucleus.GetZ_asInt(),theNucleus.GetA_asInt());
126  G4double init_E=aPart.Get4Momentum().e();
127  // residual nucleus
128 
129  const std::vector<G4Nucleon> & thy = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
130 
131  G4int resZ(0),resA(0);
132  G4double delta_m(0);
133  for(size_t them=0; them<thy.size(); them++)
134  {
135  if(thy[them].AreYouHit()) {
136  ++resA;
137  if ( thy[them].GetDefinition() == G4Proton::Proton() ) {
138  ++resZ;
139  delta_m +=G4Proton::Proton()->GetPDGMass();
140  } else {
141  delta_m +=G4Neutron::Neutron()->GetPDGMass();
142  }
143  }
144  }
145 
146  G4double final_mass(0);
147  if ( theNucleus.GetA_asInt() ) {
148  final_mass=ionTable->GetIonMass(theNucleus.GetZ_asInt()-resZ,theNucleus.GetA_asInt()- resA);
149  }
150  G4double E_excit=init_mass + init_E - final_mass - E_out;
151  G4cout << " Corrected delta mass " << init_mass - final_mass - delta_m << G4endl;
152  G4cout << "initial E, mass = " << init_E << ", " << init_mass << G4endl;
153  G4cout << " final E, mass = " << E_out <<", " << final_mass << " excitation_E " << E_excit << G4endl;
154  #endif
155 
156  G4ReactionProductVector * theTransportResult = NULL;
157 
158 // Uzhi Nov. 2012
159  G4V3DNucleus* theProjectileNucleus = theHighEnergyGenerator->GetProjectileNucleus();
160 if(theProjectileNucleus == 0) // Uzhi Nov. 2012
161 { // Uzhi Nov. 2012
162 
163  G4int hitCount = 0;
164  const std::vector<G4Nucleon>& they = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
165  for(size_t them=0; them<they.size(); them++)
166  {
167  if(they[them].AreYouHit()) hitCount ++;
168  }
170  {
171  theTransport->SetPrimaryProjectile(thePrimary); // For Bertini Cascade
172  theTransportResult =
174  if ( !theTransportResult ) {
175  G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
176  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
177  }
178  }
179  else
180  {
181  theTransportResult = theDecay.Propagate(theInitialResult, theHighEnergyGenerator->GetWoundedNucleus());
182  if ( !theTransportResult ) {
183  G4cout << "G4TheoFSGenerator: null ptr from decay propagate " << G4endl;
184  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from decay propagate");
185  }
186  }
187 
188 } else // Uzhi Nov. 2012
189 { // Uzhi Nov. 2012
190  theTransport->SetPrimaryProjectile(thePrimary);
191  theTransportResult =
192  theTransport->PropagateNuclNucl(theInitialResult,
195  if ( !theTransportResult ) {
196  G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
197  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
198  }
199 } // Uzhi Nov. 2012
200 
201  // Fill particle change
202  for(auto i=theTransportResult->begin(); i!=theTransportResult->end(); i++)
203  {
204  G4DynamicParticle * aNewDP =
205  new G4DynamicParticle((*i)->GetDefinition(),
206  (*i)->GetTotalEnergy(),
207  (*i)->GetMomentum());
208  G4HadSecondary aNew = G4HadSecondary(aNewDP);
209  G4double time=(*i)->GetFormationTime();
210  if(time < 0.0) { time = 0.0; }
211  aNew.SetTime(timePrimary + time);
212  aNew.SetCreatorModelType((*i)->GetCreatorModel());
214  delete (*i);
215  }
216 
217  // some garbage collection
218  delete theTransportResult;
219 
220  // Done
221  return theParticleChange;
222 }
223 
224 std::pair<G4double, G4double> G4TheoFSGenerator::GetEnergyMomentumCheckLevels() const
225 {
226  if ( theHighEnergyGenerator ) {
228  } else {
229  return std::pair<G4double, G4double>(DBL_MAX, DBL_MAX);
230  }
231 }