Changeset 1347 for trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc

Timestamp:

Dec 22, 2010, 3:52:27 PM (13 years ago)

Author:

garnier

Message:

geant4 tag 9.4

File:

• trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc (modified) (22 diffs)

trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclCascadeInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclCascadeInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -34 +34 @@
 
 #include "G4InclCascadeInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
 #include "CLHEP/Random/Random.h"
 
-G4InclCascadeInterface::G4InclCascadeInterface()
+G4InclCascadeInterface::G4InclCascadeInterface(const G4String& nam)
+  :G4VIntraNuclearTransportModel(nam)
 {
   hazard = new G4Hazard();
@@ -45 +47 @@
 
   varntp = new G4VarNtp();
-  inclInput = 0;
+  calincl = 0;
   ws = new G4Ws();
   mat = new G4Mat();
-  incl = new G4Incl(hazard, inclInput, ws, mat, varntp);
+  incl = new G4Incl(hazard, calincl, ws, mat, varntp);
+
+  theExcitationHandler = new G4ExcitationHandler;
+  thePrecoModel = new G4PreCompoundModel(theExcitationHandler);
+
+  if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
+    incl->setUseFermiBreakUp(true);
+  }
 
   verboseLevel = 0;
@@ -55 +64 @@
 G4InclCascadeInterface::~G4InclCascadeInterface()
 {
+  delete thePrecoModel;
+  delete theExcitationHandler;
+
   delete hazard;
   delete varntp;
-  delete inclInput;
   delete ws;
   delete mat;
@@ -68 +79 @@
 
   G4int particleI;
+  G4int bulletType = 0;
 
   // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
@@ -91 +103 @@
   G4DynamicParticle *cascadeParticle = 0;
   G4ParticleDefinition *aParticleDefinition = 0;
+  
+  G4ReactionProductVector *thePrecoResult = 0;
+  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
   // INCL assumes the projectile particle is going in the direction of
@@ -102 +117 @@
   G4LorentzRotation toLabFrame = toZ.inverse();
 
-  inclInput = new G4InclInput(aTrack, theNucleus, false);
-
   theResult.Clear(); // Make sure the output data structure is clean.
+
+  calincl = new G4InclInput(aTrack, theNucleus, false);
+  incl->setInput(calincl);
+
+  //  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
+  //  calincl->printInfo();
 
 #ifdef DEBUGINCL
@@ -117 +136 @@
   G4double eKinSum = bulletE;
   G4LorentzVector labv = G4LorentzVector(0.0, 0.0, std::sqrt(bulletE*(bulletE + 2.*mass)), bulletE + mass + amass);
+  G4LorentzVector labvA = G4LorentzVector(0.0, 0.0, 0.0, 0.0);
   G4cout <<"Energy in the beginning = " << labv.e() / MeV << G4endl;
 #endif
 
   // Check wheter the input is acceptable.
-  if((inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) {
+  if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -130 +150 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(inclInput->targetA());
-    mat->zmat[0] = int(inclInput->targetZ());
+    mat->amat[0] = int(calincl->targetA());
+    mat->zmat[0] = int(calincl->targetZ());
 
     incl->initIncl(true);
@@ -140 +160 @@
         G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventIncl(inclInput);
+      incl->processEventIncl(calincl);
 
       if(verboseLevel > 1) {
@@ -151 +171 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;
-
-      G4cout <<"G4InclCascadeInterface: Target A:  " << inclInput->targetA() << G4endl;
-      G4cout <<"G4InclCascadeInterface: Target Z:  " << inclInput->targetZ() << G4endl;
+      G4cout <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
+
+      G4cout <<"G4InclCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+      G4cout <<"G4InclCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;
-        diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclCascadeInterface: Target A:  " << inclInput->targetA() << G4endl;
-        diagdata <<"G4InclCascadeInterface: Target Z:  " << inclInput->targetZ() << G4endl;
-      }
-
-      // for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-      //   G4cout << n                       << " " << inclInput->f[6]             << " " << inclInput->f[2] << " ";
-      //   G4cout << varntp->massini         << " " << varntp->mzini             << " ";
-      //   G4cout << varntp->exini           << " " << varntp->mulncasc          << " " << varntp->mulnevap          << " " << varntp->mulntot << " ";
-      //   G4cout << varntp->bimpact         << " " << varntp->jremn             << " " << varntp->kfis              << " " << varntp->estfis << " ";
-      //   G4cout << varntp->izfis           << " " << varntp->iafis             << " " << varntp->ntrack            << " " << varntp->itypcasc[particleI] << " ";
-      //   G4cout << varntp->avv[particleI]  << " " << varntp->zvv[particleI]    << " " << varntp->enerj[particleI]  << " ";
-      //   G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
-      //   // For diagnostic output
-      //   if(verboseLevel > 3) {
-      //     diagdata << n                       << " " << inclInput->f[6]             << " " << inclInput->f[2] << " ";
-      //     diagdata << varntp->massini         << " " << varntp->mzini             << " ";
-      //     diagdata << varntp->exini           << " " << varntp->mulncasc          << " " << varntp->mulnevap          << " " << varntp->mulntot << " ";
-      //     diagdata << varntp->bimpact         << " " << varntp->jremn             << " " << varntp->kfis              << " " << varntp->estfis << " ";
-      //     diagdata << varntp->izfis           << " " << varntp->iafis             << " " << varntp->ntrack            << " ";
-      //          diagdata                                                                       << varntp->itypcasc[particleI] << " ";
-      //     diagdata << varntp->avv[particleI]  << " " << varntp->zvv[particleI]    << " " << varntp->enerj[particleI]  << " ";
-      //     diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
-      //   }
-      // }
+        diagdata <<"G4InclCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+        diagdata <<"G4InclCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
+      }
+
     }
 
@@ -196 +196 @@
 
       theResult.SetStatusChange(stopAndKill);
-      aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType());
-
-      cascadeParticle = new G4DynamicParticle();
-      cascadeParticle->SetDefinition(aParticleDefinition);
-      cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
-      theResult.AddSecondary(cascadeParticle);
+      
+      G4int bulletType = calincl->bulletType();
+      aParticleDefinition = G4InclInput::getParticleDefinition(bulletType);
+
+      if(aParticleDefinition != 0) {
+        cascadeParticle = new G4DynamicParticle();
+        cascadeParticle->SetDefinition(aParticleDefinition);
+        cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
+        theResult.AddSecondary(cascadeParticle);
+      }
     }
 
