source: trunk/source/processes/hadronic/models/incl/src/G4AblaEvaporation.cc @ 1358

Last change on this file since 1358 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 12.1 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: G4AblaEvaporation.cc,v 1.6 2010/10/26 02:47:59 kaitanie Exp $
27//
28#include <numeric>
29// #include "G4IonTable.hh"
30// #include "globals.hh"
31// #include "G4V3DNucleus.hh"
32// #include "G4DynamicParticleVector.hh"
33// #include "G4EvaporationInuclCollider.hh"
34// #include "G4InuclEvaporation.hh"
35// #include "G4InuclNuclei.hh"
36// #include "G4Track.hh"
37// #include "G4Nucleus.hh"
38// #include "G4Nucleon.hh"
39// #include "G4NucleiModel.hh"
40#include "G4HadronicException.hh"
41// #include "G4LorentzVector.hh"
42// #include "G4EquilibriumEvaporator.hh"
43// #include "G4Fissioner.hh"
44// #include "G4BigBanger.hh"
45// #include "G4InuclElementaryParticle.hh"
46// #include "G4InuclParticle.hh"
47// #include "G4CollisionOutput.hh"
48
49#include "G4AblaEvaporation.hh"
50
51#include "G4PionPlus.hh"
52#include "G4PionMinus.hh"
53#include "G4PionZero.hh"
54
55G4AblaEvaporation::G4AblaEvaporation() {
56  verboseLevel=0;
57 hazard = new G4Hazard();
58  // set initial values:
59  // First random seed:
60  // (Premiere graine)
61  //  hazard->ial = 38035;
62  hazard->ial = 979678188;
63  // other seeds:
64  hazard->igraine[0] = 3997;
65  hazard->igraine[1] = 15573;
66  hazard->igraine[2] = 9971;
67  hazard->igraine[3] = 9821; 
68  hazard->igraine[4] = 99233; 
69  hazard->igraine[5] = 11167; 
70  hazard->igraine[6] = 12399;
71  hazard->igraine[7] = 11321; 
72  hazard->igraine[8] = 9825;
73  hazard->igraine[9] = 2587; 
74  hazard->igraine[10] = 1775;
75  hazard->igraine[11] = 56799; 
76  hazard->igraine[12] = 1156;
77  //  hazard->igraine[13] = 11207;
78  hazard->igraine[13] = 38957; 
79  hazard->igraine[14] = 35779; 
80  hazard->igraine[15] = 10055; 
81  hazard->igraine[16] = 76533; 
82  hazard->igraine[17] = 33759;
83  hazard->igraine[18] = 13227;
84
85  G4VarNtp *evaporationResult = new G4VarNtp();
86  G4Volant *volant = new G4Volant();
87
88  // Initialize evaporation.
89  abla = new G4Abla(hazard, volant, evaporationResult);
90  abla->initEvapora();
91}
92
93G4AblaEvaporation::G4AblaEvaporation(const G4AblaEvaporation &) : G4VEvaporation() {
94  throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::copy_constructor meant to not be accessable.");
95}
96
97G4AblaEvaporation::~G4AblaEvaporation() {
98}
99
100const G4AblaEvaporation & G4AblaEvaporation::operator=(const G4AblaEvaporation &) {
101  throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::operator= meant to not be accessable.");
102  return *this;
103}
104
105G4bool G4AblaEvaporation::operator==(const G4AblaEvaporation &) const {
106  return false;
107}
108
109G4bool G4AblaEvaporation::operator!=(const G4AblaEvaporation &) const {
110  return true;
111}
112
113void G4AblaEvaporation::setVerboseLevel( const G4int verbose ) {
114  verboseLevel = verbose;
115}
116
117G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) { 
118  G4int nucleusA = theNucleus.GetA_asInt();
119  G4int nucleusZ = theNucleus.GetZ_asInt();
120  G4double nucleusMass = G4NucleiProperties::GetNuclearMass(nucleusA, nucleusZ);
121  G4double excitationEnergy = theNucleus.GetExcitationEnergy();
122  G4double angularMomentum = 0.0; // Don't know how to get this quantity... From Geant4???
123
124  G4LorentzVector tmp = theNucleus.GetMomentum();
125
126  G4ThreeVector momentum = tmp.vect();
127
128  G4double recoilEnergy = tmp.e();
129  G4double momX = momentum.x();
130  G4double momY = momentum.y();
131  G4double momZ = momentum.z();
132  //  G4double energy = tmp.e();
133  G4double exitationE = theNucleus.GetExcitationEnergy() * MeV;
134
135  varntp->ntrack = -1;
136  varntp->massini = theNucleus.GetA_asInt();
137  varntp->mzini = theNucleus.GetZ_asInt();
138
139  std::vector<G4DynamicParticle*> cascadeParticles;
140  G4FragmentVector * theResult = new G4FragmentVector;
141  if (theNucleus.GetExcitationEnergy() <= 0.0) { // Check that Excitation Energy > 0
142    theResult->push_back(new G4Fragment(theNucleus));
143    return theResult;
144  }
145
146  //  G4double mTar  = G4NucleiProperties::GetNuclearMass(A, Z); // Mass of the target nucleus
147  varntp->exini = exitationE;
148
149  G4int particleI, n = 0;
150
151  // Print diagnostic messages. 0 = silent, 1 and 2 = verbose
152  //  verboseLevel = 3;
153
154  // Increase the event number:
155  eventNumber++;
156
157  G4DynamicParticle *cascadeParticle = 0;
158  //  G4ParticleDefinition *aParticleDefinition = 0;
159
160  // Map Geant4 particle types to corresponding INCL4 types.
161  enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 
162                           pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9};
163
164  // Check wheter the input is acceptable. This will contain more tests in the future.
165
166//   void breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
167//                     G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ)
168  G4cout <<"Calling the actual ABLA model..." << G4endl;
169  G4cout <<"Excitation energy: " << excitationEnergy << G4endl;
170  abla->breakItUp(nucleusA, nucleusZ, nucleusMass, excitationEnergy, angularMomentum, recoilEnergy, momX, momY, momZ,
171                  eventNumber);
172  G4cout <<"Done." << G4endl;
173
174  if(verboseLevel > 0) {
175    // Diagnostic output
176    G4cout <<"G4AblaEvaporation: Target A:  " << nucleusA << G4endl;
177    G4cout <<"G4AblaEvaporation: Target Z:  " << nucleusZ << G4endl;
178
179    for(particleI = 0; particleI < varntp->ntrack; particleI++) {
180      G4cout << n << " ";
181      G4cout << varntp->massini << " " << varntp->mzini << " ";
182      G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";
183      G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";
184      G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " ";
185      G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";
186      G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
187    }
188  }
189
190  // Loop through the INCL4+ABLA output.
191  G4double momx, momy, momz; // Momentum components of the outcoming particles.
192  G4double eKin;
193  G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
194  for(particleI = 0; particleI < varntp->ntrack; particleI++) {
195    // Get energy/momentum and construct momentum vector:
196    // In INCL4 coordinates!
197    momx = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
198    momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
199    momz = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*MeV;
200
201    eKin = varntp->enerj[particleI] * MeV;
202
203    if(verboseLevel > 1) {
204      //      G4cout <<"Momentum direction: (x ,y,z)";
205      //      G4cout << "(" << momx <<"," << momy << "," << momz << ")" << G4endl;
206    }
207
208    // This vector tells the direction of the particle.
209    G4ThreeVector momDirection(momx, momy, momz);
210    momDirection = momDirection.unit();
211       
212    // Identify the particle/nucleus:
213    G4int particleIdentified = 0;
214
215    // Proton
216    if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) {
217      cascadeParticle = 
218        new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
219      particleIdentified++;
220    }
221
222    // Neutron
223    if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) {
224      cascadeParticle = 
225        new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
226      particleIdentified++;
227    }
228
229    // PionPlus
230    if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) {
231      cascadeParticle = 
232        new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
233      particleIdentified++;
234    }
235
236    // PionZero
237    if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) {
238      cascadeParticle = 
239        new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
240      particleIdentified++;
241    }
242
243    // PionMinus
244    if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) {
245      cascadeParticle = 
246        new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
247      particleIdentified++;
248    }
249
250    // Nuclei fragment
251    if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) {
252      G4ParticleDefinition * aIonDef = 0;
253      G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
254
255      G4int A = G4int(varntp->avv[particleI]);
256      G4int Z = G4int(varntp->zvv[particleI]);
257      aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
258       
259      cascadeParticle = 
260        new G4DynamicParticle(aIonDef, momDirection, eKin);
261      particleIdentified++;
262    }
263
264    // Check that the particle was identified properly.
265    if(particleIdentified == 1) {
266      // Put data into G4HadFinalState:
267      cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum());
268      cascadeParticles.push_back(cascadeParticle);
269      //      theResult.AddSecondary(cascadeParticle);
270    }
271    // Particle identification failed. Checking why...
272    else {
273      // Particle was identified as more than one particle type.
274      if(particleIdentified > 1) {
275        G4cout <<"G4InclCascadeInterface: One outcoming particle was identified as";
276        G4cout <<"more than one particle type. This is probably due to a bug in the interface." << G4endl;
277        G4cout <<"Particle A:" << varntp->avv[particleI] << "Z: " << varntp->zvv[particleI] << G4endl;
278        G4cout << "(particleIdentified =" << particleIdentified << ")"  << G4endl;
279      }
280    }
281  }
282
283  // End of conversion
284
285  // Clean up: Clean up the number of generated particles in the
286  // common block VARNTP_ for the processing of the next event.
287  varntp->ntrack = 0;
288  // End of cleanup.
289
290// Free allocated memory
291  delete varntp;
292  delete abla;
293 
294  fillResult(cascadeParticles, theResult);
295  return theResult;
296}
297
298void G4AblaEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
299                                     G4FragmentVector * aResult )
300{
301  // Fill the vector pParticleChange with secondary particles stored in vector.
302  G4cout <<"Size of the secondary particle vector = " << secondaryParticleVector.size() << G4endl;
303  for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ ) {
304      G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
305      G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
306      G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
307      if(aA>0) {
308        aResult->push_back( new G4Fragment(aA, aZ, aMomentum) ); 
309      } else {
310        aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) ); 
311      }
312    }
313  return;
314}
Note: See TracBrowser for help on using the repository browser.