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