00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include <numeric>
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "G4HadronicException.hh"
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include "G4AblaEvaporation.hh"
00050
00051 #include "G4PionPlus.hh"
00052 #include "G4PionMinus.hh"
00053 #include "G4PionZero.hh"
00054
00055 G4AblaEvaporation::G4AblaEvaporation() {
00056 verboseLevel=0;
00057 hazard = new G4Hazard();
00058
00059
00060
00061
00062 hazard->ial = 979678188;
00063
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
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 }
00085
00086 G4AblaEvaporation::G4AblaEvaporation(const G4AblaEvaporation &) : G4VEvaporation() {
00087 throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::copy_constructor meant to not be accessable.");
00088 }
00089
00090 G4AblaEvaporation::~G4AblaEvaporation() {
00091 }
00092
00093 const G4AblaEvaporation & G4AblaEvaporation::operator=(const G4AblaEvaporation &) {
00094 throw G4HadronicException(__FILE__, __LINE__, "G4AblaEvaporation::operator= meant to not be accessable.");
00095 return *this;
00096 }
00097
00098 G4bool G4AblaEvaporation::operator==(const G4AblaEvaporation &) const {
00099 return false;
00100 }
00101
00102 G4bool G4AblaEvaporation::operator!=(const G4AblaEvaporation &) const {
00103 return true;
00104 }
00105
00106 void G4AblaEvaporation::setVerboseLevel( const G4int verbose ) {
00107 verboseLevel = verbose;
00108 }
00109
00110 G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) {
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;
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
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) {
00145 theResult->push_back(new G4Fragment(theNucleus));
00146 return theResult;
00147 }
00148
00149
00150 varntp->exini = exitationE;
00151
00152 G4int particleI, n = 0;
00153
00154
00155
00156
00157
00158 eventNumber++;
00159
00160 G4DynamicParticle *cascadeParticle = 0;
00161
00162
00163
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
00168
00169
00170
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
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
00194 G4double momx, momy, momz;
00195 G4double eKin;
00196 G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl;
00197 for(particleI = 0; particleI < varntp->ntrack; particleI++) {
00198
00199
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
00208
00209 }
00210
00211
00212 G4ThreeVector momDirection(momx, momy, momz);
00213 momDirection = momDirection.unit();
00214
00215
00216 G4int particleIdentified = 0;
00217
00218
00219 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 1)) {
00220 cascadeParticle =
00221 new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
00222 particleIdentified++;
00223 }
00224
00225
00226 if((varntp->avv[particleI] == 1) && (varntp->zvv[particleI] == 0)) {
00227 cascadeParticle =
00228 new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
00229 particleIdentified++;
00230 }
00231
00232
00233 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 1)) {
00234 cascadeParticle =
00235 new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
00236 particleIdentified++;
00237 }
00238
00239
00240 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == 0)) {
00241 cascadeParticle =
00242 new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
00243 particleIdentified++;
00244 }
00245
00246
00247 if((varntp->avv[particleI] == -1) && (varntp->zvv[particleI] == -1)) {
00248 cascadeParticle =
00249 new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
00250 particleIdentified++;
00251 }
00252
00253
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
00268 if(particleIdentified == 1) {
00269
00270 cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum());
00271 cascadeParticles.push_back(cascadeParticle);
00272
00273 }
00274
00275 else {
00276
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
00287
00288
00289
00290 varntp->ntrack = 0;
00291
00292
00293
00294 delete varntp;
00295 delete abla;
00296
00297 fillResult(cascadeParticles, theResult);
00298 return theResult;
00299 }
00300
00301 void G4AblaEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
00302 G4FragmentVector * aResult )
00303 {
00304
00305 G4cout <<"Size of the secondary particle vector = " << secondaryParticleVector.size() << G4endl;
00306 for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ ) {
00307 G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
00308 G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
00309 G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
00310 if(aA>0) {
00311 aResult->push_back( new G4Fragment(aA, aZ, aMomentum) );
00312 } else {
00313 aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) );
00314 }
00315 }
00316 return;
00317 }