| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th <<
6 // * the Geant4 Collaboration. It is provided <<
7 // * conditions of the Geant4 Software License <<
8 // * LICENSE and available at http://cern.ch/ <<
9 // * include a list of copyright holders. <<
10 // * <<
11 // * Neither the authors of this software syst <<
12 // * institutes,nor the agencies providing fin <<
13 // * work make any representation or warran <<
14 // * regarding this software system or assum <<
15 // * use. Please see the license in the file <<
16 // * for the full disclaimer and the limitatio <<
17 // * <<
18 // * This code implementation is the result <<
19 // * technical work of the GEANT4 collaboratio <<
20 // * By using, copying, modifying or distri <<
21 // * any work based on the software) you ag <<
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // 7 //
26 // G4PrimaryTransformer class implementation << 8 // $Id: G4PrimaryTransformer.cc,v 1.8 2001/02/09 00:11:31 asaim Exp $
>> 9 // GEANT4 tag $Name: geant4-03-01 $
27 // 10 //
28 // Author: Makoto Asai, 1999 <<
29 // ------------------------------------------- <<
30 11
31 #include "G4PrimaryTransformer.hh" 12 #include "G4PrimaryTransformer.hh"
32 #include "G4SystemOfUnits.hh" <<
33 #include "G4Event.hh" 13 #include "G4Event.hh"
34 #include "G4PrimaryVertex.hh" 14 #include "G4PrimaryVertex.hh"
35 #include "G4ParticleDefinition.hh" 15 #include "G4ParticleDefinition.hh"
36 #include "G4DynamicParticle.hh" 16 #include "G4DynamicParticle.hh"
37 #include "G4Track.hh" 17 #include "G4Track.hh"
38 #include "G4ThreeVector.hh" 18 #include "G4ThreeVector.hh"
39 #include "G4DecayProducts.hh" 19 #include "G4DecayProducts.hh"
40 #include "G4UnitsTable.hh" <<
41 #include "G4ios.hh" 20 #include "G4ios.hh"
42 #include "Randomize.hh" <<
43 <<
44 G4double G4PrimaryTransformer::kETolerance = 1 <<
45 G4ExceptionSeverity G4PrimaryTransformer::kETS <<
46 21
47 G4PrimaryTransformer::G4PrimaryTransformer() 22 G4PrimaryTransformer::G4PrimaryTransformer()
>> 23 :verboseLevel(0)
48 { 24 {
49 particleTable = G4ParticleTable::GetParticle 25 particleTable = G4ParticleTable::GetParticleTable();
50 CheckUnknown(); <<
51 } 26 }
52 27
53 void G4PrimaryTransformer::CheckUnknown() << 28 G4PrimaryTransformer::~G4PrimaryTransformer()
54 { << 29 {;}
55 unknown = particleTable->FindParticle("unkno << 30
56 unknownParticleDefined = unknown != nullptr; << 31 G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent)
57 chargedunknown = particleTable->FindParticle << 32 {
58 chargedUnknownParticleDefined = chargedunkno << 33 //TV.clearAndDestroy();
59 opticalphoton = particleTable->FindParticle( << 34 for( int ii=0; ii<TV.size();ii++)
60 opticalphotonDefined = opticalphoton != null << 35 { delete TV[ii]; }
61 } <<
62 <<
63 G4TrackVector* <<
64 G4PrimaryTransformer::GimmePrimaries(G4Event* <<
65 { <<
66 trackID = trackIDCounter; <<
67 <<
68 for(auto tr : TV) delete tr; <<
69 TV.clear(); 36 TV.clear();
70 << 37 G4int n_vertex = anEvent->GetNumberOfPrimaryVertex();
71 // Loop over vertices << 38 if(n_vertex==0) return NULL;
72 // << 39 for( int i=0; i<n_vertex; i++ )
73 G4PrimaryVertex* nextVertex = anEvent->GetPr << 40 { GenerateTracks( anEvent->GetPrimaryVertex(i) ); }
74 while(nextVertex != nullptr) // Loop checkin <<
75 { <<
76 GenerateTracks(nextVertex); <<
77 nextVertex = nextVertex->GetNext(); <<
78 } <<
79 return &TV; 41 return &TV;
80 } 42 }
81 43
82 void G4PrimaryTransformer::GenerateTracks(G4Pr 44 void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex)
83 { 45 {
84 G4double X0 = primaryVertex->GetX0(); 46 G4double X0 = primaryVertex->GetX0();
85 G4double Y0 = primaryVertex->GetY0(); 47 G4double Y0 = primaryVertex->GetY0();
86 G4double Z0 = primaryVertex->GetZ0(); 48 G4double Z0 = primaryVertex->GetZ0();
87 G4double T0 = primaryVertex->GetT0(); 49 G4double T0 = primaryVertex->GetT0();
88 G4double WV = primaryVertex->GetWeight(); 50 G4double WV = primaryVertex->GetWeight();
89 51
90 #ifdef G4VERBOSE 52 #ifdef G4VERBOSE
91 if(verboseLevel>2) << 53 if(verboseLevel>1)
92 { 54 {
93 primaryVertex->Print(); <<
94 } <<
95 else if (verboseLevel==1) <<
96 { <<
97 G4cout << "G4PrimaryTransformer::PrimaryVe 55 G4cout << "G4PrimaryTransformer::PrimaryVertex ("
98 << X0 / mm << "(mm)," 56 << X0 / mm << "(mm),"
99 << Y0 / mm << "(mm)," 57 << Y0 / mm << "(mm),"
100 << Z0 / mm << "(mm)," 58 << Z0 / mm << "(mm),"
101 << T0 / nanosecond << "(nsec))" << 59 << T0 / nanosecond << "(nsec))" << G4endl;
102 } 60 }
103 #endif 61 #endif
104 62
105 G4PrimaryParticle* primaryParticle = primary 63 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary();
106 while( primaryParticle != nullptr ) // Loop << 64 while( primaryParticle != NULL )
107 { 65 {
108 GenerateSingleTrack( primaryParticle, X0, 66 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV );
109 primaryParticle = primaryParticle->GetNext 67 primaryParticle = primaryParticle->GetNext();
110 } 68 }
111 } 69 }
112 70
113 void G4PrimaryTransformer:: << 71 void G4PrimaryTransformer::GenerateSingleTrack
114 GenerateSingleTrack( G4PrimaryParticle* primar << 72 (G4PrimaryParticle* primaryParticle,
115 G4double x0, G4double y0, << 73 G4double x0,G4double y0,G4double z0,G4double t0,G4double wv)
116 G4double t0, G4double wv) <<
117 { 74 {
118 G4ParticleDefinition* partDef = GetDefinitio 75 G4ParticleDefinition* partDef = GetDefinition(primaryParticle);
119 if(!IsGoodForTrack(partDef)) << 76 if((!partDef)||partDef->IsShortLived())
120 { // The particle cannot be converted to G4 << 77 // The particle is not defined in GEANT4, check daughters
>> 78 {
121 #ifdef G4VERBOSE 79 #ifdef G4VERBOSE
122 if(verboseLevel>1) << 80 if(verboseLevel>2)
123 { 81 {
124 G4cout << "Primary particle (PDGcode " < 82 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode()
125 << ") --- Ignored" << G4endl; 83 << ") --- Ignored" << G4endl;
126 } 84 }
127 #endif 85 #endif
128 G4PrimaryParticle* daughter = primaryParti 86 G4PrimaryParticle* daughter = primaryParticle->GetDaughter();
129 while(daughter != nullptr) // Loop checkin << 87 while(daughter)
130 { 88 {
131 GenerateSingleTrack(daughter,x0,y0,z0,t0 89 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv);
132 daughter = daughter->GetNext(); 90 daughter = daughter->GetNext();
133 } 91 }
134 } 92 }
135 else // The particle is defined in GEANT4 << 93
>> 94 // The particle is defined in GEANT4
>> 95 else
136 { 96 {
137 // Create G4DynamicParticle object 97 // Create G4DynamicParticle object
138 #ifdef G4VERBOSE 98 #ifdef G4VERBOSE
139 if(verboseLevel>1) 99 if(verboseLevel>1)
140 { 100 {
141 G4cout << "Primary particle (" << partDe 101 G4cout << "Primary particle (" << partDef->GetParticleName()
142 << ") --- Transferred with moment << 102 << ") --- Transfered with momentum " << primaryParticle->GetMomentum()
143 << primaryParticle->GetMomentum() <<
144 << G4endl; 103 << G4endl;
145 } 104 }
146 #endif 105 #endif
147 // Check the mass of the "real" particle << 106 G4DynamicParticle* DP =
148 // N.B. PDG code 0 is used for artificial << 107 new G4DynamicParticle(partDef,primaryParticle->GetMomentum());
149 if(primaryParticle->GetPDGcode() != 0 && k << 108 DP->SetPolarization(primaryParticle->GetPolX(),
150 { << 109 primaryParticle->GetPolY(),
151 G4double pmas = partDef->GetPDGMass(); << 110 primaryParticle->GetPolZ());
152 if(std::abs(pmas - primaryParticle->GetM << 111 if(primaryParticle->GetProperTime()>0.0)
153 G4ExceptionDescription ed; << 112 { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); }
154 ed << "Primary particle PDG=" << prima << 113 // Set Charge
155 << " deltaMass(MeV)=" << (std::abs( << 114 if (abs(primaryParticle->GetCharge()-DP->GetCharge())>eplus) {
156 << " is larger than the tolerance(M << 115 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 } 116 }
227 // Set decay products to the DynamicPartic 117 // Set decay products to the DynamicParticle
228 // <<
229 SetDecayProducts( primaryParticle, DP ); 118 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 119 // Create G4Track object
251 // << 120 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 <<
255 // <<
256 ++trackID; <<
257 track->SetTrackID(trackID); <<
258 primaryParticle->SetTrackID(trackID); <<
259 <<
260 // Set parentID to 0 as a primary particle 121 // Set parentID to 0 as a primary particle
261 // <<
262 track->SetParentID(0); 122 track->SetParentID(0);
263 <<
264 // Set weight ( vertex weight * particle w 123 // Set weight ( vertex weight * particle weight )
265 // <<
266 track->SetWeight(wv*(primaryParticle->GetW 124 track->SetWeight(wv*(primaryParticle->GetWeight()));
267 <<
268 // Store it to G4TrackVector 125 // Store it to G4TrackVector
269 // <<
270 TV.push_back( track ); 126 TV.push_back( track );
271 } 127 }
272 } 128 }
273 129
274 void G4PrimaryTransformer:: << 130 void G4PrimaryTransformer::SetDecayProducts
275 SetDecayProducts(G4PrimaryParticle* mother, G4 << 131 (G4PrimaryParticle* mother, G4DynamicParticle* motherDP)
276 { 132 {
277 G4PrimaryParticle* daughter = mother->GetDau 133 G4PrimaryParticle* daughter = mother->GetDaughter();
278 if(daughter == nullptr) return; << 134 if(!daughter) return;
279 auto* decayProducts << 135 G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() );
280 = (G4DecayProducts*)(motherDP->GetPreAssig << 136 if(!decayProducts)
281 if(decayProducts == nullptr) <<
282 { 137 {
283 decayProducts = new G4DecayProducts(*mothe 138 decayProducts = new G4DecayProducts(*motherDP);
284 motherDP->SetPreAssignedDecayProducts(deca 139 motherDP->SetPreAssignedDecayProducts(decayProducts);
285 } 140 }
286 while(daughter != nullptr) << 141 while(daughter)
287 { 142 {
288 G4ParticleDefinition* partDef = GetDefinit 143 G4ParticleDefinition* partDef = GetDefinition(daughter);
289 if(!IsGoodForTrack(partDef)) << 144 if(!partDef)
290 { 145 {
291 #ifdef G4VERBOSE 146 #ifdef G4VERBOSE
292 if(verboseLevel>2) 147 if(verboseLevel>2)
293 { 148 {
294 G4cout << " >> Decay product (PDGcode 149 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode()
295 << ") --- Ignored" << G4endl; 150 << ") --- Ignored" << G4endl;
296 } 151 }
297 #endif 152 #endif
298 SetDecayProducts(daughter,motherDP); 153 SetDecayProducts(daughter,motherDP);
299 } 154 }
300 else 155 else
301 { 156 {
302 #ifdef G4VERBOSE 157 #ifdef G4VERBOSE
303 if(verboseLevel>1) 158 if(verboseLevel>1)
304 { 159 {
305 G4cout << " >> Decay product (" << par 160 G4cout << " >> Decay product (" << partDef->GetParticleName()
306 << ") --- Attached with momentu 161 << ") --- Attached with momentum " << daughter->GetMomentum()
307 << G4endl; 162 << G4endl;
308 } 163 }
309 #endif 164 #endif
310 auto* DP << 165 G4DynamicParticle*DP
311 = new G4DynamicParticle(partDef,daught 166 = 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); 167 decayProducts->PushProducts(DP);
339 SetDecayProducts(daughter,DP); 168 SetDecayProducts(daughter,DP);
340 <<
341 // Check the particle is properly constr <<
342 // <<
343 if(!CheckDynamicParticle(DP)) <<
344 { <<
345 delete DP; <<
346 return; <<
347 } <<
348 } 169 }
349 daughter = daughter->GetNext(); 170 daughter = daughter->GetNext();
350 } 171 }
351 } 172 }
352 173
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 174
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 175
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 176
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 177