Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4DecayProducts class implementation << 27 // 23 // 28 // Author: H.Kurashige, 12 July 1996 << 24 // $Id: G4DecayProducts.cc,v 1.9 2001/07/11 10:01:59 gunter Exp $ >> 25 // GEANT4 tag $Name: geant4-04-01 $ >> 26 // >> 27 // >> 28 // ------------------------------------------------------------ >> 29 // GEANT 4 class implementation file >> 30 // >> 31 // History: first implementation, based on object model of >> 32 // 10 July 1996 H.Kurashige >> 33 // 21 Oct 1996 H.Kurashige >> 34 // 12 Dec 1997 H.Kurashige 29 // ------------------------------------------- 35 // ------------------------------------------------------------ 30 36 >> 37 #include "G4ios.hh" >> 38 #include "globals.hh" 31 #include "G4DecayProducts.hh" 39 #include "G4DecayProducts.hh" 32 40 33 #include "G4LorentzRotation.hh" << 34 #include "G4LorentzVector.hh" 41 #include "G4LorentzVector.hh" 35 #include "G4PhysicalConstants.hh" << 42 #include "G4LorentzRotation.hh" 36 #include "G4SystemOfUnits.hh" << 43 37 #include "G4ios.hh" << 44 G4Allocator<G4DecayProducts> aDecayProductsAllocator; 38 #include "globals.hh" << 45 39 46 40 G4DecayProducts::G4DecayProducts() 47 G4DecayProducts::G4DecayProducts() 41 { << 48 :numberOfProducts(0),theParentParticle(0) 42 theProductVector = new G4DecayProductVector( << 49 { >> 50 43 } 51 } 44 52 45 G4DecayProducts::G4DecayProducts(const G4Dynam << 53 G4DecayProducts::G4DecayProducts(const G4DynamicParticle &aParticle) >> 54 :numberOfProducts(0),theParentParticle(0) 46 { 55 { 47 theParentParticle = new G4DynamicParticle(aP 56 theParentParticle = new G4DynamicParticle(aParticle); 48 theProductVector = new G4DecayProductVector( << 49 } 57 } 50 58 51 G4DecayProducts::G4DecayProducts(const G4Decay << 59 G4DecayProducts::G4DecayProducts(const G4DecayProducts &right) >> 60 :numberOfProducts(0) 52 { 61 { 53 theProductVector = new G4DecayProductVector( << 54 << 55 // copy parent (Deep Copy) 62 // copy parent (Deep Copy) 56 theParentParticle = new G4DynamicParticle(*r 63 theParentParticle = new G4DynamicParticle(*right.theParentParticle); 57 64 58 // copy daughters (Deep Copy) << 65 //copy daughters (Deep Copy) 59 for (G4int index = 0; index < right.numberOf << 66 for (G4int index=0; index < right.numberOfProducts; index++) 60 G4DynamicParticle* daughter = right.thePro << 67 { 61 auto pDaughter = new G4DynamicParticle(*da << 68 G4DynamicParticle* daughter = right.theProductVector[index]; 62 << 63 G4double properTime = daughter->GetPreAssi << 64 if (properTime > 0.0) pDaughter->SetPreAss << 65 << 66 const G4DecayProducts* pPreAssigned = daug 69 const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts(); 67 if (pPreAssigned != nullptr) { << 70 G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter); 68 auto pPA = new G4DecayProducts(*pPreAssi << 71 >> 72 if (pPreAssigned) { >> 73 G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned); 69 pDaughter->SetPreAssignedDecayProducts(p 74 pDaughter->SetPreAssignedDecayProducts(pPA); 70 } 75 } 71 theProductVector->push_back(pDaughter); << 76 >> 77 PushProducts( pDaughter ); 72 } 78 } 73 numberOfProducts = right.numberOfProducts; << 74 } 79 } 75 80 76 G4DecayProducts& G4DecayProducts::operator=(co << 81 G4DecayProducts & G4DecayProducts::operator=(const G4DecayProducts &right) 77 { 82 { 78 G4int index; 83 G4int index; 79 84 80 if (this != &right) { << 85 if (this != &right) >> 86 { 81 // recreate parent 87 // recreate parent 82 delete theParentParticle; << 88 if (theParentParticle != 0) delete theParentParticle; 83 theParentParticle = new G4DynamicParticle( 89 theParentParticle = new G4DynamicParticle(*right.theParentParticle); 84 90 85 // delete G4DynamicParticle objects 91 // delete G4DynamicParticle objects 86 for (index = 0; index < numberOfProducts; << 92 for (index=0; index < numberOfProducts; index++) 87 delete theProductVector->at(index); << 93 { >> 94 delete theProductVector[index]; 88 } 95 } 89 theProductVector->clear(); << 90 96 91 // copy daughters (Deep Copy) << 97 //copy daughters (Deep Copy) 92 for (index = 0; index < right.numberOfProd << 98 for (index=0; index < right.numberOfProducts; index++) { 93 G4DynamicParticle* daughter = right.theP << 99 PushProducts( new G4DynamicParticle(*right.theProductVector[index]) ); 94 auto pDaughter = new G4DynamicParticle(* << 100 } 95 << 96 G4double properTime = daughter->GetPreAs << 97 if (properTime > 0.0) pDaughter->SetPreA << 98 << 99 const G4DecayProducts* pPreAssigned = da << 100 if (pPreAssigned != nullptr) { << 101 auto pPA = new G4DecayProducts(*pPreAs << 102 pDaughter->SetPreAssignedDecayProducts << 103 } << 104 theProductVector->push_back(pDaughter); << 105 } << 106 numberOfProducts = right.numberOfProducts; 101 numberOfProducts = right.numberOfProducts; >> 102 107 } 103 } 108 return *this; 104 return *this; 109 } 105 } 110 106 111 G4DecayProducts::~G4DecayProducts() 107 G4DecayProducts::~G4DecayProducts() 112 { 108 { 113 // delete parent << 109 //delete parent 114 delete theParentParticle; << 110 if (theParentParticle != 0) delete theParentParticle; 115 theParentParticle = nullptr; << 111 116 << 117 // delete G4DynamicParticle object 112 // delete G4DynamicParticle object 118 for (G4int index = 0; index < numberOfProduc << 113 for (G4int index=0; index < numberOfProducts; index++) 119 delete theProductVector->at(index); << 114 { >> 115 delete theProductVector[index]; 120 } 116 } 121 theProductVector->clear(); << 117 numberOfProducts = 0; 122 numberOfProducts = 0; << 123 delete theProductVector; << 124 theProductVector = nullptr; << 125 } 118 } 126 119 127 G4DynamicParticle* G4DecayProducts::PopProduct 120 G4DynamicParticle* G4DecayProducts::PopProducts() 128 { 121 { 129 if (numberOfProducts > 0) { << 122 if ( numberOfProducts >0 ) { 130 numberOfProducts -= 1; << 123 numberOfProducts -= 1; 131 G4DynamicParticle* part = theProductVector << 124 return theProductVector[numberOfProducts]; 132 theProductVector->pop_back(); << 125 } else { 133 return part; << 126 return 0; 134 } << 127 } 135 << 136 return nullptr; << 137 } 128 } 138 129 139 G4int G4DecayProducts::PushProducts(G4DynamicP << 130 G4int G4DecayProducts::PushProducts(G4DynamicParticle *aParticle) 140 { << 131 { 141 theProductVector->push_back(aParticle); << 132 if ( numberOfProducts < MaxNumberOfProducts) { 142 numberOfProducts += 1; << 133 theProductVector[numberOfProducts]= aParticle; 143 return numberOfProducts; << 134 numberOfProducts += 1; >> 135 } else { >> 136 #ifdef G4VERBOSE >> 137 G4cout << "G4DecayProducts::PushProducts "; >> 138 G4cout << " exceeds MaxNumberOfProducts(=" >> 139 << G4int(MaxNumberOfProducts) << ")"; >> 140 G4cout << G4endl; >> 141 #endif >> 142 } >> 143 return numberOfProducts; 144 } 144 } 145 145 146 G4DynamicParticle* G4DecayProducts::operator[] 146 G4DynamicParticle* G4DecayProducts::operator[](G4int anIndex) const 147 { 147 { 148 if ((numberOfProducts > anIndex) && (anIndex << 148 if ((numberOfProducts > anIndex) && (anIndex >=0) ) { 149 return theProductVector->at(anIndex); << 149 return theProductVector[anIndex]; 150 } << 150 } else { 151 << 151 return 0; 152 return nullptr; << 152 } 153 } 153 } 154 154 155 void G4DecayProducts::SetParentParticle(const << 155 void G4DecayProducts::SetParentParticle(const G4DynamicParticle &aParticle) 156 { 156 { 157 delete theParentParticle; << 157 if (theParentParticle != 0) delete theParentParticle; 158 theParentParticle = new G4DynamicParticle(aP 158 theParentParticle = new G4DynamicParticle(aParticle); 159 } 159 } 160 160 161 void G4DecayProducts::Boost(G4double totalEner << 161 >> 162 void G4DecayProducts::Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection) 162 { 163 { 163 // calculate new beta << 164 // calcurate new beta 164 G4double mass = theParentParticle->GetMass() << 165 G4double mass = theParentParticle->GetMass(); 165 G4double totalMomentum(0); << 166 G4double totalMomentum = sqrt( (totalEnergy - mass)*(totalEnergy + mass) ); 166 if (totalEnergy > mass) { << 167 G4double betax = momentumDirection.x()*totalMomentum/totalEnergy; 167 totalMomentum = std::sqrt((totalEnergy - m << 168 G4double betay = momentumDirection.y()*totalMomentum/totalEnergy; 168 } << 169 G4double betaz = momentumDirection.z()*totalMomentum/totalEnergy; 169 G4double betax = momentumDirection.x() * tot << 170 this->Boost(betax, betay, betaz); 170 G4double betay = momentumDirection.y() * tot << 171 G4double betaz = momentumDirection.z() * tot << 172 Boost(betax, betay, betaz); << 173 } 171 } 174 172 175 void G4DecayProducts::Boost(G4double newbetax, 173 void G4DecayProducts::Boost(G4double newbetax, G4double newbetay, G4double newbetaz) 176 { << 174 { 177 G4double mass = theParentParticle->GetMass() << 175 G4double mass = theParentParticle->GetMass(); 178 G4double energy = theParentParticle->GetTota << 176 G4double energy = theParentParticle->GetTotalEnergy(); 179 G4double momentum = 0.0; << 177 G4double momentum = 0.0; 180 178 181 G4ThreeVector direction(0.0, 0.0, 1.0); << 179 G4ThreeVector direction(0.0,0.0,1.0); 182 G4LorentzVector p4; 180 G4LorentzVector p4; 183 181 184 if (energy - mass > DBL_MIN) { 182 if (energy - mass > DBL_MIN) { 185 // calcurate beta of initial state 183 // calcurate beta of initial state 186 momentum = theParentParticle->GetTotalMome << 184 momentum = theParentParticle->GetTotalMomentum(); 187 direction = theParentParticle->GetMomentum 185 direction = theParentParticle->GetMomentumDirection(); 188 G4double betax = -1.0 * direction.x() * mo << 186 G4double betax = -1.0*direction.x()*momentum/energy; 189 G4double betay = -1.0 * direction.y() * mo << 187 G4double betay = -1.0*direction.y()*momentum/energy; 190 G4double betaz = -1.0 * direction.z() * mo << 188 G4double betaz = -1.0*direction.z()*momentum/energy; 191 << 189 192 for (G4int index = 0; index < numberOfProd << 190 for (G4int index=0; index < numberOfProducts; index++) { 193 // make G4LorentzVector for secondaries << 191 // make G4LorentzVector for secondaries 194 p4 = (theProductVector->at(index))->Get4 << 192 p4 = theProductVector[index]->Get4Momentum(); 195 << 193 196 // boost secondaries to theParentParticl << 194 // boost secondaries to theParentParticle's rest frame 197 p4.boost(betax, betay, betaz); << 195 p4.boost(betax, betay, betaz); >> 196 >> 197 // boost secondaries to new frame >> 198 p4.boost(newbetax, newbetay, newbetaz); >> 199 >> 200 // change energy/momentum >> 201 theProductVector[index]->Set4Momentum(p4); >> 202 } >> 203 } else { >> 204 for (G4int index=0; index < numberOfProducts; index++) { >> 205 // make G4LorentzVector for secondaries >> 206 p4 = theProductVector[index]->Get4Momentum(); 198 207 199 // boost secondaries to new frame << 208 // boost secondaries to new frame 200 p4.boost(newbetax, newbetay, newbetaz); << 209 p4.boost(newbetax, newbetay, newbetaz); 201 << 202 // change energy/momentum << 203 (theProductVector->at(index))->Set4Momen << 204 } << 205 } << 206 else { << 207 for (G4int index = 0; index < numberOfProd << 208 // make G4LorentzVector for secondaries << 209 p4 = (theProductVector->at(index))->Get4 << 210 210 211 // boost secondaries to new frame << 211 // change energy/momentum 212 p4.boost(newbetax, newbetay, newbetaz); << 212 theProductVector[index]->Set4Momentum(p4); 213 << 213 } 214 // change energy/momentum << 214 } 215 (theProductVector->at(index))->Set4Momen << 215 // make G4LorentzVector for parent in its rest frame 216 } << 216 mass = theParentParticle->GetMass(); 217 } << 217 G4LorentzVector parent4( 0.0, 0.0, 0.0, mass); 218 << 219 // make G4LorentzVector for parent in its re << 220 mass = theParentParticle->GetMass(); << 221 G4LorentzVector parent4(0.0, 0.0, 0.0, mass) << 222 218 223 // boost parent to new frame << 219 // boost parent to new frame 224 parent4.boost(newbetax, newbetay, newbetaz); << 220 parent4.boost(newbetax, newbetay, newbetaz); 225 221 226 // change energy/momentum << 222 // change energy/momentum 227 theParentParticle->Set4Momentum(parent4); << 223 theParentParticle->Set4Momentum(parent4); 228 } 224 } 229 225 230 G4bool G4DecayProducts::IsChecked() const 226 G4bool G4DecayProducts::IsChecked() const 231 { 227 { 232 G4bool returnValue = true; 228 G4bool returnValue = true; 233 << 229 // check parent 234 // check parent << 230 // energy/momentum 235 // energy/momentum << 231 G4double parent_energy = theParentParticle->GetTotalEnergy(); 236 G4double parent_energy = theParentParticle-> << 237 G4ThreeVector direction = theParentParticle- 232 G4ThreeVector direction = theParentParticle->GetMomentumDirection(); 238 G4ThreeVector parent_momentum = direction * << 233 G4ThreeVector parent_momentum = direction*(theParentParticle->GetTotalMomentum()); 239 << 234 // check momentum dirction is a unit vector 240 // check momentum direction is a unit vector << 235 if ( (parent_momentum.mag() >0.0) && (abs(direction.mag()-1.0) >1.0e-6 ) ) { 241 if ((parent_momentum.mag() > 0.0) && (std::f << 242 #ifdef G4VERBOSE 236 #ifdef G4VERBOSE 243 G4cout << "G4DecayProducts::IsChecked():: << 237 G4cout << " Momentum Direction Vector of Parent is not normalized "; 244 << " Momentum Direction Vector of P << 238 G4cout << " (=" << direction.mag() << ")" << G4endl; 245 << " (=" << direction.mag() << ")" << 246 #endif 239 #endif 247 returnValue = false; 240 returnValue = false; 248 parent_momentum = parent_momentum * (1. / << 241 parent_momentum = parent_momentum * (1./direction.mag()); 249 } 242 } 250 243 251 // daughters << 244 //daughters 252 G4double mass, energy; << 245 G4double mass, energy; 253 G4ThreeVector momentum; 246 G4ThreeVector momentum; 254 G4double total_energy = parent_energy; << 247 G4double total_energy = parent_energy; 255 G4ThreeVector total_momentum = parent_moment << 248 G4ThreeVector total_momentum = parent_momentum; 256 << 249 for (G4int index=0; index < numberOfProducts; index++) 257 for (G4int index = 0; index < numberOfProduc << 250 { 258 G4DynamicParticle* part = theProductVector << 251 mass = theProductVector[index]->GetMass(); 259 mass = part->GetMass(); << 252 energy = theProductVector[index]->GetTotalEnergy(); 260 energy = part->GetTotalEnergy(); << 253 direction = theProductVector[index]->GetMomentumDirection(); 261 direction = part->GetMomentumDirection(); << 254 momentum = direction*(theProductVector[index]->GetTotalMomentum()); 262 momentum = direction * (part->GetTotalMome << 255 // check momentum dirction is a unit vector 263 << 256 if ( (momentum.mag()>0.0) && (abs(direction.mag()-1.0) > 1.0e-6)) { 264 // check momentum direction is a unit vect << 265 if ((momentum.mag() > 0.0) && (std::fabs(d << 266 #ifdef G4VERBOSE 257 #ifdef G4VERBOSE 267 G4cout << "G4DecayProducts::IsChecked(): << 258 G4cout << " Momentum Direction Vector of Daughter [" << index; 268 << " Momentum Direction Vector of << 259 G4cout << "] is not normalized (=" << direction.mag() << ")" << G4endl; 269 << "] is not normalized (=" << d << 270 #endif 260 #endif 271 returnValue = false; 261 returnValue = false; 272 momentum = momentum * (1. / direction.ma << 262 momentum = momentum * (1./direction.mag()); 273 } 263 } 274 // whether daughter stops or not 264 // whether daughter stops or not 275 if (energy - mass < DBL_MIN) { << 265 if (energy - mass < DBL_MIN ) { 276 #ifdef G4VERBOSE 266 #ifdef G4VERBOSE 277 G4cout << "G4DecayProducts::IsChecked(): << 267 G4cout << "Daughter [" << index << "] has no kinetic energy "<< G4endl; 278 << " Daughter [" << index << "] << 279 #endif 268 #endif 280 returnValue = false; 269 returnValue = false; 281 } 270 } 282 total_energy -= energy; << 271 total_energy -= energy; 283 total_momentum -= momentum; 272 total_momentum -= momentum; 284 } 273 } 285 // check energy/momentum conservation 274 // check energy/momentum conservation 286 if ((std::fabs(total_energy) > 1.0e-9 * MeV) << 275 if ( (abs(total_energy) >1.0e-6) || (total_momentum.mag() >1.0e-6 ) ){ 287 #ifdef G4VERBOSE 276 #ifdef G4VERBOSE 288 G4cout << "G4DecayProducts::IsChecked():: << 277 G4cout << " Energy/Momentum is not conserved "; 289 << " Energy/Momentum is not conserv << 278 G4cout << total_energy << " " << total_momentum << G4endl; 290 G4cout << " difference between parent ener << 291 << "[MeV] " << G4endl; << 292 G4cout << " difference between parent mome << 293 << " x:" << total_momentum.getX() / << 294 << " z:" << total_momentum.getZ() / << 295 #endif 279 #endif 296 returnValue = false; 280 returnValue = false; 297 } 281 } 298 return returnValue; 282 return returnValue; 299 } 283 } 300 284 301 void G4DecayProducts::DumpInfo() const 285 void G4DecayProducts::DumpInfo() const 302 { 286 { 303 G4cout << " ----- List of DecayProducts --- << 287 G4cout << " ----- List of DecayProducts -----" << G4endl; 304 G4cout << " ------ Parent Particle --------- << 288 G4cout << " ------ Parent Particle ----------" << G4endl; 305 if (theParentParticle != nullptr) theParentP << 289 if (theParentParticle != 0) theParentParticle->DumpInfo(); 306 G4cout << " ------ Daughter Particles ----- << 290 G4cout << " ------ Daughter Particles ------" << G4endl; 307 for (G4int index = 0; index < numberOfProduc << 291 for (G4int index=0; index < numberOfProducts; index++) 308 G4cout << " ----------" << index + 1 << " << 292 { 309 (theProductVector->at(index))->DumpInfo(); << 293 G4cout << " ----------" << index+1 << " -------------" << G4endl; 310 } << 294 theProductVector[index]-> DumpInfo(); 311 G4cout << " ----- End List of DecayProducts << 295 } 312 G4cout << G4endl; << 296 G4cout << " ----- End List of DecayProducts -----" << G4endl; 313 } << 297 G4cout << G4endl; >> 298 } 314 299