source: trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc@ 1314

Last change on this file since 1314 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 7.9 KB
Line 
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#include "G4GeneratorPrecompoundInterface.hh"
27#include "G4DynamicParticleVector.hh"
28#include "G4IonTable.hh"
29
30//
31// HPW, 10DEC 98, the decay part originally written by Gunter Folger in his FTF-test-program.
32//
33
34 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface()
35 : CaptureThreshold(80*MeV)
36 {}
37
38 G4HadFinalState* G4GeneratorPrecompoundInterface::
39 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
40 {
41 std::cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."<< G4endl;
42 std::cout << "This class is only a mediator between generator and precompound"<<G4endl;
43 std::cout << "Please remove from your physics list."<<G4endl;
44 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
45 return new G4HadFinalState;
46 }
47
48 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
49 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
50 {
51 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
52
53 // decay the strong resonances
54 G4KineticTrackVector *result1, *secondaries, *result;
55 result1=theSecondaries;
56 result=new G4KineticTrackVector();
57
58 for (unsigned int aResult=0; aResult < result1->size(); aResult++)
59 {
60 G4ParticleDefinition * pdef;
61 pdef=result1->operator[](aResult)->GetDefinition();
62 secondaries=NULL;
63 if ( pdef->IsShortLived() )
64 {
65 secondaries = result1->operator[](aResult)->Decay();
66 }
67 if ( secondaries == NULL )
68 {
69 result->push_back(result1->operator[](aResult));
70 result1->operator[](aResult)=NULL; //protect for clearAndDestroy
71 }
72 else
73 {
74 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
75 {
76 result1->push_back(secondaries->operator[](aSecondary));
77 }
78 delete secondaries;
79 }
80 }
81 std::for_each(result1->begin(), result1->end(), DeleteKineticTrack());
82 delete result1;
83
84
85
86 // prepare the fragment
87 G4Fragment anInitialState;
88 G4int anA=theNucleus->GetMassNumber();
89 G4int aZ=theNucleus->GetCharge();
90 G4int numberOfEx = 0;
91 G4int numberOfCh = 0;
92 G4int numberOfHoles = 0;
93 G4double exEnergy = 0;
94 G4ThreeVector exciton3Momentum(0,0,0);
95 // loop over secondaries
96 for(unsigned int list=0; list < result->size(); list++)
97 {
98 G4KineticTrack *aTrack = result->operator[](list);
99 if(aTrack->GetDefinition() != G4Proton::Proton() &&
100 aTrack->GetDefinition() != G4Neutron::Neutron())
101 {
102 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
103 theNew->SetMomentum(aTrack->Get4Momentum().vect());
104 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
105 theTotalResult->push_back(theNew);
106 }
107 else if(aTrack->Get4Momentum().t() - aTrack->Get4Momentum().mag()>CaptureThreshold)
108 {
109 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
110 theNew->SetMomentum(aTrack->Get4Momentum().vect());
111 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
112 theTotalResult->push_back(theNew);
113 }
114 else if(aTrack->GetPosition().mag() > theNucleus->GetNuclearRadius())
115 {
116 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
117 theNew->SetMomentum(aTrack->Get4Momentum().vect());
118 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
119 theTotalResult->push_back(theNew);
120 }
121 else
122 {
123 // within the nucleus, neutron or proton
124 // now calculate A, Z of the fragment, momentum, number of exciton states
125 anA++;;
126 numberOfEx++;
127 aZ += G4int(aTrack->GetDefinition()->GetPDGCharge());
128 numberOfCh += G4int(aTrack->GetDefinition()->GetPDGCharge());
129 exciton3Momentum += aTrack->Get4Momentum().vect();
130 exEnergy += (aTrack->Get4Momentum().t()-aTrack->Get4Momentum().m());
131 }
132 }
133
134 // loop over wounded nucleus
135 G4Nucleon * theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : NULL;
136 while(theCurrentNucleon != NULL)
137 {
138 if(theCurrentNucleon->AreYouHit())
139 {
140 numberOfHoles++;
141 numberOfEx++;
142 anA--;
143 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge());
144 exciton3Momentum -= theCurrentNucleon->Get4Momentum().vect();
145 exEnergy+=theCurrentNucleon->GetBindingEnergy();
146 }
147 theCurrentNucleon = theNucleus->GetNextNucleon();
148 }
149
150 if(!theDeExcitation)
151 {
152 // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
153 }
154 else if(0!=anA && 0!=aZ)
155 {
156 G4double residualMass =
157 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(aZ ,anA);
158 residualMass += exEnergy;
159
160 G4LorentzVector exciton4Momentum(exciton3Momentum,
161 std::sqrt(exciton3Momentum.mag2()+residualMass*residualMass));
162
163 anInitialState.SetA(anA);
164 anInitialState.SetZ(aZ);
165 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
166 anInitialState.SetNumberOfCharged(numberOfCh);
167 anInitialState.SetNumberOfHoles(numberOfHoles);
168 anInitialState.SetMomentum(exciton4Momentum);
169// anInitialState.SetExcitationEnergy(exEnergy); // now a redundant call.
170
171 // call pre-compound
172 const G4Fragment aFragment(anInitialState);
173 G4ReactionProductVector * aPreResult = theDeExcitation->DeExcite(aFragment);
174
175 // fill pre-compound part into the result, and return
176 for(unsigned int ll=0; ll<aPreResult->size(); ll++)
177 {
178 theTotalResult->push_back(aPreResult->operator[](ll));
179 }
180 delete aPreResult;
181 }
182 else
183 {
184 // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
185 }
186 // now return
187
188 std::for_each(result->begin(), result->end(), DeleteKineticTrack());
189 delete result;
190 return theTotalResult;
191 }
192
193G4double G4GeneratorPrecompoundInterface::SetCaptureThreshold(G4double value)
194{
195 G4double old=CaptureThreshold;
196 CaptureThreshold=value;
197 return old;
198
199}
Note: See TracBrowser for help on using the repository browser.