Changeset 1315 for trunk/source/processes/hadronic/models/cascade/cascade/src/G4InuclEvaporation.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/G4InuclEvaporation.cc (modified) (8 diffs)

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

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: G4InuclEvaporation.cc,v 1.18 2010/05/21 18:26:16 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
 //
-// $Id: G4InuclEvaporation.cc,v 1.8 2008/10/24 20:49:07 dennis Exp $
-//
+// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
+// 20100405  M. Kelsey -- Pass const-ref std::vector<>
+// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide(), use
+//              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.  Make
+//              G4EvaporationInuclCollider a data member.
+// 20100520  M. Kelsey -- Simplify collision loop, move momentum rotations to
+//              G4CollisionOutput, copy G4DynamicParticle directly from
+//              G4InuclParticle, no switch-block required.  Fix scaling factors.
+
 #include <numeric>
 #include "G4IonTable.hh"
@@ -41 +53 @@
 #include "G4LorentzVector.hh"
 #include "G4EquilibriumEvaporator.hh"
-#include "G4Fissioner.hh"
-#include "G4BigBanger.hh"
 #include "G4InuclElementaryParticle.hh"
 #include "G4InuclParticle.hh"
 #include "G4CollisionOutput.hh"
+#include "G4InuclParticleNames.hh"
 
-typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
-typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
+using namespace G4InuclParticleNames;
+
+typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
+typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
 
 
-G4InuclEvaporation::G4InuclEvaporation() {
-  verboseLevel=0;
-}
+G4InuclEvaporation::G4InuclEvaporation() 
+  : verboseLevel(0),   evaporator(new G4EvaporationInuclCollider) {}
 
 G4InuclEvaporation::G4InuclEvaporation(const G4InuclEvaporation &) : G4VEvaporation() {
@@ -59 +71 @@
 }
 
-
 G4InuclEvaporation::~G4InuclEvaporation() {
+  delete evaporator;
 }
 
@@ -81 +93 @@
 }
 
