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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DecayProducts.hh"
00043 #include "G4VDecayChannel.hh"
00044 #include "G4PhaseSpaceDecayChannel.hh"
00045 #include "Randomize.hh"
00046 #include "G4LorentzVector.hh"
00047 #include "G4LorentzRotation.hh"
00048
00049
00050 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
00051 :G4VDecayChannel("Phase Space", Verbose)
00052 {
00053
00054 }
00055
00056 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
00057 const G4String& theParentName,
00058 G4double theBR,
00059 G4int theNumberOfDaughters,
00060 const G4String& theDaughterName1,
00061 const G4String& theDaughterName2,
00062 const G4String& theDaughterName3,
00063 const G4String& theDaughterName4 )
00064 :G4VDecayChannel("Phase Space",
00065 theParentName,theBR,
00066 theNumberOfDaughters,
00067 theDaughterName1,
00068 theDaughterName2,
00069 theDaughterName3,
00070 theDaughterName4),
00071 current_parent_mass(0.0)
00072 {
00073
00074 }
00075
00076 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
00077 {
00078 }
00079
00080 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
00081 {
00082 #ifdef G4VERBOSE
00083 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00084 #endif
00085
00086 G4DecayProducts * products = 0;
00087
00088 if (parent == 0) FillParent();
00089 if (daughters == 0) FillDaughters();
00090
00091 if (parentMass >0.0) current_parent_mass = parentMass;
00092 else current_parent_mass = parent_mass;
00093
00094 switch (numberOfDaughters){
00095 case 0:
00096 #ifdef G4VERBOSE
00097 if (GetVerboseLevel()>0) {
00098 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00099 G4cout << " daughters not defined " <<G4endl;
00100 }
00101 #endif
00102 break;
00103 case 1:
00104 products = OneBodyDecayIt();
00105 break;
00106 case 2:
00107 products = TwoBodyDecayIt();
00108 break;
00109 case 3:
00110 products = ThreeBodyDecayIt();
00111 break;
00112 default:
00113 products = ManyBodyDecayIt();
00114 break;
00115 }
00116 #ifdef G4VERBOSE
00117 if ((products == 0) && (GetVerboseLevel()>0)) {
00118 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00119 G4cout << *parent_name << " can not decay " << G4endl;
00120 DumpInfo();
00121 }
00122 #endif
00123 return products;
00124 }
00125
00126 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
00127 {
00128 #ifdef G4VERBOSE
00129 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
00130 #endif
00131
00132
00133 G4ThreeVector dummy;
00134 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00135
00136 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00137 delete parentparticle;
00138
00139
00140 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
00141 products->PushProducts(daughterparticle);
00142
00143 #ifdef G4VERBOSE
00144 if (GetVerboseLevel()>1) {
00145 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
00146 G4cout << " create decay products in rest frame " <<G4endl;
00147 products->DumpInfo();
00148 }
00149 #endif
00150 return products;
00151 }
00152
00153 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
00154 {
00155 #ifdef G4VERBOSE
00156 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
00157 #endif
00158
00159 G4double parentmass = current_parent_mass;
00160
00161
00162 G4double daughtermass[2];
00163 G4double daughtermomentum;
00164 daughtermass[0] = daughters_mass[0];
00165 daughtermass[1] = daughters_mass[1];
00166
00167
00168 G4ThreeVector dummy;
00169 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00170
00171 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00172 delete parentparticle;
00173
00174
00175 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
00176 if (daughtermomentum <0.0) {
00177 #ifdef G4VERBOSE
00178 if (GetVerboseLevel()>0) {
00179 G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
00180 << "sum of daughter mass is larger than parent mass" << G4endl;
00181 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
00182 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
00183 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
00184 }
00185 #endif
00186 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
00187 "PART112", JustWarning,
00188 "Can not create decay products: sum of daughter mass is larger than parent mass");
00189 return products;
00190 }
00191
00192 G4double costheta = 2.*G4UniformRand()-1.0;
00193 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00194 G4double phi = twopi*G4UniformRand()*rad;
00195 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00196
00197
00198 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
00199 products->PushProducts(daughterparticle);
00200 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
00201 products->PushProducts(daughterparticle);
00202
00203 #ifdef G4VERBOSE
00204 if (GetVerboseLevel()>1) {
00205 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
00206 G4cout << " create decay products in rest frame " <<G4endl;
00207 products->DumpInfo();
00208 }
00209 #endif
00210 return products;
00211 }
00212
00213 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
00214
00215 {
00216 #ifdef G4VERBOSE
00217 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
00218 #endif
00219
00220 G4double parentmass = current_parent_mass;
00221
00222 G4double daughtermass[3];
00223 G4double sumofdaughtermass = 0.0;
00224 for (G4int index=0; index<3; index++){
00225 daughtermass[index] = daughters_mass[index];
00226 sumofdaughtermass += daughtermass[index];
00227 }
00228
00229
00230 G4ThreeVector dummy;
00231 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00232
00233
00234
00235 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00236 delete parentparticle;
00237
00238 if (sumofdaughtermass >parentmass) {
00239 #ifdef G4VERBOSE
00240 if (GetVerboseLevel()>0) {
00241 G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
00242 << "sum of daughter mass is larger than parent mass" << G4endl;
00243 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
00244 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
00245 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
00246 G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
00247 }
00248 #endif
00249 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
00250 "PART112", JustWarning,
00251 "Can not create decay products: sum of daughter mass is larger than parent mass");
00252 return products;
00253 }
00254
00255
00256
00257
00258
00259 G4double rd1, rd2, rd;
00260 G4double daughtermomentum[3];
00261 G4double momentummax=0.0, momentumsum = 0.0;
00262 G4double energy;
00263
00264 do {
00265 rd1 = G4UniformRand();
00266 rd2 = G4UniformRand();
00267 if (rd2 > rd1) {
00268 rd = rd1;
00269 rd1 = rd2;
00270 rd2 = rd;
00271 }
00272 momentummax = 0.0;
00273 momentumsum = 0.0;
00274
00275 energy = rd2*(parentmass - sumofdaughtermass);
00276 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
00277 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
00278 momentumsum += daughtermomentum[0];
00279
00280 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
00281 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
00282 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
00283 momentumsum += daughtermomentum[1];
00284
00285 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
00286 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
00287 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
00288 momentumsum += daughtermomentum[2];
00289 } while (momentummax > momentumsum - momentummax );
00290
00291 #ifdef G4VERBOSE
00292 if (GetVerboseLevel()>1) {
00293 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
00294 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
00295 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
00296 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
00297 }
00298 #endif
00299
00300
00301 G4double costheta, sintheta, phi, sinphi, cosphi;
00302 G4double costhetan, sinthetan, phin, sinphin, cosphin;
00303 costheta = 2.*G4UniformRand()-1.0;
00304 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00305 phi = twopi*G4UniformRand()*rad;
00306 sinphi = std::sin(phi);
00307 cosphi = std::cos(phi);
00308 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
00309 G4DynamicParticle * daughterparticle
00310 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
00311 products->PushProducts(daughterparticle);
00312
00313 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
00314 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00315 phin = twopi*G4UniformRand()*rad;
00316 sinphin = std::sin(phin);
00317 cosphin = std::cos(phin);
00318 G4ThreeVector direction2;
00319 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
00320 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
00321 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
00322 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
00323 products->PushProducts(daughterparticle);
00324
00325 daughterparticle =
00326 new G4DynamicParticle(
00327 daughters[1],
00328 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)
00329 );
00330 products->PushProducts(daughterparticle);
00331
00332 #ifdef G4VERBOSE
00333 if (GetVerboseLevel()>1) {
00334 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
00335 G4cout << " create decay products in rest frame " <<G4endl;
00336 products->DumpInfo();
00337 }
00338 #endif
00339 return products;
00340 }
00341
00342 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 {
00353 G4int index, index2;
00354
00355
00356 G4DecayProducts *products;
00357
00358 #ifdef G4VERBOSE
00359 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
00360 #endif
00361
00362 G4double parentmass = current_parent_mass;
00363
00364 G4double *daughtermass = new G4double[numberOfDaughters];
00365 G4double sumofdaughtermass = 0.0;
00366 for (index=0; index<numberOfDaughters; index++){
00367 daughtermass[index] = daughters_mass[index];
00368 sumofdaughtermass += daughtermass[index];
00369 }
00370
00371
00372 G4double *daughtermomentum = new G4double[numberOfDaughters];
00373 G4ThreeVector direction;
00374 G4DynamicParticle **daughterparticle;
00375 G4double *sm = new G4double[numberOfDaughters];
00376 G4double tmas;
00377 G4double weight = 1.0;
00378 G4int numberOfTry = 0;
00379
00380 do {
00381
00382 G4double temp;
00383 G4double *rd = new G4double[numberOfDaughters];
00384 rd[0] = 1.0;
00385 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
00386 rd[ numberOfDaughters -1] = 0.0;
00387 for(index =1; index < numberOfDaughters -1; index++) {
00388 for(index2 = index+1; index2 < numberOfDaughters; index2++) {
00389 if (rd[index] < rd[index2]){
00390 temp = rd[index];
00391 rd[index] = rd[index2];
00392 rd[index2] = temp;
00393 }
00394 }
00395 }
00396
00397
00398 tmas = parentmass - sumofdaughtermass;
00399 temp = sumofdaughtermass;
00400 for(index =0; index < numberOfDaughters; index++) {
00401 sm[index] = rd[index]*tmas + temp;
00402 temp -= daughtermass[index];
00403 if (GetVerboseLevel()>1) {
00404 G4cout << index << " rundom number:" << rd[index];
00405 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
00406 }
00407 }
00408 delete [] rd;
00409
00410
00411 weight = 1.0;
00412 index =numberOfDaughters-1;
00413 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
00414 #ifdef G4VERBOSE
00415 if (GetVerboseLevel()>1) {
00416 G4cout << " daughter " << index << ":" << *daughters_name[index];
00417 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00418 }
00419 #endif
00420 for(index =numberOfDaughters-2; index>=0; index--) {
00421
00422 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
00423 if(daughtermomentum[index] < 0.0) {
00424
00425 #ifdef G4VERBOSE
00426 if (GetVerboseLevel()>0) {
00427 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
00428 G4cout << " can not calculate daughter momentum " <<G4endl;
00429 G4cout << " parent:" << *parent_name;
00430 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
00431 G4cout << " daughter " << index << ":" << *daughters_name[index];
00432 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
00433 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00434 }
00435 #endif
00436 delete [] sm;
00437 delete [] daughtermass;
00438 delete [] daughtermomentum;
00439
00440 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
00441 "PART112", JustWarning,
00442 "Can not create decay products: sum of daughter mass is larger than parent mass");
00443
00444 return 0;
00445
00446 } else {
00447
00448 weight *= daughtermomentum[index]/sm[index];
00449 #ifdef G4VERBOSE
00450 if (GetVerboseLevel()>1) {
00451 G4cout << " daughter " << index << ":" << *daughters_name[index];
00452 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00453 }
00454 #endif
00455 }
00456 }
00457
00458 #ifdef G4VERBOSE
00459 if (GetVerboseLevel()>1) {
00460 G4cout << " weight: " << weight <<G4endl;
00461 }
00462 #endif
00463
00464
00465 if (numberOfTry++ >100) {
00466 #ifdef G4VERBOSE
00467 if (GetVerboseLevel()>0) {
00468 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
00469 G4cout << " can not determine Decay Kinematics " << G4endl;
00470 G4cout << "parent : " << *parent_name << G4endl;
00471 G4cout << "daughters : ";
00472 for(index=0; index<numberOfDaughters; index++) {
00473 G4cout << *daughters_name[index] << " , ";
00474 }
00475 G4cout << G4endl;
00476 }
00477 #endif
00478 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
00479 "PART113", JustWarning,
00480 " Cannot decay : Decay Kinematics cannot be calculated ");
00481
00482 delete [] sm;
00483 delete [] daughtermass;
00484 delete [] daughtermomentum;
00485 return 0;
00486 }
00487 } while ( weight > G4UniformRand());
00488
00489 #ifdef G4VERBOSE
00490 if (GetVerboseLevel()>1) {
00491 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
00492 }
00493 #endif
00494
00495 G4double costheta, sintheta, phi;
00496 G4double beta;
00497 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
00498
00499 index = numberOfDaughters -2;
00500 costheta = 2.*G4UniformRand()-1.0;
00501 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00502 phi = twopi*G4UniformRand()*rad;
00503 direction.setZ(costheta);
00504 direction.setY(sintheta*std::sin(phi));
00505 direction.setX(sintheta*std::cos(phi));
00506 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
00507 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
00508
00509 for (index = numberOfDaughters -3; index >= 0; index--) {
00510
00511 costheta = 2.*G4UniformRand()-1.0;
00512 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00513 phi = twopi*G4UniformRand()*rad;
00514 direction.setZ(costheta);
00515 direction.setY(sintheta*std::sin(phi));
00516 direction.setX(sintheta*std::cos(phi));
00517
00518
00519 beta = daughtermomentum[index];
00520 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
00521 for (index2 = index+1; index2<numberOfDaughters; index2++) {
00522 G4LorentzVector p4;
00523
00524 p4 = daughterparticle[index2]->Get4Momentum();
00525
00526
00527 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
00528
00529
00530 daughterparticle[index2]->Set4Momentum(p4);
00531 }
00532
00533 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
00534 }
00535
00536
00537 G4DynamicParticle *parentparticle;
00538 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
00539 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
00540 products = new G4DecayProducts(*parentparticle);
00541 delete parentparticle;
00542 for (index = 0; index<numberOfDaughters; index++) {
00543 products->PushProducts(daughterparticle[index]);
00544 }
00545
00546 #ifdef G4VERBOSE
00547 if (GetVerboseLevel()>1) {
00548 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
00549 G4cout << " create decay products in rest frame " << G4endl;
00550 products->DumpInfo();
00551 }
00552 #endif
00553
00554 delete [] daughterparticle;
00555 delete [] daughtermomentum;
00556 delete [] daughtermass;
00557 delete [] sm;
00558
00559 return products;
00560 }
00561
00562
00563 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
00564 {
00565
00566 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
00567 if (ppp>0) return std::sqrt(ppp);
00568 else return -1.;
00569 }
00570
00571
00572