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

    r1196 r1315  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
     25// $Id: G4CascadeInterface.cc,v 1.77 2010/05/21 18:07:30 mkelsey Exp $
     26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
    2527//
     28// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
     29// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
     30// 20100414  M. Kelsey -- Check for K0L/K0S before using G4InuclElemPart::type
     31// 20100418  M. Kelsey -- Reference output particle lists via const-ref, use
     32//              const_iterator for both.
     33// 20100428  M. Kelsey -- Use G4InuclParticleNames enum
     34// 20100429  M. Kelsey -- Change "case gamma:" to "case photon:"
     35// 20100517  M. Kelsey -- Follow new ctors for G4*Collider family.
     36// 20100520  M. Kelsey -- Simplify collision loop, move momentum rotations to
     37//              G4CollisionOutput, copy G4DynamicParticle directly from
     38//              G4InuclParticle, no switch-block required.
    2639
    2740#include "G4CascadeInterface.hh"
    2841#include "globals.hh"
    29 #include "G4DynamicParticleVector.hh"
    30 #include "G4IonTable.hh"
     42#include "G4CollisionOutput.hh"
     43#include "G4DynamicParticle.hh"
    3144#include "G4InuclCollider.hh"
    32 // #include "G4IntraNucleiCascader.hh"
    3345#include "G4InuclElementaryParticle.hh"
    3446#include "G4InuclNuclei.hh"
    3547#include "G4InuclParticle.hh"
    36 #include "G4CollisionOutput.hh"
     48#include "G4InuclParticleNames.hh"
     49#include "G4KaonZeroShort.hh"
     50#include "G4KaonZeroLong.hh"
     51#include "G4LorentzRotation.hh"
     52#include "G4Nucleus.hh"
     53#include "G4ParticleDefinition.hh"
     54#include "G4Track.hh"
    3755#include "G4V3DNucleus.hh"
    38 #include "G4Track.hh"
    39 #include "G4Nucleus.hh"
    40 #include "G4NucleiModel.hh"
    41 #include "G4LorentzRotation.hh"
    42 
     56
     57using namespace G4InuclParticleNames;
    4358
    4459//#define BERTDEV 1  // A flag to activate a development version of Bertini cascade
    4560
    46 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
    47 typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
     61typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
     62typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
    4863
    4964G4CascadeInterface::G4CascadeInterface(const G4String& nam)
    … …  
    89104  // NOTE: Geant4 units are MeV = 1 and GeV = 1000. Cascade code by default use GeV = 1.
    90105
    91   enum particleType { nuclei = 0, proton = 1, neutron = 2, pionPlus = 3,
    92                       pionMinus = 5, pionZero = 7, photon = 10,
    93                       kaonPlus = 11, kaonMinus = 13, kaonZero = 15,
    94                       kaonZeroBar = 17, lambda = 21, sigmaPlus = 23,
    95                       sigmaZero = 25, sigmaMinus = 27, xiZero = 29, xiMinus = 31 };
    96 
    97   G4int bulletType = 0;
    98 
    99   // Coding particles
    100   if (aTrack.GetDefinition() ==    G4Proton::Proton()    ) bulletType = proton;
    101   if (aTrack.GetDefinition() ==   G4Neutron::Neutron()   ) bulletType = neutron;
    102   if (aTrack.GetDefinition() ==  G4PionPlus::PionPlus()  ) bulletType = pionPlus;
    103   if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;
    104   if (aTrack.GetDefinition() ==  G4PionZero::PionZero()  ) bulletType = pionZero;
    105   if (aTrack.GetDefinition() ==     G4Gamma::Gamma()     ) bulletType = photon;
    106   if (aTrack.GetDefinition() == G4KaonPlus::KaonPlus()     ) bulletType = kaonPlus;
    107   if (aTrack.GetDefinition() == G4KaonMinus::KaonMinus()   ) bulletType = kaonMinus;
    108   if (aTrack.GetDefinition() == G4Lambda::Lambda()         ) bulletType = lambda;
    109   if (aTrack.GetDefinition() == G4SigmaPlus::SigmaPlus()   ) bulletType = sigmaPlus;
    110   if (aTrack.GetDefinition() == G4SigmaZero::SigmaZero()   ) bulletType = sigmaZero;
    111   if (aTrack.GetDefinition() == G4SigmaMinus::SigmaMinus() ) bulletType = sigmaMinus;
    112   if (aTrack.GetDefinition() == G4XiZero::XiZero()         ) bulletType = xiZero;
    113   if (aTrack.GetDefinition() == G4XiMinus::XiMinus()       ) bulletType = xiMinus; 
    114 
     106  G4int bulletType;
    115107  if (aTrack.GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
    116       aTrack.GetDefinition() == G4KaonZeroShort::KaonZeroShort() ) {
    117     if (G4UniformRand() > 0.5) {
    118       bulletType = kaonZero;
    119     } else {
    120       bulletType = kaonZeroBar;
    121     }
    122   }
     108      aTrack.GetDefinition() == G4KaonZeroShort::KaonZeroShort() )
     109    bulletType = (G4UniformRand() > 0.5) ? kaonZero : kaonZeroBar;
     110  else
     111    bulletType = G4InuclElementaryParticle::type(aTrack.GetDefinition());
    123112
    124113  // Code momentum and energy.
    125   G4double px,py,pz;
    126   px=aTrack.Get4Momentum().px() / GeV;
    127   py=aTrack.Get4Momentum().py() / GeV;
    128   pz=aTrack.Get4Momentum().pz() / GeV;
    129 
    130114  G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
    131115  G4LorentzRotation toZ;
    … …  
    134118  G4LorentzRotation toLabFrame = toZ.inverse();
    135119
    136   G4CascadeMomentum momentumBullet;
    137   momentumBullet[0] =0.;
    138   momentumBullet[1] =0;
    139   momentumBullet[2] =0;
    140   momentumBullet[3] =std::sqrt(px*px+py*py+pz*pz);
    141 
    142   G4InuclElementaryParticle *  bullet = new G4InuclElementaryParticle(momentumBullet, bulletType);
     120  G4LorentzVector momentumBullet(0., 0., aTrack.GetTotalMomentum()/GeV,
     121                                 aTrack.GetTotalEnergy()/GeV);
     122
     123  G4InuclElementaryParticle* bullet =
     124    new G4InuclElementaryParticle(momentumBullet, bulletType);
    143125
    144126  sumEnergy = bullet->getKineticEnergy(); // In GeV
    145 
    146   if (bulletType == proton || bulletType == neutron || bulletType == lambda ||
    147       bulletType == sigmaPlus || bulletType == sigmaZero || bulletType == sigmaMinus ||
    148       bulletType == xiZero || bulletType == xiMinus) {
    149 
    150     sumBaryon += 1;
    151   }
     127  sumBaryon += bullet->baryon();        // Returns baryon number (0, 1, or 2)
    152128
    153129  // Set target
    154130  G4InuclNuclei*   target  = 0;
    155131  G4InuclParticle* targetH = 0;
    156   // and outcoming particles
    157   G4DynamicParticle* cascadeParticle = 0;
    158 
    159   G4CascadeMomentum targetMomentum;
    160132
    161133  G4double theNucleusA = theNucleus.GetN();
    162134
    163   if ( !(G4int(theNucleusA) == 1) ) {
    164     target  = new G4InuclNuclei(targetMomentum,
    165                                 theNucleusA,
    166                                 theNucleus.GetZ());
    167     target->setEnergy();
    168 
    169     const G4CascadeMomentum& bmom = bullet->getMomentum();
    170     eInit = std::sqrt(bmom[0] * bmom[0]);
    171     const G4CascadeMomentum& tmom = target->getMomentum();
    172     eInit += std::sqrt(tmom[0] * tmom[0]);
     135  if ( G4int(theNucleusA) != 1 ) {
     136    target  = new G4InuclNuclei(theNucleusA, theNucleus.GetZ());
     137    eInit = bullet->getEnergy() + target->getEnergy();
    173138
    174139    sumBaryon += theNucleusA;
    … …  
    187152
    188153  // Colliders initialisation
    189   //  G4ElementaryParticleCollider*   colep = new G4ElementaryParticleCollider;
    190   //  G4IntraNucleiCascader*            inc = new G4IntraNucleiCascader; // the actual cascade
    191   //  inc->setInteractionCase(1); // Interaction type is particle with nuclei.
    192   inc.setInteractionCase(1); // Interaction type is particle with nuclei.
    193 
    194   //  G4NonEquilibriumEvaporator*     noneq = new G4NonEquilibriumEvaporator;
    195   //  G4EquilibriumEvaporator*         eqil = new G4EquilibriumEvaporator;
    196   //  G4Fissioner*                     fiss = new G4Fissioner;
    197   //  G4BigBanger*                     bigb = new G4BigBanger;
    198   G4InuclCollider* collider = new G4InuclCollider(&colep, &inc, &noneq, &eqil, &fiss, &bigb);
    199 
    200       G4int  maxTries = 100; // maximum tries for inelastic collision to avoid infinite loop
    201       G4int  nTries   = 0;  // try counter
     154  G4InuclCollider* collider = new G4InuclCollider;
     155
     156  G4int  maxTries = 100; // maximum tries for inelastic collision to avoid infinite loop
     157  G4int  nTries   = 0;  // try counter
    202158
    203159#ifdef BERTDEV
    204       G4int coulombOK =0;  // flag for correct Coulomb barrier
    205 #endif
    206       if (G4int(theNucleusA) == 1) { // special treatment for target H(1,1) (proton)
    207 
    208         targetH = new G4InuclElementaryParticle(targetMomentum, 1);
    209 
    210         G4float cutElastic[32];
    211 
    212         cutElastic[proton   ] = 1.0; // GeV
    213         cutElastic[neutron  ] = 1.0;
    214         cutElastic[pionPlus ] = 0.6;
    215         cutElastic[pionMinus] = 0.2;
    216         cutElastic[pionZero ] = 0.2;
    217         cutElastic[kaonPlus ] = 0.5;
    218         cutElastic[kaonMinus] = 0.5;
    219         cutElastic[kaonMinus] = 0.5;
    220         cutElastic[kaonZero] = 0.5;
    221         cutElastic[kaonZeroBar] = 0.5;
    222         cutElastic[lambda] = 1.0;
    223         cutElastic[sigmaPlus] = 1.0;
    224         cutElastic[sigmaZero] = 1.0;
    225         cutElastic[sigmaMinus] = 1.0;
    226         cutElastic[xiZero] = 1.0;
    227         cutElastic[xiMinus] = 1.0;
    228 
    229         if (momentumBullet[3] > cutElastic[bulletType]) { // inelastic collision possible
    230 
    231           do {   // we try to create inelastic interaction
    232             output = collider->collide(bullet, targetH);
    233             nTries++;
    234           } while(
    235                   (nTries < maxTries)                                           &&
    236                   (output.getOutgoingParticles().size() == 2                    && // elastic: bullet + p = H(1,1) coming out
    237                    (output.getOutgoingParticles().begin()->type() == bulletType ||
    238                     output.getOutgoingParticles().begin()->type() == proton)
    239                    )
    240                   );
    241         } else { // only elastic collision is energetically possible
    242           output = collider->collide(bullet, targetH);
     160  G4int coulombOK =0;  // flag for correct Coulomb barrier
     161#endif
     162  if (G4int(theNucleusA) == 1) { // special treatment for target H(1,1) (proton)
     163   
     164    targetH = new G4InuclElementaryParticle(proton);
     165   
     166    G4float cutElastic[32];
     167   
     168    cutElastic[proton   ] = 1.0; // GeV
     169    cutElastic[neutron  ] = 1.0;
     170    cutElastic[pionPlus ] = 0.6;
     171    cutElastic[pionMinus] = 0.2;
     172    cutElastic[pionZero ] = 0.2;
     173    cutElastic[kaonPlus ] = 0.5;
     174    cutElastic[kaonMinus] = 0.5;
     175    cutElastic[kaonZero] = 0.5;
     176    cutElastic[kaonZeroBar] = 0.5;
     177    cutElastic[lambda] = 1.0;
     178    cutElastic[sigmaPlus] = 1.0;
     179    cutElastic[sigmaZero] = 1.0;
     180    cutElastic[sigmaMinus] = 1.0;
     181    cutElastic[xiZero] = 1.0;
     182    cutElastic[xiMinus] = 1.0;
     183   
     184    if (momentumBullet.z() > cutElastic[bulletType]) { // inelastic collision possible
     185     
     186      do {   // we try to create inelastic interaction
     187        output.reset();
     188        collider->collide(bullet, targetH, output);
     189        nTries++;
     190      } while(
     191              (nTries < maxTries) &&
     192              (output.getOutgoingParticles().size() == 2 && // elastic: bullet + p = H(1,1) coming out
     193               (output.getOutgoingParticles().begin()->type() == bulletType ||
     194                output.getOutgoingParticles().begin()->type() == proton)
     195               )
     196              );
     197    } else { // only elastic collision is energetically possible
     198      collider->collide(bullet, targetH, output);
     199    }
     200   
     201    sumBaryon += 1;
     202   
     203    eInit = bullet->getEnergy() + targetH->getEnergy();
     204   
     205    if (verboseLevel > 2) {
     206      G4cout << "Target:  " << G4endl;
     207      targetH->printParticle();
     208    }
     209   
     210  } else {  // treat all other targets excepet H(1,1)
     211   
     212    do { // we try to create inelastic interaction
     213
     214#ifdef BERTDEV
     215      coulombOK=0;  // by default coulomb analysis is OK
     216#endif
     217      output.reset();
     218      collider->collide(bullet, target, output);
     219      nTries++;
     220     
     221#ifdef BERTDEV
     222      G4double coulumbBarrier = 8.7 * MeV;
     223      const std::vector<G4InuclElementaryParticle>& p= output.getOutgoingParticles();
     224      if(!p.empty()) {
     225        for(    particleIterator ipart = p.begin(); ipart != p.end(); ipart++) {
     226          if (ipart->type() == proton) {
     227            G4double e = ipart->getKineticEnergy()*GeV;
     228            if (e < coulumbBarrier) coulombOK= 1; // If event with coulomb barrier violation detected -> retry
     229            //            G4cout << "///AH "<< e << "" << coulumbBarrier <<G4endl;
     230          }
    243231        }
    244 
    245         sumBaryon += 1;
    246 
    247         const G4CascadeMomentum& bmom = bullet->getMomentum();
    248         eInit = std::sqrt(bmom[0] * bmom[0]);
    249         const G4CascadeMomentum& tmom = targetH->getMomentum();
    250         eInit += std::sqrt(tmom[0] * tmom[0]);
    251 
    252         if (verboseLevel > 2) {
    253           G4cout << "Target:  " << G4endl;
    254           targetH->printParticle();
    255         }
    256 
    257       } else {  // treat all other targets excepet H(1,1)
    258 
    259         do  // we try to create inelastic interaction
    260           {
     232      }
     233#endif
     234    } while(
     235            (nTries < maxTries) &&  // conditions for next try
     236            (output.getOutgoingParticles().size()!=0) &&
    261237#ifdef BERTDEV
    262             coulombOK=0;  // by default coulomb analysis is OK
    263 #endif
    264             output = collider->collide(bullet, target );
    265             nTries++;
    266 
    267 #ifdef BERTDEV
    268             G4double coulumbBarrier = 8.7 * MeV;
    269             std::vector<G4InuclElementaryParticle> p= output.getOutgoingParticles();
    270             if(!p.empty()) {
    271               for(    particleIterator ipart = p.begin(); ipart != p.end(); ipart++) {
    272                 if (ipart->type() == proton) {
    273                   G4double e = ipart->getKineticEnergy()*GeV;
    274                   if (e < coulumbBarrier) coulombOK= 1; // If event with coulomb barrier violation detected -> retry
    275                   //              G4cout << "///AH "<< e << "" << coulumbBarrier <<G4endl;
    276                 }
    277               }
    278             }
    279           } while(
    280                   (nTries < maxTries) &&  // conditions for next try
    281                   (coulombOK==1) &&
    282                   ((output.getOutgoingParticles().size() + output.getNucleiFragments().size()) > 2.5) && 
    283                   (output.getOutgoingParticles().size()!=0)                                       
    284                   );
     238            (coulombOK==1) &&
     239            ((output.getOutgoingParticles().size() + output.getNucleiFragments().size()) > 2.5)
    285240#else
    286 
    287           } while(
    288                    (nTries < maxTries)                                                               &&
    289                    (output.getOutgoingParticles().size() + output.getNucleiFragments().size() < 2.5) &&
    290                    (output.getOutgoingParticles().size()!=0) &&
    291                    (output.getOutgoingParticles().begin()->type()==bullet->type())
    292                    );
    293 #endif
    294      }
    295 
    296   if (verboseLevel > 1)
    297     {
    298       G4cout << " Cascade output: " << G4endl;
    299       output.printCollisionOutput();
    300     }
     241            ((output.getOutgoingParticles().size() + output.getNucleiFragments().size()) < 2.5) && 
     242            (output.getOutgoingParticles().begin()->type()==bullet->type())
     243#endif
     244             );
     245  }
     246
     247  if (verboseLevel > 1) {
     248    G4cout << " Cascade output: " << G4endl;
     249    output.printCollisionOutput();
     250  }
    301251 
     252  // Rotate event to put Z axis along original projectile direction
     253  output.rotateEvent(toLabFrame);
     254
    302255  // Convert cascade data to use hadronics interface
    303   std::vector<G4InuclNuclei>            nucleiFragments = output.getNucleiFragments();
    304   std::vector<G4InuclElementaryParticle> particles =      output.getOutgoingParticles();
     256  const std::vector<G4InuclNuclei>& nucleiFragments = output.getNucleiFragments();
     257  const std::vector<G4InuclElementaryParticle>& particles = output.getOutgoingParticles();
    305258
    306259  theResult.SetStatusChange(stopAndKill);
    307260
     261  // Get outcoming particles
     262  G4DynamicParticle* cascadeParticle = 0;
    308263  if (!particles.empty()) {
    309     particleIterator ipart;
    310     G4int outgoingParticle;
    311 
    312     for (ipart = particles.begin(); ipart != particles.end(); ipart++) {
    313       outgoingParticle = ipart->type();
    314       const G4CascadeMomentum& mom = ipart->getMomentum();
    315       eTot   += std::sqrt(mom[0] * mom[0]);
    316 
    317       G4double ekin = ipart->getKineticEnergy() * GeV;
    318       G4ThreeVector aMom(mom[1], mom[2], mom[3]);
    319       aMom = aMom.unit();
    320 
    321       if (ipart->baryon() ) {
    322         sumBaryon -= 1;
     264    particleIterator ipart = particles.begin();
     265    for (; ipart != particles.end(); ipart++) {
     266      G4int outgoingType = ipart->type();
     267
     268      eTot += ipart->getEnergy();
     269      sumBaryon -= ipart->baryon();
     270      sumEnergy -= ipart->getKineticEnergy();
     271
     272      if (!ipart->valid() || ipart->quasi_deutron()) {
     273        G4cerr << " ERROR: G4CascadeInterface::Propagate incompatible"
     274               << " particle type " << outgoingType << G4endl;
     275        continue;
    323276      }
    324277
    325       sumEnergy -= ekin / GeV;
    326 
    327       switch(outgoingParticle) {
    328 
    329       case proton:
    330 #ifdef debug_G4CascadeInterface
    331         G4cerr << "proton " << counter << " " << aMom << " " << ekin << G4endl;
    332 #endif
    333         cascadeParticle =
    334           new G4DynamicParticle(G4Proton::ProtonDefinition(), aMom, ekin);
    335         break;
    336 
    337       case neutron:
    338 
    339 #ifdef debug_G4CascadeInterface
    340         G4cerr << "neutron "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
    341 #endif
    342         cascadeParticle =
    343           new G4DynamicParticle(G4Neutron::NeutronDefinition(), aMom, ekin);
    344         break;
    345 
    346       case pionPlus:
    347         cascadeParticle =
    348           new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), aMom, ekin);
    349 
    350 #ifdef debug_G4CascadeInterface
    351         G4cerr << "pionPlus "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
    352 #endif
    353         break;
    354 
    355       case pionMinus:
    356         cascadeParticle =
    357           new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), aMom, ekin);
    358 
    359 #ifdef debug_G4CascadeInterface
    360         G4cerr << "pionMinus "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
    361 #endif
    362         break;
    363 
    364       case pionZero:
    365         cascadeParticle =
    366           new G4DynamicParticle(G4PionZero::PionZeroDefinition(), aMom, ekin);
    367 
    368 #ifdef debug_G4CascadeInterface
    369         G4cerr << "pionZero "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
    370 #endif
    371         break;
    372 
    373       case photon:
    374         cascadeParticle =
    375           new G4DynamicParticle(G4Gamma::Gamma(), aMom, ekin);
    376 
    377 #ifdef debug_G4CascadeInterface
    378         G4cerr << "photon "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
    379 #endif
    380         break;
    381 
    382 
    383       case kaonPlus:
    384         cascadeParticle =
    385           new G4DynamicParticle(G4KaonPlus::KaonPlusDefinition(), aMom, ekin);
    386         break;
    387 
    388       case kaonMinus:
    389         cascadeParticle =
    390          new G4DynamicParticle(G4KaonMinus::KaonMinusDefinition(), aMom, ekin);
    391         break;
    392 
    393       case kaonZero:
    394         if (G4UniformRand() > 0.5) {
    395           cascadeParticle = new G4DynamicParticle(
    396                               G4KaonZeroLong::KaonZeroLongDefinition(),
    397                               aMom, ekin);
    398         } else {
    399           cascadeParticle = new G4DynamicParticle(
    400                               G4KaonZeroShort::KaonZeroShortDefinition(),
    401                               aMom, ekin);
    402         }
    403         break;
    404 
    405       case kaonZeroBar:
    406         if (G4UniformRand() > 0.5) {
    407           cascadeParticle = new G4DynamicParticle(
    408                               G4KaonZeroLong::KaonZeroLongDefinition(),
    409                               aMom, ekin);
    410         } else {
    411           cascadeParticle = new G4DynamicParticle(
    412                               G4KaonZeroShort::KaonZeroShortDefinition(),
    413                               aMom, ekin);
    414         }
    415         break;
    416 
    417       case lambda:
    418         cascadeParticle =
    419          new G4DynamicParticle(G4Lambda::LambdaDefinition(), aMom, ekin);
    420         break;
    421 
    422       case sigmaPlus:
    423         cascadeParticle =
    424          new G4DynamicParticle(G4SigmaPlus::SigmaPlusDefinition(), aMom, ekin);
    425         break;
    426 
    427       case sigmaZero:
    428         cascadeParticle =
    429          new G4DynamicParticle(G4SigmaZero::SigmaZeroDefinition(), aMom, ekin);
    430         break;
    431 
    432       case sigmaMinus:
    433         cascadeParticle =
    434          new G4DynamicParticle(G4SigmaMinus::SigmaMinusDefinition(), aMom, ekin);
    435         break;
    436 
    437       case xiZero:
    438         cascadeParticle =
    439          new G4DynamicParticle(G4XiZero::XiZeroDefinition(), aMom, ekin);
    440         break;
    441 
    442       case xiMinus:
    443         cascadeParticle =
    444          new G4DynamicParticle(G4XiMinus::XiMinusDefinition(), aMom, ekin);
    445         break;
    446 
    447       default:
    448         G4cout << " ERROR: G4CascadeInterface::Propagate undefined particle type"
    449                << G4endl;
     278      // Copy local G4DynPart to public output (handle kaon mixing specially)
     279      if (outgoingType == kaonZero || outgoingType == kaonZeroBar) {
     280        G4ThreeVector momDir = ipart->getMomentum().vect().unit();
     281        G4double ekin = ipart->getKineticEnergy();
     282
     283        G4ParticleDefinition* pd = G4KaonZeroShort::Definition();
     284        if (G4UniformRand() > 0.5) pd = G4KaonZeroLong::Definition();
     285
     286        cascadeParticle = new G4DynamicParticle(pd, momDir, ekin);
     287      } else {
     288        cascadeParticle = new G4DynamicParticle(ipart->getDynamicParticle());
    450289      }
    451290
    452       cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
    453291      theResult.AddSecondary(cascadeParticle);
    454292    }
    … …  
    457295  // get nuclei fragments
    458296  G4DynamicParticle * aFragment = 0;
    459   G4ParticleDefinition * aIonDef = 0;
    460   G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
    461 
    462297  if (!nucleiFragments.empty()) {
    463     nucleiIterator ifrag;
    464 
    465     for (ifrag = nucleiFragments.begin(); ifrag != nucleiFragments.end(); ifrag++)
    466       {
    467         G4double eKin = ifrag->getKineticEnergy() * GeV;
    468         const G4CascadeMomentum& mom = ifrag->getMomentum();
    469         eTot   += std::sqrt(mom[0] * mom[0]);
    470 
    471         G4ThreeVector aMom(mom[1], mom[2], mom[3]);
    472         aMom = aMom.unit();
    473 
    474         // hpw @@@ ==> Should be zero: G4double fragmentExitation = ifrag->getExitationEnergyInGeV();
    475 
    476         if (verboseLevel > 2) {
    477           G4cout << " Nuclei fragment: " << G4endl;
    478           ifrag->printParticle();
    479         }
    480 
    481         G4int A = G4int(ifrag->getA());
    482         G4int Z = G4int(ifrag->getZ());
    483         aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
    484      
    485         aFragment =  new G4DynamicParticle(aIonDef, aMom, eKin);
    486 
    487         sumBaryon -= A;
    488         sumEnergy -= eKin / GeV;
    489 
    490         aFragment->Set4Momentum(aFragment->Get4Momentum()*=toLabFrame);
    491         theResult.AddSecondary(aFragment);
     298    nucleiIterator ifrag = nucleiFragments.begin();
     299    for (; ifrag != nucleiFragments.end(); ifrag++) {
     300      eTot += ifrag->getEnergy();
     301      sumBaryon -= ifrag->getA();
     302      sumEnergy -= ifrag->getKineticEnergy();
     303     
     304      // hpw @@@ ==> Should be zero: G4double fragmentExitation = ifrag->getExitationEnergyInGeV();
     305     
     306      if (verboseLevel > 2) {
     307        G4cout << " Nuclei fragment: " << G4endl;
     308        ifrag->printParticle();
    492309      }
    493   }
    494 
     310     
     311      // Copy local G4DynPart to public output
     312      aFragment =  new G4DynamicParticle(ifrag->getDynamicParticle());
     313      theResult.AddSecondary(aFragment);
     314    }
     315  }
     316
     317  // Report violations of energy, baryon conservation
    495318  if (verboseLevel > 2) {
    496319    if (sumBaryon != 0) {
    … …  
    514337
    515338  delete bullet;
    516 //  delete inc;
    517339  delete collider;
    518340
    519341  if(target != 0) delete target;
    520342  if(targetH != 0) delete targetH;
    521  // if(cascadeParticle != 0) delete cascadeParticle;
    522  // if(aFragment != 0) delete aFragment;
    523343
    524344  return &theResult;
Note: See TracChangeset for help on using the changeset viewer.