Geant4  10.01.p03
G4PrimaryTransformer.cc
Go to the documentation of this file.
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 72252 2013-07-12 09:04:11Z gcosmo $
28 //
29 
30 #include "G4PrimaryTransformer.hh"
31 #include "G4SystemOfUnits.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 
44 :verboseLevel(0),trackID(0),
45  unknown(0),unknownParticleDefined(false),
46  opticalphoton(0),opticalphotonDefined(false),
47  nWarn(0)
48 {
50  CheckUnknown();
51 }
52 
54 {;}
55 
57 {
58  unknown = particleTable->FindParticle("unknown");
59  if(unknown)
60  { unknownParticleDefined = true; }
61  else
62  { unknownParticleDefined = false; }
63  opticalphoton = particleTable->FindParticle("opticalphoton");
64  if(opticalphoton)
65  { opticalphotonDefined = true; }
66  else
67  { opticalphotonDefined = false; }
68 }
69 
71 {
72  trackID = trackIDCounter;
73 
74  //TV.clearAndDestroy();
75  for( size_t ii=0; ii<TV.size();ii++)
76  { delete TV[ii]; }
77  TV.clear();
78  G4PrimaryVertex* nextVertex = anEvent->GetPrimaryVertex();
79  while(nextVertex)
80  {
81  GenerateTracks(nextVertex);
82  nextVertex = nextVertex->GetNext();
83  }
84  return &TV;
85 }
86 
88 {
89  G4double X0 = primaryVertex->GetX0();
90  G4double Y0 = primaryVertex->GetY0();
91  G4double Z0 = primaryVertex->GetZ0();
92  G4double T0 = primaryVertex->GetT0();
93  G4double WV = primaryVertex->GetWeight();
94 
95 #ifdef G4VERBOSE
96  if(verboseLevel>2) {
97  primaryVertex->Print();
98  } else if (verboseLevel==1) {
99  G4cout << "G4PrimaryTransformer::PrimaryVertex ("
100  << X0 / mm << "(mm),"
101  << Y0 / mm << "(mm),"
102  << Z0 / mm << "(mm),"
103  << T0 / nanosecond << "(nsec))" << G4endl;
104  }
105 #endif
106 
107  G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary();
108  while( primaryParticle != 0 )
109  {
110  GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV );
111  primaryParticle = primaryParticle->GetNext();
112  }
113 }
114 
116  (G4PrimaryParticle* primaryParticle,
118 {
119  G4ParticleDefinition* partDef = GetDefinition(primaryParticle);
120  if(!IsGoodForTrack(partDef))
121  // The particle cannot be converted to G4Track, check daughters
122  {
123 #ifdef G4VERBOSE
124  if(verboseLevel>2)
125  {
126  G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode()
127  << ") --- Ignored" << G4endl;
128  }
129 #endif
130  G4PrimaryParticle* daughter = primaryParticle->GetDaughter();
131  while(daughter)
132  {
133  GenerateSingleTrack(daughter,x0,y0,z0,t0,wv);
134  daughter = daughter->GetNext();
135  }
136  }
137 
138  // The particle is defined in GEANT4
139  else
140  {
141  // Create G4DynamicParticle object
142 #ifdef G4VERBOSE
143  if(verboseLevel>1)
144  {
145  G4cout << "Primary particle (" << partDef->GetParticleName()
146  << ") --- Transfered with momentum " << primaryParticle->GetMomentum()
147  << G4endl;
148  }
149 #endif
150  G4DynamicParticle* DP =
151  new G4DynamicParticle(partDef,
152  primaryParticle->GetMomentumDirection(),
153  primaryParticle->GetKineticEnergy());
154  if(opticalphotonDefined && partDef==opticalphoton && primaryParticle->GetPolarization().mag2()==0.)
155  {
156  if(nWarn<10)
157  {
158  G4Exception("G4PrimaryTransformer::GenerateSingleTrack","ZeroPolarization",JustWarning,
159  "Polarization of the optical photon is null. Random polarization is assumed.");
160  G4cerr << "This warning message is issued up to 10 times." << G4endl;
161  nWarn++;
162  }
163 
164  G4double angle = G4UniformRand() * 360.0*deg;
165  G4ThreeVector normal (1., 0., 0.);
166  G4ThreeVector kphoton = DP->GetMomentumDirection();
167  G4ThreeVector product = normal.cross(kphoton);
168  G4double modul2 = product*product;
169 
170  G4ThreeVector e_perpend (0., 0., 1.);
171  if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
172  G4ThreeVector e_paralle = e_perpend.cross(kphoton);
173 
174  G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
175  DP->SetPolarization(polar.x(),polar.y(),polar.z());
176  }
177  else
178  {
179  DP->SetPolarization(primaryParticle->GetPolX(),
180  primaryParticle->GetPolY(),
181  primaryParticle->GetPolZ());
182  }
183  if(primaryParticle->GetProperTime()>0.0)
184  { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); }
185 
186  // Set Mass if it is specified
187  G4double pmas = primaryParticle->GetMass();
188  if(pmas>=0.)
189  { DP->SetMass(pmas); }
190 
191  // Set Charge if it is specified
192  if (primaryParticle->GetCharge()<DBL_MAX) {
193  if (partDef->GetAtomicNumber() <0) {
194  DP->SetCharge(primaryParticle->GetCharge());
195  } else {
196  // ions
197  G4int iz = partDef->GetAtomicNumber();
198  G4int iq = static_cast<int>(primaryParticle->GetCharge()/eplus);
199  G4int n_e = iz - iq;
200  if (n_e>0) DP->AddElectron(0,n_e);
201  }
202  }
203  // Set decay products to the DynamicParticle
204  SetDecayProducts( primaryParticle, DP );
205  // Set primary particle
206  DP->SetPrimaryParticle(primaryParticle);
207  // Set PDG code if it is different from G4ParticleDefinition
208  if(partDef->GetPDGEncoding()==0 && primaryParticle->GetPDGcode()!=0)
209  {
210  DP->SetPDGcode(primaryParticle->GetPDGcode());
211  }
212  // Check the particle is properly constructed
213  if(!CheckDynamicParticle(DP))
214  {
215  delete DP;
216  return;
217  }
218  // Create G4Track object
219  G4Track* track = new G4Track(DP,t0,G4ThreeVector(x0,y0,z0));
220  // Set trackID and let primary particle know it
221  trackID++;
222  track->SetTrackID(trackID);
223  primaryParticle->SetTrackID(trackID);
224  // Set parentID to 0 as a primary particle
225  track->SetParentID(0);
226  // Set weight ( vertex weight * particle weight )
227  track->SetWeight(wv*(primaryParticle->GetWeight()));
228  // Store it to G4TrackVector
229  TV.push_back( track );
230 
231  }
232 }
233 
236 {
237  G4PrimaryParticle* daughter = mother->GetDaughter();
238  if(!daughter) return;
239  G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() );
240  if(!decayProducts)
241  {
242  decayProducts = new G4DecayProducts(*motherDP);
243  motherDP->SetPreAssignedDecayProducts(decayProducts);
244  }
245  while(daughter)
246  {
247  G4ParticleDefinition* partDef = GetDefinition(daughter);
248  if(!IsGoodForTrack(partDef))
249  {
250 #ifdef G4VERBOSE
251  if(verboseLevel>2)
252  {
253  G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode()
254  << ") --- Ignored" << G4endl;
255  }
256 #endif
257  SetDecayProducts(daughter,motherDP);
258  }
259  else
260  {
261 #ifdef G4VERBOSE
262  if(verboseLevel>1)
263  {
264  G4cout << " >> Decay product (" << partDef->GetParticleName()
265  << ") --- Attached with momentum " << daughter->GetMomentum()
266  << G4endl;
267  }
268 #endif
270  = new G4DynamicParticle(partDef,daughter->GetMomentum());
271  DP->SetPrimaryParticle(daughter);
272  // Decay proper time for daughter
273  if(daughter->GetProperTime()>0.0)
274  { DP->SetPreAssignedDecayProperTime(daughter->GetProperTime()); }
275  // Set Charge is specified
276  if (daughter->GetCharge()<DBL_MAX) {
277  DP->SetCharge(daughter->GetCharge());
278  }
279  G4double pmas = daughter->GetMass();
280  if(pmas>=0.)
281  { DP->SetMass(pmas); }
282  decayProducts->PushProducts(DP);
283  SetDecayProducts(daughter,DP);
284  // Check the particle is properly constructed
285  if(!CheckDynamicParticle(DP))
286  {
287  delete DP;
288  return;
289  }
290  }
291  daughter = daughter->GetNext();
292  }
293 }
294 
296 {
299  { G4cerr << "unknownParticleDefined cannot be set true because G4UnknownParticle is not defined in the physics list."
300  << G4endl << "Command ignored." << G4endl;
301  unknownParticleDefined = false;
302  }
303 }
304 
306 {
307  if(IsGoodForTrack(DP->GetDefinition())) return true;
309  if(decayProducts && decayProducts->entries()>0) return true;
310  G4cerr << G4endl
311  << "G4PrimaryTransformer: a shortlived primary particle is found" << G4endl
312  << " without any valid decay table nor pre-assigned decay mode." << G4endl;
313  G4Exception("G4PrimaryTransformer","InvalidPrimary",JustWarning,
314  "This primary particle will be ignored.");
315  return false;
316 }
317 
319 {
320  G4ParticleDefinition* partDef = pp->GetG4code();
321  if(!partDef) partDef = particleTable->FindParticle(pp->GetPDGcode());
322  if(unknownParticleDefined && ((!partDef)||partDef->IsShortLived())) partDef = unknown;
323  return partDef;
324 }
325 
327 {
328  if(!pd)
329  { return false; }
330  else if(!(pd->IsShortLived()))
331  { return true; }
332 // Following two lines should be removed if the user does not want to make shortlived
333 // primary particle with proper decay table to be converted into a track.
334  else if(pd->GetDecayTable())
335  { return true; }
336  return false;
337 }
338 
G4double GetX0() const
G4double GetPolX() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4ThreeVector GetMomentum() const
G4TrackVector * GimmePrimaries(G4Event *anEvent, G4int trackIDCounter=0)
void SetPreAssignedDecayProperTime(G4double)
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
G4double GetKineticEnergy() const
void Print() const
G4ThreeVector GetPolarization() const
const G4DecayProducts * GetPreAssignedDecayProducts() const
G4double GetPolZ() const
G4bool CheckDynamicParticle(G4DynamicParticle *DP)
G4double GetT0() const
G4ParticleDefinition * GetDefinition() const
G4double GetPolY() const
int G4int
Definition: G4Types.hh:78
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * opticalphoton
const G4String & GetParticleName() const
G4ParticleDefinition * GetG4code() const
G4int GetAtomicNumber() const
void SetWeight(G4double aValue)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double GetWeight() const
G4DecayTable * GetDecayTable() const
G4PrimaryParticle * GetPrimary(G4int i=0) const
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
void SetUnknnownParticleDefined(G4bool vl)
static const double deg
Definition: G4SIunits.hh:133
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double iz
Definition: TRTMaterials.hh:39
G4PrimaryParticle * GetDaughter() const
void SetTrackID(G4int id)
G4double GetZ0() const
void SetPrimaryParticle(G4PrimaryParticle *p)
G4PrimaryVertex * GetPrimaryVertex(G4int i=0) const
Definition: G4Event.hh:167
void SetPolarization(G4double polX, G4double polY, G4double polZ)
G4PrimaryVertex * GetNext() const
void GenerateTracks(G4PrimaryVertex *primaryVertex)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetProperTime() const
virtual G4bool IsGoodForTrack(G4ParticleDefinition *pd)
G4int GetPDGcode() const
G4double GetWeight() const
void AddElectron(G4int orbit, G4int number=1)
static const double nanosecond
Definition: G4SIunits.hh:139
static G4ParticleTable * GetParticleTable()
G4PrimaryParticle * GetNext() const
void GenerateSingleTrack(G4PrimaryParticle *primaryParticle, G4double x0, G4double y0, G4double z0, G4double t0, G4double wv)
std::vector< G4Track * > G4TrackVector
void SetParentID(const G4int aValue)
G4double GetCharge() const
G4double GetY0() const
G4ParticleTable * particleTable
G4ParticleDefinition * unknown
#define G4endl
Definition: G4ios.hh:61
virtual G4ParticleDefinition * GetDefinition(G4PrimaryParticle *pp)
G4int entries() const
void SetPDGcode(G4int c)
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
void SetPreAssignedDecayProducts(G4DecayProducts *aDecayProducts)
void SetCharge(G4double charge)
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:102
void SetMass(G4double mass)
void SetTrackID(const G4int aValue)
G4GLOB_DLL std::ostream G4cerr
void SetDecayProducts(G4PrimaryParticle *mother, G4DynamicParticle *motherDP)