TGeant4/TG4ExtDecayer.cxx

   1 // $Id$
   2 // Category: physics
   3 //
   4 // See the class description in the header file.
   5
   6 #include "AliDecayer.h"
   7
   8 #include "TG4ExtDecayer.h"
   9 #include "TG4ParticlesManager.h"
  10 #include "TG4G3Units.h"
  11
  12 #include <G4DynamicParticle.hh>
  13 #include <G4DecayProducts.hh>
  14 #include <G4DecayTable.hh>
  15 #include <G4ParticleTable.hh>
  16 #include <G4DynamicParticle.hh>
  17 #include <G4Track.hh>
  18
  19 #include <TParticle.h>
  20 #include <TLorentzVector.h>
  21 #include <TClonesArray.h>
  22
  23 #include <math.h>
  24
  25 //_____________________________________________________________________________
  26 TG4ExtDecayer::TG4ExtDecayer(AliDecayer* externalDecayer)
  27   : G4VExtDecayer("TG4ExtDecayer"),
  28     fExternalDecayer(externalDecayer),
  29     fVerboseLevel(0) {
  30 //
  31   fDecayProductsArray = new  TClonesArray("TParticle", 1000);
  32   fParticlesManager = TG4ParticlesManager::Instance();
  33 }
  34
  35 //_____________________________________________________________________________
  36 TG4ExtDecayer::TG4ExtDecayer(const TG4ExtDecayer& right) {
  37 //
  38   TG4Globals::Exception(
  39     "TG4ExtDecayer is protected from copying.");
  40 }
  41
  42 //_____________________________________________________________________________
  43 TG4ExtDecayer::~TG4ExtDecayer() {
  44 //
  45   delete fDecayProductsArray;
  46 }
  47
  48
  49 // operators
  50
  51 //_____________________________________________________________________________
  52 TG4ExtDecayer& TG4ExtDecayer::operator=(const TG4ExtDecayer& right)
  53 {
  54   // check assignement to self
  55   if (this == &right) return *this;
  56
  57   TG4Globals::Exception(
  58     "TG4ExtDecayer is protected from assigning.");
  59
  60   return *this;
  61 }
  62
  63 // public methods
  64
  65 //_____________________________________________________________________________
  66 G4DecayProducts* TG4ExtDecayer::ImportDecayProducts(const G4Track& track)
  67 {
  68   // check if external decayer is defined
  69   if (!fExternalDecayer) {
  70 #ifdef G4VERBOSE
  71      G4cerr << "TG4ExtDecayer::ImportDecayProducts: " << G4endl
  72             << " No fExternalDecayer is defined." << G4endl;
  73 #endif
  74     return 0;
  75   }
  76
  77   // get particle momentum
  78   G4ThreeVector momentum = track.GetMomentum();
  79   G4double mag = momentum.mag();
  80   TLorentzVector p;
  81   p[0] = momentum.x()/mag;
  82   p[1] = momentum.x()/mag;
  83   p[2] = momentum.x()/mag;
  84   p[3] = mag/TG4G3Units::Energy();
  85
  86   // get particle PDG
  87   // ask TG4ParticlesManager to get PDG encoding
  88   // (in order to get PDG from extended TDatabasePDG
  89   // in case the standard PDG code is not defined)
  90   G4ParticleDefinition* particleDef = track.GetDefinition();
  91   G4int pdgEncoding = fParticlesManager->GetPDGEncodingFast(particleDef);
  92
  93   // let AliDecayer decay the particle
  94   // and import the decay products
  95   fExternalDecayer->Decay(pdgEncoding, &p);
  96   G4int nofParticles
  97     = fExternalDecayer->ImportParticles(fDecayProductsArray);
  98
  99   if (fVerboseLevel>0)
 100     G4cout << "nofParticles: " <<  nofParticles << G4endl;
 101
 102   // convert decay products TParticle type
 103   // to G4DecayProducts
 104   G4DecayProducts* decayProducts
 105     = new G4DecayProducts(*(track.GetDynamicParticle()));
 106
 107   G4int counter = 0;
 108   for (G4int i=0; i<nofParticles; i++) {
 109
 110     // get particle from TClonesArray
 111     TParticle* particle
 112       = fParticlesManager->GetParticle(fDecayProductsArray, i);
 113
 114     G4int status = particle->GetStatusCode();
 115     G4int pdg = particle->GetPdgCode();
 116     if ( status>0 && status<11 &&
 117          abs(pdg)!=12 && abs(pdg)!=14 && abs(pdg)!=16 ) {
 118       // pass to tracking final particles only;
 119       // skip neutrinos
 120
 121       if (fVerboseLevel>0)
 122         G4cout << "  " << i << "th particle PDG: " << pdg << "   ";
 123
 124       // create G4DynamicParticle
 125       G4DynamicParticle* dynamicParticle
 126         = fParticlesManager->CreateDynamicParticle(particle);
 127
 128       if (dynamicParticle) {
 129
 130         if (fVerboseLevel>0)
 131           G4cout << "  G4 particle name: "
 132                  << dynamicParticle->GetDefinition()->GetParticleName()
 133                  << G4endl;
 134
 135         // add dynamicParticle to decayProducts
 136         decayProducts->PushProducts(dynamicParticle);
 137
 138         counter++;
 139       }
 140     }
 141   }
 142   if (fVerboseLevel>0)
 143     G4cout << "nofParticles for tracking: " <<  counter << G4endl;
 144
 145   return decayProducts;
 146 }
 147