@@ -209 +213 @@
 
 #ifdef DEBUGINCL
-    G4cout << "E [MeV]" << std::setw(12) << " Ekin [MeV]" << std::setw(12) << " E* [MeV]" << std::setw(12) << "Px [MeV]" << std::setw(12) << " Py [MeV]" << std::setw(12) << "Pz [MeV]" << std::setw(12) << "Pt [MeV]" << std::setw(12) << "A" << std::setw(12) << "Z" << G4endl;  
+    G4cout << "E [MeV]" << std::setw(12)
+           << " Ekin [MeV]" << std::setw(12)
+           << "Px [MeV]" << std::setw(12)
+           << " Py [MeV]" << std::setw(12)
+           << "Pz [MeV]" << std::setw(12)
+           << "Pt [MeV]" << std::setw(12)
+           << "A" << std::setw(12)
+           << "Z" << G4endl;
 #endif
 
@@ -276 +287 @@
       if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
         G4ParticleDefinition * aIonDef = 0;
-        G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
         G4int A = G4int(varntp->avv[particleI]);
@@ -315 +325 @@
         G4double p = mom.mag();
         labv -= fm;
-        G4double px = mom.x() * MeV;
-        G4double py = mom.y() * MeV;
-        G4double pz = mom.z() * MeV;
+        if(varntp->avv[particleI] > 1) {
+          labvA += fm;
+        }
+        G4double px = mom.x() * MeV;
+        G4double py = mom.y() * MeV;
+        G4double pz = mom.z() * MeV;
         G4double pt = std::sqrt(px*px+py*py);
         G4double e  = fm.e();
@@ -330 +343 @@
         G4cout << fm.e() / MeV
                << std::setw(12) << cascadeParticle->GetKineticEnergy() / MeV
-               << std::setw(12) << exE / MeV
                << std::setw(12) << mom.x() / MeV
                << std::setw(12) << mom.y() / MeV
@@ -351 +363 @@
       }
     }
