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

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
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-cand-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.