Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $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  {
46 }
47 
49 {
50  delete theParticleChange;
51 }
52 
53 void G4TheoFSGenerator::ModelDescription(std::ostream& outFile) const
54 {
55  outFile << GetModelName() <<" consists of a " << theHighEnergyGenerator->GetModelName()
56  << " string model and a stage to de-excite the excited nuclear fragment."
57  << "\n<p>"
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  << "<p>The string model:\n";
64  outFile <<"\n<p>";
66 }
67 
68 
70 {
71  // init particle change
74  G4double timePrimary=thePrimary.GetGlobalTime();
75 
76  // check if models have been registered, and use default, in case this is not true @@
77 
78  const G4DynamicParticle aPart(thePrimary.GetDefinition(),thePrimary.Get4Momentum().vect());
79 
80  if ( theQuasielastic ) {
81 
82  if ( theQuasielastic->GetFraction(theNucleus, aPart) > G4UniformRand() )
83  {
84  //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl;
85  G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, aPart);
86  //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl;
87  if (result)
88  {
89  for(unsigned int i=0; i<result->size(); i++)
90  {
91  G4DynamicParticle * aNew =
92  new G4DynamicParticle(result->operator[](i)->GetDefinition(),
93  result->operator[](i)->Get4Momentum().e(),
94  result->operator[](i)->Get4Momentum().vect());
96  delete result->operator[](i);
97  }
98  delete result;
99 
100  } else
101  {
105 
106  }
107  return theParticleChange;
108  }
109  }
110 
111  // get result from high energy model
112 
113  G4KineticTrackVector * theInitialResult =
114  theHighEnergyGenerator->Scatter(theNucleus, aPart);
115 
116 //#define DEBUG_initial_result
117  #ifdef DEBUG_initial_result
118  G4double E_out(0);
120  std::vector<G4KineticTrack *>::iterator ir_iter;
121  for(ir_iter=theInitialResult->begin(); ir_iter!=theInitialResult->end(); ir_iter++)
122  {
123  //G4cout << "TheoFS secondary, mom " << (*ir_iter)->GetDefinition()->GetParticleName() << " " << (*ir_iter)->Get4Momentum() << G4endl;
124  E_out += (*ir_iter)->Get4Momentum().e();
125  }
126  G4double init_mass= ionTable->GetIonMass(theNucleus.GetZ_asInt(),theNucleus.GetA_asInt());
127  G4double init_E=aPart.Get4Momentum().e();
128  // residual nucleus
129 
130  const std::vector<G4Nucleon> & thy = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
131 
132  G4int resZ(0),resA(0);
133  G4double delta_m(0);
134  for(size_t them=0; them<thy.size(); them++)
135  {
136  if(thy[them].AreYouHit()) {
137  ++resA;
138  if ( thy[them].GetDefinition() == G4Proton::Proton() ) {
139  ++resZ;
140  delta_m +=G4Proton::Proton()->GetPDGMass();
141  } else {
142  delta_m +=G4Neutron::Neutron()->GetPDGMass();
143  }
144  }
145  }
146 
147  G4double final_mass(0);
148  if ( theNucleus.GetA_asInt() ) {
149  final_mass=ionTable->GetIonMass(theNucleus.GetZ_asInt()-resZ,theNucleus.GetA_asInt()- resA);
150  }
151  G4double E_excit=init_mass + init_E - final_mass - E_out;
152  G4cout << " Corrected delta mass " << init_mass - final_mass - delta_m << G4endl;
153  G4cout << "initial E, mass = " << init_E << ", " << init_mass << G4endl;
154  G4cout << " final E, mass = " << E_out <<", " << final_mass << " excitation_E " << E_excit << G4endl;
155  #endif
156 
157  G4ReactionProductVector * theTransportResult = NULL;
158 
159 // Uzhi Nov. 2012
160  G4V3DNucleus* theProjectileNucleus = theHighEnergyGenerator->GetProjectileNucleus();
161 if(theProjectileNucleus == 0) // Uzhi Nov. 2012
162 { // Uzhi Nov. 2012
163 
164  G4int hitCount = 0;
165  const std::vector<G4Nucleon>& they = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
166  for(size_t them=0; them<they.size(); them++)
167  {
168  if(they[them].AreYouHit()) hitCount ++;
169  }
171  {
172  theTransport->SetPrimaryProjectile(thePrimary); // For Bertini Cascade
173  theTransportResult =
175  if ( !theTransportResult ) {
176  G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
177  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
178  }
179  }
180  else
181  {
182  theTransportResult = theDecay.Propagate(theInitialResult, theHighEnergyGenerator->GetWoundedNucleus());
183  if ( !theTransportResult ) {
184  G4cout << "G4TheoFSGenerator: null ptr from decay propagate " << G4endl;
185  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from decay propagate");
186  }
187  }
188 
189 } else // Uzhi Nov. 2012
190 { // Uzhi Nov. 2012
191  theTransport->SetPrimaryProjectile(thePrimary);
192  theTransportResult =
193  theTransport->PropagateNuclNucl(theInitialResult,
196  if ( !theTransportResult ) {
197  G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
198  throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
199  }
200 } // Uzhi Nov. 2012
201 
202  // Fill particle change
203  for(auto i=theTransportResult->begin(); i!=theTransportResult->end(); i++)
204  {
205  G4DynamicParticle * aNewDP =
206  new G4DynamicParticle((*i)->GetDefinition(),
207  (*i)->GetTotalEnergy(),
208  (*i)->GetMomentum());
209  G4HadSecondary aNew = G4HadSecondary(aNewDP);
210  G4double time=(*i)->GetFormationTime();
211  if(time < 0.0) { time = 0.0; }
212  aNew.SetTime(timePrimary + time);
213  aNew.SetCreatorModelType((*i)->GetCreatorModel());
215  delete (*i);
216  }
217 
218  // some garbage collection
219  delete theTransportResult;
220 
221  // Done
222  return theParticleChange;
223 }
224 
225 std::pair<G4double, G4double> G4TheoFSGenerator::GetEnergyMomentumCheckLevels() const
226 {
227  if ( theHighEnergyGenerator ) {
229  } else {
230  return std::pair<G4double, G4double>(DBL_MAX, DBL_MAX);
231  }
232 }
std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
const XML_Char * name
Definition: expat.h:151
void ModelDescription(std::ostream &outFile) const
static G4ParticleTable * GetParticleTable()
virtual const std::vector< G4Nucleon > & GetNucleons()=0
void SetMomentumChange(const G4ThreeVector &aV)
#define G4endl
Definition: G4ios.hh:61
G4KineticTrackVector * Scatter(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
virtual G4ReactionProductVector * PropagateNuclNucl(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4V3DNucleus *theProjectileNucleus)
G4IonTable * GetIonTable() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4VHighEnergyGenerator * theHighEnergyGenerator
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4TheoFSGenerator(const G4String &name="TheoFSGenerator")
G4double GetPDGMass() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1336
void SetEnergyChange(G4double anEnergy)
double G4double
Definition: G4Types.hh:76
G4double GetGlobalTime() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
virtual G4String GetModelName() const
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetFraction(G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)
void SetPrimaryProjectile(const G4HadProjectile &aPrimary)
std::vector< G4ReactionProduct * > G4ReactionProductVector
const G4LorentzVector & Get4Momentum() const
virtual G4int GetMassNumber()=0
Hep3Vector unit() const
G4double G4ParticleHPJENDLHEData::G4double result
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *)
G4double GetKineticEnergy() const
const G4ParticleDefinition * GetDefinition() const
virtual G4KineticTrackVector * Scatter(const G4Nucleus &theNucleus, const G4DynamicParticle &thePrimary)=0
int G4int
Definition: G4Types.hh:78
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
virtual void ModelDescription(std::ostream &) const
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
virtual G4V3DNucleus * GetWoundedNucleus() const =0
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
G4GLOB_DLL std::ostream G4cout
virtual void PropagateModelDescription(std::ostream &outFile) const
Hep3Vector vect() const
virtual G4V3DNucleus * GetProjectileNucleus() const
G4DecayStrongResonances theDecay
G4QuasiElasticChannel * theQuasielastic
#define DBL_MAX
Definition: templates.hh:83
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)=0
G4VIntraNuclearTransportModel * theTransport
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void SetStatusChange(G4HadFinalStateStatus aS)
G4HadFinalState * theParticleChange
const G4String & GetModelName() const