ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PrimaryTransformer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PrimaryTransformer.cc
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 //
28 
29 #include "G4PrimaryTransformer.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4Event.hh"
32 #include "G4PrimaryVertex.hh"
33 #include "G4ParticleDefinition.hh"
34 #include "G4DynamicParticle.hh"
35 #include "G4Track.hh"
36 #include "G4ThreeVector.hh"
37 #include "G4DecayProducts.hh"
38 #include "G4UnitsTable.hh"
39 #include "G4ios.hh"
40 #include "Randomize.hh"
41 
43 :verboseLevel(0),trackID(0),
44  unknown(nullptr),unknownParticleDefined(false),
45  opticalphoton(nullptr),opticalphotonDefined(false),
46  nWarn(0)
47 {
49  CheckUnknown();
50 }
51 
53 {;}
54 
56 {
57  unknown = particleTable->FindParticle("unknown");
58  if(unknown)
59  { unknownParticleDefined = true; }
60  else
61  { unknownParticleDefined = false; }
62  opticalphoton = particleTable->FindParticle("opticalphoton");
63  if(opticalphoton)
64  { opticalphotonDefined = true; }
65  else
66  { opticalphotonDefined = false; }
67 }
68 
70 {
71  trackID = trackIDCounter;
72 
73  //TV.clearAndDestroy();
74  for(auto tr : TV) delete tr;
75  TV.clear();
76 
77  //Loop over vertices
78  G4PrimaryVertex* nextVertex = anEvent->GetPrimaryVertex();
79  while(nextVertex) // Loop checking 12.28.2015 M.Asai
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 ) // Loop checking 12.28.2015 M.Asai
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) // Loop checking 12.28.2015 M.Asai
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  DP->SetPolarization(daughter->GetPolX(),daughter->GetPolY(),daughter->GetPolZ());
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