| Geant4 Cross Reference |
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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /*
39 * G4INCLNucleus.hh
40 *
41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #ifndef G4INCLNUCLEUS_HH_
46 #define G4INCLNUCLEUS_HH_
47
48 #include <list>
49 #include <string>
50
51 #include "G4INCLParticle.hh"
52 #include "G4INCLEventInfo.hh"
53 #include "G4INCLCluster.hh"
54 #include "G4INCLFinalState.hh"
55 #include "G4INCLStore.hh"
56 #include "G4INCLGlobals.hh"
57 #include "G4INCLParticleTable.hh"
58 #include "G4INCLConfig.hh"
59 #include "G4INCLConfigEnums.hh"
60 #include "G4INCLCluster.hh"
61 #include "G4INCLProjectileRemnant.hh"
62
63 namespace G4INCL {
64
65 enum AnnihilationType {Def=0, PType, NType, PTypeInFlight, NTypeInFlight, NbarPTypeInFlight, NbarNTypeInFlight};
66
67 class Nucleus : public Cluster {
68 public:
69 Nucleus(G4int mass, G4int charge, G4int strangess, Config const * const conf, const G4double universeRadius=-1., AnnihilationType AType=Def);
70 virtual ~Nucleus();
71
72 /// \brief Dummy copy constructor to silence Coverity warning
73 Nucleus(const Nucleus &rhs);
74
75 /// \brief Dummy assignment operator to silence Coverity warning
76 Nucleus &operator=(const Nucleus &rhs);
77
78 AnnihilationType getAType() const;
79 void setAType(AnnihilationType type);
80
81 /**
82 * Call the Cluster method to generate the initial distribution of
83 * particles. At the beginning all particles are assigned as spectators.
84 */
85 void initializeParticles();
86
87 /// \brief Insert a new particle (e.g. a projectile) in the nucleus.
88 void insertParticle(Particle *p) {
89 theZ += p->getZ();
90 theA += p->getA();
91 theS += p->getS();
92 theStore->particleHasEntered(p);
93 if(p->isNucleon()) {
94 theNpInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
95 theNnInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
96 }
97 if(p->isPion()) {
98 theNpionplusInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
99 theNpionminusInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
100 }
101 if(p->isKaon() || p->isAntiKaon()) {
102 theNkaonplusInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
103 theNkaonminusInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType()));
104 }
105 if(p->isAntiNucleon()) {
106 theNantiprotonInitial += Math::heaviside(ParticleTable::getIsospin(p->getType()));
107 }
108 if(!p->isTargetSpectator()) theStore->getBook().incrementCascading();
109 };
110
111 /**
112 * Apply reaction final state information to the nucleus.
113 */
114 void applyFinalState(FinalState *);
115
116 G4int getInitialA() const { return theInitialA; };
117 G4int getInitialZ() const { return theInitialZ; };
118 G4int getInitialS() const { return theInitialS; };
119
120 /**
121 * Propagate the particles one time step.
122 *
123 * @param step length of the time step
124 */
125 void propagateParticles(G4double step);
126
127 G4int getNumberOfEnteringProtons() const { return theNpInitial; };
128 G4int getNumberOfEnteringNeutrons() const { return theNnInitial; };
129 G4int getNumberOfEnteringPions() const { return theNpionplusInitial+theNpionminusInitial; };
130 G4int getNumberOfEnteringKaons() const { return theNkaonplusInitial+theNkaonminusInitial; };
131 G4int getNumberOfEnteringantiProtons() const { return theNantiprotonInitial; };
132
133 /** \brief Outgoing - incoming separation energies.
134 *
135 * Used by CDPP.
136 */
137 G4double computeSeparationEnergyBalance() const {
138 G4double S = 0.0;
139 ParticleList const &outgoing = theStore->getOutgoingParticles();
140 for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
141 const ParticleType t = (*i)->getType();
142 switch(t) {
143 case Proton:
144 case Neutron:
145 case DeltaPlusPlus:
146 case DeltaPlus:
147 case DeltaZero:
148 case DeltaMinus:
149 case Lambda:
150 case PiPlus:
151 case PiMinus:
152 case KPlus:
153 case KMinus:
154 case KZero:
155 case KZeroBar:
156 case KShort:
157 case KLong:
158 case SigmaPlus:
159 case SigmaZero:
160 case SigmaMinus:
161 case antiProton:
162 //case antiNeutron:
163 //case antiLambda:
164 S += thePotential->getSeparationEnergy(*i);
165 break;
166 case Composite:
167 S += (*i)->getZ() * thePotential->getSeparationEnergy(Proton)
168 + ((*i)->getA() + (*i)->getS() - (*i)->getZ()) * thePotential->getSeparationEnergy(Neutron)
169 - (*i)->getS() * thePotential->getSeparationEnergy(Lambda);
170 break;
171 default:
172 break;
173 }
174 }
175
176 S -= theNpInitial * thePotential->getSeparationEnergy(Proton);
177 S -= theNnInitial * thePotential->getSeparationEnergy(Neutron);
178 S -= theNpionplusInitial*thePotential->getSeparationEnergy(PiPlus);;
179 S -= theNkaonplusInitial*thePotential->getSeparationEnergy(KPlus);
180 S -= theNpionminusInitial*thePotential->getSeparationEnergy(PiMinus);
181 S -= theNkaonminusInitial*thePotential->getSeparationEnergy(KMinus);
182 S -= theNantiprotonInitial*thePotential->getSeparationEnergy(antiProton);
183 return S;
184 }
185
186 /** \brief Force the decay of outgoing deltas.
187 *
188 * \return true if any delta was forced to decay.
189 */
190 G4bool decayOutgoingDeltas();
191
192 /** \brief Force the decay of deltas inside the nucleus.
193 *
194 * \return true if any delta was forced to decay.
195 */
196 G4bool decayInsideDeltas();
197
198 /** \brief Force the transformation of strange particles into a Lambda;
199 *
200 * \return true if any strange particles was forced to absorb.
201 */
202 G4bool decayInsideStrangeParticles();
203
204 /** \brief Force the decay of outgoing PionResonances (eta/omega).
205 *
206 * \return true if any eta was forced to decay.
207 */
208 G4bool decayOutgoingPionResonances(G4double timeThreshold);
209
210 /** \brief Force the decay of outgoing Neutral Sigma.
211 *
212 * \return true if any Sigma was forced to decay.
213 */
214 G4bool decayOutgoingSigmaZero(G4double timeThreshold);
215
216 /** \brief Force the transformation of outgoing Neutral Kaon into propation eigenstate.
217 *
218 * \return true if any kaon was forced to decay.
219 */
220 G4bool decayOutgoingNeutralKaon();
221
222 /** \brief Force the decay of unstable outgoing clusters.
223 *
224 * \return true if any cluster was forced to decay.
225 */
226 G4bool decayOutgoingClusters();
227
228 /** \brief Force the phase-space decay of the Nucleus.
229 *
230 * Only applied if Z==0 or N==0.
231 *
232 * \return true if the nucleus was forced to decay.
233 */
234 G4bool decayMe();
235
236 /// \brief Force emission of all pions inside the nucleus.
237 void emitInsidePions();
238
239 /// \brief Force emission of all strange particles inside the nucleus.
240 void emitInsideStrangeParticles();
241
242 /// \brief Force emission of all Lambda (desexitation code with strangeness not implanted yet)
243 G4int emitInsideLambda();
244
245 /// \brief Force emission of all Kaon inside the nucleus
246 G4bool emitInsideKaon();
247
248 /** \brief Compute the recoil momentum and spin of the nucleus. */
249 void computeRecoilKinematics();
250
251 /** \brief Compute the current center-of-mass position.
252 *
253 * \return the center-of-mass position vector [fm].
254 */
255 ThreeVector computeCenterOfMass() const;
256
257 /** \brief Compute the current total energy.
258 *
259 * \return the total energy [MeV]
260 */
261 G4double computeTotalEnergy() const;
262
263 /** \brief Compute the current excitation energy.
264 *
265 * \return the excitation energy [MeV]
266 */
267 G4double computeExcitationEnergy() const;
268
269 /** \brief Set the incoming angular-momentum vector. */
270 void setIncomingAngularMomentum(const ThreeVector &j) {
271 incomingAngularMomentum = j;
272 }
273
274 /** \brief Get the incoming angular-momentum vector. */
275 const ThreeVector &getIncomingAngularMomentum() const { return incomingAngularMomentum; }
276
277 /** \brief Set the incoming momentum vector. */
278 void setIncomingMomentum(const ThreeVector &p) {
279 incomingMomentum = p;
280 }
281
282 /** \brief Get the incoming momentum vector. */
283 const ThreeVector &getIncomingMomentum() const {
284 return incomingMomentum;
285 }
286
287 /** \brief Set the initial energy. */
288 void setInitialEnergy(const G4double e) { initialEnergy = e; }
289
290 /** \brief Get the initial energy. */
291 G4double getInitialEnergy() const { return initialEnergy; }
292
293 /** \brief Get the excitation energy of the nucleus.
294 *
295 * Method computeRecoilKinematics() should be called first.
296 */
297 G4double getExcitationEnergy() const { return theExcitationEnergy; }
298
299 ///\brief Returns true if the nucleus contains any deltas.
300 inline G4bool containsDeltas() {
301 ParticleList const &inside = theStore->getParticles();
302 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
303 if((*i)->isDelta()) return true;
304 return false;
305 }
306
307 ///\brief Returns true if the nucleus contains any anti Kaons.
308 inline G4bool containsAntiKaon() {
309 ParticleList const &inside = theStore->getParticles();
310 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
311 if((*i)->isAntiKaon()) return true;
312 return false;
313 }
314
315 ///\brief Returns true if the nucleus contains any Lambda.
316 inline G4bool containsLambda() {
317 ParticleList const &inside = theStore->getParticles();
318 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
319 if((*i)->isLambda()) return true;
320 return false;
321 }
322
323 ///\brief Returns true if the nucleus contains any Sigma.
324 inline G4bool containsSigma() {
325 ParticleList const &inside = theStore->getParticles();
326 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
327 if((*i)->isSigma()) return true;
328 return false;
329 }
330
331 ///\brief Returns true if the nucleus contains any Kaons.
332 inline G4bool containsKaon() {
333 ParticleList const &inside = theStore->getParticles();
334 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
335 if((*i)->isKaon()) return true;
336 return false;
337 }
338
339 ///\brief Returns true if the nucleus contains any etas.
340 inline G4bool containsEtas() {
341 ParticleList const &inside = theStore->getParticles();
342 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
343 if((*i)->isEta()) return true;
344 return false;
345 }
346
347 ///\brief Returns true if the nucleus contains any omegas.
348 inline G4bool containsOmegas() {
349 ParticleList const &inside = theStore->getParticles();
350 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
351 if((*i)->isOmega()) return true;
352 return false;
353 }
354
355
356 /**
357 * Print the nucleus info
358 */
359 std::string print();
360
361 Store* getStore() const {return theStore; };
362 void setStore(Store *str) {
363 delete theStore;
364 theStore = str;
365 };
366
367 G4double getInitialInternalEnergy() const { return initialInternalEnergy; };
368
369 /** \brief Is the event transparent?
370 *
371 * To be called at the end of the cascade.
372 **/
373 G4bool isEventTransparent() const;
374
375 /** \brief Does the nucleus give a cascade remnant?
376 *
377 * To be called after computeRecoilKinematics().
378 **/
379 G4bool hasRemnant() const { return remnant; }
380
381 /**
382 * Fill the event info which contains INCL output data
383 */
384 void fillEventInfo(EventInfo *eventInfo);
385
386 G4bool getTryCompoundNucleus() { return tryCN; }
387
388 /// \brief Get the transmission barrier
389 G4double getTransmissionBarrier(Particle const * const p) {
390 const G4double theTransmissionRadius = theDensity->getTransmissionRadius(p);
391 const G4double theParticleZ = p->getZ();
392 return PhysicalConstants::eSquared*(theZ-theParticleZ)*theParticleZ/theTransmissionRadius;
393 }
394
395 /// \brief Struct for conservation laws
396 struct ConservationBalance {
397 ThreeVector momentum;
398 G4double energy;
399 G4int Z, A, S;
400 };
401
402 /// \brief Compute charge, mass, energy and momentum balance
403 ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const;
404
405 /// \brief Adjust the kinematics for complete-fusion events
406 void useFusionKinematics();
407
408 /** \brief Get the maximum allowed radius for a given particle.
409 *
410 * Calls the NuclearDensity::getMaxRFromP() method for nucleons and deltas,
411 * and the NuclearDensity::getTrasmissionRadius() method for pions.
412 *
413 * \param particle pointer to a particle
414 * \return surface radius
415 */
416 G4double getSurfaceRadius(Particle const * const particle) const {
417 if(particle->isNucleon() || particle->isLambda() || particle->isResonance()){
418 const G4double pr = particle->getReflectionMomentum()/thePotential->getFermiMomentum(particle);
419 if(pr>=1.)
420 return getUniverseRadius();
421 else
422 return theDensity->getMaxRFromP(particle->getType(), pr);
423 }
424 else {
425 // Temporarily set RPION = RMAX
426 return getUniverseRadius();
427 //return 0.5*(theDensity->getTransmissionRadius(particle)+getUniverseRadius());
428 }
429 }
430
431 /// \brief Getter for theUniverseRadius.
432 G4double getUniverseRadius() const { return theUniverseRadius; }
433
434 /// \brief Setter for theUniverseRadius.
435 void setUniverseRadius(const G4double universeRadius) { theUniverseRadius=universeRadius; }
436
437 /// \brief Is it a nucleus-nucleus collision?
438 G4bool isNucleusNucleusCollision() const { return isNucleusNucleus; }
439
440 /// \brief Set a nucleus-nucleus collision
441 void setNucleusNucleusCollision() { isNucleusNucleus=true; }
442
443 /// \brief Set a particle-nucleus collision
444 void setParticleNucleusCollision() { isNucleusNucleus=false; }
445
446 /// \brief Set the projectile remnant
447 void setProjectileRemnant(ProjectileRemnant * const c) {
448 delete theProjectileRemnant;
449 theProjectileRemnant = c;
450 }
451
452 /// \brief Get the projectile remnant
453 ProjectileRemnant *getProjectileRemnant() const { return theProjectileRemnant; }
454
455 /// \brief Delete the projectile remnant
456 void deleteProjectileRemnant() {
457 delete theProjectileRemnant;
458 theProjectileRemnant = NULL;
459 }
460
461 /** \brief Finalise the projectile remnant
462 *
463 * Complete the treatment of the projectile remnant. If it contains
464 * nucleons, assign its excitation energy and spin. Move stuff to the
465 * outgoing list, if appropriate.
466 *
467 * \param emissionTime the emission time of the projectile remnant
468 */
469 void finalizeProjectileRemnant(const G4double emissionTime);
470
471 /// \brief Update the particle potential energy.
472 inline void updatePotentialEnergy(Particle *p) const {
473 p->setPotentialEnergy(thePotential->computePotentialEnergy(p));
474 }
475
476 /// \brief Setter for theDensity
477 void setDensity(NuclearDensity const * const d) {
478 theDensity=d;
479 if(theParticleSampler)
480 theParticleSampler->setDensity(theDensity);
481 };
482
483 /// \brief Getter for theDensity
484 NuclearDensity const *getDensity() const { return theDensity; };
485
486 /// \brief Getter for thePotential
487 NuclearPotential::INuclearPotential const *getPotential() const { return thePotential; };
488
489 /// \brief Getter for theAnnihilationType
490 AnnihilationType getAnnihilationType() const { return theAType; }; //D
491
492 /// \brief Setter for theAnnihilationType
493 void setAnnihilationType(const AnnihilationType at){
494 theAType = at;
495 }; //D
496
497 private:
498 /** \brief Compute the recoil kinematics for a 1-nucleon remnant.
499 *
500 * Puts the remnant nucleon on mass shell and tries to enforce approximate
501 * energy conservation by modifying the masses of the outgoing particles.
502 */
503 void computeOneNucleonRecoilKinematics();
504
505 private:
506
507 G4int theInitialZ, theInitialA, theInitialS;
508 /// \brief The number of entering protons
509 G4int theNpInitial;
510 /// \brief The number of entering neutrons
511 G4int theNnInitial;
512 /// \brief The number of entering pions
513 G4int theNpionplusInitial;
514 G4int theNpionminusInitial;
515 /// \brief The number of entering kaons
516 G4int theNkaonplusInitial;
517 G4int theNkaonminusInitial;
518 /// \brief The number of entering antiprotons
519 G4int theNantiprotonInitial;
520
521 G4double initialInternalEnergy;
522 ThreeVector incomingAngularMomentum, incomingMomentum;
523 ThreeVector initialCenterOfMass;
524 G4bool remnant;
525
526 G4double initialEnergy;
527 Store *theStore;
528 G4bool tryCN;
529
530 /// \brief The radius of the universe
531 G4double theUniverseRadius;
532
533 /** \brief true if running a nucleus-nucleus collision
534 *
535 * Tells INCL whether to make a projectile-like pre-fragment or not.
536 */
537 G4bool isNucleusNucleus;
538
539 /** \brief Pointer to the quasi-projectile
540 *
541 * Owned by the Nucleus object.
542 */
543 ProjectileRemnant *theProjectileRemnant;
544
545 /// \brief Pointer to the NuclearDensity object
546 NuclearDensity const *theDensity;
547
548 /// \brief Pointer to the NuclearPotential object
549 NuclearPotential::INuclearPotential const *thePotential;
550
551 AnnihilationType theAType; //D same order as in the cc
552
553 INCL_DECLARE_ALLOCATION_POOL(Nucleus)
554 };
555
556 }
557
558 #endif /* G4INCLNUCLEUS_HH_ */
559