source: trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc

Last change on this file was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 30.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// $Id: G4InclAblaLightIonInterface.cc,v 1.16 2010/11/17 20:19:09 kaitanie Exp $
27// Translation of INCL4.2/ABLA V3
28// Pekka Kaitaniemi, HIP (translation)
29// Christelle Schmidt, IPNL (fission code)
30// Alain Boudard, CEA (contact person INCL/ABLA)
31// Aatos Heikkinen, HIP (project coordination)
32
33#include <vector>
34
35#include "G4InclAblaLightIonInterface.hh"
36#include "G4FermiBreakUp.hh"
37#include "math.h"
38#include "G4GenericIon.hh"
39#include "CLHEP/Random/Random.h"
40
41G4InclAblaLightIonInterface::G4InclAblaLightIonInterface()
42{
43  hazard = new G4Hazard();
44
45  const G4long* table_entry = CLHEP::HepRandom::getTheSeeds(); // Get random seed from CLHEP.
46  hazard->ial = (*table_entry);
47
48  varntp = new G4VarNtp();
49  calincl = 0;
50  ws = new G4Ws();
51  mat = new G4Mat();
52  incl = new G4Incl(hazard, calincl, ws, mat, varntp);
53  useProjectileSpectator = true;
54  useFermiBreakup = true;
55  incl->setUseProjectileSpectators(useProjectileSpectator);
56  if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
57    incl->setUseFermiBreakUp(true);
58    useFermiBreakup = true;
59  }
60  verboseLevel = 0;
61  if(getenv("G4INCLVERBOSE")) {
62    verboseLevel = 1;
63  }
64}
65
66G4InclAblaLightIonInterface::~G4InclAblaLightIonInterface()
67{
68  delete hazard;
69  delete varntp;
70  delete calincl;
71  delete ws;
72  delete mat;
73  delete incl;
74}
75
76G4HadFinalState* G4InclAblaLightIonInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
77{
78  //  const G4bool useFermiBreakup = false;
79  G4int maxTries = 200;
80
81  G4int particleI;
82
83  G4int baryonNumberBalanceInINCL = 0;
84  G4int chargeNumberBalanceInINCL = 0;
85
86  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
87
88  // Increase the event number:
89  eventNumber++;
90
91  // Clean up the INCL input
92  if(calincl != 0) {
93    delete calincl;
94    calincl = 0;
95  }
96
97  if (verboseLevel > 1) {
98    G4cout << " >>> G4InclAblaLightIonInterface::ApplyYourself called" << G4endl;
99  }
100
101  if(verboseLevel > 1) {
102    G4cout <<"G4InclAblaLightIonInterface: Now processing INCL4 event number:" << eventNumber << G4endl;
103  }
104
105  // Inverse kinematics for targets with Z = 1 and A = 1
106  //  if(false) {
107  G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector();
108
109  if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4InclInput::canUseInverseKinematics(aTrack, theNucleus)) {
110    G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0);
111    const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition();
112
113    if(oldProjectileDef != 0 && oldTargetDef != 0) {
114    G4int oldTargetA = oldTargetDef->GetAtomicMass();
115    G4int newTargetA = oldProjectileDef->GetAtomicMass();
116    G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
117
118    if(newTargetA > 0 && newTargetZ > 0) {
119      G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber());
120
121      //      G4cout <<"Original projectile kinE = " << aTrack.GetKineticEnergy() / MeV << G4endl;
122
123      // We need the same energy/nucleon.
124      G4double projectileE = ((aTrack.GetKineticEnergy() / MeV) / newTargetA) * oldTargetA * MeV;
125
126      //    G4cout <<"projectileE = " << projectileE << G4endl;
127      G4DynamicParticle swappedProjectileParticle(oldTargetDef, G4ThreeVector(0.0, 0.0, 1.0), projectileE);
128      const G4LorentzVector swapped4Momentum = (swappedProjectileParticle.Get4Momentum()*=toBreit);
129      swappedProjectileParticle.Set4Momentum(swapped4Momentum);
130      const G4HadProjectile swappedProjectile(swappedProjectileParticle);
131      //  G4cout <<"New projectile kinE = " << swappedProjectile.GetKineticEnergy() / MeV << G4endl;
132      calincl = new G4InclInput(swappedProjectile, swappedTarget, true);
133    } else {
134      G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl;
135      calincl = new G4InclInput(aTrack, theNucleus, false);
136    }
137    }
138  } else {
139    calincl = new G4InclInput(aTrack, theNucleus, false);
140  }
141
142  G4double eKin;
143  G4double momx = 0.0, momy = 0.0, momz = 0.0;
144  G4DynamicParticle *cascadeParticle = 0;
145  G4ParticleDefinition *aParticleDefinition = 0;
146
147  // INCL assumes the projectile particle is going in the direction of
148  // the Z-axis. Here we construct proper rotation to convert the
149  // momentum vectors of the outcoming particles to the original
150  // coordinate system.
151  G4LorentzVector projectileMomentum = aTrack.Get4Momentum();
152  G4LorentzRotation toZ;
153  toZ.rotateZ(-projectileMomentum.phi());
154  toZ.rotateY(-projectileMomentum.theta());
155  G4LorentzRotation toLabFrame = toZ.inverse();
156
157  /*
158  G4cout <<"Projectile theta = " << projectileMomentum.theta() << " phi = " << projectileMomentum.phi() << G4endl;
159  G4cout <<"Projectile momentum "
160         << "(px = " << projectileMomentum.px()
161         << ", py = " << projectileMomentum.py()
162         << ", pz = " << projectileMomentum.pz() << ")" << G4endl;
163  G4cout << "Projectile energy = " << bulletE << " MeV" << G4endl;
164  */
165
166  G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp();
167  G4FragmentVector *theSpectatorFermiBreakupResult = 0;
168  G4FragmentVector *theFermiBreakupResult = 0;
169
170  theResult.Clear(); // Make sure the output data structure is clean.
171
172  std::vector<G4DynamicParticle*> result; // Temporary list for the results
173
174  // Map Geant4 particle types to corresponding INCL4 types.
175  enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 
176                           pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9,
177                           c12 = -12}; // Carbon beam support.
178
179  G4int bulletType = calincl->bulletType();
180  chargeNumberBalanceInINCL = calincl->targetZ();
181  baryonNumberBalanceInINCL = calincl->targetA();
182
183  //  G4cout <<"Type of the projectile (INCL projectile code): " << bulletType << G4endl;
184
185  if(bulletType == proton) {
186    chargeNumberBalanceInINCL += 1;
187    baryonNumberBalanceInINCL += 1;
188  } else if(bulletType == neutron) {
189    baryonNumberBalanceInINCL += 1;
190  } else if(bulletType == pionPlus) { //Note: positive pion doesn't contribute to the baryon and charge number counters
191    chargeNumberBalanceInINCL += 1;
192  } else if(bulletType == pionMinus) {
193    chargeNumberBalanceInINCL -= 1;
194  } else if(bulletType == deuteron) {
195    chargeNumberBalanceInINCL += 1;
196    baryonNumberBalanceInINCL += 2;
197  } else if(bulletType == triton) {
198    chargeNumberBalanceInINCL += 1;
199    baryonNumberBalanceInINCL += 3;
200  } else if(bulletType == he3) {
201    chargeNumberBalanceInINCL += 2;
202    baryonNumberBalanceInINCL += 3;
203  } else if(bulletType == he4) {
204    chargeNumberBalanceInINCL += 2;
205    baryonNumberBalanceInINCL += 4;
206  } if(bulletType == c12) {
207    chargeNumberBalanceInINCL += 6;
208    baryonNumberBalanceInINCL += 12;
209  } if(bulletType == -666) {
210    chargeNumberBalanceInINCL += calincl->extendedProjectileZ();
211    baryonNumberBalanceInINCL += calincl->extendedProjectileA();
212  }
213
214  // Check wheter the input is acceptable.
215  if((bulletType != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
216    ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon
217    ws->xfoisa = 8;
218    ws->npaulstr = 0;
219
220    int nTries = 0;
221    varntp->ntrack = 0;
222
223    mat->nbmat = 1;
224    mat->amat[0] = int(calincl->targetA());
225    mat->zmat[0] = int(calincl->targetA());
226
227    incl->setInput(calincl);
228    incl->initIncl(true);
229
230    while((varntp->ntrack <= 0) && (nTries < maxTries)) { // Loop until we produce real cascade
231      nTries++;
232      if(verboseLevel > 1) {
233        G4cout <<"G4InclAblaLightIonInterface: Try number = " << nTries << G4endl; 
234      }
235      incl->processEventInclAbla(calincl, eventNumber);
236
237      if(verboseLevel > 1) {
238        G4cout <<"G4InclAblaLightIonInterface: number of tracks = " <<  varntp->ntrack <<G4endl;
239      }
240    }
241
242    if(verboseLevel > 1) {
243      /**
244       * Diagnostic output
245       */
246      G4cout <<"G4InclAblaLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl;
247      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
248      if(bulletType == -666) {
249        G4cout <<"   Extended projectile: A = " << calincl->extendedProjectileA()
250               <<" Z = " << calincl->extendedProjectileZ() << G4endl;
251      }
252
253      G4cout <<"G4InclAblaLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
254      G4cout <<"G4InclAblaLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
255
256      if(verboseLevel > 3) {
257        diagdata <<"G4InclAblaLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl;
258        diagdata <<"G4InclAblaLightIonInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
259       
260        diagdata <<"G4InclAblaLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
261        diagdata <<"G4InclAblaLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
262      }
263    }
264
265    // Check whether a valid cascade was produced.
266    // If not return the original bullet particle with the same momentum.
267    if(varntp->ntrack <= 0) {
268      if(verboseLevel > 1) {
269        G4cout <<"WARNING G4InclAblaLightIonInterface: No cascade. Returning original particle with original momentum." << G4endl;
270        G4cout <<"\t Reached maximum trials of 200 to produce inelastic scattering." << G4endl;
271      }
272
273      theResult.SetStatusChange(stopAndKill);
274     
275      if(bulletType == proton) {
276        aParticleDefinition = G4Proton::ProtonDefinition();
277      } else if(bulletType == neutron) {
278        aParticleDefinition = G4Neutron::NeutronDefinition();
279      } else if(bulletType == pionPlus) {
280        aParticleDefinition = G4PionPlus::PionPlusDefinition();
281      } else if(bulletType == pionZero) {
282        aParticleDefinition = G4PionZero::PionZeroDefinition();
283      } else if(bulletType == pionMinus) {
284        aParticleDefinition = G4PionMinus::PionMinusDefinition();
285      } else if(bulletType == deuteron) {
286        aParticleDefinition = G4Deuteron::DeuteronDefinition();
287      } else if(bulletType == triton) {
288        aParticleDefinition = G4Triton::TritonDefinition();
289      } else if(bulletType == he3) {
290        aParticleDefinition = G4He3::He3Definition();
291      } else if(bulletType == he4) {
292        aParticleDefinition = G4Alpha::AlphaDefinition();
293      } else { // Particle was not recognized. Probably an unsupported particle was given as input
294        aParticleDefinition = 0;
295      }
296
297      if(aParticleDefinition != 0) {
298        cascadeParticle = new G4DynamicParticle();
299        cascadeParticle->SetDefinition(aParticleDefinition);
300        cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
301        result.push_back(cascadeParticle);
302      }
303    }
304
305    // Convert INCL4 output to Geant4 compatible data structures.
306    // Elementary particles are converted to G4DynamicParticle.
307    theResult.SetStatusChange(stopAndKill);
308   
309    for(particleI = 0; particleI <= varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.
310      // Get energy/momentum and construct momentum vector in INCL4 coordinates.
311      //      if(varntp->itypcasc[particleI] == -1) continue; // Avoid nucleons that are part of the spectator
312      if(varntp->avv[particleI] == 0 && varntp->zvv[particleI] == 0) continue;
313      momx = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::cos(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
314      momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
315      momz = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*MeV;
316
317      eKin = varntp->enerj[particleI] * MeV;
318
319      G4ThreeVector momDirection(momx, momy, momz); // Direction of the particle.
320      momDirection = momDirection.unit();
321      if(verboseLevel > 2) {
322        G4cout <<"G4InclAblaLightIonInterface: " << G4endl;
323        G4cout <<"A    = " << varntp->avv[particleI] << " Z = "  << varntp->zvv[particleI] << G4endl;
324        G4cout <<"eKin = " << eKin                   << " MeV"   << G4endl;
325        G4cout <<"px   = " << momDirection.x()       << " py = " << momDirection.y()       <<" pz = " << momDirection.z() << G4endl;
326      }
327
328      G4int particleIdentified = 0; // Check particle ID.
329
330      if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) { // Proton
331        cascadeParticle = 
332          new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
333        particleIdentified++;
334        baryonNumberBalanceInINCL -= 1;
335        chargeNumberBalanceInINCL -= 1;
336      }
337
338      if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) { // Neutron
339        cascadeParticle = 
340          new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
341        particleIdentified++;
342        baryonNumberBalanceInINCL -= 1;
343      }
344
345      if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) { // PionPlus
346        cascadeParticle = 
347          new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
348        particleIdentified++;
349        chargeNumberBalanceInINCL -= 1;
350      }
351
352      if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) { // PionZero
353        cascadeParticle = 
354          new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
355        particleIdentified++;
356        chargeNumberBalanceInINCL -= 0;
357      }
358
359      if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) { // PionMinus
360        cascadeParticle = 
361          new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
362        particleIdentified++;
363        chargeNumberBalanceInINCL -= -1;
364      }
365
366      if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
367        G4ParticleDefinition * aIonDef = 0;
368
369        G4int A = G4int(varntp->avv[particleI]);
370        G4int Z = G4int(varntp->zvv[particleI]);
371        G4double excitationE = G4double(varntp->exini) * MeV;
372
373        if(verboseLevel > 1) {
374          G4cout <<"Finding ion: A = " << A << " Z = " << Z << " E* = " << excitationE/MeV << G4endl;
375        }
376        aIonDef = theTableOfParticles->GetIon(Z, A, excitationE);
377       
378        if(aIonDef == 0) {
379          if(verboseLevel > 1) {
380            G4cout <<"G4InclAblaLightIonInterface: " << G4endl;
381            G4cout <<"FATAL ERROR: aIonDef = 0" << G4endl;
382            G4cout <<"A = " << A << " Z = " << Z << " E* = " << excitationE << G4endl;
383          }
384        }
385
386        if(aIonDef != 0) { // If the ion was identified add it to output.
387          cascadeParticle =
388            new G4DynamicParticle(aIonDef, momDirection, eKin);
389          particleIdentified++;
390          baryonNumberBalanceInINCL -= A;
391          chargeNumberBalanceInINCL -= Z;
392        }
393      }
394       
395      if(particleIdentified == 1) { // Particle identified properly.
396        cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
397        result.push_back(cascadeParticle);
398      }
399      else { // Particle identification failed.
400        if(particleIdentified > 1) { // Particle was identified as more than one particle type.
401          if(verboseLevel > 1) {
402            G4cout <<"G4InclAblaLightIonInterface: One outcoming particle was identified as";
403            G4cout <<"more than one particle type. This is probably due to a bug in the interface." << G4endl;
404            G4cout <<"Particle A:" << varntp->avv[particleI] << "Z: " << varntp->zvv[particleI] << G4endl;
405            G4cout << "(particleIdentified =" << particleIdentified << ")"  << G4endl;
406          }
407        }
408      }
409    }
410
411    // Spectator nucleus Fermi break-up
412    if(useFermiBreakup && useProjectileSpectator && varntp->masp > 1) {
413      baryonNumberBalanceInINCL -= G4int(varntp->masp);
414      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->masp), G4int(varntp->mzsp)) + varntp->exsp * MeV;
415      // Use momentum scaling to compensate for different masses in G4 and INCL:
416      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->masp),
417                                                                        G4int(varntp->mzsp),
418                                                                        varntp->exsp,
419                                                                        varntp->spectatorT,
420                                                                        varntp->spectatorP1,
421                                                                        varntp->spectatorP2,
422                                                                        varntp->spectatorP3);
423      G4LorentzVector p4(momentumScaling * varntp->spectatorP1 * MeV, momentumScaling * varntp->spectatorP2 * MeV,
424                         momentumScaling * varntp->spectatorP3 * MeV,
425                         varntp->spectatorT * MeV + nuclearMass);
426      // Four-momentum, baryon number and charge balance:
427      G4LorentzVector fourMomentumBalance = p4;
428      G4int baryonNumberBalance = G4int(varntp->masp);
429      chargeNumberBalanceInINCL -= G4int(varntp->mzsp);
430      G4int chargeBalance = G4int(varntp->mzsp);
431
432      G4LorentzRotation toFragmentZ;
433      // Assume that Fermi breakup uses Z as the direction of the projectile
434      toFragmentZ.rotateZ(-p4.theta());
435      toFragmentZ.rotateY(-p4.phi());
436      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
437      //      p4 *= toFragmentZ;
438     
439      G4LorentzVector p4rest = p4;
440      //      p4rest.boost(-p4.boostVector());
441      if(verboseLevel > 0) {
442        G4cout <<"Spectator nucleus:" << G4endl;
443        G4cout <<"p4: " << G4endl;
444        G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
445        G4cout <<" E = " << p4.e() << G4endl;
446        G4cout <<"p4rest: " << G4endl;
447        G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
448        G4cout <<" E = " << p4rest.e() << G4endl;
449      }
450      G4Fragment theSpectatorNucleus(G4int(varntp->masp), G4int(varntp->mzsp), p4rest);
451      theSpectatorFermiBreakupResult = fermiBreakUp->BreakItUp(theSpectatorNucleus);
452      if(theSpectatorFermiBreakupResult != 0) {
453      G4FragmentVector::iterator fragment;
454      for(fragment = theSpectatorFermiBreakupResult->begin(); fragment != theSpectatorFermiBreakupResult->end(); fragment++) {
455        G4ParticleDefinition *theFragmentDefinition = 0;
456        if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron
457          theFragmentDefinition = G4Neutron::NeutronDefinition();
458        } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) {
459          theFragmentDefinition = G4Proton::ProtonDefinition();
460        } else {
461          theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
462        }
463        if(theFragmentDefinition != 0) {
464          G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
465          G4LorentzVector labMomentum = theFragment->Get4Momentum();
466          //      labMomentum.boost(p4.boostVector());
467          //      labMomentum *= toFragmentLab;
468          //      labMomentum *= toLabFrame;
469          theFragment->Set4Momentum(labMomentum);
470          fourMomentumBalance -= theFragment->Get4Momentum();
471          baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
472          chargeBalance -= theFragmentDefinition->GetAtomicNumber();
473          if(verboseLevel > 0) {
474            G4cout <<"Resulting fragment: " << G4endl;
475            G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
476            G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
477          }
478          result.push_back(theFragment);
479        } else {
480          G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." 
481                 << G4endl;
482          G4cout <<"Resulting fragment: " << G4endl;
483          G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
484          G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
485          G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
486          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
487        }
488      }
489      delete theSpectatorFermiBreakupResult;
490      theSpectatorFermiBreakupResult = 0;
491
492      if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
493        G4cout <<"Four-momentum balance after spectator nucleus Fermi break-up:" << G4endl;
494        G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
495        G4cout <<"Vector components (px, py, pz, E) = ("
496               << fourMomentumBalance.px() << ", "
497               << fourMomentumBalance.py() << ", "
498               << fourMomentumBalance.pz() << ", "
499               << fourMomentumBalance.e() << ")" << G4endl;
500      }
501      if(baryonNumberBalance != 0) {
502        G4cout <<"Event " << eventNumber << ": Baryon number balance after spectator nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
503      }
504      if(chargeBalance != 0) {
505        G4cout <<"Event " << eventNumber <<": Charge balance after spectator nucleus Fermi break-up: " << chargeBalance << G4endl;
506      }
507    }
508  }
509
510    // Finally do Fermi break-up if needed
511    if(varntp->needsFermiBreakup && varntp->massini > 0) {
512      baryonNumberBalanceInINCL -= G4int(varntp->massini);
513      chargeNumberBalanceInINCL -= G4int(varntp->mzini);
514      // Call Fermi Break-up
515      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV;
516      G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV,
517                                       varntp->erecrem * MeV + nuclearMass);
518      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini),
519                                                                        varntp->exini,
520                                                                        varntp->erecrem,
521                                                                        varntp->pxrem,
522                                                                        varntp->pyrem,
523                                                                        varntp->pzrem);
524      G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV,
525                         momentumScaling * varntp->pzrem * MeV,
526                         varntp->erecrem + nuclearMass);
527
528      // For four-momentum, baryon number and charge conservation check:
529      G4LorentzVector fourMomentumBalance = p4;
530      G4int baryonNumberBalance = G4int(varntp->massini);
531      G4int chargeBalance = G4int(varntp->mzini);
532
533      G4LorentzRotation toFragmentZ;
534      toFragmentZ.rotateZ(-p4.theta());
535      toFragmentZ.rotateY(-p4.phi());
536      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
537      //      p4 *= toFragmentZ;
538
539      G4LorentzVector p4rest = p4;
540      //      p4rest.boost(-p4.boostVector());
541      if(verboseLevel > 0) {
542        G4cout <<"Cascade remnant nucleus:" << G4endl;
543        G4cout <<"p4: " << G4endl;
544        G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
545        G4cout <<" E = " << p4.e() << G4endl;
546
547        G4cout <<"p4rest: " << G4endl;
548        G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
549        G4cout <<" E = " << p4rest.e() << G4endl;
550      }
551
552      G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest);
553      theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant);
554      if(theFermiBreakupResult != 0) {
555      G4FragmentVector::iterator fragment;
556      for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) {
557        G4ParticleDefinition *theFragmentDefinition = 0;
558        if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron
559          theFragmentDefinition = G4Neutron::NeutronDefinition();
560        } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) {
561          theFragmentDefinition = G4Proton::ProtonDefinition();
562        } else {
563          theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
564        }
565
566        if(theFragmentDefinition != 0) {
567          G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
568          G4LorentzVector labMomentum = theFragment->Get4Momentum();
569          //      labMomentum.boost(p4.boostVector());
570          //      labMomentum *= toFragmentLab;
571          //      labMomentum *= toLabFrame;
572          theFragment->Set4Momentum(labMomentum);
573          fourMomentumBalance -= theFragment->Get4Momentum();
574          baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
575          chargeBalance -= theFragmentDefinition->GetAtomicNumber();
576          if(verboseLevel > 0) {
577            G4cout <<"Resulting fragment: " << G4endl;
578            G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
579            G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
580          }
581          result.push_back(theFragment);
582        } else {
583          G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
584          G4cout <<"Resulting fragment: " << G4endl;
585          G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
586          G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
587          G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
588          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
589        }
590      }
591      delete theFermiBreakupResult;
592      theFermiBreakupResult = 0;
593
594      if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
595        G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
596        G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
597        G4cout <<"Vector components (px, py, pz, E) = ("
598               << fourMomentumBalance.px() << ", "
599               << fourMomentumBalance.py() << ", "
600               << fourMomentumBalance.pz() << ", "
601               << fourMomentumBalance.e() << ")" << G4endl;
602      }
603      if(baryonNumberBalance != 0) {
604        G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
605      }
606      if(chargeBalance != 0) {
607        G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
608      }
609    }
610    }
611
612    varntp->ntrack = 0; // Clean up the number of generated particles in the event.
613
614    if(baryonNumberBalanceInINCL != 0 && verboseLevel > 1) {
615      G4cout <<"Event " << eventNumber <<": G4InclAblaLightIonInterface: Baryon number conservation problem in INCL detected!" << G4endl;
616      G4cout <<"Baryon number balance: " << baryonNumberBalanceInINCL << G4endl;
617      if(baryonNumberBalanceInINCL < 0) {
618        G4cout <<"Event " << eventNumber <<": Too many outcoming baryons!" << G4endl;
619      } else if(baryonNumberBalanceInINCL > 0) {
620        G4cout <<"Event " << eventNumber <<": Too few outcoming baryons!" << G4endl;
621      }
622    }
623
624    if(chargeNumberBalanceInINCL != 0 && verboseLevel > 1) {
625      G4cout <<"Event " << eventNumber <<": G4InclAblaLightIonInterface: Charge number conservation problem in INCL detected!" << G4endl;
626      G4cout <<"Event " << eventNumber <<": Charge number balance: " << chargeNumberBalanceInINCL << G4endl;
627    }
628  }
629  /**
630   * Report unsupported features.
631   * (Check bullet, target, energy range)
632   */
633  else { // If the bullet type was not recognized by the interface, it will be returned back without any interaction.
634    theResult.SetStatusChange(stopAndKill);
635
636    G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
637    cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
638
639    result.push_back(cascadeParticle);
640
641    if(verboseLevel > 1) {
642      G4cout <<"G4InclAblaLightIonInterface: Error processing event number (internal) " << eventNumber << G4endl;
643    }
644    if(verboseLevel > 3) {
645      diagdata <<"G4InclAblaLightIonInterface: Error processing event number (internal) " << eventNumber << G4endl;
646    }
647
648    if(bulletType == 0) {
649      if(verboseLevel > 1) {
650        G4cout <<"G4InclAblaLightIonInterface: Unknown bullet type" << G4endl;     
651        G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
652      }
653      if(verboseLevel > 3) {
654        diagdata <<"G4InclAblaLightIonInterface: Unknown bullet type" << G4endl;     
655        diagdata <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
656      }
657    }
658
659    if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
660      if(verboseLevel > 1) {
661        G4cout <<"Unsupported target: " << G4endl;
662        G4cout <<"Target A: " << calincl->targetA() << G4endl;
663        G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
664      }
665      if(verboseLevel > 3) {
666        diagdata <<"Unsupported target: " << G4endl;
667        diagdata <<"Target A: " << calincl->targetA() << G4endl;
668        diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
669      }
670    }
671
672    if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
673      if(verboseLevel > 1) {
674        G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
675        G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
676      }
677      if(verboseLevel > 3) {
678        diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
679      }
680    }
681
682    if(verboseLevel > 3) {
683      diagdata <<"WARNING: returning the original bullet with original energy back to Geant4." << G4endl;
684    }
685  }
686
687  // Finally copy the accumulated secondaries into the result collection:
688  G4ThreeVector boostVector = aTrack.Get4Momentum().boostVector();
689  G4LorentzRotation boostBack = toBreit.inverse();
690
691  for(std::vector<G4DynamicParticle*>::iterator i = result.begin(); i != result.end(); ++i) {
692    // If the calculation was performed in inverse kinematics we have to
693    // convert the result back...
694    if(calincl->isInverseKinematics()) {
695      G4LorentzVector mom = (*i)->Get4Momentum();
696      mom.setPz(-1.0 * mom.pz()); // Reverse the z-component of the momentum vector
697      mom *= boostBack;
698      (*i)->Set4Momentum(mom);
699    }
700    theResult.AddSecondary((*i));
701  }
702
703  delete fermiBreakUp;
704  delete calincl;
705  calincl = 0;
706  return &theResult;
707} 
708
709G4ReactionProductVector* G4InclAblaLightIonInterface::Propagate(G4KineticTrackVector* , G4V3DNucleus* ) {
710  return 0;
711}
712
713
Note: See TracBrowser for help on using the repository browser.