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

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

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: G4IntraNucleiCascader.cc,v 1.35 2010/05/21 17:56:34 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
 //
-
+// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
+// 20100307  M. Kelsey -- Bug fix: momentum_out[0] should be momentum_out.e()
+// 20100309  M. Kelsey -- Eliminate some unnecessary std::pow()
+// 20100407  M. Kelsey -- Pass "all_particles" as argument to initializeCascad,
+//              following recent change to G4NucleiModel.
+// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
+// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
+//              simple data members
+
+#include "G4IntraNucleiCascader.hh"
 #define RUN
 
-#include "G4IntraNucleiCascader.hh"
+#include "G4CascadParticle.hh"
+#include "G4CollisionOutput.hh"
+#include "G4ElementaryParticleCollider.hh"
+#include "G4HadTmpUtil.hh"
 #include "G4InuclElementaryParticle.hh"
 #include "G4InuclNuclei.hh"
-#include "G4LorentzConvertor.hh"
+#include "G4InuclSpecialFunctions.hh"
+#include "G4NucleiModel.hh"
+#include "G4NucleiProperties.hh"
 #include "G4ParticleLargerEkin.hh"
-#include "G4NucleiModel.hh"
-#include "G4CascadParticle.hh"
 #include "Randomize.hh"
 #include <algorithm>
 
+using namespace G4InuclSpecialFunctions;
+
+
 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
 
 G4IntraNucleiCascader::G4IntraNucleiCascader()
-  : verboseLevel(1) {
+  : G4VCascadeCollider("G4IntraNucleiCascader"),
+    theElementaryParticleCollider(new G4ElementaryParticleCollider) {}
+
+G4IntraNucleiCascader::~G4IntraNucleiCascader() {
+  delete theElementaryParticleCollider;
+}
+
+
+void G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
+                                    G4InuclParticle* target,
+                                    G4CollisionOutput& output) {
 
   if (verboseLevel > 3) {
-    G4cout << " >>> G4IntraNucleiCascader::G4IntraNucleiCascader" << G4endl;
-  }
-
-}
-
-G4CollisionOutput G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
-                                                 G4InuclParticle* target) {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4IntraNucleiCascader::collide" << G4endl;
+    G4cout << " >>> G4IntraNucleiCascader::collide inter_case " << inter_case 
+           << G4endl;
   }
 
@@ -63 +82 @@
     target->printParticle();
   }
-
-  G4CollisionOutput output;
 
 #ifdef RUN
@@ -73 +90 @@
   G4NucleiModel model(tnuclei);
   G4double coulombBarrier = 0.00126*tnuclei->getZ()/
-                                      (1.+std::pow(tnuclei->getA(),0.333));
-
-  G4CascadeMomentum momentum_in = bullet->getMomentum();
-
-  momentum_in[0] += tnuclei->getMass();
+                                      (1.+G4cbrt(tnuclei->getA()));
+
+  G4LorentzVector momentum_in = bullet->getMomentum();
+
+  momentum_in.setE(momentum_in.e()+tnuclei->getMass());
 
   G4double ekin_in; 
@@ -83 +100 @@
   if (verboseLevel > 3) {
     model.printModel();
-    G4cout << " intitial momentum  E " << momentum_in[0] << " Px " << momentum_in[1] 
-           << " Py " << momentum_in[2] << " Pz " << momentum_in[3] << G4endl;
+    G4cout << " intitial momentum  E " << momentum_in.e() << " Px "
+           << momentum_in.x() << " Py " << momentum_in.y() << " Pz "
+           << momentum_in.z() << G4endl;
   }
 
@@ -117 +135 @@
       zfin += zb;
 
-      std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > 
-        all_particles = model.initializeCascad(bnuclei, tnuclei);
+      G4NucleiModel::modelLists all_particles;    // Buffer to receive lists
+      model.initializeCascad(bnuclei, tnuclei, all_particles);
 
       cascad_particles = all_particles.first;
@@ -139 +157 @@
 
         for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
-
         for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
       };