+
+      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV;
+      G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV,
+                                       varntp->erecrem * MeV + nuclearMass);
+      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini),
+                                                                        varntp->exini,
+                                                                        varntp->erecrem,
+                                                                        varntp->pxrem,
+                                                                        varntp->pyrem,
+                                                                        varntp->pzrem);
+      G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV,
+                         momentumScaling * varntp->pzrem * MeV,
+                         varntp->erecrem + nuclearMass);
+
+      // For four-momentum, baryon number and charge conservation check:
+      G4LorentzVector fourMomentumBalance = p4;
+      G4int baryonNumberBalance = G4int(varntp->massini);
+      G4int chargeBalance = G4int(varntp->mzini);
+
+      G4LorentzRotation toFragmentZ;
+      toFragmentZ.rotateZ(-p4.theta());
+      toFragmentZ.rotateY(-p4.phi());
+      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
+      p4 *= toFragmentZ;
+
+      G4LorentzVector p4rest = p4;
+      p4rest.boost(-p4.boostVector());
+      if(verboseLevel > 0) {
+        G4cout <<"Cascade remnant nucleus:" << G4endl;
+        G4cout <<"p4: " << G4endl;
+        G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
+        G4cout <<" E = " << p4.e() << G4endl;
+
+        G4cout <<"p4rest: " << G4endl;
+        G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
+        G4cout <<" E = " << p4rest.e() << G4endl;
+      }
+
+      G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest);
+      thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant);
+      if(thePrecoResult != 0) {
+      G4ReactionProductVector::iterator fragment;
+      for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) {
+        G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition();
+
+        if(theFragmentDefinition != 0) {
+          G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
+          G4LorentzVector labMomentum = theFragment->Get4Momentum();
+          labMomentum.boost(p4.boostVector());
+          labMomentum *= toFragmentLab;
+          labMomentum *= toLabFrame;
+          theFragment->Set4Momentum(labMomentum);
+          fourMomentumBalance -= theFragment->Get4Momentum();
+          baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
+          chargeBalance -= theFragmentDefinition->GetAtomicNumber();
+          if(verboseLevel > 0) {
+            G4cout <<"Resulting fragment: " << G4endl;
+            G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
+            G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
+          }
+          theResult.AddSecondary(theFragment);
+        } else {
+          G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+          G4cout <<"Resulting fragment: " << G4endl;
+          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+        }
+      }
+      delete thePrecoResult;
+      thePrecoResult = 0;
+
+      if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
+        G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
+        G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
+        G4cout <<"Vector components (px, py, pz, E) = ("
+               << fourMomentumBalance.px() << ", "
+               << fourMomentumBalance.py() << ", "
+               << fourMomentumBalance.pz() << ", "
+               << fourMomentumBalance.e() << ")" << G4endl;
+      }
+      if(baryonNumberBalance != 0 && verboseLevel > 1) {
+        G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0 && verboseLevel > 1) {
+        G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+      }
+      //    } // if(needsFermiBreakUp)
+
 #ifdef DEBUGINCL
     G4cout <<"--------------------------------------------------------------------------------" << G4endl;
     G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2));
-    G4cout << labv.e() / MeV << std::setw(12) << eKinSum / MeV << std::setw(12) << labv.x() << std::setw(12) << labv.y() << std::setw(12) << labv.z() << std::setw(12) <<  pt / MeV << std::setw(12) << baryonNumber << std::setw(12) << chargeNumber << " totals" << G4endl;
+    G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2));
+    G4cout << labv.e() / MeV << std::setw(12)
+           << eKinSum  / MeV << std::setw(12)
+           << labv.x() / MeV << std::setw(12)
+           << labv.y() / MeV << std::setw(12)
+           << labv.z() / MeV << std::setw(12)
+           << pt       / MeV << std::setw(12)
+           << baryonNumber << std::setw(12)
+           << chargeNumber << " totals" << G4endl;
+    G4cout << " - " << std::setw(12)
+           << " - " << std::setw(12)
+           << labvA.x() / MeV << std::setw(12)
+           << labvA.y() / MeV << std::setw(12)
+           << labvA.z() / MeV << std::setw(12)
+           << ptA       / MeV << std::setw(12)
+           << " - " << std::setw(12) << " - " << " totals ABLA" << G4endl;
     G4cout << G4endl;
-
+ 
     if(verboseLevel > 3) {
       if(baryonNumber != 0) {
@@ -369 +484 @@
     }
 #endif
-    
+   
     varntp->ntrack = 0; // Clean up the number of generated particles in the event.
   }
@@ -379 +494 @@
     theResult.SetStatusChange(stopAndKill);
 
-    G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
@@ -388 +502 @@
     }
     if(verboseLevel > 3) {
-      diagdata <<"ERROR G4InclCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;
-    }
-
-    if(inclInput->bulletType() == 0) {
+      diagdata <<"ERROR G4InclCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl;
+    }
+
+    if(bulletType == 0) {
       if(verboseLevel > 1) {
         G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl;      
@@ -402 +516 @@
     }
 
-    if((inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target
+    if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
       if(verboseLevel > 1) {
         G4cout <<"Unsupported target: " << G4endl;
-        G4cout <<"Target A: " << inclInput->targetA() << G4endl;
-        G4cout <<"TargetZ: " << inclInput->targetZ() << G4endl;
+        G4cout <<"Target A: " << calincl->targetA() << G4endl;
+        G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
       }
       if(verboseLevel > 3) {
         diagdata <<"Unsupported target: " << G4endl;
-        diagdata <<"Target A: " << inclInput->targetA() << G4endl;
-        diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl;
-      }
-    }
-
-    if(inclInput->bulletE() < 100) { // INCL does not support E < 100 MeV.
+        diagdata <<"Target A: " << calincl->targetA() << G4endl;
+        diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
+      }
+    }
+
+    if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
       if(verboseLevel > 1) {
-        G4cout <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
+        G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
         G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
       }
       if(verboseLevel > 3) {
-        diagdata <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
+        diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
       }
     }
@@ -430 +544 @@
   }
 
+  delete calincl;
+  calincl = 0;
   return &theResult;
 }