88{
89 static G4int eventcounter=0;
90 eventcounter++;
91 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction number starts ######### "<<eventcounter<<
G4endl;
92 G4ping debug(
"debug_G4BinaryLightIonReaction");
95 tA=targetNucleus.GetA_asInt();
96 tZ=targetNucleus.GetZ_asInt();
97
99
101
102 G4bool swapped=SetLighterAsProjectile(mom, toBreit);
103
106
107
108
109
110
111
112
113
114
115
116 if( (mom.t()-mom.mag())/pA < 50*MeV )
117 {
118
119
120 cascaders=FuseNucleiAndPrompound(mom);
121 }
122 else
123 {
124 result=Interact(mom,toBreit);
125
126 if(! result )
127 {
128 {
129
130
131 G4cerr <<
"G4BinaryLightIonReaction no final state for: " <<
G4endl;
137 G4cerr <<
" Target nucleus (A,Z)=("
138 << (swapped?pA:tA) << ","
139 << (swapped?pZ:tZ) <<
")" <<
G4endl;
140 G4cerr <<
" if frequent, please submit above information as bug report"
142
147 return &theResult;
148
149 }
150 }
151
152
153 G4double theStatisticalExEnergy = GetProjectileExcitation();
154
155
156 pInitialState = mom;
157
158 pInitialState.
setT(pInitialState.
getT() +
160
161
162 delete target3dNucleus;target3dNucleus=0;
163 delete projectile3dNucleus;projectile3dNucleus=0;
164
166
168
170
171
172 std::vector<G4ReactionProduct *>::iterator iter;
173
174
175
176
177
178
179
180 delete result;
181 result=0;
184
185 while (std::abs(momentum.e()-pspectators.
e()) > 10*MeV)
186 {
188
189
190 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
191 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
192 {
193 G4cout <<
"Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" <<
G4endl;
194 }
196 unsigned int i;
197 for(i=0; i<cascaders->size(); i++)
198 {
199 pFinalState +=
G4LorentzVector( (*cascaders)[i]->GetMomentum(), (*cascaders)[i]->GetTotalEnergy() );
200 }
201 momentum=pInitialState-pFinalState;
202 if (++loopcount > 10 )
203 {
204 if ( momentum.vect().mag() - momentum.e()> 10*keV )
205 {
206 G4cerr <<
"G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" <<
G4endl;
208 } else {
209 break;
210 }
211
212 }
213 }
214
215 if (spectatorA > 0 )
216 {
217
218 if ( momentum.vect().mag() - momentum.e()> 10*keV )
219 {
220
221 G4ReactionProductVector::iterator ispectator;
222 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
223 {
224 delete *ispectator;
225 }
226 delete spectators;
227
228 G4cout <<
"G4BinaryLightIonReaction.cc: mom check: " << momentum
229 <<
" 3.mag "<< momentum.vect().mag() <<
G4endl
230 << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
231 << pFinalState <<
" " << pspectators <<
G4endl
232 <<
" .. A,Z " << spectatorA <<
" "<< spectatorZ <<
G4endl;
233 G4cout <<
"G4BinaryLightIonReaction invalid final state for: " <<
G4endl;
239 G4cout <<
" Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt()
240 <<
"," << targetNucleus.GetZ_asInt() <<
")" <<
G4endl;
241 G4cout <<
" if frequent, please submit above information as bug report"
243
248 return &theResult;
249 }
250
251
252 DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
253 }
254 }
255
260
261
262
266 size_t i=0;
267 for(i=0; i<cascaders->size(); i++)
268 {
269 if((*cascaders)[i]->GetNewlyAdded())
270 {
273 (*cascaders)[i]->GetTotalEnergy(),
274 (*cascaders)[i]->GetMomentum() );
276 if(swapped)
277 {
278 tmp*=toBreit.inverse();
280 }
281 tmp *= toLab;
283
286
287
288
289
290 }
291 delete (*cascaders)[i];
292 }
293 delete cascaders;
294
295#ifdef debug_BLIR_result
297 G4cout <<
" Energy conservation initial/primary/nucleus/final/delta(init-final) "
300#endif
301
302 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction number ends ######### "<<eventcounter<<
G4endl;
303
304 return &theResult;
305}
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
void setVect(const Hep3Vector &)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
void SetStatusChange(G4HadFinalStateStatus aS)
G4int GetNumberOfSecondaries() const
void AddSecondary(G4DynamicParticle *aP)
void SetEnergyChange(G4double anEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4double GetTotalEnergy() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0) const
!! Only ground states are supported now
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
static G4ParticleTable * GetParticleTable()
G4IonTable * GetIonTable()