Changeset 1315 for trunk/source/processes/hadronic/models/cascade/cascade/src/G4PreCompoundCascadeInterface.cc

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:

• trunk/source/processes/hadronic/models/cascade/cascade/src/G4PreCompoundCascadeInterface.cc (modified) (11 diffs)

trunk/source/processes/hadronic/models/cascade/cascade/src/G4PreCompoundCascadeInterface.cc

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: G4PreCompoundCascadeInterface.cc,v 1.13 2010/05/21 18:07:30 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
 //
+// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
+// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
+// 20100414  M. Kelsey -- Check for K0L/K0S before using G4InuclElemPart::type
+// 20100419  M. Kelsey -- Access G4CollisionOutput lists by const-ref, and 
+//              const_iterator
+// 20100428  M. Kelsey -- Use G4InuclParticleNames enum
+// 20100429  M. Kelsey -- Change "case gamma:" to "case photon:"
+// 20100517  M. Kelsey -- Follow new ctors for G4*Collider family.
+// 20100520  M. Kelsey -- Add missing name string to ctor, follow code changes
+//              from G4CascadeInterface.
 
 #include "G4PreCompoundCascadeInterface.hh"
 #include "globals.hh"
-#include "G4DynamicParticleVector.hh"
-#include "G4IonTable.hh"
-#include "G4PreCompoundInuclCollider.hh"
-#include "G4IntraNucleiCascader.hh"
-#include "G4ElementaryParticleCollider.hh"
-#include "G4NonEquilibriumEvaporator.hh"
-#include "G4BigBanger.hh"
+#include "G4CollisionOutput.hh"
+#include "G4DynamicParticle.hh"
 #include "G4InuclElementaryParticle.hh"
 #include "G4InuclNuclei.hh"
 #include "G4InuclParticle.hh"
-#include "G4CollisionOutput.hh"
+#include "G4InuclParticleNames.hh"
+#include "G4KaonZeroShort.hh"
+#include "G4KaonZeroLong.hh"
+#include "G4LorentzRotation.hh"
+#include "G4Nucleus.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4PreCompoundInuclCollider.hh"
+#include "G4Track.hh"
 #include "G4V3DNucleus.hh"
