TGeant4/TG4ExtDecayer.cxx

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