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.18 2003/10/01 13:17:34 asaim Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-06-00-patch-01 $ >> 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" 40 #include "G4UnitsTable.hh" 36 #include "G4UnitsTable.hh" 41 #include "G4ios.hh" 37 #include "G4ios.hh" 42 #include "Randomize.hh" << 43 << 44 G4double G4PrimaryTransformer::kETolerance = 1 << 45 G4ExceptionSeverity G4PrimaryTransformer::kETS << 46 38 47 G4PrimaryTransformer::G4PrimaryTransformer() 39 G4PrimaryTransformer::G4PrimaryTransformer() >> 40 :verboseLevel(0),trackID(0) 48 { 41 { 49 particleTable = G4ParticleTable::GetParticle 42 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 } 43 } 62 44 63 G4TrackVector* << 45 G4PrimaryTransformer::~G4PrimaryTransformer() 64 G4PrimaryTransformer::GimmePrimaries(G4Event* << 46 {;} >> 47 >> 48 G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent,G4int trackIDCounter) 65 { 49 { 66 trackID = trackIDCounter; 50 trackID = trackIDCounter; 67 51 68 for(auto tr : TV) delete tr; << 52 //TV.clearAndDestroy(); >> 53 for( size_t ii=0; ii<TV.size();ii++) >> 54 { delete TV[ii]; } 69 TV.clear(); 55 TV.clear(); 70 << 56 G4int n_vertex = anEvent->GetNumberOfPrimaryVertex(); 71 // Loop over vertices << 57 if(n_vertex==0) return 0; 72 // << 58 for( G4int i=0; i<n_vertex; i++ ) 73 G4PrimaryVertex* nextVertex = anEvent->GetPr << 59 { GenerateTracks( anEvent->GetPrimaryVertex(i) ); } 74 while(nextVertex != nullptr) // Loop checkin << 75 { << 76 GenerateTracks(nextVertex); << 77 nextVertex = nextVertex->GetNext(); << 78 } << 79 return &TV; 60 return &TV; 80 } 61 } 81 62 82 void G4PrimaryTransformer::GenerateTracks(G4Pr 63 void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex) 83 { 64 { 84 G4double X0 = primaryVertex->GetX0(); 65 G4double X0 = primaryVertex->GetX0(); 85 G4double Y0 = primaryVertex->GetY0(); 66 G4double Y0 = primaryVertex->GetY0(); 86 G4double Z0 = primaryVertex->GetZ0(); 67 G4double Z0 = primaryVertex->GetZ0(); 87 G4double T0 = primaryVertex->GetT0(); 68 G4double T0 = primaryVertex->GetT0(); 88 G4double WV = primaryVertex->GetWeight(); 69 G4double WV = primaryVertex->GetWeight(); 89 70 90 #ifdef G4VERBOSE 71 #ifdef G4VERBOSE 91 if(verboseLevel>2) << 72 if(verboseLevel>1) 92 { 73 { 93 primaryVertex->Print(); << 94 } << 95 else if (verboseLevel==1) << 96 { << 97 G4cout << "G4PrimaryTransformer::PrimaryVe 74 G4cout << "G4PrimaryTransformer::PrimaryVertex (" 98 << X0 / mm << "(mm)," 75 << X0 / mm << "(mm)," 99 << Y0 / mm << "(mm)," 76 << Y0 / mm << "(mm)," 100 << Z0 / mm << "(mm)," 77 << Z0 / mm << "(mm)," 101 << T0 / nanosecond << "(nsec))" << 78 << T0 / nanosecond << "(nsec))" << G4endl; 102 } 79 } 103 #endif 80 #endif 104 81 105 G4PrimaryParticle* primaryParticle = primary 82 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary(); 106 while( primaryParticle != nullptr ) // Loop << 83 while( primaryParticle != 0 ) 107 { 84 { 108 GenerateSingleTrack( primaryParticle, X0, 85 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV ); 109 primaryParticle = primaryParticle->GetNext 86 primaryParticle = primaryParticle->GetNext(); 110 } 87 } 111 } 88 } 112 89 113 void G4PrimaryTransformer:: << 90 void G4PrimaryTransformer::GenerateSingleTrack 114 GenerateSingleTrack( G4PrimaryParticle* primar << 91 (G4PrimaryParticle* primaryParticle, 115 G4double x0, G4double y0, << 92 G4double x0,G4double y0,G4double z0,G4double t0,G4double wv) 116 G4double t0, G4double wv) << 117 { 93 { 118 G4ParticleDefinition* partDef = GetDefinitio 94 G4ParticleDefinition* partDef = GetDefinition(primaryParticle); 119 if(!IsGoodForTrack(partDef)) << 95 if((!partDef)||partDef->IsShortLived()) 120 { // The particle cannot be converted to G4 << 96 // The particle is not defined in GEANT4, check daughters >> 97 { 121 #ifdef G4VERBOSE 98 #ifdef G4VERBOSE 122 if(verboseLevel>1) << 99 if(verboseLevel>2) 123 { 100 { 124 G4cout << "Primary particle (PDGcode " < 101 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode() 125 << ") --- Ignored" << G4endl; 102 << ") --- Ignored" << G4endl; 126 } 103 } 127 #endif 104 #endif 128 G4PrimaryParticle* daughter = primaryParti 105 G4PrimaryParticle* daughter = primaryParticle->GetDaughter(); 129 while(daughter != nullptr) // Loop checkin << 106 while(daughter) 130 { 107 { 131 GenerateSingleTrack(daughter,x0,y0,z0,t0 108 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv); 132 daughter = daughter->GetNext(); 109 daughter = daughter->GetNext(); 133 } 110 } 134 } 111 } 135 else // The particle is defined in GEANT4 << 112 >> 113 // The particle is defined in GEANT4 >> 114 else 136 { 115 { 137 // Create G4DynamicParticle object 116 // Create G4DynamicParticle object 138 #ifdef G4VERBOSE 117 #ifdef G4VERBOSE 139 if(verboseLevel>1) 118 if(verboseLevel>1) 140 { 119 { 141 G4cout << "Primary particle (" << partDe 120 G4cout << "Primary particle (" << partDef->GetParticleName() 142 << ") --- Transferred with moment << 121 << ") --- Transfered with momentum " << primaryParticle->GetMomentum() 143 << primaryParticle->GetMomentum() << 144 << G4endl; 122 << G4endl; 145 } 123 } 146 #endif 124 #endif 147 // Check the mass of the "real" particle << 125 G4DynamicParticle* DP = 148 // N.B. PDG code 0 is used for artificial << 126 new G4DynamicParticle(partDef,primaryParticle->GetMomentum()); 149 if(primaryParticle->GetPDGcode() != 0 && k << 127 DP->SetPolarization(primaryParticle->GetPolX(), 150 { << 128 primaryParticle->GetPolY(), 151 G4double pmas = partDef->GetPDGMass(); << 129 primaryParticle->GetPolZ()); 152 if(std::abs(pmas - primaryParticle->GetM << 130 if(primaryParticle->GetProperTime()>0.0) 153 G4ExceptionDescription ed; << 131 { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); } 154 ed << "Primary particle PDG=" << prima << 132 // Set Charge is specified 155 << " deltaMass(MeV)=" << (std::abs( << 133 if (primaryParticle->GetCharge()<DBL_MAX) { 156 << " is larger than the tolerance(M << 134 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 } 135 } 227 // Set decay products to the DynamicPartic 136 // Set decay products to the DynamicParticle 228 // << 229 SetDecayProducts( primaryParticle, DP ); 137 SetDecayProducts( primaryParticle, DP ); 230 << 231 // Set primary particle 138 // Set primary particle 232 // << 233 DP->SetPrimaryParticle(primaryParticle); 139 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 140 // Create G4Track object 251 // << 141 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 142 // Set trackID and let primary particle know it 255 // << 143 trackID++; 256 ++trackID; << 257 track->SetTrackID(trackID); 144 track->SetTrackID(trackID); 258 primaryParticle->SetTrackID(trackID); 145 primaryParticle->SetTrackID(trackID); 259 << 260 // Set parentID to 0 as a primary particle 146 // Set parentID to 0 as a primary particle 261 // << 262 track->SetParentID(0); 147 track->SetParentID(0); 263 << 264 // Set weight ( vertex weight * particle w 148 // Set weight ( vertex weight * particle weight ) 265 // << 266 track->SetWeight(wv*(primaryParticle->GetW 149 track->SetWeight(wv*(primaryParticle->GetWeight())); 267 << 268 // Store it to G4TrackVector 150 // Store it to G4TrackVector 269 // << 270 TV.push_back( track ); 151 TV.push_back( track ); 271 } 152 } 272 } 153 } 273 154 274 void G4PrimaryTransformer:: << 155 void G4PrimaryTransformer::SetDecayProducts 275 SetDecayProducts(G4PrimaryParticle* mother, G4 << 156 (G4PrimaryParticle* mother, G4DynamicParticle* motherDP) 276 { 157 { 277 G4PrimaryParticle* daughter = mother->GetDau 158 G4PrimaryParticle* daughter = mother->GetDaughter(); 278 if(daughter == nullptr) return; << 159 if(!daughter) return; 279 auto* decayProducts << 160 G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() ); 280 = (G4DecayProducts*)(motherDP->GetPreAssig << 161 if(!decayProducts) 281 if(decayProducts == nullptr) << 282 { 162 { 283 decayProducts = new G4DecayProducts(*mothe 163 decayProducts = new G4DecayProducts(*motherDP); 284 motherDP->SetPreAssignedDecayProducts(deca 164 motherDP->SetPreAssignedDecayProducts(decayProducts); 285 } 165 } 286 while(daughter != nullptr) << 166 while(daughter) 287 { 167 { 288 G4ParticleDefinition* partDef = GetDefinit 168 G4ParticleDefinition* partDef = GetDefinition(daughter); 289 if(!IsGoodForTrack(partDef)) << 169 if(!partDef) 290 { 170 { 291 #ifdef G4VERBOSE 171 #ifdef G4VERBOSE 292 if(verboseLevel>2) 172 if(verboseLevel>2) 293 { 173 { 294 G4cout << " >> Decay product (PDGcode 174 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode() 295 << ") --- Ignored" << G4endl; 175 << ") --- Ignored" << G4endl; 296 } 176 } 297 #endif 177 #endif 298 SetDecayProducts(daughter,motherDP); 178 SetDecayProducts(daughter,motherDP); 299 } 179 } 300 else 180 else 301 { 181 { 302 #ifdef G4VERBOSE 182 #ifdef G4VERBOSE 303 if(verboseLevel>1) 183 if(verboseLevel>1) 304 { 184 { 305 G4cout << " >> Decay product (" << par 185 G4cout << " >> Decay product (" << partDef->GetParticleName() 306 << ") --- Attached with momentu 186 << ") --- Attached with momentum " << daughter->GetMomentum() 307 << G4endl; 187 << G4endl; 308 } 188 } 309 #endif 189 #endif 310 auto* DP << 190 G4DynamicParticle*DP 311 = new G4DynamicParticle(partDef,daught 191 = new G4DynamicParticle(partDef,daughter->GetMomentum()); 312 DP->SetPrimaryParticle(daughter); 192 DP->SetPrimaryParticle(daughter); 313 << 314 // Decay proper time for daughter 193 // Decay proper time for daughter 315 // << 194 if(daughter->GetProperTime()>0.0) 316 if(daughter->GetProperTime()>=0.0) << 195 { DP->SetPreAssignedDecayProperTime(daughter->GetProperTime()); } 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); 196 decayProducts->PushProducts(DP); 339 SetDecayProducts(daughter,DP); 197 SetDecayProducts(daughter,DP); 340 << 341 // Check the particle is properly constr << 342 // << 343 if(!CheckDynamicParticle(DP)) << 344 { << 345 delete DP; << 346 return; << 347 } << 348 } 198 } 349 daughter = daughter->GetNext(); 199 daughter = daughter->GetNext(); 350 } 200 } 351 } 201 } 352 202 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 203 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 204 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 205 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 << 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 << 429 return false; << 430 } << 431 206