-#include "G4Track.hh"
-#include "G4Nucleus.hh"
-#include "G4NucleiModel.hh"
-#include "G4LorentzRotation.hh"
-
-
-typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
-typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
-
-G4PreCompoundCascadeInterface::G4PreCompoundCascadeInterface()
-  :verboseLevel(0)  {
+
+using namespace G4InuclParticleNames;
+
+
+typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
+typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
+
+G4PreCompoundCascadeInterface::G4PreCompoundCascadeInterface(const G4String& nam)
+  :G4VIntraNuclearTransportModel(nam), verboseLevel(0)  {
 
   if (verboseLevel > 3) {
@@ -85 +97 @@
   // NOTE: Geant4 units are MeV = 1 and GeV = 1000. Cascade code by default use GeV = 1.
 
-  enum particleType { nuclei      = 0,  proton     = 1,  neutron   = 2,  pionPlus = 3,
-                      pionMinus   = 5,  pionZero   = 7,  photon    = 10,
-                      kaonPlus    = 11, kaonMinus  = 13, kaonZero  = 15,
-                      kaonZeroBar = 17, lambda     = 21, sigmaPlus = 23,
-                      sigmaZero   = 25, sigmaMinus = 27, xiZero    = 29, xiMinus  = 31 };
-
-  G4int bulletType = 0;
-
-  // Coding particles 
-  if (aTrack.GetDefinition() == G4Proton::Proton()         ) bulletType = proton;
-  if (aTrack.GetDefinition() == G4Neutron::Neutron()       ) bulletType = neutron;
-  if (aTrack.GetDefinition() == G4PionPlus::PionPlus()     ) bulletType = pionPlus;
-  if (aTrack.GetDefinition() == G4PionMinus::PionMinus()   ) bulletType = pionMinus;
-  if (aTrack.GetDefinition() == G4PionZero::PionZero()     ) bulletType = pionZero;
-  if (aTrack.GetDefinition() == G4Gamma::Gamma()           ) bulletType = photon;
-  if (aTrack.GetDefinition() == G4KaonPlus::KaonPlus()     ) bulletType = kaonPlus;
-  if (aTrack.GetDefinition() == G4KaonMinus::KaonMinus()   ) bulletType = kaonMinus;
-  if (aTrack.GetDefinition() == G4Lambda::Lambda()         ) bulletType = lambda;
-  if (aTrack.GetDefinition() == G4SigmaPlus::SigmaPlus()   ) bulletType = sigmaPlus;
-  if (aTrack.GetDefinition() == G4SigmaZero::SigmaZero()   ) bulletType = sigmaZero;
-  if (aTrack.GetDefinition() == G4SigmaMinus::SigmaMinus() ) bulletType = sigmaMinus;
-  if (aTrack.GetDefinition() == G4XiZero::XiZero()         ) bulletType = xiZero;
-  if (aTrack.GetDefinition() == G4XiMinus::XiMinus()       ) bulletType = xiMinus;  
-
+  G4int bulletType;
   if (aTrack.GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
-      aTrack.GetDefinition() == G4KaonZeroShort::KaonZeroShort() ) {
-    if (G4UniformRand() > 0.5) {
-      bulletType = kaonZero;
-    } else {
-      bulletType = kaonZeroBar;
-    }
-  }
+      aTrack.GetDefinition() == G4KaonZeroShort::KaonZeroShort() )
+    bulletType = (G4UniformRand() > 0.5) ? kaonZero : kaonZeroBar;
+  else 
+    bulletType = G4InuclElementaryParticle::type(aTrack.GetDefinition());
 
   // Code momentum and energy.
-  G4double px,py,pz;
-  px=aTrack.Get4Momentum().px() / GeV;
-  py=aTrack.Get4Momentum().py() / GeV;
-  pz=aTrack.Get4Momentum().pz() / GeV;
-
   G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
   G4LorentzRotation toZ;
@@ -130 +111 @@
   G4LorentzRotation toLabFrame = toZ.inverse();
 
-  G4CascadeMomentum momentumBullet;
-  momentumBullet[0] =0.;
-  momentumBullet[1] =0;
-  momentumBullet[2] =0;
-  momentumBullet[3] =std::sqrt(px*px+py*py+pz*pz);
-
-  G4InuclElementaryParticle *  bullet = new G4InuclElementaryParticle(momentumBullet, bulletType); 
+  G4LorentzVector momentumBullet(0., 0., aTrack.GetTotalMomentum()/GeV,
+                                 aTrack.GetTotalEnergy()/GeV);
+
+  G4InuclElementaryParticle *  bullet =
+    new G4InuclElementaryParticle(momentumBullet, bulletType); 
 
   sumEnergy = bullet->getKineticEnergy(); // In GeV 
-
-  if (bulletType == proton || bulletType == neutron || bulletType == lambda ||
-      bulletType == sigmaPlus || bulletType == sigmaZero || bulletType == sigmaMinus ||
-      bulletType == xiZero || bulletType == xiMinus) {
-
-    sumBaryon += 1;
-  } 
+  sumBaryon += bullet->baryon();
 
   // Set target
   G4InuclNuclei*   target  = 0;
   G4InuclParticle* targetH = 0;
-  // and outcoming particles
-  G4DynamicParticle* cascadeParticle = 0;
-
-  G4CascadeMomentum targetMomentum;
 
   G4double theNucleusA = theNucleus.GetN();
 
   if ( !(G4int(theNucleusA) == 1) ) {
-    target  = new G4InuclNuclei(targetMomentum, 
-                                theNucleusA, 
-                                theNucleus.GetZ());
-    target->setEnergy();
-
-    const G4CascadeMomentum& bmom = bullet->getMomentum();
-    eInit = std::sqrt(bmom[0] * bmom[0]);
-    const G4CascadeMomentum& tmom = target->getMomentum();
-    eInit += std::sqrt(tmom[0] * tmom[0]);
+    target  = new G4InuclNuclei(theNucleusA, theNucleus.GetZ());
+    eInit = bullet->getEnergy() + target->getEnergy();
 
     sumBaryon += theNucleusA;
@@ -183 +145 @@
 
   // Colliders initialisation
-  G4ElementaryParticleCollider*   colep = new G4ElementaryParticleCollider;
-
-  G4IntraNucleiCascader*            inc = new G4IntraNucleiCascader; // the actual cascade
-  inc->setInteractionCase(1); // Interaction type is particle with nuclei.
-
-  G4NonEquilibriumEvaporator*     noneq = new G4NonEquilibriumEvaporator;
-  G4BigBanger*                     bigb = new G4BigBanger;
-  G4PreCompoundInuclCollider*  collider = new G4PreCompoundInuclCollider(colep, inc, noneq, bigb);
+  G4PreCompoundInuclCollider*  collider = new G4PreCompoundInuclCollider;
 
   G4int  maxTries = 10; // maximum tries for inelastic collision to avoid infinite loop
@@ -197 +152 @@
   if (G4int(theNucleusA) == 1) { // special treatment for target H(1,1) (proton)
 
-    targetH = new G4InuclElementaryParticle(targetMomentum, 1);
+    targetH = new G4InuclElementaryParticle(1);
 
     G4float cutElastic[32];
@@ -214 +169 @@
     cutElastic[kaonPlus ]   = 0.5; // 0.5 GeV
     cutElastic[kaonMinus]   = 0.5;
-    cutElastic[kaonMinus]   = 0.5;
     cutElastic[kaonZero]    = 0.5;
     cutElastic[kaonZeroBar] = 0.5;
@@ -222 +176 @@
 
 
-    if (momentumBullet[3] > cutElastic[bulletType]) { // inelastic collision possible
+    if (momentumBullet.z() > cutElastic[bulletType]) { // inelastic collision possible
 
       do {   // we try to create inelastic interaction
-        output = collider->collide(bullet, targetH);
+        output.reset();
+        collider->collide(bullet, targetH, output);
         nTries++;
       } while(
@@ -235 +190 @@
 
     } else { // only elastic collision is energetically possible
-      output = collider->collide(bullet, targetH);
+      collider->collide(bullet, targetH, output);
     }
 
     sumBaryon += 1;
 
-    const G4CascadeMomentum& bmom = bullet->getMomentum();
-    eInit = std::sqrt(bmom[0] * bmom[0]);
-    const G4CascadeMomentum& tmom = targetH->getMomentum();
-    eInit += std::sqrt(tmom[0] * tmom[0]);
+    eInit = bullet->getEnergy() + target->getEnergy();
 
     if (verboseLevel > 2) {
@@ -254 +206 @@
     do  // we try to create inelastic interaction
       {
-        output = collider->collide(bullet, target );
+        output.reset();
+        collider->collide(bullet, target, output);
         nTries++;
       } while (
-               //             (nTries < maxTries)                                                               &&
-               //(output.getOutgoingParticles().size() + output.getNucleiFragments().size() < 2.5) &&
-               //(output.getOutgoingParticles().size()!=0)                                         &&
-               //(output.getOutgoingParticles().begin()->type()==bullet->type())
-               //);
-
                (nTries < maxTries) &&
                output.getOutgoingParticles().size() == 1 &&     // we retry when elastic collision happened
@@ -272 +219 @@
   }
  
-  if (verboseLevel > 1) 
-    {
+  if (verboseLevel > 1) {
       G4cout << " Cascade output: " << G4endl;
       output.printCollisionOutput();
-    }
+  }
   
+  // Rotate event to put Z axis along original projectile direction
+  output.rotateEvent(toLabFrame);
+
   // Convert cascade data to use hadronics interface
-  std::vector<G4InuclNuclei>             nucleiFragments = output.getNucleiFragments();
-  std::vector<G4InuclElementaryParticle> particles =       output.getOutgoingParticles();
+  const std::vector<G4InuclNuclei>& nucleiFragments = output.getNucleiFragments();
+  const std::vector<G4InuclElementaryParticle>& particles = output.getOutgoingParticles();
 
   theResult.SetStatusChange(stopAndKill);
 
+  // Get outcoming particles
+  G4DynamicParticle* cascadeParticle = 0;
   if (!particles.empty()) { 
-    particleIterator ipart;
-    G4int outgoingParticle;
-
-    for (ipart = particles.begin(); ipart != particles.end(); ipart++) {
-      outgoingParticle = ipart->type();
-      const G4CascadeMomentum& mom = ipart->getMomentum();
-      eTot   += std::sqrt(mom[0] * mom[0]);
-
-      G4double ekin = ipart->getKineticEnergy() * GeV;
-      G4ThreeVector aMom(mom[1], mom[2], mom[3]);
-      aMom = aMom.unit();
-
-      if (ipart->baryon() ) {
-        sumBaryon -= 1;
+    particleIterator ipart = particles.begin();
+    for (; ipart != particles.end(); ipart++) {
+      G4int outgoingType = ipart->type();
+
+      eTot += ipart->getEnergy();
+      sumBaryon -= ipart->baryon();
+      sumEnergy -= ipart->getKineticEnergy();
+
+      if (!ipart->valid() || ipart->quasi_deutron()) {
+        G4cerr << " ERROR: G4PreCompoundCascadeInterface::Propagate incompatible"
+               << " particle type " << ipart->type() << G4endl;
+        continue;
       }
 
-      sumEnergy -= ekin / GeV;
-
-      switch(outgoingParticle) {
-
-      case proton: 
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "proton " << counter << " " << aMom << " " << ekin << G4endl;
-#endif
-        cascadeParticle = new G4DynamicParticle(G4Proton::ProtonDefinition(), aMom, ekin);
-        break; 
-
-      case neutron: 
-
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "neutron "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
-#endif
-        cascadeParticle = new G4DynamicParticle(G4Neutron::NeutronDefinition(), aMom, ekin);
-        break;
-
-      case pionPlus: 
-        cascadeParticle = new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), aMom, ekin);
-
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "pionPlus "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
-#endif
-        break;
-
-      case pionMinus:
-        cascadeParticle = new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), aMom, ekin);
-
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "pionMinus "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
-#endif
-        break;
-
-      case pionZero: 
-        cascadeParticle = new G4DynamicParticle(G4PionZero::PionZeroDefinition(), aMom, ekin);
-
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "pionZero "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
-#endif
-        break;
-
-      case photon: 
-        cascadeParticle = new G4DynamicParticle(G4Gamma::Gamma(), aMom, ekin);
-
-#ifdef debug_G4PreCompoundCascadeInterface
-        G4cerr << "photon "<< counter<<" "<<aMom<<" "<<  ekin<<G4endl;
-#endif
-        break;
-
-
-      case kaonPlus:
-        cascadeParticle = new G4DynamicParticle(G4KaonPlus::KaonPlusDefinition(), aMom, ekin);
-        break;
-
-      case kaonMinus:
-        cascadeParticle = new G4DynamicParticle(G4KaonMinus::KaonMinusDefinition(), aMom, ekin);
-        break;
-
-      case kaonZero:
-        if (G4UniformRand() > 0.5) {
-          cascadeParticle = new G4DynamicParticle(G4KaonZeroLong::KaonZeroLongDefinition(), aMom, ekin);
-        } else {
-          cascadeParticle = new G4DynamicParticle(G4KaonZeroShort::KaonZeroShortDefinition(), aMom, ekin);
-        }
-        break;
-
-      case kaonZeroBar:
-        if (G4UniformRand() > 0.5) {
-          cascadeParticle = new G4DynamicParticle(G4KaonZeroLong::KaonZeroLongDefinition(), aMom, ekin);
-        } else {
-          cascadeParticle = new G4DynamicParticle(G4KaonZeroShort::KaonZeroShortDefinition(), aMom, ekin);
-        }
-        break;
-
-      case lambda:
-        cascadeParticle = new G4DynamicParticle(G4Lambda::LambdaDefinition(), aMom, ekin);
-        break;
-
-      case sigmaPlus:
-        cascadeParticle = new G4DynamicParticle(G4SigmaPlus::SigmaPlusDefinition(), aMom, ekin);
-        break;
-
-      case sigmaZero:
-        cascadeParticle = new G4DynamicParticle(G4SigmaZero::SigmaZeroDefinition(), aMom, ekin);
-        break;
-
-      case sigmaMinus:
-        cascadeParticle = new G4DynamicParticle(G4SigmaMinus::SigmaMinusDefinition(), aMom, ekin);
-        break;
-
-      case xiZero:
-        cascadeParticle = new G4DynamicParticle(G4XiZero::XiZeroDefinition(), aMom, ekin);
-        break;
-
-      case xiMinus:
-        cascadeParticle = new G4DynamicParticle(G4XiMinus::XiMinusDefinition(), aMom, ekin);
-        break;
-
-      default:
-        G4cout << " ERROR: G4PreCompoundCascadeInterface::Propagate undefined particle type" << G4endl;
+      // Copy local G4DynPart to public output (handle kaon mixing specially)
+      if (outgoingType == kaonZero || outgoingType == kaonZeroBar) {
+        G4ThreeVector momDir = ipart->getMomentum().vect().unit();
+        G4double ekin = ipart->getKineticEnergy();
+
+        G4ParticleDefinition* pd = G4KaonZeroShort::Definition();
+        if (G4UniformRand() > 0.5) pd = G4KaonZeroLong::Definition();
+
+        cascadeParticle = new G4DynamicParticle(pd, momDir, ekin);
+      } else {
+        cascadeParticle = new G4DynamicParticle(ipart->getDynamicParticle());
       }
-
-      cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
+ 
       theResult.AddSecondary(cascadeParticle); 
     }
@@ -412 +269 @@
   // get nuclei fragments
   G4DynamicParticle * aFragment = 0;
-  G4ParticleDefinition * aIonDef = 0;
-  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
-
   if (!nucleiFragments.empty()) { 
-    nucleiIterator ifrag;
-
-    for (ifrag = nucleiFragments.begin(); ifrag != nucleiFragments.end(); ifrag++) 
-      {
-        G4double eKin = ifrag->getKineticEnergy() * GeV;
-        const G4CascadeMomentum& mom = ifrag->getMomentum();
-        eTot   += std::sqrt(mom[0] * mom[0]);
-
-        G4ThreeVector aMom(mom[1], mom[2], mom[3]);
-        aMom = aMom.unit();
-
-        // hpw @@@ ==> Should be zero: G4double fragmentExitation = ifrag->getExitationEnergyInGeV();
-
-        if (verboseLevel > 2) {
-          G4cout << " Nuclei fragment: " << G4endl;
-          ifrag->printParticle();
-        }
-
-        G4int A = G4int(ifrag->getA());
-        G4int Z = G4int(ifrag->getZ());
-        aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
-      
-        aFragment =  new G4DynamicParticle(aIonDef, aMom, eKin);
-
-        sumBaryon -= A;
-        sumEnergy -= eKin / GeV;
-
-        aFragment->Set4Momentum(aFragment->Get4Momentum()*=toLabFrame);
-        theResult.AddSecondary(aFragment); 
+    nucleiIterator ifrag = nucleiFragments.begin();
+    for (; ifrag != nucleiFragments.end(); ifrag++) {
+      eTot += ifrag->getEnergy();
+      sumBaryon -= ifrag->getA();
+      sumEnergy -= ifrag->getKineticEnergy();
+
+      if (verboseLevel > 2) {
+        G4cout << " Nuclei fragment: " << G4endl;
+        ifrag->printParticle();
       }
-  }
-
+
+      // Copy local G4DynPart to public output 
+      aFragment =  new G4DynamicParticle(ifrag->getDynamicParticle());
+      theResult.AddSecondary(aFragment); 
+    }
+  }
+
+  // Report violations of energy, baryon conservation
   if (verboseLevel > 2) {
     if (sumBaryon != 0) {
-      G4cout << "ERROR: no baryon number conservation, sum of baryons = " << sumBaryon << G4endl;
+      G4cout << "ERROR: no baryon number conservation, sum of baryons = "
+             << sumBaryon << G4endl;
     }
 
     if (sumEnergy > 0.01 ) {
-      G4cout << "Kinetic energy conservation violated by " << sumEnergy << " GeV" << G4endl;
+      G4cout << "Kinetic energy conservation violated by "
+             << sumEnergy << " GeV" << G4endl;
     }
      
-    G4cout << "Total energy conservation at level ~" << (eInit - eTot) * GeV << " MeV" << G4endl;
+    G4cout << "Total energy conservation at level ~"
+           << (eInit - eTot) * GeV << " MeV" << G4endl;
     
     if (sumEnergy < -5.0e-5 ) { // 0.05 MeV
-      G4cout << "FATAL ERROR: energy created  " << sumEnergy * GeV << " MeV" << G4endl;
+      G4cout << "FATAL ERROR: energy created  "
+             << sumEnergy * GeV << " MeV" << G4endl;
     }
   }
 
   delete bullet;
-  delete colep;
-  delete inc;
-  delete noneq; 
-  delete bigb;
   delete collider;
 
   if(target != 0) delete target;
   if(targetH != 0) delete targetH;
-  // if(cascadeParticle != 0) delete cascadeParticle;
-  // if(aFragment != 0) delete aFragment;
 
   return &theResult;