Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4InuclCollider.cc,v 1.30.2.1 2010/06/25 09:44:36 gunter Exp $
    27 // Geant4 tag: $Name: geant4-09-04-beta-01 $
     26// $Id: G4InuclCollider.cc,v 1.51 2010/10/19 19:48:39 mkelsey Exp $
     27// Geant4 tag: $Name: hadr-casc-V09-03-85 $
    2828//
    2929// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
    … …  
    3434// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
    3535//              simple data members, consolidate code
     36// 20100620  M. Kelsey -- Reorganize top level if-blocks to reduce nesting,
     37//              use new four-vector conservation check.
     38// 20100701  M. Kelsey -- Bug fix energy-conservation after equilibrium evap,
     39//              pass verbosity through to G4CollisionOutput
     40// 20100714  M. Kelsey -- Move conservation checking to base class, report
     41//              number of iterations at end
     42// 20100715  M. Kelsey -- Remove all the "Before xxx" and "After xxx"
     43//              conservation checks, as the colliders now all do so.  Move
     44//              local buffers outside while() loop, use new "combined add()"
     45//              interface for copying local buffers to global.
     46// 20100716  M. Kelsey -- Drop G4IntraNucleiCascader::setInteractionCase()
     47// 20100720  M. Kelsey -- Make all the collders pointer members (to reducde
     48//              external compile dependences).
     49// 20100915  M. Kelsey -- Move post-cascade colliders to G4CascadeDeexcitation,
     50//              simplify operational code somewhat
     51// 20100922  M. Kelsey -- Add functions to select de-excitation method;
     52//              default is G4CascadeDeexcitation (i.e., built-in modules)
     53// 20100924  M. Kelsey -- Migrate to integer A and Z
     54// 20101019  M. Kelsey -- CoVerity report: check dynamic_cast<> for null
    3655
    3756#include "G4InuclCollider.hh"
    38 #include "G4BigBanger.hh"
     57#include "G4CascadeDeexcitation.hh"
     58#include "G4PreCompoundDeexcitation.hh"
    3959#include "G4CollisionOutput.hh"
    4060#include "G4ElementaryParticleCollider.hh"
    41 #include "G4EquilibriumEvaporator.hh"
    4261#include "G4IntraNucleiCascader.hh"
    4362#include "G4InuclElementaryParticle.hh"
     63#include "G4InuclNuclei.hh"
    4464#include "G4LorentzConvertor.hh"
    45 #include "G4NonEquilibriumEvaporator.hh"
    4665
    4766
    4867G4InuclCollider::G4InuclCollider()
    49   : G4VCascadeCollider("G4InuclCollider"),
     68  : G4CascadeColliderBase("G4InuclCollider"),
    5069    theElementaryParticleCollider(new G4ElementaryParticleCollider),
    5170    theIntraNucleiCascader(new G4IntraNucleiCascader),
    52     theNonEquilibriumEvaporator(new G4NonEquilibriumEvaporator),
    53     theEquilibriumEvaporator(new G4EquilibriumEvaporator),
    54     theBigBanger(new G4BigBanger) {}
     71    theDeexcitation(new G4CascadeDeexcitation) {}
    5572
    5673G4InuclCollider::~G4InuclCollider() {
    5774  delete theElementaryParticleCollider;
    5875  delete theIntraNucleiCascader;
    59   delete theNonEquilibriumEvaporator;
    60   delete theEquilibriumEvaporator;
    61   delete theBigBanger;
    62 }
    63 
     76  delete theDeexcitation;
     77}
     78
     79
     80// Select post-cascade processing (default will be CascadeDeexcitation)
     81
     82void G4InuclCollider::useCascadeDeexcitation() {
     83  delete theDeexcitation;
     84  theDeexcitation = new G4CascadeDeexcitation;
     85}
     86
     87void G4InuclCollider::usePreCompoundDeexcitation() {
     88  delete theDeexcitation;
     89  theDeexcitation = new G4PreCompoundDeexcitation;
     90}
     91
     92
     93// Main action
    6494
    6595void G4InuclCollider::collide(G4InuclParticle* bullet, G4InuclParticle* target,
    6696                              G4CollisionOutput& globalOutput) {
     97  if (verboseLevel) G4cout << " >>> G4InuclCollider::collide" << G4endl;
     98
     99  // Initialize colliders verbosity
     100  theElementaryParticleCollider->setVerboseLevel(verboseLevel);
     101  theIntraNucleiCascader->setVerboseLevel(verboseLevel);
     102  theDeexcitation->setVerboseLevel(verboseLevel);
     103
     104  output.setVerboseLevel(verboseLevel);
     105  DEXoutput.setVerboseLevel(verboseLevel);
     106
     107  const G4int itry_max = 1000;
     108
     109  // Particle-on-particle collision; no nucleus involved
     110  if (useEPCollider(bullet,target)) {
     111    if (verboseLevel > 2)
     112      G4cout << " InuclCollider -> particle on particle collision" << G4endl;
     113 
     114    theElementaryParticleCollider->collide(bullet, target, globalOutput);
     115    return;
     116  }
     117 
     118  interCase.set(bullet,target);         // Classify collision type
     119  if (verboseLevel > 2) {
     120    G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
     121  }
     122
     123  if (!interCase.valid()) {
     124    if (verboseLevel > 1)
     125      G4cerr << " InuclCollider -> no collision possible " << G4endl;
     126
     127    globalOutput.trivialise(bullet, target);
     128    return;
     129  }
     130
     131  // Target must be a nucleus
     132  G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
     133  if (!ntarget) {
     134    G4cerr << " InuclCollider -> ERROR target is not a nucleus " << G4endl;
     135
     136    globalOutput.trivialise(bullet, target);
     137    return;
     138  }
     139
     140  G4int btype = 0;
     141  G4int ab = 0;
     142  G4int zb = 0;
     143 
     144  if (interCase.hadNucleus()) {         // particle with nuclei
     145    G4InuclElementaryParticle* pbullet =
     146      dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
     147    if (!pbullet) {
     148      G4cerr << " InuclCollider -> ERROR bullet is not a hadron " << G4endl;
     149      globalOutput.trivialise(bullet, target);
     150      return;
     151    } else if (pbullet->isPhoton()) {
     152      G4cerr << " InuclCollider -> can not collide with photon " << G4endl;
     153      globalOutput.trivialise(bullet, target);
     154      return;
     155    } else {
     156      btype = pbullet->type();
     157    }
     158  } else {                              // nuclei with nuclei
     159    G4InuclNuclei* nbullet =
     160      dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
     161    if (!nbullet) {
     162      G4cerr << " InuclCollider -> ERROR bullet is not a nucleus " << G4endl;
     163      globalOutput.trivialise(bullet, target);
     164      return;
     165    }
     166   
     167    ab = nbullet->getA();
     168    zb = nbullet->getZ();
     169  }
     170
     171  G4LorentzConvertor convertToTargetRestFrame(bullet, ntarget);
     172  G4double ekin = convertToTargetRestFrame.getKinEnergyInTheTRS();
     173 
     174  if (verboseLevel > 3) G4cout << " ekin in trs " << ekin << G4endl;
     175
     176  if (!inelasticInteractionPossible(bullet, target, ekin)) {
     177    if (verboseLevel > 3)
     178      G4cout << " InuclCollider -> inelastic interaction is impossible\n"
     179             << " due to the coulomb barirer " << G4endl;
     180
     181    globalOutput.trivialise(bullet, target);
     182    return;
     183  }
     184
     185  // Generate interaction secondaries in rest frame of target nucleus
     186  convertToTargetRestFrame.toTheTargetRestFrame();
    67187  if (verboseLevel > 3) {
    68     G4cout << " >>> G4InuclCollider::collide" << G4endl;
    69   }
    70 
    71   const G4int itry_max = 1000;
    72                      
    73   if (useEPCollider(bullet,target)) {
    74     if (verboseLevel > 2) {
    75       bullet->printParticle();
    76       target->printParticle();
     188    G4cout << " degenerated? " << convertToTargetRestFrame.trivial()
     189           << G4endl;
     190  }
     191 
     192  G4LorentzVector bmom;                 // Bullet is along local Z
     193  bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
     194
     195  // Need to make copy of bullet with momentum realigned
     196  G4InuclParticle* zbullet = 0;
     197  if (interCase.hadNucleus())
     198    zbullet = new G4InuclElementaryParticle(bmom, btype);
     199  else
     200    zbullet = new G4InuclNuclei(bmom, ab, zb);
     201
     202  G4int itry = 0;
     203  while (itry < itry_max) {
     204    itry++;
     205    if (verboseLevel > 2)
     206      G4cout << " IntraNucleiCascader itry " << itry << G4endl;
     207
     208    globalOutput.reset();               // Clear buffers for this attempt
     209    output.reset();     
     210    DEXoutput.reset();
     211
     212    theIntraNucleiCascader->collide(zbullet, target, output);
     213   
     214    if (verboseLevel > 1) G4cout << " After Cascade " << G4endl;
     215
     216    // FIXME:  The code below still does too much copying!  Would rather
     217    //         remove initial fragment from list (or get it a different way)
     218    DEXoutput.addOutgoingParticles(output.getOutgoingParticles());
     219   
     220    if (output.numberOfOutgoingNuclei() == 1) { // Residual fragment
     221      // FIXME:  Making a copy here because of constness issues
     222      G4InuclNuclei recoil_nucleus = output.getOutgoingNuclei()[0];
     223      theDeexcitation->collide(0, &recoil_nucleus, DEXoutput);
    77224    }
    78225
    79     theElementaryParticleCollider->collide(bullet, target, globalOutput);
    80   } else { // needs to call all machinery       
    81     G4LorentzConvertor convertToTargetRestFrame;
    82 
    83     interCase.set(bullet,target);
    84     if (interCase.valid()) { // ok
    85       G4InuclNuclei* ntarget =
    86         dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
    87 
    88       convertToTargetRestFrame.setTarget(ntarget);
    89       G4int btype = 0;
    90       G4double ab = 0.0;
    91       G4double zb = 0.0;
    92       G4double at = ntarget->getA();
    93       G4double zt = ntarget->getZ();
    94        
    95       if (interCase.hadNucleus()) { // particle with nuclei
    96         G4InuclElementaryParticle* pbullet =
    97           dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
    98          
    99         if (pbullet->isPhoton()) {
    100           G4cerr << " InuclCollider -> can not collide with photon " << G4endl;
    101 
    102           globalOutput.trivialise(bullet, target);
    103           return;
    104         } else {
    105           convertToTargetRestFrame.setBullet(pbullet);   
    106           btype = pbullet->type();
    107         };
    108 
    109       } else { // nuclei with nuclei
    110         G4InuclNuclei* nbullet =
    111           dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
    112 
    113         convertToTargetRestFrame.setBullet(nbullet);   
    114         ab = nbullet->getA();
    115         zb = nbullet->getZ();
    116       };
    117        
    118       G4double ekin = convertToTargetRestFrame.getKinEnergyInTheTRS();
    119 
    120       if (verboseLevel > 3) {
    121         G4cout << " ekin in trs " << ekin << G4endl;
    122       }
    123 
    124       if (inelasticInteractionPossible(bullet, target, ekin)) {
    125         convertToTargetRestFrame.toTheTargetRestFrame();
    126 
    127         if (verboseLevel > 3) {
    128           G4cout << " degenerated? " << convertToTargetRestFrame.trivial() << G4endl;
    129         }
    130 
    131         G4LorentzVector bmom;
    132         bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
    133 
    134         G4InuclNuclei ntarget(at, zt);          // Default is at rest
    135 
    136         theIntraNucleiCascader->setInteractionCase(interCase.code());
    137          
    138         G4bool bad = true;
    139         G4int itry = 0;
    140          
    141         G4CollisionOutput TRFoutput;
    142         G4CollisionOutput output;
    143         while (bad && itry < itry_max) {
    144           itry++;
    145 
    146           output.reset();       // Clear buffers for this attempt
    147           TRFoutput.reset();
    148 
    149           if (interCase.hadNucleus()) {
    150             G4InuclElementaryParticle pbullet(bmom, btype);
    151 
    152             theIntraNucleiCascader->collide(&pbullet, &ntarget, output);
    153           } else {
    154             G4InuclNuclei nbullet(bmom, ab, zb);
    155             theIntraNucleiCascader->collide(&nbullet, &ntarget, output);
    156           };   
    157 
    158           if (verboseLevel > 3) {
    159             G4cout << " After Cascade " << G4endl;
    160             output.printCollisionOutput();
    161           }
    162          
    163           // the rest, if any
    164           // FIXME:  The code below still does too much copying!
    165           TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
    166 
    167           if (output.numberOfNucleiFragments() == 1) { // there is smth. after
    168             G4InuclNuclei cascad_rec_nuclei = output.getNucleiFragments()[0];
    169             if (explosion(&cascad_rec_nuclei)) {
    170               if (verboseLevel > 3) {
    171                 G4cout << " big bang after cascade " << G4endl;
    172               };
    173 
    174               theBigBanger->collide(0,&cascad_rec_nuclei, TRFoutput);
    175             } else {
    176               output.reset();
    177               theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei, output);
    178 
    179               if (verboseLevel > 3) {
    180                 G4cout << " After NonEquilibriumEvaporator " << G4endl;
    181                 output.printCollisionOutput();
    182               };
    183 
    184               TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
    185               G4InuclNuclei exiton_rec_nuclei = output.getNucleiFragments()[0];
    186 
    187               output.reset();
    188               theEquilibriumEvaporator->collide(0, &exiton_rec_nuclei, output);
    189 
    190               if (verboseLevel > 3) {
    191                 G4cout << " After EquilibriumEvaporator " << G4endl;
    192                 output.printCollisionOutput();
    193               };
    194 
    195               TRFoutput.addOutgoingParticles(output.getOutgoingParticles()); 
    196               TRFoutput.addTargetFragments(output.getNucleiFragments());
    197             };
    198           };
    199          
    200           // convert to the LAB
    201           TRFoutput.boostToLabFrame(convertToTargetRestFrame);
    202 
    203           globalOutput.addOutgoingParticles(TRFoutput.getOutgoingParticles());
    204           globalOutput.addTargetFragments(TRFoutput.getNucleiFragments());
    205           globalOutput.setOnShell(bullet, target);
    206           if (globalOutput.acceptable()) return;
    207 
    208           globalOutput.reset();         // Clear and try again
    209         };
    210 
    211         if (verboseLevel > 3) {
    212           G4cout << " InuclCollider -> can not generate acceptable inter. after "
    213                  << itry_max << " attempts " << G4endl;
    214         }
    215       } else {
    216         if (verboseLevel > 3) {
    217           G4cout << " InuclCollider -> inelastic interaction is impossible " << G4endl
    218                  << " due to the coulomb barirer " << G4endl;
    219         }
    220       }
    221 
    222       globalOutput.trivialise(bullet, target);
    223       return;
    224     } else {
    225       if (verboseLevel > 3) {
    226         G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
    227       };
    228     };       
    229   };
    230 
     226    if (verboseLevel > 2)
     227      G4cout << " itry " << itry << " finished, moving to lab frame" << G4endl;
     228
     229    // convert to the LAB frame and add to final result
     230    DEXoutput.boostToLabFrame(convertToTargetRestFrame);
     231    globalOutput.add(DEXoutput);
     232
     233    // Adjust final state particles to balance momentum and energy
     234    // FIXME:  This should no longer be necessary!
     235    globalOutput.setOnShell(bullet, target);
     236    if (globalOutput.acceptable()) {
     237      if (verboseLevel)
     238        G4cout << " InuclCollider output after trials " << itry << G4endl;
     239      return;
     240    }
     241  }     // while (itry < itry_max)
     242 
     243  if (verboseLevel) {
     244    G4cout << " InuclCollider -> can not generate acceptable inter. after "
     245           << itry_max << " attempts " << G4endl;
     246  }
     247 
     248  globalOutput.trivialise(bullet, target);
    231249  return;
    232250}
Note: See TracChangeset for help on using the changeset viewer.