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