| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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 /* \file G4INCLInteractionAvatar.cc
39 * \brief Virtual class for interaction avatar
40 *
41 * This class is inherited by decay and collis
42 * provide a uniform treatment of common physi
43 * enforcement of energy conservation, etc.
44 *
45 * \date Mar 1st, 2011
46 * \author Davide Mancusi
47 */
48
49 #include "G4INCLInteractionAvatar.hh"
50 #include "G4INCLKinematicsUtils.hh"
51 #include "G4INCLCrossSections.hh"
52 #include "G4INCLPauliBlocking.hh"
53 #include "G4INCLRootFinder.hh"
54 #include "G4INCLLogger.hh"
55 #include "G4INCLConfigEnums.hh"
56 #include "G4INCLConfig.hh"
57 // #include <cassert>
58
59 namespace G4INCL {
60
61 const G4double InteractionAvatar::locEAccura
62 const G4int InteractionAvatar::maxIterLocE =
63 G4ThreadLocal Particle *InteractionAvatar::b
64 G4ThreadLocal Particle *InteractionAvatar::b
65
66 InteractionAvatar::InteractionAvatar(G4doubl
67 : IAvatar(time), theNucleus(n),
68 particle1(p1), particle2(NULL),
69 isPiN(false),
70 weight(1.),
71 violationEFunctor(NULL)
72 {
73 }
74
75 InteractionAvatar::InteractionAvatar(G4doubl
76 G4INCL::Particle *p2)
77 : IAvatar(time), theNucleus(n),
78 particle1(p1), particle2(p2),
79 isPiN((p1->isPion() && p2->isNucleon()) ||
80 weight(1.),
81 violationEFunctor(NULL)
82 {
83 }
84
85 InteractionAvatar::~InteractionAvatar() {
86 }
87
88 void InteractionAvatar::deleteBackupParticle
89 delete backupParticle1;
90 if(backupParticle2)
91 delete backupParticle2;
92 backupParticle1 = NULL;
93 backupParticle2 = NULL;
94 }
95
96 void InteractionAvatar::preInteractionBlocki
97 if(backupParticle1)
98 (*backupParticle1) = (*particle1);
99 else
100 backupParticle1 = new Particle(*particle
101
102 if(particle2) {
103 if(backupParticle2)
104 (*backupParticle2) = (*particle2);
105 else
106 backupParticle2 = new Particle(*partic
107
108 oldTotalEnergy = particle1->getEnergy()
109 - particle1->getPotentialEnergy() - pa
110 oldXSec = CrossSections::total(particle1
111 } else {
112 oldTotalEnergy = particle1->getEnergy()
113 }
114 }
115
116 void InteractionAvatar::preInteractionLocalE
117 if(!theNucleus || p->isMeson() || p->isPho
118
119 if(shouldUseLocalEnergy())
120 KinematicsUtils::transformToLocalEnergyF
121 }
122
123 void InteractionAvatar::preInteraction() {
124 preInteractionBlocking();
125
126 preInteractionLocalEnergy(particle1);
127
128 if(particle2) {
129 preInteractionLocalEnergy(particle2);
130 boostVector = KinematicsUtils::makeBoost
131 particle2->boost(boostVector);
132 } else {
133 boostVector = particle1->getMomentum()/p
134 }
135 particle1->boost(boostVector);
136 }
137
138 G4bool InteractionAvatar::bringParticleInsid
139 if(!theNucleus)
140 return false;
141
142 ThreeVector pos = p->getPosition();
143 p->rpCorrelate();
144 G4double pos2 = pos.mag2();
145 const G4double r = theNucleus->getSurfaceR
146 short iterations=0;
147 const short maxIterations=50;
148
149 if(pos2 < r*r) return true;
150
151 while( pos2 >= r*r && iterations<maxIterat
152 {
153 pos *= std::sqrt(r*r*0.9801/pos2); // 0.
154 pos2 = pos.mag2();
155 iterations++;
156 }
157 if( iterations < maxIterations)
158 {
159 INCL_DEBUG("Particle position vector len
160 p->setPosition(pos);
161 return true;
162 }
163 else
164 return false;
165 }
166
167 void InteractionAvatar::postInteraction(Fina
168 INCL_DEBUG("postInteraction: final state:
169 modified = fs->getModifiedParticles();
170 created = fs->getCreatedParticles();
171 Destroyed = fs->getDestroyedParticles();
172 modifiedAndCreated = modified;
173 modifiedAndCreated.insert(modifiedAndCreat
174 ModifiedAndDestroyed = modified;
175 ModifiedAndDestroyed.insert(ModifiedAndDes
176
177 // Boost back to lab
178 modifiedAndCreated.boost(-boostVector);
179
180 // If there is no Nucleus, just return
181 if(!theNucleus) return;
182
183 // Mark pions and kaons that have been cre
184 // to be emitted later).
185 for(ParticleIter i=created.begin(), e=crea
186 if(((*i)->isPion() || (*i)->isKaon() ||
187 (*i)->makeParticipant();
188 (*i)->setOutOfWell();
189 fs->addOutgoingParticle(*i);
190 INCL_DEBUG("Pion was created outside i
191 << (*i)->print());
192 }
193
194 // Try to enforce energy conservation
195 fs->setTotalEnergyBeforeInteraction(oldTot
196 G4bool success = enforceEnergyConservation
197 if(!success) {
198 INCL_DEBUG("Enforcing energy conservatio
199
200 // Restore the state of the initial part
201 restoreParticles();
202
203 // Delete newly created particles
204 for(ParticleIter i=created.begin(), e=cr
205 delete *i;
206
207 fs->reset();
208 fs->makeNoEnergyConservation();
209 fs->setTotalEnergyBeforeInteraction(0.0)
210
211 return; // Interaction is blocked. Retur
212 }
213 INCL_DEBUG("Enforcing energy conservation:
214
215 INCL_DEBUG("postInteraction after energy c
216
217 // Check that outgoing delta resonances ca
218 for(ParticleIter i=modified.begin(), e=mod
219 if((*i)->isDelta() &&
220 (*i)->getMass() < ParticleTable::min
221 INCL_DEBUG("Mass of the produced delta
222 (*i)->getMass() << '\n');
223
224 // Restore the state of the initial pa
225 restoreParticles();
226
227 // Delete newly created particles
228 for(ParticleIter j=created.begin(), en
229 delete *j;
230
231 fs->reset();
232 fs->makeNoEnergyConservation();
233 fs->setTotalEnergyBeforeInteraction(0.
234
235 return; // Interaction is blocked. Ret
236 }
237
238 INCL_DEBUG("Random seeds before Pauli bloc
239 // Test Pauli blocking
240 G4bool isBlocked = Pauli::isBlocked(modifi
241
242 if(isBlocked) {
243 INCL_DEBUG("Pauli: Blocked!" << '\n');
244
245 // Restore the state of the initial part
246 restoreParticles();
247
248 // Delete newly created particles
249 for(ParticleIter i=created.begin(), e=cr
250 delete *i;
251
252 fs->reset();
253 fs->makePauliBlocked();
254 fs->setTotalEnergyBeforeInteraction(0.0)
255
256 return; // Interaction is blocked. Retur
257 }
258 INCL_DEBUG("Pauli: Allowed!" << '\n');
259
260 // Test CDPP blocking
261 G4bool isCDPPBlocked = Pauli::isCDPPBlocke
262
263 if(isCDPPBlocked) {
264 INCL_DEBUG("CDPP: Blocked!" << '\n');
265
266 // Restore the state of the initial part
267 restoreParticles();
268
269 // Delete newly created particles
270 for(ParticleIter i=created.begin(), e=cr
271 delete *i;
272
273 fs->reset();
274 fs->makePauliBlocked();
275 fs->setTotalEnergyBeforeInteraction(0.0)
276
277 return; // Interaction is blocked. Retur
278 }
279 INCL_DEBUG("CDPP: Allowed!" << '\n');
280
281 // If all went well, try to bring particle
282 for(ParticleIter i=modifiedAndCreated.begi
283 {
284 // ...except for pions beyond their surf
285 if((*i)->isOutOfWell()) continue;
286
287 const G4bool successBringParticlesInside
288 if( !successBringParticlesInside ) {
289 INCL_ERROR("Failed to bring particle i
290 }
291 }
292
293 // Collision accepted!
294 // Biasing of the final state
295 std::vector<G4int> newBiasCollisionVector;
296 newBiasCollisionVector = ModifiedAndDestro
297 if(std::fabs(weight-1.) > 1E-6){
298 newBiasCollisionVector.push_back(Particl
299 Particle::FillINCLBiasVector(1./weight);
300 weight = 1.; // useless?
301 }
302 for(ParticleIter i=modifiedAndCreated.begi
303 (*i)->setBiasCollisionVector(newBiasColl
304 if(!(*i)->isOutOfWell()) {
305 // Decide if the particle should be ma
306 // (Back to spectator)
307 G4bool goesBackToSpectator = false;
308 if((*i)->isNucleon() && theNucleus->ge
309 G4double threshold = (*i)->getPotent
310 if((*i)->getType()==Proton)
311 threshold += Math::twoThirds*theNu
312 if((*i)->getKineticEnergy() < thresh
313 goesBackToSpectator = true;
314 }
315 // Thaw the particle propagation
316 (*i)->thawPropagation();
317
318 // Increment or decrement the particip
319 if(goesBackToSpectator) {
320 INCL_DEBUG("The following particle g
321 << (*i)->print() << '\n');
322 if(!(*i)->isTargetSpectator()) {
323 theNucleus->getStore()->getBook().
324 }
325 (*i)->makeTargetSpectator();
326 } else {
327 if((*i)->isTargetSpectator()) {
328 theNucleus->getStore()->getBook().
329 }
330 (*i)->makeParticipant();
331 }
332 }
333 }
334 ParticleList destroyed = fs->getDestroyedP
335 for(ParticleIter i=destroyed.begin(), e=de
336 if(!(*i)->isTargetSpectator())
337 theNucleus->getStore()->getBook().decr
338 return;
339 }
340
341 void InteractionAvatar::restoreParticles() c
342 (*particle1) = (*backupParticle1);
343 if(particle2)
344 (*particle2) = (*backupParticle2);
345 }
346
347 G4bool InteractionAvatar::shouldUseLocalEner
348 if(!theNucleus) return false;
349 LocalEnergyType theLocalEnergyType;
350 if(theNucleus->getStore()->getConfig()->ge
351 theNucleus->getStore()->getConfig()->getP
352 return false;
353 }
354 if(getType()==DecayAvatarType || isPiN)
355 theLocalEnergyType = theNucleus->getStor
356 else
357 theLocalEnergyType = theNucleus->getStor
358
359 const G4bool firstAvatar = (theNucleus->ge
360 return ((theLocalEnergyType == FirstCollis
361 theLocalEnergyType == AlwaysLocalE
362 }
363
364 G4bool InteractionAvatar::enforceEnergyConse
365 // Set up the violationE calculation
366 const G4bool manyBodyFinalState = (modifie
367
368 if(manyBodyFinalState)
369 violationEFunctor = new ViolationEMoment
370 else {
371 if (modified.empty()) {
372 Particle * const p1 = created.front();
373 // The following condition is necessa
374 // correctly. A similar condition exi
375 if(p1->getMass() < ParticleTable::minD
376 return false;
377 violationEFunctor = new ViolationEEner
378 }
379 else{
380 Particle * const p2 = modified.front()
381 // The following condition is necessa
382 // correctly. A similar condition exi
383 if(p2->getMass() < ParticleTable::min
384 return false;
385 violationEFunctor = new ViolationEEne
386 }
387 }
388
389 // Apply the root-finding algorithm
390 const RootFinder::Solution theSolution = R
391 if(theSolution.success) { // Apply the sol
392 (*violationEFunctor)(theSolution.x);
393 } else if(theNucleus){
394 INCL_DEBUG("Couldn't enforce energy cons
395 theNucleus->getStore()->getBook().increm
396 }
397 delete violationEFunctor;
398 violationEFunctor = NULL;
399 return theSolution.success;
400 }
401
402 /* ***
403 * *** InteractionAvatar::ViolationEMomentum
404 * ***
405
406 InteractionAvatar::ViolationEMomentumFunctor
407 RootFunctor(0., 1E6),
408 finalParticles(modAndCre),
409 initialEnergy(totalEnergyBeforeInteraction
410 theNucleus(nucleus),
411 boostVector(boost),
412 shouldUseLocalEnergy(localE)
413 {
414 // Store the particle momenta (necessary f
415 // scaleParticleMomenta() to work)
416 for(ParticleIter i=finalParticles.begin(),
417 (*i)->boost(boostVector);
418 particleMomenta.push_back((*i)->getMomen
419 }
420 }
421
422 InteractionAvatar::ViolationEMomentumFunctor
423 particleMomenta.clear();
424 }
425
426 G4double InteractionAvatar::ViolationEMoment
427 scaleParticleMomenta(alpha);
428
429 G4double deltaE = 0.0;
430 for(ParticleIter i=finalParticles.begin(),
431 deltaE += (*i)->getEnergy() - (*i)->getP
432 deltaE -= initialEnergy;
433 return deltaE;
434 }
435
436 void InteractionAvatar::ViolationEMomentumFu
437
438 std::vector<ThreeVector>::const_iterator i
439 for(ParticleIter i=finalParticles.begin(),
440 (*i)->setMomentum((*iP)*alpha);
441 (*i)->adjustEnergyFromMomentum();
442 (*i)->rpCorrelate();
443 (*i)->boost(-boostVector);
444 if(theNucleus){
445 theNucleus->updatePotentialEnergy(*i);
446 } else {
447 (*i)->setPotentialEnergy(0.);
448 }
449
450 if(shouldUseLocalEnergy && !(*i)->isPion
451 !(*i)->isKaon() && !(*i)->isAntiKaon(
452 // assert(theNucleus); // Local energy without
453 const G4double energy = (*i)->getEner
454 G4double locE = KinematicsUtils::getL
455 G4double locEOld;
456 G4double deltaLocE = InteractionAvata
457 for(G4int iterLocE=0;
458 deltaLocE>InteractionAvatar::locEA
459 ++iterLocE) {
460 locEOld = locE;
461 (*i)->setEnergy(energy + locE); // U
462 (*i)->adjustMomentumFromEnergy();
463 theNucleus->updatePotentialEnergy(*i
464 locE = KinematicsUtils::getLocalEner
465 deltaLocE = std::abs(locE-locEOld);
466 }
467 }
468
469 //jlrs For lambdas and nuclei with masses hig
470 if(shouldUseLocalEnergy && (*i)->isLam
471 // assert(theNucleus); // Local energy without
472 const G4double energy = (*i)->getEner
473 G4double locE = KinematicsUtils::getL
474 G4double locEOld;
475 G4double deltaLocE = InteractionAvata
476 for(G4int iterLocE=0;
477 deltaLocE>InteractionAvatar::locEA
478 ++iterLocE) {
479 locEOld = locE;
480 (*i)->setEnergy(energy + locE); // U
481 (*i)->adjustMomentumFromEnergy();
482 theNucleus->updatePotentialEnergy(*i
483 locE = KinematicsUtils::getLocalEner
484 deltaLocE = std::abs(locE-locEOld);
485 }
486 }
487 }
488 }
489
490 void InteractionAvatar::ViolationEMomentumFu
491 if(!success)
492 scaleParticleMomenta(1.);
493 }
494
495 /* ***
496 * *** InteractionAvatar::ViolationEEnergyFu
497 * ***
498
499 InteractionAvatar::ViolationEEnergyFunctor::
500 RootFunctor(0., 1E6),
501 initialEnergy(totalEnergyBeforeInteraction
502 theNucleus(nucleus),
503 theParticle(aParticle),
504 theEnergy(theParticle->getEnergy()),
505 theMomentum(theParticle->getMomentum()),
506 energyThreshold(KinematicsUtils::energy(th
507 shouldUseLocalEnergy(localE)
508 {
509 // assert(theParticle->isDelta());
510 }
511
512 G4double InteractionAvatar::ViolationEEnergy
513 setParticleEnergy(alpha);
514 return theParticle->getEnergy() - theParti
515 }
516
517 void InteractionAvatar::ViolationEEnergyFunc
518
519 G4double locE;
520 if(shouldUseLocalEnergy) {
521 // assert(theNucleus); // Local energy without
522 locE = KinematicsUtils::getLocalEnergy(t
523 } else
524 locE = 0.;
525 G4double locEOld;
526 G4double deltaLocE = InteractionAvatar::lo
527 for(G4int iterLocE=0;
528 deltaLocE>InteractionAvatar::locEAccur
529 ++iterLocE) {
530 locEOld = locE;
531 G4double particleEnergy = energyThreshol
532 const G4double theMass2 = std::pow(parti
533 G4double theMass;
534 if(theMass2>ParticleTable::minDeltaMass2
535 theMass = std::sqrt(theMass2);
536 else {
537 theMass = ParticleTable::minDeltaMass;
538 particleEnergy = energyThreshold;
539 }
540 theParticle->setMass(theMass);
541 theParticle->setEnergy(particleEnergy);
542 if(theNucleus) {
543 theNucleus->updatePotentialEnergy(theP
544 if(shouldUseLocalEnergy)
545 locE = KinematicsUtils::getLocalEner
546 else
547 locE = 0.;
548 } else
549 locE = 0.;
550 deltaLocE = std::abs(locE-locEOld);
551 }
552
553 }
554
555 void InteractionAvatar::ViolationEEnergyFunc
556 if(!success)
557 setParticleEnergy(1.);
558 }
559
560 }
561