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

Last change on this file since 1345 was 1337, checked in by garnier, 15 years ago

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