source: trunk/source/event/src/G4PrimaryTransformer.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 10.7 KB
Line 
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//
27// $Id: G4PrimaryTransformer.cc,v 1.27 2006/06/29 18:09:58 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30
31#include "G4PrimaryTransformer.hh"
32#include "G4Event.hh"
33#include "G4PrimaryVertex.hh"
34#include "G4ParticleDefinition.hh"
35#include "G4DynamicParticle.hh"
36#include "G4Track.hh"
37#include "G4ThreeVector.hh"
38#include "G4DecayProducts.hh"
39#include "G4UnitsTable.hh"
40#include "G4ios.hh"
41#include "Randomize.hh"
42
43G4PrimaryTransformer::G4PrimaryTransformer()
44:verboseLevel(0),trackID(0),unknown(0),unknownParticleDefined(false)
45{
46 particleTable = G4ParticleTable::GetParticleTable();
47 CheckUnknown();
48}
49
50G4PrimaryTransformer::~G4PrimaryTransformer()
51{;}
52
53void G4PrimaryTransformer::CheckUnknown()
54{
55 unknown = particleTable->FindParticle("unknown");
56 if(unknown)
57 { unknownParticleDefined = true; }
58 else
59 { unknownParticleDefined = false; }
60}
61
62G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent,G4int trackIDCounter)
63{
64 trackID = trackIDCounter;
65
66 //TV.clearAndDestroy();
67 for( size_t ii=0; ii<TV.size();ii++)
68 { delete TV[ii]; }
69 TV.clear();
70 G4int n_vertex = anEvent->GetNumberOfPrimaryVertex();
71 if(n_vertex==0) return 0;
72 for( G4int i=0; i<n_vertex; i++ )
73 { GenerateTracks( anEvent->GetPrimaryVertex(i) ); }
74 return &TV;
75}
76
77void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex)
78{
79 G4double X0 = primaryVertex->GetX0();
80 G4double Y0 = primaryVertex->GetY0();
81 G4double Z0 = primaryVertex->GetZ0();
82 G4double T0 = primaryVertex->GetT0();
83 G4double WV = primaryVertex->GetWeight();
84
85#ifdef G4VERBOSE
86 if(verboseLevel>1)
87 {
88 G4cout << "G4PrimaryTransformer::PrimaryVertex ("
89 << X0 / mm << "(mm),"
90 << Y0 / mm << "(mm),"
91 << Z0 / mm << "(mm),"
92 << T0 / nanosecond << "(nsec))" << G4endl;
93 }
94#endif
95
96 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary();
97 while( primaryParticle != 0 )
98 {
99 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV );
100 primaryParticle = primaryParticle->GetNext();
101 }
102}
103
104void G4PrimaryTransformer::GenerateSingleTrack
105 (G4PrimaryParticle* primaryParticle,
106 G4double x0,G4double y0,G4double z0,G4double t0,G4double wv)
107{
108 static G4ParticleDefinition* optPhoton = 0;
109 static G4int nWarn = 0;
110 if(!optPhoton) optPhoton = particleTable->FindParticle("opticalphoton");
111
112 G4ParticleDefinition* partDef = GetDefinition(primaryParticle);
113 if(!IsGoodForTrack(partDef))
114 // The particle cannot be converted to G4Track, check daughters
115 {
116#ifdef G4VERBOSE
117 if(verboseLevel>2)
118 {
119 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode()
120 << ") --- Ignored" << G4endl;
121 }
122#endif
123 G4PrimaryParticle* daughter = primaryParticle->GetDaughter();
124 while(daughter)
125 {
126 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv);
127 daughter = daughter->GetNext();
128 }
129 }
130
131 // The particle is defined in GEANT4
132 else
133 {
134 // Create G4DynamicParticle object
135#ifdef G4VERBOSE
136 if(verboseLevel>1)
137 {
138 G4cout << "Primary particle (" << partDef->GetParticleName()
139 << ") --- Transfered with momentum " << primaryParticle->GetMomentum()
140 << G4endl;
141 }
142#endif
143 G4DynamicParticle* DP =
144 new G4DynamicParticle(partDef,primaryParticle->GetMomentum());
145 if(partDef==optPhoton && primaryParticle->GetPolarization().mag2()==0.)
146 {
147 if(nWarn<10)
148 {
149 G4Exception("G4PrimaryTransformer::GenerateSingleTrack","ZeroPolarization",JustWarning,
150 "Polarization of the optical photon is null. Random polarization is assumed.");
151 G4cerr << "This warning message is issued up to 10 times." << G4endl;
152 nWarn++;
153 }
154
155 G4double angle = G4UniformRand() * 360.0*deg;
156 G4ThreeVector normal (1., 0., 0.);
157 G4ThreeVector kphoton = DP->GetMomentumDirection();
158 G4ThreeVector product = normal.cross(kphoton);
159 G4double modul2 = product*product;
160
161 G4ThreeVector e_perpend (0., 0., 1.);
162 if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
163 G4ThreeVector e_paralle = e_perpend.cross(kphoton);
164
165 G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
166 DP->SetPolarization(polar.x(),polar.y(),polar.z());
167 }
168 else
169 {
170 DP->SetPolarization(primaryParticle->GetPolX(),
171 primaryParticle->GetPolY(),
172 primaryParticle->GetPolZ());
173 }
174 if(primaryParticle->GetProperTime()>0.0)
175 { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); }
176 // Set Charge if it is specified
177 if (primaryParticle->GetCharge()<DBL_MAX) {
178 DP->SetCharge(primaryParticle->GetCharge());
179 }
180 // Set Mass if it is specified
181 G4double pmas = primaryParticle->GetMass();
182 if(pmas>=0.)
183 { DP->SetMass(pmas); }
184 // Set decay products to the DynamicParticle
185 SetDecayProducts( primaryParticle, DP );
186 // Set primary particle
187 DP->SetPrimaryParticle(primaryParticle);
188 // Set PDG code if it is different from G4ParticleDefinition
189 if(partDef->GetPDGEncoding()==0 && primaryParticle->GetPDGcode()!=0)
190 {
191 DP->SetPDGcode(primaryParticle->GetPDGcode());
192 }
193 // Check the particle is properly constructed
194 if(!CheckDynamicParticle(DP))
195 {
196 delete DP;
197 return;
198 }
199 // Create G4Track object
200 G4Track* track = new G4Track(DP,t0,G4ThreeVector(x0,y0,z0));
201 // Set trackID and let primary particle know it
202 trackID++;
203 track->SetTrackID(trackID);
204 primaryParticle->SetTrackID(trackID);
205 // Set parentID to 0 as a primary particle
206 track->SetParentID(0);
207 // Set weight ( vertex weight * particle weight )
208 track->SetWeight(wv*(primaryParticle->GetWeight()));
209 // Store it to G4TrackVector
210 TV.push_back( track );
211 }
212}
213
214void G4PrimaryTransformer::SetDecayProducts
215 (G4PrimaryParticle* mother, G4DynamicParticle* motherDP)
216{
217 G4PrimaryParticle* daughter = mother->GetDaughter();
218 if(!daughter) return;
219 G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() );
220 if(!decayProducts)
221 {
222 decayProducts = new G4DecayProducts(*motherDP);
223 motherDP->SetPreAssignedDecayProducts(decayProducts);
224 }
225 while(daughter)
226 {
227 G4ParticleDefinition* partDef = GetDefinition(daughter);
228 if(!IsGoodForTrack(partDef))
229 {
230#ifdef G4VERBOSE
231 if(verboseLevel>2)
232 {
233 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode()
234 << ") --- Ignored" << G4endl;
235 }
236#endif
237 SetDecayProducts(daughter,motherDP);
238 }
239 else
240 {
241#ifdef G4VERBOSE
242 if(verboseLevel>1)
243 {
244 G4cout << " >> Decay product (" << partDef->GetParticleName()
245 << ") --- Attached with momentum " << daughter->GetMomentum()
246 << G4endl;
247 }
248#endif
249 G4DynamicParticle*DP
250 = new G4DynamicParticle(partDef,daughter->GetMomentum());
251 DP->SetPrimaryParticle(daughter);
252 // Decay proper time for daughter
253 if(daughter->GetProperTime()>0.0)
254 { DP->SetPreAssignedDecayProperTime(daughter->GetProperTime()); }
255 // Set Charge is specified
256 if (daughter->GetCharge()<DBL_MAX) {
257 DP->SetCharge(daughter->GetCharge());
258 }
259 G4double pmas = daughter->GetMass();
260 if(pmas>=0.)
261 { DP->SetMass(pmas); }
262 decayProducts->PushProducts(DP);
263 SetDecayProducts(daughter,DP);
264 // Check the particle is properly constructed
265 if(!CheckDynamicParticle(DP))
266 {
267 delete DP;
268 return;
269 }
270 }
271 daughter = daughter->GetNext();
272 }
273}
274
275void G4PrimaryTransformer::SetUnknnownParticleDefined(G4bool vl)
276{
277 unknownParticleDefined = vl;
278 if(unknownParticleDefined && !unknown)
279 { G4cerr << "unknownParticleDefined cannot be set true because G4UnknownParticle is not defined in the physics list."
280 << G4endl << "Command ignored." << G4endl;
281 unknownParticleDefined = false;
282 }
283}
284
285G4bool G4PrimaryTransformer::CheckDynamicParticle(G4DynamicParticle*DP)
286{
287 if(IsGoodForTrack(DP->GetDefinition())) return true;
288 G4DecayProducts* decayProducts = (G4DecayProducts*)(DP->GetPreAssignedDecayProducts());
289 if(decayProducts && decayProducts->entries()>0) return true;
290 G4cerr << G4endl
291 << "G4PrimaryTransformer: a shortlived primary particle is found" << G4endl
292 << " without any valid decay table nor pre-assigned decay mode." << G4endl;
293 G4Exception("G4PrimaryTransformer","InvalidPrimary",JustWarning,
294 "This primary particle will be ignored.");
295 return false;
296}
297
298G4ParticleDefinition* G4PrimaryTransformer::GetDefinition(G4PrimaryParticle*pp)
299{
300 G4ParticleDefinition* partDef = pp->GetG4code();
301 if(!partDef) partDef = particleTable->FindParticle(pp->GetPDGcode());
302 if(unknownParticleDefined && ((!partDef)||partDef->IsShortLived())) partDef = unknown;
303 return partDef;
304}
305
306G4bool G4PrimaryTransformer::IsGoodForTrack(G4ParticleDefinition* pd)
307{
308 if(!pd)
309 { return false; }
310 else if(!(pd->IsShortLived()))
311 { return true; }
312// Following two lines should be removed if the user does not want to make shortlived
313// primary particle with proper decay table to be converted into a track.
314 else if(pd->GetDecayTable())
315 { return true; }
316 return false;
317}
318
Note: See TracBrowser for help on using the repository browser.