source: trunk/source/processes/hadronic/models/cascade/evaporation/src/G4BertiniEvaporation.cc @ 1196

Last change on this file since 1196 was 962, checked in by garnier, 15 years ago

update processes
File size: 15.9 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//
27// Implementation of the HETC88 code into Geant4.
28// Evaporation and De-excitation parts
29// T. Lampen, Helsinki Institute of Physics, May-2000
30
31#include "globals.hh"
32#include "G4BertiniEvaporation.hh"
33#include "G4BENeutronChannel.hh"
34#include "G4BEProtonChannel.hh"
35#include "G4BEDeuteronChannel.hh"
36#include "G4BETritonChannel.hh"
37#include "G4BEHe3Channel.hh"
38#include "G4BEHe4Channel.hh"
39#include "G4BEGammaDeexcitation.hh"
40
41// verboseLevels : 4 inform about basic probabilities and branchings
42//                 6 some details about calculations
43//                 8 inform about final values
44//                10 inform everything
45
46G4BertiniEvaporation::G4BertiniEvaporation()
47{
48    G4cout << "Info G4BertiniEvaporation: This is still very fresh code."<< G4endl;
49    G4cout << "     G4BertiniEvaporation: feed-back for improvement is very wellcome."<< G4endl;
50    verboseLevel = 0;
51 
52    channelVector.push_back( new G4BENeutronChannel );
53    channelVector.push_back( new G4BEProtonChannel  );
54    channelVector.push_back( new G4BEDeuteronChannel);
55    channelVector.push_back( new G4BETritonChannel  );
56    channelVector.push_back( new G4BEHe3Channel     );
57    channelVector.push_back( new G4BEHe4Channel     );
58}
59
60
61G4BertiniEvaporation::~G4BertiniEvaporation()
62{
63  while ( !channelVector.empty() )
64    {
65      delete channelVector.back();
66      channelVector.pop_back();
67    }
68}
69
70
71void G4BertiniEvaporation::setVerboseLevel( const G4int verbose )
72{
73  verboseLevel = verbose;
74
75  // Update verbose level to all evaporation channels.
76  std::vector< G4BertiniEvaporationChannel * >::iterator iChannel = channelVector.begin(); 
77  for ( ; iChannel != channelVector.end() ; *iChannel++ )                 
78    ( *iChannel )->setVerboseLevel( verboseLevel );
79}
80
81
82G4FragmentVector * G4BertiniEvaporation::BreakItUp( G4LayeredNucleus & nucleus )
83{
84  G4int nucleusA;
85  G4int nucleusZ;
86  G4int i;
87  G4double totalProbability;
88  G4double excE;
89  G4double newExcitation;
90  G4double nucleusTotalMomentum;
91  G4double nucleusKineticEnergy;
92  G4double mRes; // Mass of residual nucleus.
93  G4ThreeVector nucleusMomentumVector;
94  G4DynamicParticle *pEmittedParticle;
95  std::vector< G4DynamicParticle * > secondaryParticleVector;
96  G4FragmentVector * result = new G4FragmentVector;
97 
98  // Read properties of the nucleus.
99  nucleusA = ( G4int ) nucleus.GetN(); // GetN should in fact get GetA
100  nucleusZ = ( G4int ) nucleus.GetZ();
101  excE                  = nucleus.GetEnergyDeposit();
102  nucleusMomentumVector = nucleus.GetMomentum();
103
104  // Move to CMS frame, save initial velocity of the nucleus to boostToLab vector.
105  G4ThreeVector boostToLab( ( 1/G4NucleiProperties::GetNuclearMass( nucleusA, nucleusZ ) ) 
106                            * nucleusMomentumVector ); // xx mass ok?
107
108  if ( verboseLevel >= 10 )
109    G4cout << " G4BertiniEvaporation : initial kinematics : boostToLab vector = " << boostToLab << G4endl
110           << "                     excitation energy  : " << excE<< G4endl;
111
112  if ( nucleusA == 8 && nucleusZ == 4 )
113    {
114      splitBe8( excE, boostToLab, secondaryParticleVector );
115      fillResult( secondaryParticleVector, result);
116      return result;
117    }
118 
119  // Initialize evaporation channels and calculate sum of emission
120  // probabilities.
121  std::vector< G4BertiniEvaporationChannel * >::iterator iChannel = channelVector.begin(); 
122  totalProbability = 0;
123  for ( ; iChannel != channelVector.end() ; *iChannel++ )                 
124     {
125       ( *iChannel )->setNucleusA( nucleusA );
126       ( *iChannel )->setNucleusZ( nucleusZ );
127       ( *iChannel )->setExcitationEnergy( excE );
128       ( *iChannel )->setPairingCorrection( 1 );
129       ( *iChannel )->calculateProbability();
130       totalProbability += ( *iChannel )->getProbability();
131     }
132
133  // If no evaporation process is possible, try without pairing energy.
134  if ( totalProbability == 0 )
135    {
136      if ( verboseLevel >= 4 )
137        G4cout << "G4BertiniEvaporation : no emission possible with pairing correction, trying without it" << G4endl;
138      for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ )
139        {
140          ( *iChannel )->setPairingCorrection(0);
141          ( *iChannel )->calculateProbability();
142          totalProbability += ( *iChannel )->getProbability();
143        }
144      if ( verboseLevel >= 4 )
145        G4cout << "                       probability without correction " << totalProbability << G4endl;
146
147    }
148
149  // Normal evaporation cycle follows.
150  // Particles are evaporated till all probabilities are zero.
151  while ( totalProbability > 0 )
152    {
153      G4BertiniEvaporationChannel *pSelectedChannel;
154
155      // Sample active emission channel.
156      G4double sampledProbability = G4UniformRand() * totalProbability;
157      if ( verboseLevel >= 10 ) G4cout << "          RandomProb : " << sampledProbability << G4endl;
158      for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ )
159        {
160          sampledProbability -= ( *iChannel )->getProbability();
161          if ( sampledProbability < 0 ) break;
162        }
163      pSelectedChannel = ( *iChannel );
164      if ( iChannel == channelVector.end() )
165        throw G4HadronicException(__FILE__, __LINE__,  "G4BertiniEvaporation : Error while sampling evaporation particle" );
166
167      if ( verboseLevel >= 4 ) 
168        G4cout << "G4BertiniEvaporation : particle " << pSelectedChannel->getName() << " selected " << G4endl;
169
170      // Max 10 tries to get a physically acceptable particle energy
171      // in this emission channel.
172      i = 0;
173
174      do
175        {
176          pEmittedParticle = pSelectedChannel->emit();
177          // This loop checks that particle with too large energy is not emitted.
178          // CMS frame is considered in this loop. Nonrelativistic treatment. xxx
179
180
181          const G4int zRes = nucleusZ - pSelectedChannel->getParticleZ(); 
182          const G4int aRes  = nucleusA - pSelectedChannel->getParticleA(); 
183          // const G4double eBind = G4NucleiProperties::GetBindingEnergy( aRes, zRes );  // Binding energy of the nucleus.
184          mRes  = G4NucleiProperties::GetNuclearMass( aRes, zRes ); // Mass of the target nucleus
185          //      In HETC88:
186          //      eBind = Z * (-0.78244) + A * 8.36755 - cameron ( A , Z );
187          //      mRes =  zRes * 938.79304 + ( aRes - zRes ) * 939.57548 - eBind;
188          //      where cameron ( A, Z ) is the mass excess by Cameron, see Canadian Journal of Physics,
189          //      vol. 35, 1957, p.1022
190
191          nucleusTotalMomentum = pEmittedParticle->GetTotalMomentum(); // CMS frame
192          nucleusKineticEnergy = std::pow( nucleusTotalMomentum, 2 ) / ( 2 * mRes );
193          newExcitation = excE - pEmittedParticle->GetKineticEnergy() - nucleusKineticEnergy - pSelectedChannel->getQ();
194
195          if ( verboseLevel >= 10)
196            G4cout << "G4BertiniEvaporation : Kinematics " << G4endl
197                   << "                    part kinetic E in CMS = " 
198                   << pEmittedParticle->GetKineticEnergy() << G4endl
199                   << "                    new excitation E      =  " 
200                   << newExcitation << G4endl;
201
202          i++;
203          if ( !( newExcitation > 0 ) && i < 10) delete pEmittedParticle;
204        } while ( !( newExcitation > 0 )  &&  i < 10 );
205
206      if ( i >= 10 ) 
207        {
208          // No appropriate particle energy found.
209          // Set probability of this channel to zero
210          // and try to sample another particle on the
211          // next round.
212          delete pEmittedParticle;
213
214          if ( verboseLevel >= 4 )
215            G4cout << "G4BertiniEvaporation : No appropriate energy for particle " 
216                   << pSelectedChannel->getName() << " found." << G4endl;
217
218          pSelectedChannel->setProbability( 0 );
219
220          totalProbability = 0;
221          for ( ; iChannel != channelVector.end() ; *iChannel++ )                 
222            totalProbability += ( *iChannel )->getProbability();
223        } // Treatment of physically unacceptable particle ends.
224      else
225        {
226          // Good particle found.
227
228          // Transform particle properties to lab frame and save it.
229          G4LorentzVector particleFourVector = pEmittedParticle->Get4Momentum();
230          G4LorentzVector nucleusFourVector(  - pEmittedParticle->GetMomentum(), // CMS Frame.
231                                              nucleusKineticEnergy + mRes ); // Total energy.
232          particleFourVector.boost( boostToLab );
233          nucleusFourVector.boost(  boostToLab );
234          G4DynamicParticle *pPartLab  = new G4DynamicParticle( pEmittedParticle->GetDefinition(),
235                                                                particleFourVector.e(), // Total energy
236                                                                particleFourVector.vect() ); // momentum vector
237          secondaryParticleVector.push_back( pPartLab );
238
239          // Update residual nucleus and boostToLab vector.
240          nucleusA = nucleusA - pSelectedChannel->getParticleA();
241          nucleusZ = nucleusZ - pSelectedChannel->getParticleZ();
242          excE = newExcitation;
243          boostToLab = G4ThreeVector ( ( 1/mRes ) * nucleusFourVector.vect() );
244
245          if ( verboseLevel >= 10 )
246            G4cout << "  particle mom in cms " << pEmittedParticle->GetMomentum() << G4endl
247                   << "  particle mom in cms " << pEmittedParticle->GetTotalMomentum() << G4endl
248                   << "          Etot in cms " << pEmittedParticle->GetTotalEnergy() << G4endl
249                   << "          Ekin in cms " << pEmittedParticle->GetKineticEnergy() << G4endl
250                   << "        particle mass " << pEmittedParticle->GetMass() << G4endl
251                   << "          boost  vect " << boostToLab << G4endl
252                   << "          boosted 4v  " << particleFourVector << G4endl
253                   << "          nucleus 4v  " << nucleusFourVector << G4endl
254                   << "          nucl cm mom " << nucleusTotalMomentum << G4endl
255                   << " particle k.e. in lab " << pPartLab->GetKineticEnergy() << G4endl
256                   << "     new boost vector " << boostToLab << G4endl;
257
258          delete pEmittedParticle;
259
260          // If the residual nucleus is Be8, split it.
261          if ( nucleusA == 8 && nucleusZ == 4 )
262            {
263              splitBe8( excE, boostToLab, secondaryParticleVector );
264              fillResult( secondaryParticleVector, result);
265              return result;
266            }
267
268          // Update evaporation channels and
269          // total emission probability.
270          totalProbability = 0;
271          for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ ) 
272            {
273              ( *iChannel )->setNucleusA( nucleusA ); 
274              ( *iChannel )->setNucleusZ( nucleusZ );
275              ( *iChannel )->setExcitationEnergy( excE ); 
276              ( *iChannel )->calculateProbability();
277              totalProbability += ( *iChannel )->getProbability();
278            }
279        } // Treatment of physically acceptable particle ends.
280    } // End of evaporation cycle.
281
282  // Gamma de-excitation.
283  G4BEGammaDeexcitation *pGammaChannel = new G4BEGammaDeexcitation;
284  pGammaChannel->setNucleusA( nucleusA );
285  pGammaChannel->setNucleusZ( nucleusZ );
286  pGammaChannel->setVerboseLevel( verboseLevel );
287  pGammaChannel->setExcitationEnergy( excE );
288
289  // Emit equally sampled photons until all the excitation energy is
290  // used; the last photon gets the remaining de-excitation energy.
291  // Change of momentum of the nucleus is neglected.
292  G4double gammaEnergy;
293  G4double totalDeExcEnergy = 0;
294
295  while ( excE > 0 )
296    {
297      pEmittedParticle = pGammaChannel->emit();
298      gammaEnergy = pEmittedParticle->GetKineticEnergy();
299     
300      if ( ( totalDeExcEnergy + gammaEnergy ) > excE ) 
301        {
302          // All the remaining energy is used here.
303          gammaEnergy = excE - totalDeExcEnergy;
304          excE = 0;
305        }
306
307      totalDeExcEnergy += gammaEnergy; 
308
309      if ( verboseLevel >= 10 )
310        G4cout << " G4BertiniEvaporation : gamma de-excitation, g of " << gammaEnergy << " MeV " << G4endl;
311     
312      G4LorentzVector gammaFourVector( pEmittedParticle->GetMomentum(),
313                                       pEmittedParticle->GetKineticEnergy() );
314      gammaFourVector.boost( boostToLab );
315      pEmittedParticle->SetKineticEnergy( gammaFourVector.e() );
316      pEmittedParticle->SetMomentumDirection( gammaFourVector.vect().unit() );
317
318      secondaryParticleVector.push_back( pEmittedParticle );
319
320      } 
321
322  delete pGammaChannel;
323
324  // Residual nucleus is not returned.
325
326  fillResult( secondaryParticleVector, result);
327 
328  return result;
329}
330
331
332void G4BertiniEvaporation::splitBe8( const G4double E,
333                                     const G4ThreeVector boostToLab,
334                                     std::vector< G4DynamicParticle * > & secondaryParticleVector )
335{
336  G4double kineticEnergy; 
337  G4double u;
338  G4double v;
339  G4double w;
340  const G4double Be8DecayEnergy = 0.093 * MeV; 
341
342  if ( E <= 0 ) throw G4HadronicException(__FILE__, __LINE__,  "G4BertiniEvaporation : excitation energy < 0 " );
343  if ( verboseLevel >= 4 ) G4cout << "     Be8 split to 2 x He4" << G4endl;
344
345  kineticEnergy = 0.5 * ( E + Be8DecayEnergy ); 
346 
347  // Create two alpha particles in CMS frame.
348  isotropicCosines( u , v , w );
349  G4ThreeVector momentumDirectionCMS( u, v, w ); 
350
351  G4DynamicParticle *pP1Cms = new G4DynamicParticle( G4Alpha::Alpha(),
352                                                     momentumDirectionCMS,
353                                                     kineticEnergy );
354  G4DynamicParticle *pP2Cms = new G4DynamicParticle( G4Alpha::Alpha(),
355                                                     -momentumDirectionCMS,
356                                                     kineticEnergy );
357  G4LorentzVector fourVector1( pP1Cms->GetMomentum(), 
358                               pP1Cms->GetTotalEnergy() );
359  G4LorentzVector fourVector2( pP2Cms->GetMomentum(), 
360                               pP2Cms->GetTotalEnergy() );
361
362  // Transform into lab frame by transforming the four vectors. Then
363  // add to the vector of secondary particles.
364  fourVector1.boost( boostToLab );
365  fourVector2.boost( boostToLab );
366  G4DynamicParticle * pP1Lab  = new G4DynamicParticle( G4Alpha::Alpha(),
367                                                      fourVector1.vect(),
368                                                      fourVector1.e() );
369  G4DynamicParticle * pP2Lab  = new G4DynamicParticle( G4Alpha::Alpha(),
370                                                      fourVector2.vect(),
371                                                      fourVector2.e() );
372  secondaryParticleVector.push_back( pP1Lab );
373  secondaryParticleVector.push_back( pP2Lab );
374
375  delete pP1Cms;
376  delete pP2Cms;
377
378  return;
379}
380
381
382void G4BertiniEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
383                                      G4FragmentVector * aResult )
384{
385  // Fill the vector pParticleChange with secondary particles stored in vector.
386  for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ )
387  {
388    G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
389    G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
390    G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
391    if(aA>0) 
392    {
393      aResult->push_back( new G4Fragment(aA, aZ, aMomentum) ); 
394    }
395    else 
396    {
397      aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) ); 
398    }
399  }
400  return;
401}
402
403
404void G4BertiniEvaporation::isotropicCosines( G4double & u, G4double & v, G4double & w )
405{
406  // Samples isotropic random direction cosines.
407  G4double CosTheta = 1.0 - 2.0 * G4UniformRand();
408  G4double SinTheta = std::sqrt( 1.0 - CosTheta * CosTheta );
409  G4double Phi = twopi * G4UniformRand();
410
411  u = std::cos( Phi ) * SinTheta;
412  v = std::cos( Phi ) * CosTheta,
413  w = std::sin( Phi );
414
415  return;
416}
Note: See TracBrowser for help on using the repository browser.