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 // G4PrimaryTransformer class implementation << 27 // 23 // 28 // Author: Makoto Asai, 1999 << 24 // $Id: G4PrimaryTransformer.cc,v 1.13 2001/11/22 18:53:13 asaim Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-00 $ >> 26 // 30 27 31 #include "G4PrimaryTransformer.hh" 28 #include "G4PrimaryTransformer.hh" 32 #include "G4SystemOfUnits.hh" << 33 #include "G4Event.hh" 29 #include "G4Event.hh" 34 #include "G4PrimaryVertex.hh" 30 #include "G4PrimaryVertex.hh" 35 #include "G4ParticleDefinition.hh" 31 #include "G4ParticleDefinition.hh" 36 #include "G4DynamicParticle.hh" 32 #include "G4DynamicParticle.hh" 37 #include "G4Track.hh" 33 #include "G4Track.hh" 38 #include "G4ThreeVector.hh" 34 #include "G4ThreeVector.hh" 39 #include "G4DecayProducts.hh" 35 #include "G4DecayProducts.hh" >> 36 #include "G4TransportationManager.hh" >> 37 #include "G4Navigator.hh" >> 38 #include "G4VPhysicalVolume.hh" >> 39 #include "G4LogicalVolume.hh" >> 40 #include "G4VSolid.hh" 40 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh" 41 #include "G4ios.hh" 42 #include "G4ios.hh" 42 #include "Randomize.hh" << 43 << 44 G4double G4PrimaryTransformer::kETolerance = 1 << 45 G4ExceptionSeverity G4PrimaryTransformer::kETS << 46 43 47 G4PrimaryTransformer::G4PrimaryTransformer() 44 G4PrimaryTransformer::G4PrimaryTransformer() >> 45 :verboseLevel(0),trackID(0) 48 { 46 { 49 particleTable = G4ParticleTable::GetParticle 47 particleTable = G4ParticleTable::GetParticleTable(); 50 CheckUnknown(); << 51 } << 52 << 53 void G4PrimaryTransformer::CheckUnknown() << 54 { << 55 unknown = particleTable->FindParticle("unkno << 56 unknownParticleDefined = unknown != nullptr; << 57 chargedunknown = particleTable->FindParticle << 58 chargedUnknownParticleDefined = chargedunkno << 59 opticalphoton = particleTable->FindParticle( << 60 opticalphotonDefined = opticalphoton != null << 61 } 48 } 62 49 63 G4TrackVector* << 50 G4PrimaryTransformer::~G4PrimaryTransformer() 64 G4PrimaryTransformer::GimmePrimaries(G4Event* << 51 {;} >> 52 >> 53 G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent,G4int trackIDCounter) 65 { 54 { 66 trackID = trackIDCounter; 55 trackID = trackIDCounter; 67 56 68 for(auto tr : TV) delete tr; << 57 //TV.clearAndDestroy(); >> 58 for( size_t ii=0; ii<TV.size();ii++) >> 59 { delete TV[ii]; } 69 TV.clear(); 60 TV.clear(); 70 << 61 G4int n_vertex = anEvent->GetNumberOfPrimaryVertex(); 71 // Loop over vertices << 62 if(n_vertex==0) return 0; 72 // << 63 for( G4int i=0; i<n_vertex; i++ ) 73 G4PrimaryVertex* nextVertex = anEvent->GetPr << 64 { GenerateTracks( anEvent->GetPrimaryVertex(i) ); } 74 while(nextVertex != nullptr) // Loop checkin << 75 { << 76 GenerateTracks(nextVertex); << 77 nextVertex = nextVertex->GetNext(); << 78 } << 79 return &TV; 65 return &TV; 80 } 66 } 81 67 82 void G4PrimaryTransformer::GenerateTracks(G4Pr 68 void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex) 83 { 69 { 84 G4double X0 = primaryVertex->GetX0(); 70 G4double X0 = primaryVertex->GetX0(); 85 G4double Y0 = primaryVertex->GetY0(); 71 G4double Y0 = primaryVertex->GetY0(); 86 G4double Z0 = primaryVertex->GetZ0(); 72 G4double Z0 = primaryVertex->GetZ0(); 87 G4double T0 = primaryVertex->GetT0(); 73 G4double T0 = primaryVertex->GetT0(); 88 G4double WV = primaryVertex->GetWeight(); 74 G4double WV = primaryVertex->GetWeight(); 89 75 90 #ifdef G4VERBOSE 76 #ifdef G4VERBOSE 91 if(verboseLevel>2) << 77 if(verboseLevel>1) 92 { 78 { 93 primaryVertex->Print(); << 94 } << 95 else if (verboseLevel==1) << 96 { << 97 G4cout << "G4PrimaryTransformer::PrimaryVe 79 G4cout << "G4PrimaryTransformer::PrimaryVertex (" 98 << X0 / mm << "(mm)," 80 << X0 / mm << "(mm)," 99 << Y0 / mm << "(mm)," 81 << Y0 / mm << "(mm)," 100 << Z0 / mm << "(mm)," 82 << Z0 / mm << "(mm)," 101 << T0 / nanosecond << "(nsec))" << 83 << T0 / nanosecond << "(nsec))" << G4endl; 102 } 84 } 103 #endif 85 #endif 104 86 >> 87 //Check whether the vertex is inside the world volume >> 88 G4Navigator* navi = >> 89 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); >> 90 G4VSolid* theWorldSolid = navi->GetWorldVolume()->GetLogicalVolume()->GetSolid(); >> 91 if(theWorldSolid->Inside(G4ThreeVector(X0,Y0,Z0))==kOutside) >> 92 { >> 93 G4cerr << "Primary vertex " >> 94 << G4BestUnit(G4ThreeVector(X0,Y0,Z0),G4String("Length")) >> 95 << " is outside of the world volume." << G4endl; >> 96 G4Exception("G4PrimaryTransformer::Primary vertex outside of the world"); >> 97 } >> 98 105 G4PrimaryParticle* primaryParticle = primary 99 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary(); 106 while( primaryParticle != nullptr ) // Loop << 100 while( primaryParticle != 0 ) 107 { 101 { 108 GenerateSingleTrack( primaryParticle, X0, 102 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV ); 109 primaryParticle = primaryParticle->GetNext 103 primaryParticle = primaryParticle->GetNext(); 110 } 104 } 111 } 105 } 112 106 113 void G4PrimaryTransformer:: << 107 void G4PrimaryTransformer::GenerateSingleTrack 114 GenerateSingleTrack( G4PrimaryParticle* primar << 108 (G4PrimaryParticle* primaryParticle, 115 G4double x0, G4double y0, << 109 G4double x0,G4double y0,G4double z0,G4double t0,G4double wv) 116 G4double t0, G4double wv) << 117 { 110 { 118 G4ParticleDefinition* partDef = GetDefinitio 111 G4ParticleDefinition* partDef = GetDefinition(primaryParticle); 119 if(!IsGoodForTrack(partDef)) << 112 if((!partDef)||partDef->IsShortLived()) 120 { // The particle cannot be converted to G4 << 113 // The particle is not defined in GEANT4, check daughters >> 114 { 121 #ifdef G4VERBOSE 115 #ifdef G4VERBOSE 122 if(verboseLevel>1) << 116 if(verboseLevel>2) 123 { 117 { 124 G4cout << "Primary particle (PDGcode " < 118 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode() 125 << ") --- Ignored" << G4endl; 119 << ") --- Ignored" << G4endl; 126 } 120 } 127 #endif 121 #endif 128 G4PrimaryParticle* daughter = primaryParti 122 G4PrimaryParticle* daughter = primaryParticle->GetDaughter(); 129 while(daughter != nullptr) // Loop checkin << 123 while(daughter) 130 { 124 { 131 GenerateSingleTrack(daughter,x0,y0,z0,t0 125 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv); 132 daughter = daughter->GetNext(); 126 daughter = daughter->GetNext(); 133 } 127 } 134 } 128 } 135 else // The particle is defined in GEANT4 << 129 >> 130 // The particle is defined in GEANT4 >> 131 else 136 { 132 { 137 // Create G4DynamicParticle object 133 // Create G4DynamicParticle object 138 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE 139 if(verboseLevel>1) 135 if(verboseLevel>1) 140 { 136 { 141 G4cout << "Primary particle (" << partDe 137 G4cout << "Primary particle (" << partDef->GetParticleName() 142 << ") --- Transferred with moment << 138 << ") --- Transfered with momentum " << primaryParticle->GetMomentum() 143 << primaryParticle->GetMomentum() << 144 << G4endl; 139 << G4endl; 145 } 140 } 146 #endif 141 #endif 147 // Check the mass of the "real" particle << 142 G4DynamicParticle* DP = 148 // N.B. PDG code 0 is used for artificial << 143 new G4DynamicParticle(partDef,primaryParticle->GetMomentum()); 149 if(primaryParticle->GetPDGcode() != 0 && k << 144 DP->SetPolarization(primaryParticle->GetPolX(), 150 { << 145 primaryParticle->GetPolY(), 151 G4double pmas = partDef->GetPDGMass(); << 146 primaryParticle->GetPolZ()); 152 if(std::abs(pmas - primaryParticle->GetM << 147 if(primaryParticle->GetProperTime()>0.0) 153 G4ExceptionDescription ed; << 148 { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); } 154 ed << "Primary particle PDG=" << prima << 149 // Set Charge 155 << " deltaMass(MeV)=" << (std::abs( << 150 if (abs(primaryParticle->GetCharge()-DP->GetCharge())>eplus) { 156 << " is larger than the tolerance(M << 151 DP->SetCharge(primaryParticle->GetCharge()); 157 << "\n Specified mass(MeV)=" << (pr << 158 << " while PDG mass(MEV)=" << pmas/ << 159 << "\n To change the tolerance or t << 160 << " use G4PrimaryTransformer::SetK << 161 G4Exception("G4PrimaryParticle::Set4Mo << 162 } << 163 } << 164 auto* DP = << 165 new G4DynamicParticle(partDef, << 166 primaryParticle->G << 167 primaryParticle->G << 168 if(opticalphotonDefined && partDef==optica << 169 && primaryParticle->GetPolarization().m << 170 { << 171 if(nWarn<10) << 172 { << 173 G4Exception("G4PrimaryTransformer::Gen << 174 "ZeroPolarization", JustWa << 175 "Polarization of the optic << 176 Random polarization is as << 177 G4cerr << "This warning message is iss << 178 ++nWarn; << 179 } << 180 << 181 G4double angle = G4UniformRand() * 360.0 << 182 G4ThreeVector normal (1., 0., 0.); << 183 G4ThreeVector kphoton = DP->GetMomentumD << 184 G4ThreeVector product = normal.cross(kph << 185 G4double modul2 = product*product; << 186 << 187 G4ThreeVector e_perpend (0., 0., 1.); << 188 if (modul2 > 0.) e_perpend = (1./std::sq << 189 G4ThreeVector e_paralle = e_perpend.c << 190 << 191 G4ThreeVector polar = std::cos(angle)*e_ << 192 + std::sin(angle)*e_ << 193 DP->SetPolarization(polar.x(),polar.y(), << 194 } << 195 else << 196 { << 197 DP->SetPolarization(primaryParticle->Get << 198 primaryParticle->Get << 199 primaryParticle->Get << 200 } << 201 if(primaryParticle->GetProperTime()>=0.0) << 202 { << 203 DP->SetPreAssignedDecayProperTime(primar << 204 } << 205 << 206 // Set Mass if it is specified << 207 // << 208 G4double pmas = primaryParticle->GetMass() << 209 if(pmas>=0.) { DP->SetMass(pmas); } << 210 << 211 // Set Charge if it is specified << 212 // << 213 if (primaryParticle->GetCharge()<DBL_MAX) << 214 { << 215 if (partDef->GetAtomicNumber() <0) << 216 { << 217 DP->SetCharge(primaryParticle->GetChar << 218 } << 219 else // ions << 220 { << 221 G4int iz = partDef->GetAtomicNumber(); << 222 auto iq = static_cast<G4int>(primaryP << 223 G4int n_e = iz - iq; << 224 if (n_e>0) DP->AddElectron(0,n_e); << 225 } << 226 } 152 } 227 // Set decay products to the DynamicPartic 153 // Set decay products to the DynamicParticle 228 // << 229 SetDecayProducts( primaryParticle, DP ); 154 SetDecayProducts( primaryParticle, DP ); 230 << 231 // Set primary particle << 232 // << 233 DP->SetPrimaryParticle(primaryParticle); << 234 << 235 // Set PDG code if it is different from G4 << 236 // << 237 if(partDef->GetPDGEncoding()==0 && primary << 238 { << 239 DP->SetPDGcode(primaryParticle->GetPDGco << 240 } << 241 << 242 // Check the particle is properly construc << 243 // << 244 if(!CheckDynamicParticle(DP)) << 245 { << 246 delete DP; << 247 return; << 248 } << 249 << 250 // Create G4Track object 155 // Create G4Track object 251 // << 156 G4Track* track = new G4Track(DP,t0,G4ThreeVector(x0,y0,z0)); 252 auto track = new G4Track(DP,t0,G4ThreeVec << 253 << 254 // Set trackID and let primary particle kn 157 // Set trackID and let primary particle know it 255 // << 158 trackID++; 256 ++trackID; << 257 track->SetTrackID(trackID); 159 track->SetTrackID(trackID); 258 primaryParticle->SetTrackID(trackID); 160 primaryParticle->SetTrackID(trackID); 259 << 260 // Set parentID to 0 as a primary particle 161 // Set parentID to 0 as a primary particle 261 // << 262 track->SetParentID(0); 162 track->SetParentID(0); 263 << 264 // Set weight ( vertex weight * particle w 163 // Set weight ( vertex weight * particle weight ) 265 // << 266 track->SetWeight(wv*(primaryParticle->GetW 164 track->SetWeight(wv*(primaryParticle->GetWeight())); 267 << 268 // Store it to G4TrackVector 165 // Store it to G4TrackVector 269 // << 270 TV.push_back( track ); 166 TV.push_back( track ); 271 } 167 } 272 } 168 } 273 169 274 void G4PrimaryTransformer:: << 170 void G4PrimaryTransformer::SetDecayProducts 275 SetDecayProducts(G4PrimaryParticle* mother, G4 << 171 (G4PrimaryParticle* mother, G4DynamicParticle* motherDP) 276 { 172 { 277 G4PrimaryParticle* daughter = mother->GetDau 173 G4PrimaryParticle* daughter = mother->GetDaughter(); 278 if(daughter == nullptr) return; << 174 if(!daughter) return; 279 auto* decayProducts << 175 G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() ); 280 = (G4DecayProducts*)(motherDP->GetPreAssig << 176 if(!decayProducts) 281 if(decayProducts == nullptr) << 282 { 177 { 283 decayProducts = new G4DecayProducts(*mothe 178 decayProducts = new G4DecayProducts(*motherDP); 284 motherDP->SetPreAssignedDecayProducts(deca 179 motherDP->SetPreAssignedDecayProducts(decayProducts); 285 } 180 } 286 while(daughter != nullptr) << 181 while(daughter) 287 { 182 { 288 G4ParticleDefinition* partDef = GetDefinit 183 G4ParticleDefinition* partDef = GetDefinition(daughter); 289 if(!IsGoodForTrack(partDef)) << 184 if(!partDef) 290 { 185 { 291 #ifdef G4VERBOSE 186 #ifdef G4VERBOSE 292 if(verboseLevel>2) 187 if(verboseLevel>2) 293 { 188 { 294 G4cout << " >> Decay product (PDGcode 189 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode() 295 << ") --- Ignored" << G4endl; 190 << ") --- Ignored" << G4endl; 296 } 191 } 297 #endif 192 #endif 298 SetDecayProducts(daughter,motherDP); 193 SetDecayProducts(daughter,motherDP); 299 } 194 } 300 else 195 else 301 { 196 { 302 #ifdef G4VERBOSE 197 #ifdef G4VERBOSE 303 if(verboseLevel>1) 198 if(verboseLevel>1) 304 { 199 { 305 G4cout << " >> Decay product (" << par 200 G4cout << " >> Decay product (" << partDef->GetParticleName() 306 << ") --- Attached with momentu 201 << ") --- Attached with momentum " << daughter->GetMomentum() 307 << G4endl; 202 << G4endl; 308 } 203 } 309 #endif 204 #endif 310 auto* DP << 205 G4DynamicParticle*DP 311 = new G4DynamicParticle(partDef,daught 206 = new G4DynamicParticle(partDef,daughter->GetMomentum()); 312 DP->SetPrimaryParticle(daughter); << 313 << 314 // Decay proper time for daughter << 315 // << 316 if(daughter->GetProperTime()>=0.0) << 317 { << 318 DP->SetPreAssignedDecayProperTime(daug << 319 } << 320 << 321 // Set Charge and Mass is specified << 322 // << 323 if (daughter->GetCharge()<DBL_MAX) << 324 { << 325 DP->SetCharge(daughter->GetCharge()); << 326 } << 327 G4double pmas = daughter->GetMass(); << 328 if(pmas>=0.) << 329 { << 330 DP->SetMass(pmas); << 331 } << 332 << 333 // Set Polarization << 334 // << 335 DP->SetPolarization(daughter->GetPolX(), << 336 daughter->GetPolY(), << 337 daughter->GetPolZ()) << 338 decayProducts->PushProducts(DP); 207 decayProducts->PushProducts(DP); 339 SetDecayProducts(daughter,DP); 208 SetDecayProducts(daughter,DP); 340 << 341 // Check the particle is properly constr << 342 // << 343 if(!CheckDynamicParticle(DP)) << 344 { << 345 delete DP; << 346 return; << 347 } << 348 } 209 } 349 daughter = daughter->GetNext(); 210 daughter = daughter->GetNext(); 350 } 211 } 351 } 212 } 352 213 353 void G4PrimaryTransformer::SetUnknownParticleD << 354 { << 355 unknownParticleDefined = vl; << 356 if(unknownParticleDefined && (unknown == nul << 357 { << 358 G4cerr << "unknownParticleDefined cannot b << 359 << "G4UnknownParticle is not define << 360 << "Command ignored." << G4endl; << 361 unknownParticleDefined = false; << 362 } << 363 } << 364 << 365 void G4PrimaryTransformer::SetChargedUnknownPa << 366 { << 367 chargedUnknownParticleDefined = vl; << 368 if(chargedUnknownParticleDefined && (charged << 369 { << 370 G4cerr << "chargedUnknownParticleDefined c << 371 << "G4ChargedUnknownParticle is not << 372 << "Command ignored." << G4endl; << 373 chargedUnknownParticleDefined = false; << 374 } << 375 } << 376 << 377 G4bool G4PrimaryTransformer::CheckDynamicParti << 378 { << 379 if(IsGoodForTrack(DP->GetDefinition())) retu << 380 auto* decayProducts << 381 = (G4DecayProducts*)(DP->GetPreAssignedDec << 382 if(decayProducts != nullptr && decayProducts << 383 G4cerr << G4endl << 384 << "G4PrimaryTransformer: a shortlive << 385 << G4endl << 386 << " without any valid decay table no << 387 << G4endl; << 388 G4Exception("G4PrimaryTransformer", "Invalid << 389 "This primary particle will be i << 390 return false; << 391 } << 392 214 393 G4ParticleDefinition* << 394 G4PrimaryTransformer::GetDefinition(G4PrimaryP << 395 { << 396 G4ParticleDefinition* partDef = pp->GetG4cod << 397 if(partDef == nullptr) << 398 { << 399 partDef = particleTable->FindParticle(pp-> << 400 } << 401 if((partDef == nullptr) || partDef->IsShortL << 402 { << 403 if (chargedUnknownParticleDefined && pp->G << 404 { << 405 partDef = chargedunknown; << 406 } << 407 else if (unknownParticleDefined) << 408 { << 409 partDef = unknown; << 410 } << 411 } << 412 return partDef; << 413 } << 414 215 415 G4bool G4PrimaryTransformer::IsGoodForTrack(G4 << 416 { << 417 if(pd == nullptr) << 418 { return false; } << 419 if(!(pd->IsShortLived())) << 420 { return true; } << 421 // << 422 // Following two lines should be removed if << 423 // shortlived primary particle with proper d << 424 // a track. << 425 // << 426 if(pd->GetDecayTable() != nullptr) << 427 { return true; } << 428 216 429 return false; << 430 } << 431 217