-G4FragmentVector * G4InuclEvaporation::BreakItUp(const G4Fragment &theNucleus) {
- 
-  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 };  
-
-  std::vector< G4DynamicParticle * > secondaryParticleVector;
-  G4FragmentVector * theResult = new G4FragmentVector;
+G4FragmentVector* G4InuclEvaporation::BreakItUp(const G4Fragment &theNucleus) {
+  G4FragmentVector* theResult = new G4FragmentVector;
 
   if (theNucleus.GetExcitationEnergy() <= 0.0) { // Check that Excitation Energy > 0
@@ -100 +104 @@
   G4double Z = theNucleus.GetZ();
   G4double mTar  = G4NucleiProperties::GetNuclearMass(A, Z); // Mass of the target nucleus
-  G4LorentzVector tmp =theNucleus.GetMomentum();
 
-  G4ThreeVector momentum = tmp.vect();
+  G4ThreeVector momentum =  theNucleus.GetMomentum().vect() / GeV;
   //  G4double energy = tmp.e();
   G4double exitationE = theNucleus.GetExcitationEnergy();
@@ -108 +111 @@
   // Move to CMS frame, save initial velocity of the nucleus to boostToLab vector.
   //   G4ThreeVector boostToLab( ( 1/G4NucleiProperties::GetNuclearMass( A, Z ) ) * momentum ); 
-  G4InuclNuclei* tempNuc = new G4InuclNuclei(A, Z);
-  G4double mass=tempNuc->getMass()*1000;
-  G4ThreeVector boostToLab( ( 1/mass) * momentum ); 
+
+  G4double mass = mTar / GeV;
+  G4ThreeVector boostToLab( momentum/mass ); 
 
   if ( verboseLevel > 2 )
@@ -122 +125 @@
 
   G4InuclNuclei* nucleus = new G4InuclNuclei(A, Z);
-  nucleus->setExitationEnergy(exitationE/1000);
-  G4CascadeMomentum tmom;
-  nucleus->setMomentum(tmom);
-  nucleus->setEnergy();
+  nucleus->setExitationEnergy(exitationE);
 
-  G4EquilibriumEvaporator*          eqil = new G4EquilibriumEvaporator;
-  G4Fissioner*                      fiss = new G4Fissioner;
-  G4BigBanger*                      bigb = new G4BigBanger;
-  G4EvaporationInuclCollider* evaporator = new G4EvaporationInuclCollider(eqil, fiss, bigb);
   G4CollisionOutput output;
+  evaporator->collide(0, nucleus, output);
 
-  output = evaporator->collide(0, nucleus);
+  const std::vector<G4InuclNuclei>& nucleiFragments = output.getNucleiFragments();
+  const std::vector<G4InuclElementaryParticle>& particles = output.getOutgoingParticles();
 
-  std::vector<G4InuclNuclei>             nucleiFragments = output.getNucleiFragments();
-  std::vector<G4InuclElementaryParticle> particles =       output.getOutgoingParticles();
+  G4double eTot=0.0;
+  G4int  i=1;
 
-  G4double ekin,emas;
-  G4double eTot=0.0;
-  G4DynamicParticle* cascadeParticle = 0;
-  G4int  i=1;
-  //G4cout << "# particles: " << output.getOutgoingParticles().size() << G4endl;
   if (!particles.empty()) { 
-    particleIterator ipart;
-    G4int outgoingParticle;
-
-    for (ipart = particles.begin(); ipart != particles.end(); ipart++) {
-     
-      outgoingParticle = ipart->type();
+    G4int outgoingType;
+    particleIterator ipart = particles.begin();
+    for (; ipart != particles.end(); ipart++) {
+      outgoingType = ipart->type();
 
       if (verboseLevel > 2) {
-        G4cout << "Evaporated particle:  " << i << " of type: " << outgoingParticle << G4endl;
+        G4cout << "Evaporated particle:  " << i << " of type: " << outgoingType << G4endl;
         i++;
         //              ipart->printParticle();
       }
 
-      const G4CascadeMomentum& mom = ipart->getMomentum();
-      eTot   += std::sqrt(mom[0]*1000 * mom[0]*1000);
+      eTot += ipart->getEnergy();
+      
+      G4LorentzVector vlab = ipart->getMomentum().boost(boostToLab);
 
-      ekin = ipart->getKineticEnergy()*1000;
-      emas = ipart->getMass()*1000;
-
-      G4ThreeVector aMom(mom[1]*1000, mom[2]*1000, mom[3]*1000);
-
-      G4LorentzVector v(aMom, (ekin+emas));
-      v.boost( boostToLab );
-
-      switch(outgoingParticle) {
-
-      case proton: 
-        cascadeParticle = new G4DynamicParticle(G4Proton::ProtonDefinition(), v.vect(), v.e());
-        break; 
-
-      case neutron: 
-        cascadeParticle = new G4DynamicParticle(G4Neutron::NeutronDefinition(), v.vect(), v.e());
-        break;
-
-      case photon: 
-        cascadeParticle = new G4DynamicParticle(G4Gamma::Gamma(), v.vect(), v.e());
-        break;
-      default:
-        G4cout << " ERROR: GInuclEvapration::Propagate undefined particle type" << G4endl;
-      };
-
-      secondaryParticleVector.push_back( cascadeParticle );
+      theResult->push_back( new G4Fragment(vlab, ipart->getDefinition()) );
       //      theResult.AddSecondary(cascadeParticle); 
     }  
   }
-  fillResult( secondaryParticleVector, theResult);
-  
 
   //  G4cout << "# fragments " << output.getNucleiFragments().size() << G4endl;
   i=1; 
   if (!nucleiFragments.empty()) { 
-    nucleiIterator ifrag;
-
-    for (ifrag = nucleiFragments.begin(); ifrag != nucleiFragments.end(); ifrag++) {
-
-      ekin = ifrag->getKineticEnergy()*1000;
-      emas = ifrag->getMass()*1000;
-      const G4CascadeMomentum& mom = ifrag->getMomentum();
-      eTot  += std::sqrt(mom[0]*1000 * mom[0]*1000);
-
-      G4ThreeVector aMom(mom[1]*1000, mom[2]*1000, mom[3]*1000);
-      //      aMom = aMom.unit();
-
-      G4LorentzVector v(aMom, (ekin+emas));
-      v.boost( boostToLab );
-
+    nucleiIterator ifrag = nucleiFragments.begin();
+    for (; ifrag != nucleiFragments.end(); ifrag++) {
       if (verboseLevel > 2) {
         G4cout << " Nuclei fragment: " << i << G4endl; i++;
-        //      ifrag->printParticle();
       }
 
+      eTot += ifrag->getEnergy();
+
+      G4LorentzVector vlab = ifrag->getMomentum().boost(boostToLab);
+ 
       G4int A = G4int(ifrag->getA());
       G4int Z = G4int(ifrag->getZ());
-      //        cascadeParticle = new G4DynamicParticle(G4Proton::ProtonDefinition(), v.vect(), v.e());
       if (verboseLevel > 2) {
-        G4cout << "boosted v" << v << G4endl;
+        G4cout << "boosted v" << vlab << G4endl;
       }
-      // theResult->push_back( new G4Fragment(A, Z, tmp) ); 
-      theResult->push_back( new G4Fragment(A, Z, v) ); 
+      theResult->push_back( new G4Fragment(A, Z, vlab) ); 
     }
   }
@@ -229 +182 @@
   return theResult;
 }
-
-void G4InuclEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
-                                     G4FragmentVector * aResult )
-{
-  // Fill the vector pParticleChange with secondary particles stored in vector.
-  for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ )
-    {
-      G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
-      G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
-      G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
-      if(aA>0) {
-        aResult->push_back( new G4Fragment(aA, aZ, aMomentum) ); 
-      } else {
-        aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) ); 
-      }
-    }
-  return;
-}
-
-
-