@@ -208 +225 @@
                                             (1./KE - 1./coulombBarrier)*
                                          std::sqrt(mass*(coulombBarrier-KE)) );
+
             if (G4UniformRand() < CBP) {
               output_particles.push_back(currentParticle);
@@ -241 +259 @@
     }
 
-    G4CascadeMomentum momentum_out;
+    G4LorentzVector momentum_out;
     particleIterator ipart;
 
     for (ipart = output_particles.begin(); ipart != output_particles.end(); ipart++) {
-      const G4CascadeMomentum& mom = ipart->getMomentum();
-
-      for (G4int j = 0; j < 4; j++) momentum_out[j] += mom[j];
+      momentum_out += ipart->getMomentum();
 
       zfin -= ipart->getCharge();
@@ -260 +276 @@
       G4InuclNuclei outgoing_nuclei(afin, zfin);
       G4double mass = outgoing_nuclei.getMass();
-      momentum_out[0] += mass;        
-
-      for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
+      momentum_out.setE(momentum_out.e()+mass);
+
+      momentum_out = momentum_in - momentum_out;
 
       if (verboseLevel > 3) {
-        G4cout << "  Eex + Ekin " << momentum_out[0]  <<  G4endl;
+        G4cout << "  Eex + Ekin " << momentum_out.e()  <<  G4endl;
       }
 
-      if (momentum_out[0] > 0.0) { // Eex + Ekin > 0.0
-        G4double pnuc = momentum_out[1] * momentum_out[1] + 
-          momentum_out[2] * momentum_out[2] +
-          momentum_out[3] * momentum_out[3]; 
+      if (momentum_out.e() > 0.0) { // Eex + Ekin > 0.0
+        G4double pnuc = momentum_out.vect().mag2(); 
         G4double ekin = std::sqrt(mass * mass + pnuc) - mass;
-        G4double Eex = 1000.0 * (momentum_out[0] - ekin);
+        G4double Eex = 1000.0 * (momentum_out.e() - ekin);
 
         if (verboseLevel > 3) {
@@ -283 +297 @@
           output.addOutgoingParticles(output_particles);
           outgoing_nuclei.setMomentum(momentum_out);
-          outgoing_nuclei.setEnergy();
           outgoing_nuclei.setExitationEnergy(Eex);
           outgoing_nuclei.setExitonConfiguration(theExitonConfiguration);                                         
           output.addTargetFragment(outgoing_nuclei);
 
-          return output;
+          return;
         };
       };
@@ -296 +309 @@
       if (afin == 1.0) { // recoiling nucleon
 
-        for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
+        momentum_out = momentum_in - momentum_out;
 
         G4InuclElementaryParticle  last_particle;
@@ -314 +327 @@
       output.addOutgoingParticles(output_particles);
 
-      return output;
+      return;
     }; 
   };
@@ -322 +335 @@
   // special branch to avoid the cascad generation but to get the input for evaporation etc
 
-  G4CascadeMomentum momentum_out;
+  G4LorentzVector momentum_out;
   G4InuclNuclei outgoing_nuclei(169, 69);
 
   outgoing_nuclei.setMomentum(momentum_out);
-  outgoing_nuclei.setEnergy();
   outgoing_nuclei.setExitationEnergy(150.0);
 
@@ -334 +346 @@
   output.addTargetFragment(outgoing_nuclei);
 
-  return output;
+  return;
 
   /*
@@ -342 +354 @@
     output.addOutgoingParticle(*bparticle);
     output.addTargetFragment(*tnuclei);
-    return output;
+    return;
   */
 
@@ -354 +366 @@
   output.trivialise(bullet, target);
 
-  return output;
+  return;
 }
 
@@ -374 +386 @@
   if (eexs > eexs_cut) {
     G4double eexs_max0z = 1000.0 * ein / ediv_cut;
-    G4double dm = bindingEnergy(a, z);
+    //    G4double dm = bindingEnergy(a, z);
+    G4double dm = G4NucleiProperties::GetBindingEnergy(G4lrint(a), G4lrint(z));
     G4double eexs_max = eexs_max0z > reason_cut*dm ? eexs_max0z : reason_cut * dm;