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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
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 $
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.