| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4PrimaryTransformer class implementation
27 //
28 // Author: Makoto Asai, 1999
29 // --------------------------------------------------------------------
30
31 #include "G4PrimaryTransformer.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4Event.hh"
34 #include "G4PrimaryVertex.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4DynamicParticle.hh"
37 #include "G4Track.hh"
38 #include "G4ThreeVector.hh"
39 #include "G4DecayProducts.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4ios.hh"
42 #include "Randomize.hh"
43
44 G4double G4PrimaryTransformer::kETolerance = 1.0 * CLHEP::MeV;
45 G4ExceptionSeverity G4PrimaryTransformer::kETSeverity = JustWarning;
46
47 G4PrimaryTransformer::G4PrimaryTransformer()
48 {
49 particleTable = G4ParticleTable::GetParticleTable();
50 CheckUnknown();
51 }
52
53 void G4PrimaryTransformer::CheckUnknown()
54 {
55 unknown = particleTable->FindParticle("unknown");
56 unknownParticleDefined = unknown != nullptr;
57 chargedunknown = particleTable->FindParticle("chargedunknown");
58 chargedUnknownParticleDefined = chargedunknown != nullptr;
59 opticalphoton = particleTable->FindParticle("opticalphoton");
60 opticalphotonDefined = opticalphoton != nullptr;
61 }
62
63 G4TrackVector*
64 G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent, G4int trackIDCounter)
65 {
66 trackID = trackIDCounter;
67
68 for(auto tr : TV) delete tr;
69 TV.clear();
70
71 // Loop over vertices
72 //
73 G4PrimaryVertex* nextVertex = anEvent->GetPrimaryVertex();
74 while(nextVertex != nullptr) // Loop checking 12.28.2015 M.Asai
75 {
76 GenerateTracks(nextVertex);
77 nextVertex = nextVertex->GetNext();
78 }
79 return &TV;
80 }
81
82 void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex)
83 {
84 G4double X0 = primaryVertex->GetX0();
85 G4double Y0 = primaryVertex->GetY0();
86 G4double Z0 = primaryVertex->GetZ0();
87 G4double T0 = primaryVertex->GetT0();
88 G4double WV = primaryVertex->GetWeight();
89
90 #ifdef G4VERBOSE
91 if(verboseLevel>2)
92 {
93 primaryVertex->Print();
94 }
95 else if (verboseLevel==1)
96 {
97 G4cout << "G4PrimaryTransformer::PrimaryVertex ("
98 << X0 / mm << "(mm),"
99 << Y0 / mm << "(mm),"
100 << Z0 / mm << "(mm),"
101 << T0 / nanosecond << "(nsec))" << G4endl;
102 }
103 #endif
104
105 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary();
106 while( primaryParticle != nullptr ) // Loop checking 12.28.2015 M.Asai
107 {
108 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV );
109 primaryParticle = primaryParticle->GetNext();
110 }
111 }
112
113 void G4PrimaryTransformer::
114 GenerateSingleTrack( G4PrimaryParticle* primaryParticle,
115 G4double x0, G4double y0, G4double z0,
116 G4double t0, G4double wv)
117 {
118 G4ParticleDefinition* partDef = GetDefinition(primaryParticle);
119 if(!IsGoodForTrack(partDef))
120 { // The particle cannot be converted to G4Track, check daughters
121 #ifdef G4VERBOSE
122 if(verboseLevel>1)
123 {
124 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode()
125 << ") --- Ignored" << G4endl;
126 }
127 #endif
128 G4PrimaryParticle* daughter = primaryParticle->GetDaughter();
129 while(daughter != nullptr) // Loop checking 12.28.2015 M.Asai
130 {
131 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv);
132 daughter = daughter->GetNext();
133 }
134 }
135 else // The particle is defined in GEANT4
136 {
137 // Create G4DynamicParticle object
138 #ifdef G4VERBOSE
139 if(verboseLevel>1)
140 {
141 G4cout << "Primary particle (" << partDef->GetParticleName()
142 << ") --- Transferred with momentum "
143 << primaryParticle->GetMomentum()
144 << G4endl;
145 }
146 #endif
147 // Check the mass of the "real" particle
148 // N.B. PDG code 0 is used for artificial particle such as geantino
149 if(primaryParticle->GetPDGcode() != 0 && kETSeverity != IgnoreTheIssue)
150 {
151 G4double pmas = partDef->GetPDGMass();
152 if(std::abs(pmas - primaryParticle->GetMass()) > kETolerance) {
153 G4ExceptionDescription ed;
154 ed << "Primary particle PDG=" << primaryParticle->GetPDGcode()
155 << " deltaMass(MeV)=" << (std::abs(pmas - primaryParticle->GetMass()))/CLHEP::MeV
156 << " is larger than the tolerance(MeV)=" << kETolerance/CLHEP::MeV
157 << "\n Specified mass(MeV)=" << (primaryParticle->GetMass())/CLHEP::MeV
158 << " while PDG mass(MEV)=" << pmas/CLHEP::MeV
159 << "\n To change the tolerance or the exception severity,"
160 << " use G4PrimaryTransformer::SetKETolerance() method.";
161 G4Exception("G4PrimaryParticle::Set4Momentum", "part0005", kETSeverity, ed);
162 }
163 }
164 auto* DP =
165 new G4DynamicParticle(partDef,
166 primaryParticle->GetMomentumDirection(),
167 primaryParticle->GetKineticEnergy());
168 if(opticalphotonDefined && partDef==opticalphoton
169 && primaryParticle->GetPolarization().mag2()==0.)
170 {
171 if(nWarn<10)
172 {
173 G4Exception("G4PrimaryTransformer::GenerateSingleTrack",
174 "ZeroPolarization", JustWarning,
175 "Polarization of the optical photon is null.\
176 Random polarization is assumed.");
177 G4cerr << "This warning message is issued up to 10 times." << G4endl;
178 ++nWarn;
179 }
180
181 G4double angle = G4UniformRand() * 360.0*deg;
182 G4ThreeVector normal (1., 0., 0.);
183 G4ThreeVector kphoton = DP->GetMomentumDirection();
184 G4ThreeVector product = normal.cross(kphoton);
185 G4double modul2 = product*product;
186
187 G4ThreeVector e_perpend (0., 0., 1.);
188 if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
189 G4ThreeVector e_paralle = e_perpend.cross(kphoton);
190
191 G4ThreeVector polar = std::cos(angle)*e_paralle
192 + std::sin(angle)*e_perpend;
193 DP->SetPolarization(polar.x(),polar.y(),polar.z());
194 }
195 else
196 {
197 DP->SetPolarization(primaryParticle->GetPolX(),
198 primaryParticle->GetPolY(),
199 primaryParticle->GetPolZ());
200 }
201 if(primaryParticle->GetProperTime()>=0.0)
202 {
203 DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime());
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->GetCharge());
218 }
219 else // ions
220 {
221 G4int iz = partDef->GetAtomicNumber();
222 auto iq = static_cast<G4int>(primaryParticle->GetCharge()/eplus);
223 G4int n_e = iz - iq;
224 if (n_e>0) DP->AddElectron(0,n_e);
225 }
226 }
227 // Set decay products to the DynamicParticle
228 //
229 SetDecayProducts( primaryParticle, DP );
230
231 // Set primary particle
232 //
233 DP->SetPrimaryParticle(primaryParticle);
234
235 // Set PDG code if it is different from G4ParticleDefinition
236 //
237 if(partDef->GetPDGEncoding()==0 && primaryParticle->GetPDGcode()!=0)
238 {
239 DP->SetPDGcode(primaryParticle->GetPDGcode());
240 }
241
242 // Check the particle is properly constructed
243 //
244 if(!CheckDynamicParticle(DP))
245 {
246 delete DP;
247 return;
248 }
249
250 // Create G4Track object
251 //
252 auto track = new G4Track(DP,t0,G4ThreeVector(x0,y0,z0));
253
254 // Set trackID and let primary particle know it
255 //
256 ++trackID;
257 track->SetTrackID(trackID);
258 primaryParticle->SetTrackID(trackID);
259
260 // Set parentID to 0 as a primary particle
261 //
262 track->SetParentID(0);
263
264 // Set weight ( vertex weight * particle weight )
265 //
266 track->SetWeight(wv*(primaryParticle->GetWeight()));
267
268 // Store it to G4TrackVector
269 //
270 TV.push_back( track );
271 }
272 }
273
274 void G4PrimaryTransformer::
275 SetDecayProducts(G4PrimaryParticle* mother, G4DynamicParticle* motherDP)
276 {
277 G4PrimaryParticle* daughter = mother->GetDaughter();
278 if(daughter == nullptr) return;
279 auto* decayProducts
280 = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() );
281 if(decayProducts == nullptr)
282 {
283 decayProducts = new G4DecayProducts(*motherDP);
284 motherDP->SetPreAssignedDecayProducts(decayProducts);
285 }
286 while(daughter != nullptr)
287 {
288 G4ParticleDefinition* partDef = GetDefinition(daughter);
289 if(!IsGoodForTrack(partDef))
290 {
291 #ifdef G4VERBOSE
292 if(verboseLevel>2)
293 {
294 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode()
295 << ") --- Ignored" << G4endl;
296 }
297 #endif
298 SetDecayProducts(daughter,motherDP);
299 }
300 else
301 {
302 #ifdef G4VERBOSE
303 if(verboseLevel>1)
304 {
305 G4cout << " >> Decay product (" << partDef->GetParticleName()
306 << ") --- Attached with momentum " << daughter->GetMomentum()
307 << G4endl;
308 }
309 #endif
310 auto* DP
311 = 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(daughter->GetProperTime());
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);
339 SetDecayProducts(daughter,DP);
340
341 // Check the particle is properly constructed
342 //
343 if(!CheckDynamicParticle(DP))
344 {
345 delete DP;
346 return;
347 }
348 }
349 daughter = daughter->GetNext();
350 }
351 }
352
353 void G4PrimaryTransformer::SetUnknownParticleDefined(G4bool vl)
354 {
355 unknownParticleDefined = vl;
356 if(unknownParticleDefined && (unknown == nullptr))
357 {
358 G4cerr << "unknownParticleDefined cannot be set true because" << G4endl
359 << "G4UnknownParticle is not defined in the physics list." << G4endl
360 << "Command ignored." << G4endl;
361 unknownParticleDefined = false;
362 }
363 }
364
365 void G4PrimaryTransformer::SetChargedUnknownParticleDefined(G4bool vl)
366 {
367 chargedUnknownParticleDefined = vl;
368 if(chargedUnknownParticleDefined && (chargedunknown == nullptr))
369 {
370 G4cerr << "chargedUnknownParticleDefined cannot be set true because" << G4endl
371 << "G4ChargedUnknownParticle is not defined in the physics list." << G4endl
372 << "Command ignored." << G4endl;
373 chargedUnknownParticleDefined = false;
374 }
375 }
376
377 G4bool G4PrimaryTransformer::CheckDynamicParticle(G4DynamicParticle* DP)
378 {
379 if(IsGoodForTrack(DP->GetDefinition())) return true;
380 auto* decayProducts
381 = (G4DecayProducts*)(DP->GetPreAssignedDecayProducts());
382 if(decayProducts != nullptr && decayProducts->entries()>0) return true;
383 G4cerr << G4endl
384 << "G4PrimaryTransformer: a shortlived primary particle is found"
385 << G4endl
386 << " without any valid decay table nor pre-assigned decay mode."
387 << G4endl;
388 G4Exception("G4PrimaryTransformer", "InvalidPrimary", JustWarning,
389 "This primary particle will be ignored.");
390 return false;
391 }
392
393 G4ParticleDefinition*
394 G4PrimaryTransformer::GetDefinition(G4PrimaryParticle* pp)
395 {
396 G4ParticleDefinition* partDef = pp->GetG4code();
397 if(partDef == nullptr)
398 {
399 partDef = particleTable->FindParticle(pp->GetPDGcode());
400 }
401 if((partDef == nullptr) || partDef->IsShortLived())
402 {
403 if (chargedUnknownParticleDefined && pp->GetCharge()!=0.0)
404 {
405 partDef = chargedunknown;
406 }
407 else if (unknownParticleDefined)
408 {
409 partDef = unknown;
410 }
411 }
412 return partDef;
413 }
414
415 G4bool G4PrimaryTransformer::IsGoodForTrack(G4ParticleDefinition* pd)
416 {
417 if(pd == nullptr)
418 { return false; }
419 if(!(pd->IsShortLived()))
420 { return true; }
421 //
422 // Following two lines should be removed if the user does not want to make
423 // shortlived primary particle with proper decay table to be converted into
424 // a track.
425 //
426 if(pd->GetDecayTable() != nullptr)
427 { return true; }
428
429 return false;
430 }
431