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