Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/cascade/cascade/src/G4PreCompoundInuclCollider.cc

    r962 r1315  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
     25// $Id: G4PreCompoundInuclCollider.cc,v 1.9 2010/05/21 17:56:34 mkelsey Exp $
     26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
    2527//
     28// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
     29// 20100309  M. Kelsey -- Eliminate some unnecessary std::pow()
     30// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
     31// 20100429  M. Kelsey -- Change "photon()" to "isPhoton()"
     32// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
     33//              simple data members, consolidate code
     34
    2635#include "G4PreCompoundInuclCollider.hh"
     36#include "G4BigBanger.hh"
     37#include "G4CollisionOutput.hh"
     38#include "G4ElementaryParticleCollider.hh"
     39#include "G4IntraNucleiCascader.hh"
    2740#include "G4InuclElementaryParticle.hh"
    2841#include "G4LorentzConvertor.hh"
    29 #include "G4ParticleLargerEkin.hh"
    30 #include <algorithm>
    31 
    32 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
    33 typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
     42#include "G4NonEquilibriumEvaporator.hh"
     43
    3444         
    3545G4PreCompoundInuclCollider::G4PreCompoundInuclCollider()
    36   : verboseLevel(0) {
    37 
    38   if (verboseLevel > 3) {
    39     G4cout << " >>> G4PreCompoundInuclCollider::G4PreCompoundInuclCollider" << G4endl;
    40   }
     46  : G4VCascadeCollider("G4PreCompoundInuclCollider"),
     47    theElementaryParticleCollider(new G4ElementaryParticleCollider),
     48    theIntraNucleiCascader(new G4IntraNucleiCascader),
     49    theNonEquilibriumEvaporator(new G4NonEquilibriumEvaporator),
     50    theBigBanger(new G4BigBanger) {}
     51
     52G4PreCompoundInuclCollider::~G4PreCompoundInuclCollider() {
     53  delete theElementaryParticleCollider;
     54  delete theIntraNucleiCascader;
     55  delete theNonEquilibriumEvaporator;
     56  delete theBigBanger;
    4157}
    4258
    43 G4CollisionOutput G4PreCompoundInuclCollider::collide(G4InuclParticle* bullet,
    44                                            G4InuclParticle* target) {
    45 
    46   verboseLevel = 0;
     59
     60void G4PreCompoundInuclCollider::collide(G4InuclParticle* bullet,
     61                                         G4InuclParticle* target,
     62                                         G4CollisionOutput& globalOutput) {
    4763  if (verboseLevel > 3) {
    4864    G4cout << " >>> G4PreCompoundInuclCollider::collide" << G4endl;
    … …  
    5167  const G4int itry_max = 1000;
    5268                     
    53   G4CollisionOutput globalOutput;
    54   G4InuclElementaryParticle* particle1 =
    55     dynamic_cast<G4InuclElementaryParticle*>(bullet);
    56   G4InuclElementaryParticle* particle2 =
    57     dynamic_cast<G4InuclElementaryParticle*>(target);
    58  
    59   if (particle1 && particle2) { // particle + particle (NOTE: also the h + H(1,1) treated here)
     69  if (useEPCollider(bullet,target)) {
    6070    if (verboseLevel > 2) {
    61       particle1->printParticle();
    62       particle2->printParticle();
     71      bullet->printParticle();
     72      target->printParticle();
    6373    }
    6474
    65     globalOutput = theElementaryParticleCollider->collide(bullet, target);
    66 
     75    theElementaryParticleCollider->collide(bullet, target, globalOutput);
    6776  } else { // needs to call all machinery       
    6877    G4LorentzConvertor convertToTargetRestFrame;
    69     G4InteractionCase interCase = bulletTargetSetter(bullet, target);
    70     G4int intcase = interCase.getInterCase();
     78
     79    interCase.set(bullet, target);
    7180     
    72     if (intcase > 0) { // ok
     81    if (interCase.valid()) { // ok
    7382      G4InuclNuclei* ntarget =
    7483        dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
    7584
    76       convertToTargetRestFrame.setTarget(ntarget->getMomentum(),
    77                                          ntarget->getMass());
     85      convertToTargetRestFrame.setTarget(ntarget);
    7886      G4int btype = 0;
    7987      G4double ab = 0.0;
    … …  
    8290      G4double zt = ntarget->getZ();
    8391       
    84       if (intcase == 1) { // particle with nuclei
     92      if (interCase.hadNucleus()) { // particle with nuclei
    8593        G4InuclElementaryParticle* pbullet =
    8694          dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
    8795         
    88         if (pbullet->photon()) {
     96        if (pbullet->isPhoton()) {
    8997          G4cout << " InuclCollider -> can not collide with photon " << G4endl;
    9098
    9199          globalOutput.trivialise(bullet, target);
    92 
    93           return globalOutput;
     100          return;
    94101        } else {
    95           convertToTargetRestFrame.setBullet(pbullet->getMomentum(),
    96                                              pbullet->getMass());   
     102          convertToTargetRestFrame.setBullet(pbullet);   
    97103          btype = pbullet->type();
    98104        };
    … …  
    102108          dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
    103109
    104         convertToTargetRestFrame.setBullet(nbullet->getMomentum(),
    105                                            nbullet->getMass());   
     110        convertToTargetRestFrame.setBullet(nbullet);   
    106111        ab = nbullet->getA();
    107112        zb = nbullet->getZ();
    … …  
    121126        }
    122127
    123         G4CascadeMomentum bmom;
    124 
    125         bmom[3] = convertToTargetRestFrame.getTRSMomentum();
    126 
    127         G4InuclNuclei ntarget(at, zt);
    128         G4CascadeMomentum tmom;
    129 
    130         ntarget.setMomentum(tmom);
    131         ntarget.setEnergy();
    132         theIntraNucleiCascader->setInteractionCase(intcase);
     128        G4LorentzVector bmom;
     129        bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
     130
     131        G4InuclNuclei ntarget(at, zt);          // Default is at rest
     132
     133        theIntraNucleiCascader->setInteractionCase(interCase.code());
    133134         
    134135        G4bool bad = true;
    135136        G4int itry = 0;
    136137         
     138        G4CollisionOutput TRFoutput;
     139        G4CollisionOutput output;
    137140        while (bad && itry < itry_max) {
    138           G4CollisionOutput TRFoutput;
    139           G4CollisionOutput output;
    140 
    141141          itry++;
    142           if (intcase == 1) {
     142
     143          output.reset();       // Clear buffers for this attempt
     144          TRFoutput.reset();
     145
     146          if (interCase.hadNucleus()) {
    143147            G4InuclElementaryParticle pbullet(bmom, btype);
    144148
    145             output = theIntraNucleiCascader->collide(&pbullet, &ntarget);
     149            theIntraNucleiCascader->collide(&pbullet, &ntarget, output);
    146150          } else {
    147             G4InuclNuclei nbullet(ab, zb);
    148             nbullet.setMomentum(bmom);
    149             nbullet.setEnergy();
    150             output = theIntraNucleiCascader->collide(&nbullet, &ntarget);
     151            G4InuclNuclei nbullet(bmom, ab, zb);
     152            theIntraNucleiCascader->collide(&nbullet, &ntarget, output);
    151153          };   
    152154
    153155          if (verboseLevel > 3) {
    154156            G4cout << " After Cascade " << G4endl;
    155 
    156157            output.printCollisionOutput();
    157158          }
    … …  
    160161          TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
    161162
    162           if (output.numberOfNucleiFragments() == 1) { // there is smth. after     
     163          if (output.numberOfNucleiFragments() == 1) { // there is smth. after
    163164            G4InuclNuclei cascad_rec_nuclei = output.getNucleiFragments()[0];
    164165            if (explosion(&cascad_rec_nuclei)) {
    … …  
    167168              };
    168169
    169               output = theBigBanger->collide(0,&cascad_rec_nuclei);
    170               TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
     170              theBigBanger->collide(0,&cascad_rec_nuclei, TRFoutput);
    171171            } else {
    172               output = theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei);
     172              output.reset();
     173              theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei, output);
    173174
    174175              if (verboseLevel > 3) {
    … …  
    178179
    179180              TRFoutput.addOutgoingParticles(output.getOutgoingParticles()); 
    180               TRFoutput.addTargetFragments(output.getNucleiFragments());         
     181              TRFoutput.addTargetFragments(output.getNucleiFragments());
    181182            };
    182183          };
    183184         
    184185          // convert to the LAB       
    185           G4bool withReflection = convertToTargetRestFrame.reflectionNeeded();       
    186           std::vector<G4InuclElementaryParticle> particles =
    187             TRFoutput.getOutgoingParticles();
    188 
    189           if (!particles.empty()) {
    190             particleIterator ipart;
    191             for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
    192               G4CascadeMomentum mom = ipart->getMomentum();
    193 
    194               if (withReflection) mom[3] = -mom[3];
    195               mom = convertToTargetRestFrame.rotate(mom);
    196               ipart->setMomentum(mom);
    197               mom = convertToTargetRestFrame.backToTheLab(ipart->getMomentum());
    198               ipart->setMomentum(mom);
    199             };
    200             std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
    201           };
    202            
    203           std::vector<G4InuclNuclei> nucleus = TRFoutput.getNucleiFragments();
    204 
    205           if (!nucleus.empty()) {
    206             nucleiIterator inuc;
    207 
    208             for (inuc = nucleus.begin(); inuc != nucleus.end(); inuc++) {
    209               G4CascadeMomentum mom = inuc->getMomentum();
    210 
    211               if (withReflection) mom[3] = -mom[3];
    212               mom = convertToTargetRestFrame.rotate(mom);
    213               inuc->setMomentum(mom);
    214               inuc->setEnergy();
    215               mom = convertToTargetRestFrame.backToTheLab(inuc->getMomentum());
    216               inuc->setMomentum(mom);
    217               inuc->setEnergy();
    218             };
    219           };
    220           globalOutput.addOutgoingParticles(particles);
    221           globalOutput.addTargetFragments(nucleus);
     186          TRFoutput.boostToLabFrame(convertToTargetRestFrame);
     187
     188          globalOutput.addOutgoingParticles(TRFoutput.getOutgoingParticles());
     189          globalOutput.addTargetFragments(TRFoutput.getNucleiFragments());
    222190          globalOutput.setOnShell(bullet, target);
    223           if(globalOutput.acceptable()) {
    224 
    225             return globalOutput;
    226           } else {
    227             globalOutput.reset();
    228           };
     191          if (globalOutput.acceptable()) return;
     192
     193          globalOutput.reset();         // Clear and try again
    229194        };
    230195
    … …  
    233198                 << itry_max << " attempts " << G4endl;
    234199        }
    235 
    236         globalOutput.trivialise(bullet, target);
    237 
    238         return globalOutput;       
    239200      } else {
    240 
    241201        if (verboseLevel > 3) {
    242202          G4cout << " InuclCollider -> inelastic interaction is impossible " << G4endl
    243203                 << " due to the coulomb barirer " << G4endl;
    244204        }
    245 
    246         globalOutput.trivialise(bullet, target);
    247 
    248         return globalOutput;
    249       };
     205      }
    250206       
     207      globalOutput.trivialise(bullet, target);
     208      return;
    251209    } else {
    252 
    253210      if (verboseLevel > 3) {
    254         G4cout << " InuclCollider -> inter case " << intcase << G4endl;
     211        G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
    255212      };
    256213    };       
    257214  };
    258215
    259   return globalOutput;
     216  return;
    260217}
    261                      
    262 G4bool G4PreCompoundInuclCollider::inelasticInteractionPossible(G4InuclParticle* bullet,
    263                                                      G4InuclParticle* target,
    264                                                      G4double ekin) const {
    265 
    266   if (verboseLevel > 3) {
    267     G4cout << " >>> G4PreCompoundInuclCollider::inelasticInteractionPossible" << G4endl;
    268   }
    269 
    270   const G4double coeff = 0.001 * 1.2;
    271   const G4double one_third = 1.0 / 3.0;
    272 
    273   G4bool possible = true;
    274   G4double at;
    275   G4double zt;
    276   G4double ab;
    277   G4double zb;
    278 
    279   if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
    280     at = nuclei_target->getA();
    281     zt = nuclei_target->getZ();
    282     if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
    283       ab = nuclei_bullet->getA();
    284       zb = nuclei_bullet->getZ();     
    285     } else {
    286       G4InuclElementaryParticle* particle =
    287         dynamic_cast<G4InuclElementaryParticle*>(bullet);
    288 
    289       ab = 1;
    290       zb = particle->getCharge();
    291     };
    292   } else {
    293     if(G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
    294       ab = nuclei_bullet->getA();
    295       zb = nuclei_bullet->getZ();     
    296 
    297       G4InuclElementaryParticle* particle =
    298         dynamic_cast<G4InuclElementaryParticle*>(target);
    299 
    300       at = 1;
    301       zt = particle->getCharge();   
    302     } else {
    303 
    304       return possible;
    305     }; 
    306   };
    307 
    308   // VCOL used  for testing if elastic collision possible
    309   G4double VCOL = coeff * zt * zb / (std::pow(at, one_third) + std::pow(ab, one_third));
    310 
    311   // possible = VCOL < ekin; // NOTE: inelastic collision if not true
    312   possible = true; // we force elastic
    313 
    314   if (verboseLevel > 3) {
    315     G4cout << " >>> G4PreCompoundInuclCollider::inelasticInteractionPossible" << G4endl;
    316     G4cout << " VCOL: " << VCOL << " ekin: " << ekin << " inelastic possible: " << possible << G4endl;
    317   }
    318 
    319   return possible;
    320 
    321 }
    322        
    323 G4InteractionCase G4PreCompoundInuclCollider::bulletTargetSetter(G4InuclParticle* bullet,
    324                                                       G4InuclParticle* target) const {
    325 
    326   if (verboseLevel > 3) {
    327     G4cout << " >>> G4PreCompoundInuclCollider::bulletTargetSetter" << G4endl;
    328   }
    329 
    330   G4InteractionCase interCase;
    331 
    332   if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {     
    333     if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) { // A + A         
    334       interCase.setInterCase(2);
    335       if (nuclei_target->getA() >= nuclei_bullet->getA()) {
    336         interCase.setBulletTarget(bullet, target);
    337       } else {
    338         interCase.setBulletTarget(target, bullet);
    339       };
    340     } else {
    341       interCase.setInterCase(1);
    342       interCase.setBulletTarget(bullet, target);
    343     };
    344   } else {
    345     G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet);
    346     if (nuclei_bullet) {
    347       G4InuclElementaryParticle* part =
    348         dynamic_cast<G4InuclElementaryParticle*>(target);
    349       if (part) {
    350         interCase.setInterCase(1);
    351         interCase.setBulletTarget(target, bullet);
    352       };
    353     };
    354   };
    355 
    356   return interCase;
    357 }       
    358 
    359 G4bool G4PreCompoundInuclCollider::explosion(G4InuclNuclei* target) const {
    360 
    361   if (verboseLevel > 3) {
    362     G4cout << " >>> G4PreCompoundInuclCollider::explosion" << G4endl;
    363   }
    364 
    365   const G4double a_cut = 20.0;
    366   const G4double be_cut = 3.0;
    367 
    368   G4double a = target->getA();
    369   G4double z = target->getZ();
    370   G4double eexs = target->getExitationEnergy();
    371   G4bool explo = true;
    372 
    373   if (a > a_cut) {
    374     explo = false;
    375   } else {
    376     if (eexs < be_cut * bindingEnergy(a, z)) explo = false;
    377   };   
    378 
    379   return explo;
    380 }
    381  
Note: See TracChangeset for help on using the changeset viewer.