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
00040 #include "G4ios.hh"
00041 #include "globals.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4DecayProducts.hh"
00045
00046 #include "G4LorentzVector.hh"
00047 #include "G4LorentzRotation.hh"
00048
00049
00050 G4DecayProducts::G4DecayProducts()
00051 :numberOfProducts(0),theParentParticle(0)
00052 {
00053 theProductVector = new G4DecayProductVector();
00054 }
00055
00056 G4DecayProducts::G4DecayProducts(const G4DynamicParticle &aParticle)
00057 :numberOfProducts(0),theParentParticle(0)
00058 {
00059 theParentParticle = new G4DynamicParticle(aParticle);
00060 theProductVector = new G4DecayProductVector();
00061 }
00062
00063 G4DecayProducts::G4DecayProducts(const G4DecayProducts &right)
00064 :numberOfProducts(0)
00065 {
00066 theProductVector = new G4DecayProductVector();
00067
00068
00069 theParentParticle = new G4DynamicParticle(*right.theParentParticle);
00070
00071
00072 for (G4int index=0; index < right.numberOfProducts; index++) {
00073 G4DynamicParticle* daughter = right.theProductVector->at(index);
00074 G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
00075
00076 G4double properTime = daughter->GetPreAssignedDecayProperTime();
00077 if(properTime>0.0)pDaughter->SetPreAssignedDecayProperTime(properTime);
00078
00079 const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
00080 if (pPreAssigned) {
00081 G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
00082 pDaughter->SetPreAssignedDecayProducts(pPA);
00083 }
00084
00085 theProductVector->push_back( pDaughter );
00086 }
00087 numberOfProducts = right.numberOfProducts;
00088 }
00089
00090 G4DecayProducts & G4DecayProducts::operator=(const G4DecayProducts &right)
00091 {
00092 G4int index;
00093
00094 if (this != &right)
00095 {
00096
00097 if (theParentParticle != 0) delete theParentParticle;
00098 theParentParticle = new G4DynamicParticle(*right.theParentParticle);
00099
00100
00101 for (index=0; index < numberOfProducts; index++) {
00102 delete theProductVector->at(index);
00103 }
00104 theProductVector->clear();
00105
00106
00107 for (index=0; index < right.numberOfProducts; index++) {
00108 G4DynamicParticle* daughter = right.theProductVector->at(index);
00109 G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
00110
00111 G4double properTime = daughter->GetPreAssignedDecayProperTime();
00112 if(properTime>0.0) pDaughter->SetPreAssignedDecayProperTime(properTime);
00113
00114 const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
00115 if (pPreAssigned) {
00116 G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
00117 pDaughter->SetPreAssignedDecayProducts(pPA);
00118 }
00119 theProductVector->push_back( pDaughter );
00120 }
00121 numberOfProducts = right.numberOfProducts;
00122
00123 }
00124 return *this;
00125 }
00126
00127 G4DecayProducts::~G4DecayProducts()
00128 {
00129
00130 if (theParentParticle != 0) delete theParentParticle;
00131
00132
00133 for (G4int index=0; index < numberOfProducts; index++) {
00134 delete theProductVector->at(index);
00135 }
00136 theProductVector->clear();
00137 numberOfProducts = 0;
00138 delete theProductVector;
00139 }
00140
00141 G4DynamicParticle* G4DecayProducts::PopProducts()
00142 {
00143 if ( numberOfProducts >0 ) {
00144 numberOfProducts -= 1;
00145 G4DynamicParticle* part = theProductVector->back();
00146 theProductVector->pop_back();
00147 return part;
00148 } else {
00149 return 0;
00150 }
00151 }
00152
00153 G4int G4DecayProducts::PushProducts(G4DynamicParticle *aParticle)
00154 {
00155 theProductVector->push_back(aParticle);
00156 numberOfProducts += 1;
00157 return numberOfProducts;
00158 }
00159
00160 G4DynamicParticle* G4DecayProducts::operator[](G4int anIndex) const
00161 {
00162 if ((numberOfProducts > anIndex) && (anIndex >=0) ) {
00163 return theProductVector->at(anIndex);
00164 } else {
00165 return 0;
00166 }
00167 }
00168
00169 void G4DecayProducts::SetParentParticle(const G4DynamicParticle &aParticle)
00170 {
00171 if (theParentParticle != 0) delete theParentParticle;
00172 theParentParticle = new G4DynamicParticle(aParticle);
00173 }
00174
00175
00176 void G4DecayProducts::Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
00177 {
00178
00179 G4double mass = theParentParticle->GetMass();
00180 G4double totalMomentum(0);
00181 if (totalEnergy > mass ) totalMomentum = std::sqrt( (totalEnergy - mass)*(totalEnergy + mass) );
00182 G4double betax = momentumDirection.x()*totalMomentum/totalEnergy;
00183 G4double betay = momentumDirection.y()*totalMomentum/totalEnergy;
00184 G4double betaz = momentumDirection.z()*totalMomentum/totalEnergy;
00185 this->Boost(betax, betay, betaz);
00186 }
00187
00188 void G4DecayProducts::Boost(G4double newbetax, G4double newbetay, G4double newbetaz)
00189 {
00190 G4double mass = theParentParticle->GetMass();
00191 G4double energy = theParentParticle->GetTotalEnergy();
00192 G4double momentum = 0.0;
00193
00194 G4ThreeVector direction(0.0,0.0,1.0);
00195 G4LorentzVector p4;
00196
00197 if (energy - mass > DBL_MIN) {
00198
00199 momentum = theParentParticle->GetTotalMomentum();
00200 direction = theParentParticle->GetMomentumDirection();
00201 G4double betax = -1.0*direction.x()*momentum/energy;
00202 G4double betay = -1.0*direction.y()*momentum/energy;
00203 G4double betaz = -1.0*direction.z()*momentum/energy;
00204
00205 for (G4int index=0; index < numberOfProducts; index++) {
00206
00207 p4 = (theProductVector->at(index))->Get4Momentum();
00208
00209
00210 p4.boost(betax, betay, betaz);
00211
00212
00213 p4.boost(newbetax, newbetay, newbetaz);
00214
00215
00216 (theProductVector->at(index))->Set4Momentum(p4);
00217 }
00218 } else {
00219 for (G4int index=0; index < numberOfProducts; index++) {
00220
00221 p4 = (theProductVector->at(index))->Get4Momentum();
00222
00223
00224 p4.boost(newbetax, newbetay, newbetaz);
00225
00226
00227 (theProductVector->at(index))->Set4Momentum(p4);
00228 }
00229 }
00230
00231 mass = theParentParticle->GetMass();
00232 G4LorentzVector parent4( 0.0, 0.0, 0.0, mass);
00233
00234
00235 parent4.boost(newbetax, newbetay, newbetaz);
00236
00237
00238 theParentParticle->Set4Momentum(parent4);
00239 }
00240
00241 G4bool G4DecayProducts::IsChecked() const
00242 {
00243 G4bool returnValue = true;
00244
00245
00246 G4double parent_energy = theParentParticle->GetTotalEnergy();
00247 G4ThreeVector direction = theParentParticle->GetMomentumDirection();
00248 G4ThreeVector parent_momentum = direction*(theParentParticle->GetTotalMomentum());
00249
00250 if ( (parent_momentum.mag() >0.0) && (std::fabs(direction.mag()-1.0) >1.0e-6 ) ) {
00251 #ifdef G4VERBOSE
00252 G4cerr << "G4DecayProducts::IsChecked():: "
00253 << " Momentum Direction Vector of Parent is not normalized "
00254 << " (=" << direction.mag() << ")" << G4endl;
00255 #endif
00256 returnValue = false;
00257 parent_momentum = parent_momentum * (1./direction.mag());
00258 }
00259
00260
00261 G4double mass, energy;
00262 G4ThreeVector momentum;
00263 G4double total_energy = parent_energy;
00264 G4ThreeVector total_momentum = parent_momentum;
00265 for (G4int index=0; index < numberOfProducts; index++)
00266 {
00267 G4DynamicParticle* part = theProductVector->at(index);
00268 mass = part->GetMass();
00269 energy = part->GetTotalEnergy();
00270 direction = part->GetMomentumDirection();
00271 momentum = direction*(part->GetTotalMomentum());
00272
00273 if ( (momentum.mag()>0.0) && (std::fabs(direction.mag()-1.0) > 1.0e-6)) {
00274 #ifdef G4VERBOSE
00275 G4cerr << "G4DecayProducts::IsChecked():: "
00276 << " Momentum Direction Vector of Daughter [" << index
00277 << "] is not normalized (=" << direction.mag() << ")" << G4endl;
00278 #endif
00279 returnValue = false;
00280 momentum = momentum * (1./direction.mag());
00281 }
00282
00283 if (energy - mass < DBL_MIN ) {
00284 #ifdef G4VERBOSE
00285 G4cerr << "G4DecayProducts::IsChecked():: "
00286 << " Daughter [" << index << "] has no kinetic energy "<< G4endl;
00287 #endif
00288 returnValue = false;
00289 }
00290 total_energy -= energy;
00291 total_momentum -= momentum;
00292 }
00293
00294 if ( (std::fabs(total_energy) >1.0e-9*MeV) || (total_momentum.mag() >1.0e-9*MeV ) ){
00295 #ifdef G4VERBOSE
00296 G4cerr << "G4DecayProducts::IsChecked():: "
00297 << " Energy/Momentum is not conserved "<< G4endl;
00298 G4cerr << " difference between parent energy and sum of dughters' energy : "
00299 << total_energy /MeV << "[MeV] " << G4endl;
00300 G4cerr << " difference between parent momentum and sum of dughters' momentum : "
00301 << " x:" << total_momentum.getX()/MeV
00302 << " y:" << total_momentum.getY()/MeV
00303 << " z:" << total_momentum.getZ()/MeV
00304 << G4endl;
00305 #endif
00306 returnValue = false;
00307 }
00308 return returnValue;
00309 }
00310
00311 void G4DecayProducts::DumpInfo() const
00312 {
00313 G4cout << " ----- List of DecayProducts -----" << G4endl;
00314 G4cout << " ------ Parent Particle ----------" << G4endl;
00315 if (theParentParticle != 0) theParentParticle->DumpInfo();
00316 G4cout << " ------ Daughter Particles ------" << G4endl;
00317 for (G4int index=0; index < numberOfProducts; index++)
00318 {
00319 G4cout << " ----------" << index+1 << " -------------" << G4endl;
00320 (theProductVector->at(index))-> DumpInfo();
00321 }
00322 G4cout << " ----- End List of DecayProducts -----" << G4endl;
00323 G4cout << G4endl;
00324 }