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