Geant4_10
G4TheoFSGenerator.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // G4TheoFSGenerator
29 //
30 // 20110307 M. Kelsey -- Add call to new theTransport->SetPrimaryProjectile()
31 // to provide access to full initial state (for Bertini)
32 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
33 
34 #include "G4DynamicParticle.hh"
35 #include "G4TheoFSGenerator.hh"
37 #include "G4ReactionProduct.hh"
38 #include "G4IonTable.hh"
39 
41  : G4HadronicInteraction(name)
42  , theTransport(0), theHighEnergyGenerator(0)
43  , theQuasielastic(0)
44  {
45  theParticleChange = new G4HadFinalState;
46 }
47 
49 {
50  delete theParticleChange;
51 }
52 
54 {
55  outFile << GetModelName() <<" consists of a " << theHighEnergyGenerator->GetModelName()
56  << " string model and of "
57  << ".\n"
58  << "The string model simulates the interaction of\n"
59  << "an incident hadron with a nucleus, forming \n"
60  << "excited strings, decays these strings into hadrons,\n"
61  << "and leaves an excited nucleus.\n"
62  << "The string model:\n";
63  theHighEnergyGenerator->ModelDescription(outFile);
64 //theTransport->IntraNuclearTransportDescription(outFile)
65 }
66 
67 
69 {
70  // init particle change
71  theParticleChange->Clear();
72  theParticleChange->SetStatusChange(stopAndKill);
73 
74  // check if models have been registered, and use default, in case this is not true @@
75 
76  // get result from high energy model
77  G4DynamicParticle aTemp(const_cast<G4ParticleDefinition *>(thePrimary.GetDefinition()),
78  thePrimary.Get4Momentum().vect());
79  const G4DynamicParticle * aPart = &aTemp;
80 
81  if ( theQuasielastic ) {
82 
83  if ( theQuasielastic->GetFraction(theNucleus, *aPart) > G4UniformRand() )
84  {
85  //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl;
86  G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart);
87  //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl;
88  if (result)
89  {
90  for(unsigned int i=0; i<result->size(); i++)
91  {
92  G4DynamicParticle * aNew =
93  new G4DynamicParticle(result->operator[](i)->GetDefinition(),
94  result->operator[](i)->Get4Momentum().e(),
95  result->operator[](i)->Get4Momentum().vect());
96  theParticleChange->AddSecondary(aNew);
97  delete result->operator[](i);
98  }
99  delete result;
100 
101  } else
102  {
103  theParticleChange->SetStatusChange(isAlive);
104  theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
105  theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
106 
107  }
108  return theParticleChange;
109  }
110  }
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  }
169  if(hitCount != theHighEnergyGenerator->GetWoundedNucleus()->GetMassNumber() )
170  {
171  theTransport->SetPrimaryProjectile(thePrimary); // For Bertini Cascade
172  theTransportResult =
173  theTransport->Propagate(theInitialResult, theHighEnergyGenerator->GetWoundedNucleus());
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,
193  theHighEnergyGenerator->GetWoundedNucleus(),
194  theHighEnergyGenerator->GetProjectileNucleus());
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  unsigned int i;
203  for(i=0; i<theTransportResult->size(); i++)
204  {
205  G4DynamicParticle * aNew =
206  new G4DynamicParticle(theTransportResult->operator[](i)->GetDefinition(),
207  theTransportResult->operator[](i)->GetTotalEnergy(),
208  theTransportResult->operator[](i)->GetMomentum());
209  // @@@ - overkill? G4double newTime = theParticleChange->GetGlobalTime(theTransportResult->operator[](i)->GetFormationTime());
210  theParticleChange->AddSecondary(aNew);
211  delete theTransportResult->operator[](i);
212  }
213 
214  // some garbage collection
215  delete theTransportResult;
216 
217  // Done
218  return theParticleChange;
219 }
220 
221 std::pair<G4double, G4double> G4TheoFSGenerator::GetEnergyMomentumCheckLevels() const
222 {
223  if ( theHighEnergyGenerator ) {
224  return theHighEnergyGenerator->GetEnergyMomentumCheckLevels();
225  } else {
226  return std::pair<G4double, G4double>(DBL_MAX, DBL_MAX);
227  }
228 }
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)=0
G4TheoFSGenerator(const G4String &name="TheoFSGenerator")
std::ofstream outFile
Definition: GammaRayTel.cc:68
void ModelDescription(std::ostream &outFile) const
virtual G4int GetMassNumber()=0
G4double G4NeutronHPJENDLHEData::G4double result
const XML_Char * name
Definition: expat.h:151
const G4String & GetModelName() const
int G4int
Definition: G4Types.hh:78
G4KineticTrackVector * Scatter(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
virtual G4KineticTrackVector * Scatter(const G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)=0
void SetStatusChange(G4HadFinalStateStatus aS)
void SetPrimaryProjectile(const G4HadProjectile &aPrimary)
virtual G4V3DNucleus * GetWoundedNucleus() const =0
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1232
G4double GetKineticEnergy() const
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
virtual G4V3DNucleus * GetProjectileNucleus() const
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetFraction(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
const G4LorentzVector & Get4Momentum() const
virtual G4String GetModelName() const
void SetEnergyChange(G4double anEnergy)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
virtual G4ReactionProductVector * PropagateNuclNucl(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4V3DNucleus *theProjectileNucleus)
Hep3Vector unit() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
virtual void ModelDescription(std::ostream &) const
virtual const std::vector< G4Nucleon > & GetNucleons()=0
#define G4endl
Definition: G4ios.hh:61
std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
void SetMomentumChange(const G4ThreeVector &aV)
void AddSecondary(G4DynamicParticle *aP)
G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)