#include <G4AblaEvaporation.hh>
Inheritance diagram for G4AblaEvaporation:
Public Member Functions | |
G4AblaEvaporation () | |
~G4AblaEvaporation () | |
G4FragmentVector * | BreakItUp (const G4Fragment &theNucleus) |
void | setVerboseLevel (const G4int verbose) |
Definition at line 48 of file G4AblaEvaporation.hh.
G4AblaEvaporation::G4AblaEvaporation | ( | ) |
Constructor.
Definition at line 55 of file G4AblaEvaporation.cc.
00055 { 00056 verboseLevel=0; 00057 hazard = new G4Hazard(); 00058 // set initial values: 00059 // First random seed: 00060 // (Premiere graine) 00061 // hazard->ial = 38035; 00062 hazard->ial = 979678188; 00063 // other seeds: 00064 hazard->igraine[0] = 3997; 00065 hazard->igraine[1] = 15573; 00066 hazard->igraine[2] = 9971; 00067 hazard->igraine[3] = 9821; 00068 hazard->igraine[4] = 99233; 00069 hazard->igraine[5] = 11167; 00070 hazard->igraine[6] = 12399; 00071 hazard->igraine[7] = 11321; 00072 hazard->igraine[8] = 9825; 00073 hazard->igraine[9] = 2587; 00074 hazard->igraine[10] = 1775; 00075 hazard->igraine[11] = 56799; 00076 hazard->igraine[12] = 1156; 00077 // hazard->igraine[13] = 11207; 00078 hazard->igraine[13] = 38957; 00079 hazard->igraine[14] = 35779; 00080 hazard->igraine[15] = 10055; 00081 hazard->igraine[16] = 76533; 00082 hazard->igraine[17] = 33759; 00083 hazard->igraine[18] = 13227; 00084 }
G4AblaEvaporation::~G4AblaEvaporation | ( | ) |
G4FragmentVector * G4AblaEvaporation::BreakItUp | ( | const G4Fragment & | theNucleus | ) | [virtual] |
The method for calling
Implements G4VEvaporation.
Definition at line 110 of file G4AblaEvaporation.cc.
References abla, G4InuclParticleNames::deuteron, G4ParticleTable::FindIon(), G4cout, G4endl, G4DynamicParticle::Get4Momentum(), G4Fragment::GetA(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetMomentum(), G4ParticleTable::GetParticleTable(), G4Fragment::GetZ(), CLHEP::detail::n, neutron, G4Neutron::NeutronDefinition(), G4INCL::Math::pi, G4InuclParticleNames::pionMinus, G4PionMinus::PionMinusDefinition(), G4InuclParticleNames::pionPlus, G4PionPlus::PionPlusDefinition(), G4InuclParticleNames::pionZero, G4PionZero::PionZeroDefinition(), G4InuclParticleNames::proton, G4Proton::ProtonDefinition(), G4DynamicParticle::Set4Momentum(), and G4InuclParticleNames::triton.
00110 { 00111 00112 00113 G4VarNtp *varntp = new G4VarNtp(); 00114 G4Volant *volant = new G4Volant(); 00115 00116 G4Abla *abla = new G4Abla(hazard, volant, varntp); 00117 G4cout <<"Initializing evaporation..." << G4endl; 00118 abla->initEvapora(); 00119 G4cout <<"Initialization complete!" << G4endl; 00120 00121 G4double nucleusA = theNucleus.GetA(); 00122 G4double nucleusZ = theNucleus.GetZ(); 00123 G4double nucleusMass = G4NucleiProperties::GetAtomicMass(nucleusA, nucleusZ); 00124 G4double excitationEnergy = theNucleus.GetExcitationEnergy(); 00125 G4double angularMomentum = 0.0; // Don't know how to get this quantity... From Geant4??? 00126 00127 G4LorentzVector tmp = theNucleus.GetMomentum(); 00128 00129 G4ThreeVector momentum = tmp.vect(); 00130 00131 G4double recoilEnergy = tmp.e(); 00132 G4double momX = momentum.x(); 00133 G4double momY = momentum.y(); 00134 G4double momZ = momentum.z(); 00135 // G4double energy = tmp.e(); 00136 G4double exitationE = theNucleus.GetExcitationEnergy() * MeV; 00137 00138 varntp->ntrack = -1; 00139 varntp->massini = theNucleus.GetA(); 00140 varntp->mzini = theNucleus.GetZ(); 00141 00142 std::vector<G4DynamicParticle*> cascadeParticles; 00143 G4FragmentVector * theResult = new G4FragmentVector; 00144 if (theNucleus.GetExcitationEnergy() <= 0.0) { // Check that Excitation Energy > 0 00145 theResult->push_back(new G4Fragment(theNucleus)); 00146 return theResult; 00147 } 00148 00149 // G4double mTar = G4NucleiProperties::GetAtomicMass(A, Z); // Mass of the target nucleus 00150 varntp->exini = exitationE; 00151 00152 G4int particleI, n = 0; 00153 00154 // Print diagnostic messages. 0 = silent, 1 and 2 = verbose 00155 // verboseLevel = 3; 00156 00157 // Increase the event number: 00158 eventNumber++; 00159 00160 G4DynamicParticle *cascadeParticle = 0; 00161 // G4ParticleDefinition *aParticleDefinition = 0; 00162 00163 // Map Geant4 particle types to corresponding INCL4 types. 00164 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 00165 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9}; 00166 00167 // Check wheter the input is acceptable. This will contain more tests in the future. 00168 00169 // void breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy, 00170 // G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ) 00171 G4cout <<"Calling the actual ABLA model..." << G4endl; 00172 G4cout <<"Excitation energy: " << excitationEnergy << G4endl; 00173 abla->breakItUp(nucleusA, nucleusZ, nucleusMass, excitationEnergy, angularMomentum, recoilEnergy, momX, momY, momZ, 00174 eventNumber); 00175 G4cout <<"Done." << G4endl; 00176 00177 if(verboseLevel > 0) { 00178 // Diagnostic output 00179 G4cout <<"G4AblaEvaporation: Target A: " << nucleusA << G4endl; 00180 G4cout <<"G4AblaEvaporation: Target Z: " << nucleusZ << G4endl; 00181 00182 for(particleI = 0; particleI < varntp->ntrack; particleI++) { 00183 G4cout << n << " "; 00184 G4cout << varntp->massini << " " << varntp->mzini << " "; 00185 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 00186 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 00187 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 00188 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 00189 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 00190 } 00191 } 00192 00193 // Loop through the INCL4+ABLA output. 00194 G4double momx, momy, momz; // Momentum components of the outcoming particles. 00195 G4double eKin; 00196 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl; 00197 for(particleI = 0; particleI < varntp->ntrack; particleI++) { 00198 // Get energy/momentum and construct momentum vector: 00199 // In INCL4 coordinates! 00200 momx = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; 00201 momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; 00202 momz = varntp->plab[particleI]*std::cos(varntp->tetlab[particleI]*CLHEP::pi/180.0)*MeV; 00203 00204 eKin = varntp->enerj[particleI] * MeV; 00205 00206 if(verboseLevel > 1) { 00207 // G4cout <<"Momentum direction: (x ,y,z)"; 00208 // G4cout << "(" << momx <<"," << momy << "," << momz << ")" << G4endl; 00209 } 00210 00211 // This vector tells the direction of the particle. 00212 G4ThreeVector momDirection(momx, momy, momz); 00213 momDirection = momDirection.unit(); 00214 00215 // Identify the particle/nucleus: 00216 G4int particleIdentified = 0; 00217 00218 // Proton 00219 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) { 00220 cascadeParticle = 00221 new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin); 00222 particleIdentified++; 00223 } 00224 00225 // Neutron 00226 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) { 00227 cascadeParticle = 00228 new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin); 00229 particleIdentified++; 00230 } 00231 00232 // PionPlus 00233 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) { 00234 cascadeParticle = 00235 new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin); 00236 particleIdentified++; 00237 } 00238 00239 // PionZero 00240 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) { 00241 cascadeParticle = 00242 new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin); 00243 particleIdentified++; 00244 } 00245 00246 // PionMinus 00247 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) { 00248 cascadeParticle = 00249 new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin); 00250 particleIdentified++; 00251 } 00252 00253 // Nuclei fragment 00254 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { 00255 G4ParticleDefinition * aIonDef = 0; 00256 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 00257 00258 G4int A = G4int(varntp->avv[particleI]); 00259 G4int Z = G4int(varntp->zvv[particleI]); 00260 aIonDef = theTableOfParticles->FindIon(Z, A, 0, Z); 00261 00262 cascadeParticle = 00263 new G4DynamicParticle(aIonDef, momDirection, eKin); 00264 particleIdentified++; 00265 } 00266 00267 // Check that the particle was identified properly. 00268 if(particleIdentified == 1) { 00269 // Put data into G4HadFinalState: 00270 cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()); 00271 cascadeParticles.push_back(cascadeParticle); 00272 // theResult.AddSecondary(cascadeParticle); 00273 } 00274 // Particle identification failed. Checking why... 00275 else { 00276 // Particle was identified as more than one particle type. 00277 if(particleIdentified > 1) { 00278 G4cout <<"G4InclCascadeInterface: One outcoming particle was identified as"; 00279 G4cout <<"more than one particle type. This is probably due to a bug in the interface." << G4endl; 00280 G4cout <<"Particle A:" << varntp->avv[particleI] << "Z: " << varntp->zvv[particleI] << G4endl; 00281 G4cout << "(particleIdentified =" << particleIdentified << ")" << G4endl; 00282 } 00283 } 00284 } 00285 00286 // End of conversion 00287 00288 // Clean up: Clean up the number of generated particles in the 00289 // common block VARNTP_ for the processing of the next event. 00290 varntp->ntrack = 0; 00291 // End of cleanup. 00292 00293 // Free allocated memory 00294 delete varntp; 00295 delete abla; 00296 00297 fillResult(cascadeParticles, theResult); 00298 return theResult; 00299 }
void G4AblaEvaporation::setVerboseLevel | ( | const G4int | verbose | ) |