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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

update processes
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.4 2008/10/24 21:07:40 dennis 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
86G4AblaEvaporation::G4AblaEvaporation(const G4AblaEvaporation &) : G4VEvaporation() {
87 throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::copy_constructor meant to not be accessable.");
88}
89
90G4AblaEvaporation::~G4AblaEvaporation() {
91}
92
93const G4AblaEvaporation & G4AblaEvaporation::operator=(const G4AblaEvaporation &) {
94 throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::operator= meant to not be accessable.");
95 return *this;
96}
97
98G4bool G4AblaEvaporation::operator==(const G4AblaEvaporation &) const {
99 return false;
100}
101
102G4bool G4AblaEvaporation::operator!=(const G4AblaEvaporation &) const {
103 return true;
104}
105
106void G4AblaEvaporation::setVerboseLevel( const G4int verbose ) {
107 verboseLevel = verbose;
108}
109
110G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) {
111
112
113 G4VarNtp *varntp = new G4VarNtp();
114 G4Volant *volant = new G4Volant();
115
116 G4Abla *abla = new G4Abla(hazard, volant, varntp);
117 G4cout <<"Initializing evaporation..." << G4endl;
118 abla->initEvapora();
119 G4cout <<"Initialization complete!" << G4endl;
120
121 G4double nucleusA = theNucleus.GetA();
122 G4double nucleusZ = theNucleus.GetZ();
123 G4double nucleusMass = G4NucleiProperties::GetNuclearMass(nucleusA, nucleusZ);
124 G4double excitationEnergy = theNucleus.GetExcitationEnergy();
125 G4double angularMomentum = 0.0; // Don't know how to get this quantity... From Geant4???
126
127 G4LorentzVector tmp = theNucleus.GetMomentum();
128
129 G4ThreeVector momentum = tmp.vect();
130
131 G4double recoilEnergy = tmp.e();
132 G4double momX = momentum.x();
133 G4double momY = momentum.y();
134 G4double momZ = momentum.z();
135 // G4double energy = tmp.e();
136 G4double exitationE = theNucleus.GetExcitationEnergy() * MeV;
137
138 varntp->ntrack = -1;
139 varntp->massini = theNucleus.GetA();
140 varntp->mzini = theNucleus.GetZ();
141
142 std::vector<G4DynamicParticle*> cascadeParticles;
143 G4FragmentVector * theResult = new G4FragmentVector;
144 if (theNucleus.GetExcitationEnergy() <= 0.0) { // Check that Excitation Energy > 0
145 theResult->push_back(new G4Fragment(theNucleus));
146 return theResult;
147 }
148
149 // G4double mTar = G4NucleiProperties::GetNuclearMass(A, Z); // Mass of the target nucleus
150 varntp->exini = exitationE;
151
152 G4int particleI, n = 0;
153
154 // Print diagnostic messages. 0 = silent, 1 and 2 = verbose
155 // verboseLevel = 3;
156
157 // Increase the event number:
158 eventNumber++;
159
160 G4DynamicParticle *cascadeParticle = 0;
161 // G4ParticleDefinition *aParticleDefinition = 0;
162
163 // Map Geant4 particle types to corresponding INCL4 types.
164 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4,
165 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9};
166
167 // Check wheter the input is acceptable. This will contain more tests in the future.
168
169// void breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
170// G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ)
171 G4cout <<"Calling the actual ABLA model..." << G4endl;
172 G4cout <<"Excitation energy: " << excitationEnergy << G4endl;
173 abla->breakItUp(nucleusA, nucleusZ, nucleusMass, excitationEnergy, angularMomentum, recoilEnergy, momX, momY, momZ,
174 eventNumber);
175 G4cout <<"Done." << G4endl;
176
177 if(verboseLevel > 0) {
178 // Diagnostic output
179 G4cout <<"G4AblaEvaporation: Target A: " << nucleusA << G4endl;
180 G4cout <<"G4AblaEvaporation: Target Z: " << nucleusZ << G4endl;
181
182 for(particleI = 0; particleI < varntp->ntrack; particleI++) {
183 G4cout << n << " ";
184 G4cout << varntp->massini << " " << varntp->mzini << " ";
185 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";
186 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";
187 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " ";
188 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";
189 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
190 }
191 }
192
193 // Loop through the INCL4+ABLA output.
194 G4double momx, momy, momz; // Momentum components of the outcoming particles.
195 G4double eKin;
196 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
197 for(particleI = 0; particleI < varntp->ntrack; particleI++) {
198 // Get energy/momentum and construct momentum vector:
199 // In INCL4 coordinates!
200 momx = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
201 momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
202 momz = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*MeV;
203
204 eKin = varntp->enerj[particleI] * MeV;
205
206 if(verboseLevel > 1) {
207 // G4cout <<"Momentum direction: (x ,y,z)";
208 // G4cout << "(" << momx <<"," << momy << "," << momz << ")" << G4endl;
209 }
210
211 // This vector tells the direction of the particle.
212 G4ThreeVector momDirection(momx, momy, momz);
213 momDirection = momDirection.unit();
214
215 // Identify the particle/nucleus:
216 G4int particleIdentified = 0;
217
218 // Proton
219 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) {
220 cascadeParticle =
221 new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
222 particleIdentified++;
223 }
224
225 // Neutron
226 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) {
227 cascadeParticle =
228 new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
229 particleIdentified++;
230 }
231
232 // PionPlus
233 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) {
234 cascadeParticle =
235 new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
236 particleIdentified++;
237 }
238
239 // PionZero
240 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) {
241 cascadeParticle =
242 new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
243 particleIdentified++;
244 }
245
246 // PionMinus
247 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) {
248 cascadeParticle =
249 new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
250 particleIdentified++;
251 }
252
253 // Nuclei fragment
254 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) {
255 G4ParticleDefinition * aIonDef = 0;
256 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
257
258 G4int A = G4int(varntp->avv[particleI]);
259 G4int Z = G4int(varntp->zvv[particleI]);
260 aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z);
261
262 cascadeParticle =
263 new G4DynamicParticle(aIonDef, momDirection, eKin);
264 particleIdentified++;
265 }
266
267 // Check that the particle was identified properly.
268 if(particleIdentified == 1) {
269 // Put data into G4HadFinalState:
270 cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum());
271 cascadeParticles.push_back(cascadeParticle);
272 // theResult.AddSecondary(cascadeParticle);
273 }
274 // Particle identification failed. Checking why...
275 else {
276 // Particle was identified as more than one particle type.
277 if(particleIdentified > 1) {
278 G4cout <<"G4InclCascadeInterface: One outcoming particle was identified as";
279 G4cout <<"more than one particle type. This is probably due to a bug in the interface." << G4endl;
280 G4cout <<"Particle A:" << varntp->avv[particleI] << "Z: " << varntp->zvv[particleI] << G4endl;
281 G4cout << "(particleIdentified =" << particleIdentified << ")" << G4endl;
282 }
283 }
284 }
285
286 // End of conversion
287
288 // Clean up: Clean up the number of generated particles in the
289 // common block VARNTP_ for the processing of the next event.
290 varntp->ntrack = 0;
291 // End of cleanup.
292
293// Free allocated memory
294 delete varntp;
295 delete abla;
296
297 fillResult(cascadeParticles, theResult);
298 return theResult;
299}
300
301void G4AblaEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
302 G4FragmentVector * aResult )
303{
304 // Fill the vector pParticleChange with secondary particles stored in vector.
305 G4cout <<"Size of the secondary particle vector = " << secondaryParticleVector.size() << G4endl;
306 for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ ) {
307 G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
308 G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
309 G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
310 if(aA>0) {
311 aResult->push_back( new G4Fragment(aA, aZ, aMomentum) );
312 } else {
313 aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) );
314 }
315 }
316 return;
317}
Note: See TracBrowser for help on using the repository browser.