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/G4InuclCollider.cc

    r962 r1315  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
     25// $Id: G4InuclCollider.cc,v 1.30 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// 20100418  M. Kelsey -- Move lab-frame transformation code to G4CollisonOutput
     32// 20100429  M. Kelsey -- Change "photon()" to "isPhoton()"
     33// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
     34//              simple data members, consolidate code
     35
    2636#include "G4InuclCollider.hh"
     37#include "G4BigBanger.hh"
     38#include "G4CollisionOutput.hh"
     39#include "G4ElementaryParticleCollider.hh"
     40#include "G4EquilibriumEvaporator.hh"
     41#include "G4IntraNucleiCascader.hh"
    2742#include "G4InuclElementaryParticle.hh"
    2843#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;
    34          
     44#include "G4NonEquilibriumEvaporator.hh"
     45
     46
    3547G4InuclCollider::G4InuclCollider()
    36   : verboseLevel(0) {
    37 
    38   if (verboseLevel > 3) {
    39     G4cout << " >>> G4InuclCollider::G4InuclCollider" << G4endl;
    40   }
     48  : G4VCascadeCollider("G4InuclCollider"),
     49    theElementaryParticleCollider(new G4ElementaryParticleCollider),
     50    theIntraNucleiCascader(new G4IntraNucleiCascader),
     51    theNonEquilibriumEvaporator(new G4NonEquilibriumEvaporator),
     52    theEquilibriumEvaporator(new G4EquilibriumEvaporator),
     53    theBigBanger(new G4BigBanger) {}
     54
     55G4InuclCollider::~G4InuclCollider() {
     56  delete theElementaryParticleCollider;
     57  delete theIntraNucleiCascader;
     58  delete theNonEquilibriumEvaporator;
     59  delete theEquilibriumEvaporator;
     60  delete theBigBanger;
    4161}
    4262
    43 G4CollisionOutput G4InuclCollider::collide(G4InuclParticle* bullet,
    44                                            G4InuclParticle* target) {
    45 
    46   verboseLevel = 0;
     63
     64void G4InuclCollider::collide(G4InuclParticle* bullet, G4InuclParticle* target,
     65                              G4CollisionOutput& globalOutput) {
    4766  if (verboseLevel > 3) {
    4867    G4cout << " >>> G4InuclCollider::collide" << G4endl;
    … …  
    5170  const G4int itry_max = 1000;
    5271                     
    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)
     72  if (useEPCollider(bullet,target)) {
    6073    if (verboseLevel > 2) {
    61       particle1->printParticle();
    62       particle2->printParticle();
     74      bullet->printParticle();
     75      target->printParticle();
    6376    }
    6477
    65     globalOutput = theElementaryParticleCollider->collide(bullet, target);
    66 
     78    theElementaryParticleCollider->collide(bullet, target, globalOutput);
    6779  } else { // needs to call all machinery       
    6880    G4LorentzConvertor convertToTargetRestFrame;
    69     G4InteractionCase interCase = bulletTargetSetter(bullet, target);
    70     G4int intcase = interCase.getInterCase();
    71      
    72     if (intcase > 0) { // ok
     81
     82    interCase.set(bullet,target);
     83    if (interCase.valid()) { // ok
    7384      G4InuclNuclei* ntarget =
    7485        dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
    7586
    76       convertToTargetRestFrame.setTarget(ntarget->getMomentum(),
    77                                          ntarget->getMass());
     87      convertToTargetRestFrame.setTarget(ntarget);
    7888      G4int btype = 0;
    7989      G4double ab = 0.0;
    … …  
    8292      G4double zt = ntarget->getZ();
    8393       
    84       if (intcase == 1) { // particle with nuclei
     94      if (interCase.hadNucleus()) { // particle with nuclei
    8595        G4InuclElementaryParticle* pbullet =
    8696          dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
    8797         
    88         if (pbullet->photon()) {
    89           G4cout << " InuclCollider -> can not collide with photon " << G4endl;
     98        if (pbullet->isPhoton()) {
     99          G4cerr << " InuclCollider -> can not collide with photon " << G4endl;
    90100
    91101          globalOutput.trivialise(bullet, target);
    92 
    93           return globalOutput;
     102          return;
    94103        } else {
    95           convertToTargetRestFrame.setBullet(pbullet->getMomentum(),
    96                                              pbullet->getMass());   
     104          convertToTargetRestFrame.setBullet(pbullet);   
    97105          btype = pbullet->type();
    98106        };
    … …  
    102110          dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
    103111
    104         convertToTargetRestFrame.setBullet(nbullet->getMomentum(),
    105                                            nbullet->getMass());   
     112        convertToTargetRestFrame.setBullet(nbullet);   
    106113        ab = nbullet->getA();
    107114        zb = nbullet->getZ();
    … …  
    121128        }
    122129
    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);
     130        G4LorentzVector bmom;
     131        bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
     132
     133        G4InuclNuclei ntarget(at, zt);          // Default is at rest
     134
     135        theIntraNucleiCascader->setInteractionCase(interCase.code());
    133136         
    134137        G4bool bad = true;
    135138        G4int itry = 0;
    136139         
     140        G4CollisionOutput TRFoutput;
     141        G4CollisionOutput output;
    137142        while (bad && itry < itry_max) {
    138           G4CollisionOutput TRFoutput;
    139           G4CollisionOutput output;
    140 
    141143          itry++;
    142           if (intcase == 1) {
     144
     145          output.reset();       // Clear buffers for this attempt
     146          TRFoutput.reset();
     147
     148          if (interCase.hadNucleus()) {
    143149            G4InuclElementaryParticle pbullet(bmom, btype);
    144150
    145             output = theIntraNucleiCascader->collide(&pbullet, &ntarget);
     151            theIntraNucleiCascader->collide(&pbullet, &ntarget, output);
    146152          } else {
    147             G4InuclNuclei nbullet(ab, zb);
    148             nbullet.setMomentum(bmom);
    149             nbullet.setEnergy();
    150             output = theIntraNucleiCascader->collide(&nbullet, &ntarget);
     153            G4InuclNuclei nbullet(bmom, ab, zb);
     154            theIntraNucleiCascader->collide(&nbullet, &ntarget, output);
    151155          };   
    152156
    153157          if (verboseLevel > 3) {
    154158            G4cout << " After Cascade " << G4endl;
    155 
    156159            output.printCollisionOutput();
    157160          }
    158161         
    159162          // the rest, if any
     163          // FIXME:  The code below still does too much copying!
    160164          TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
    161165
    162           if (output.numberOfNucleiFragments() == 1) { // there is smth. after     
     166          if (output.numberOfNucleiFragments() == 1) { // there is smth. after
    163167            G4InuclNuclei cascad_rec_nuclei = output.getNucleiFragments()[0];
    164168            if (explosion(&cascad_rec_nuclei)) {
    … …  
    167171              };
    168172
    169               output = theBigBanger->collide(0,&cascad_rec_nuclei);
    170               TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
     173              theBigBanger->collide(0,&cascad_rec_nuclei, TRFoutput);
    171174            } else {
    172               output = theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei);
     175              output.reset();
     176              theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei, output);
    173177
    174178              if (verboseLevel > 3) {
    … …  
    179183              TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
    180184              G4InuclNuclei exiton_rec_nuclei = output.getNucleiFragments()[0];
    181               output = theEquilibriumEvaporator->collide(0, &exiton_rec_nuclei);
     185
     186              output.reset();
     187              theEquilibriumEvaporator->collide(0, &exiton_rec_nuclei, output);
    182188
    183189              if (verboseLevel > 3) {
    184190                G4cout << " After EquilibriumEvaporator " << G4endl;
    185 
    186191                output.printCollisionOutput();
    187192              };
    188193
    189194              TRFoutput.addOutgoingParticles(output.getOutgoingParticles()); 
    190               TRFoutput.addTargetFragments(output.getNucleiFragments());         
     195              TRFoutput.addTargetFragments(output.getNucleiFragments());
    191196            };
    192197          };
    193198         
    194           // convert to the LAB       
    195           G4bool withReflection = convertToTargetRestFrame.reflectionNeeded();       
    196           std::vector<G4InuclElementaryParticle> particles =
    197             TRFoutput.getOutgoingParticles();
    198 
    199           if (!particles.empty()) {
    200             particleIterator ipart;
    201             for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
    202               G4CascadeMomentum mom = ipart->getMomentum();
    203 
    204               if (withReflection) mom[3] = -mom[3];
    205               mom = convertToTargetRestFrame.rotate(mom);
    206               ipart->setMomentum(mom);
    207               mom = convertToTargetRestFrame.backToTheLab(ipart->getMomentum());
    208               ipart->setMomentum(mom);
    209             };
    210             std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
    211           };
    212            
    213           std::vector<G4InuclNuclei> nucleus = TRFoutput.getNucleiFragments();
    214 
    215           if (!nucleus.empty()) {
    216             nucleiIterator inuc;
    217 
    218             for (inuc = nucleus.begin(); inuc != nucleus.end(); inuc++) {
    219               G4CascadeMomentum mom = inuc->getMomentum();
    220 
    221               if (withReflection) mom[3] = -mom[3];
    222               mom = convertToTargetRestFrame.rotate(mom);
    223               inuc->setMomentum(mom);
    224               inuc->setEnergy();
    225               mom = convertToTargetRestFrame.backToTheLab(inuc->getMomentum());
    226               inuc->setMomentum(mom);
    227               inuc->setEnergy();
    228             };
    229           };
    230           globalOutput.addOutgoingParticles(particles);
    231           globalOutput.addTargetFragments(nucleus);
     199          // convert to the LAB
     200          TRFoutput.boostToLabFrame(convertToTargetRestFrame);
     201
     202          globalOutput.addOutgoingParticles(TRFoutput.getOutgoingParticles());
     203          globalOutput.addTargetFragments(TRFoutput.getNucleiFragments());
    232204          globalOutput.setOnShell(bullet, target);
    233           if(globalOutput.acceptable()) {
    234 
    235             return globalOutput;
    236           } else {
    237             globalOutput.reset();
    238           };
     205          if (globalOutput.acceptable()) return;
     206
     207          globalOutput.reset();         // Clear and try again
    239208        };
    240209
    … …  
    243212                 << itry_max << " attempts " << G4endl;
    244213        }
    245 
    246         globalOutput.trivialise(bullet, target);
    247 
    248         return globalOutput;       
    249214      } else {
    250 
    251215        if (verboseLevel > 3) {
    252216          G4cout << " InuclCollider -> inelastic interaction is impossible " << G4endl
    253217                 << " due to the coulomb barirer " << G4endl;
    254218        }
    255 
    256         globalOutput.trivialise(bullet, target);
    257 
    258         return globalOutput;
    259       };
    260        
     219      }
     220
     221      globalOutput.trivialise(bullet, target);
     222      return;
    261223    } else {
    262 
    263224      if (verboseLevel > 3) {
    264         G4cout << " InuclCollider -> inter case " << intcase << G4endl;
     225        G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
    265226      };
    266227    };       
    267228  };
    268229
    269   return globalOutput;
     230  return;
    270231}
    271                      
    272 G4bool G4InuclCollider::inelasticInteractionPossible(G4InuclParticle* bullet,
    273                                                      G4InuclParticle* target,
    274                                                      G4double ekin) const {
    275 
    276   if (verboseLevel > 3) {
    277     G4cout << " >>> G4InuclCollider::inelasticInteractionPossible" << G4endl;
    278   }
    279 
    280   const G4double coeff = 0.001 * 1.2;
    281   const G4double one_third = 1.0 / 3.0;
    282 
    283   G4bool possible = true;
    284   G4double at;
    285   G4double zt;
    286   G4double ab;
    287   G4double zb;
    288 
    289   if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
    290     at = nuclei_target->getA();
    291     zt = nuclei_target->getZ();
    292     if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
    293       ab = nuclei_bullet->getA();
    294       zb = nuclei_bullet->getZ();     
    295     } else {
    296       G4InuclElementaryParticle* particle =
    297         dynamic_cast<G4InuclElementaryParticle*>(bullet);
    298 
    299       ab = 1;
    300       zb = particle->getCharge();
    301     };
    302   } else {
    303     if(G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
    304       ab = nuclei_bullet->getA();
    305       zb = nuclei_bullet->getZ();     
    306 
    307       G4InuclElementaryParticle* particle =
    308         dynamic_cast<G4InuclElementaryParticle*>(target);
    309 
    310       at = 1;
    311       zt = particle->getCharge();   
    312     } else {
    313 
    314       return possible;
    315     }; 
    316   };
    317 
    318   // VCOL used  for testing if elastic collision possible
    319   G4double VCOL = coeff * zt * zb / (std::pow(at, one_third) + std::pow(ab, one_third));
    320 
    321   // possible = VCOL < ekin; // NOTE: inelastic collision if not true
    322   possible = true; // we force elastic
    323 
    324   if (verboseLevel > 3) {
    325     G4cout << " >>> G4InuclCollider::inelasticInteractionPossible" << G4endl;
    326     G4cout << " VCOL: " << VCOL << " ekin: " << ekin << " inelastic possible: " << possible << G4endl;
    327   }
    328 
    329   return possible;
    330 
    331 }
    332        
    333 G4InteractionCase G4InuclCollider::bulletTargetSetter(G4InuclParticle* bullet,
    334                                                       G4InuclParticle* target) const {
    335 
    336   if (verboseLevel > 3) {
    337     G4cout << " >>> G4InuclCollider::bulletTargetSetter" << G4endl;
    338   }
    339 
    340   G4InteractionCase interCase;
    341 
    342   if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {     
    343     if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) { // A + A         
    344       interCase.setInterCase(2);
    345       if (nuclei_target->getA() >= nuclei_bullet->getA()) {
    346         interCase.setBulletTarget(bullet, target);
    347       } else {
    348         interCase.setBulletTarget(target, bullet);
    349       };
    350     } else {
    351       interCase.setInterCase(1);
    352       interCase.setBulletTarget(bullet, target);
    353     };
    354   } else {
    355     G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet);
    356     if (nuclei_bullet) {
    357       G4InuclElementaryParticle* part =
    358         dynamic_cast<G4InuclElementaryParticle*>(target);
    359       if (part) {
    360         interCase.setInterCase(1);
    361         interCase.setBulletTarget(target, bullet);
    362       };
    363     };
    364   };
    365 
    366   return interCase;
    367 }       
    368 
    369 G4bool G4InuclCollider::explosion(G4InuclNuclei* target) const {
    370 
    371   if (verboseLevel > 3) {
    372     G4cout << " >>> G4InuclCollider::explosion" << G4endl;
    373   }
    374 
    375   const G4double a_cut = 20.0;
    376   const G4double be_cut = 3.0;
    377 
    378   G4double a = target->getA();
    379   G4double z = target->getZ();
    380   G4double eexs = target->getExitationEnergy();
    381   G4bool explo = true;
    382 
    383   if (a > a_cut) {
    384     explo = false;
    385   } else {
    386     if (eexs < be_cut * bindingEnergy(a, z)) explo = false;
    387   };   
    388 
    389   return explo;
    390 }
    391  
Note: See TracChangeset for help on using the changeset viewer.