FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FairFastSimModel.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 // -------------------------------------------------------------------------
9 // ----- FairFastSimModel source file -----
10 // ----- Created 2019/01/30 by R. Karabowicz -----
11 // -------------------------------------------------------------------------
12 
13 #include "FairFastSimModel.h"
14 
15 #include "FairGenericStack.h"
16 #include "FairLogger.h"
17 #include "FairMCApplication.h"
18 #include "FairTrajFilter.h"
19 #include "G4Electron.hh"
20 #include "G4Gamma.hh"
21 #include "G4ParticleTable.hh"
22 #include "G4StackManager.hh"
23 #include "G4SystemOfUnits.hh"
24 #include "G4VPhysicalVolume.hh"
25 
26 #include <TClonesArray.h> // for TClonesArray
27 #include <TGeoManager.h>
28 #include <TParticle.h>
29 #include <TVector3.h> // for TVector3
30 #include <TVirtualMC.h>
31 #include <tuple> // for tie, tuple
32 
33 // I.H. make this optional
34 // #include "G4GDMLParser.hh"
35 
36 FairFastSimModel::FairFastSimModel(G4String modelName, G4Region* envelope)
37  : G4VFastSimulationModel(modelName, envelope)
38 {}
39 
41  : G4VFastSimulationModel(modelName)
42 {}
43 
45 
46 // I.H. make this optional
47 // void FairFastSimModel::WriteGeometryToGDML(
48 // G4VPhysicalVolume* physicalVolume) {
49 
50 // G4GDMLParser* parser = new G4GDMLParser();
51 // remove("garfieldGeometry.gdml");
52 // parser->Write("garfieldGeometry.gdml", physicalVolume, false);
53 // delete parser;
54 // }
55 
56 G4bool FairFastSimModel::IsApplicable([[gnu::unused]] const G4ParticleDefinition& particleType) { return true; }
57 
58 G4bool FairFastSimModel::ModelTrigger([[gnu::unused]] const G4FastTrack& fastTrack) { return true; }
59 
60 void FairFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
61 {
62 
63  LOG(debug) << "FairFastSimModel::DoIt() called";
64  // G4TouchableHandle theTouchable =
65  // fastTrack.GetPrimaryTrack()->GetTouchableHandle();
66  // G4String name = theTouchable->GetVolume()->GetName();
67 
68  int firstSecondary(0);
69  int nofSecondaries(0);
70  int movedParticleIndex(0);
71  std::tie(movedParticleIndex, firstSecondary, nofSecondaries) =
72  ((FairGenericStack*)(gMC->GetStack()))->FastSimGetMovedIndex();
73  if (movedParticleIndex == -2) {
75  std::tie(movedParticleIndex, firstSecondary, nofSecondaries) =
76  ((FairGenericStack*)(gMC->GetStack()))->FastSimGetMovedIndex();
77  }
78  TClonesArray* particles = ((FairGenericStack*)(gMC->GetStack()))->GetListOfParticles();
79  if (nofSecondaries != 0) {
80  fastStep.SetNumberOfSecondaryTracks(nofSecondaries);
81 
82  int addedParticleIndex = firstSecondary;
83  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
84  for (int ipart = 0; ipart < nofSecondaries; ipart++) {
85  if (addedParticleIndex == movedParticleIndex)
86  addedParticleIndex++; // added protection in case user adds, moves and adds again
87  TParticle* particle = (TParticle*)particles->At(addedParticleIndex);
88 
89  G4ParticleDefinition* particleDefinition = 0;
90  particleDefinition = particleTable->FindParticle(particle->GetPdgCode());
91  if (!particleDefinition)
92  LOG(fatal) << "FairFastSimModel::DoIt() PDG " << particle->GetPdgCode() << " unknown!";
93  G4ThreeVector pos(
94  particle->Vx() * 10., particle->Vy() * 10., particle->Vz() * 10.); // change from cm to mm
95  G4ThreeVector mom(particle->Px() * 1000., particle->Py() * 1000., particle->Pz() * 1000.);
96  TVector3 polVect;
97  particle->GetPolarisation(polVect);
98  G4ThreeVector pol(polVect.X(), polVect.Y(), polVect.Z()); // get polarisation
99  G4double tim = particle->T() * 10.e9; // change from ns to s
100  G4double ek = particle->Ek() * 1000.;
101  G4DynamicParticle dynParticle(particleDefinition, mom.unit(), ek);
102  G4Track* tempTrack = fastStep.CreateSecondaryTrack(dynParticle, pos, tim, false);
103  tempTrack->SetTouchableHandle(fastTrack.GetPrimaryTrack()->GetTouchableHandle());
104  addedParticleIndex++;
105  }
106  }
107 
108  LOG(debug) << "FairFastSimModel::DoIt() moving particle " << gMC->GetStack()->GetCurrentTrackNumber() << " (pos #"
109  << movedParticleIndex << ").";
110  if (movedParticleIndex != -1) {
111  TParticle* particle = (TParticle*)particles->At(movedParticleIndex);
112 
113  LOG(debug) << "FAST SIM (moving particle to <" << particle->Vx() << "," << particle->Vy() << ","
114  << particle->Vz() << "," << particle->T() << "> with p=<" << particle->Px() << "," << particle->Py()
115  << "," << particle->Pz() << "," << particle->Ek() << ">)";
116 
117  G4ThreeVector pos(particle->Vx() * 10., particle->Vy() * 10., particle->Vz() * 10.); // change from cm to mm
118  G4ThreeVector mom(particle->Px() * 1000., particle->Py() * 1000., particle->Pz() * 1000.);
119  G4ThreeVector pol(0., 0., 0.); // get polarisation
120  G4double tim = particle->T() * 10.e9; // change from ns to s
121  G4double len = 0.;
122  G4double ek = particle->Ek() * 1000.;
123 
124  fastStep.SetPrimaryTrackFinalPosition(pos, false);
125  fastStep.SetPrimaryTrackFinalTime(tim);
126  fastStep.SetPrimaryTrackFinalProperTime(tim);
127  fastStep.SetPrimaryTrackFinalMomentum(mom, false);
128  fastStep.SetPrimaryTrackFinalKineticEnergy(ek);
129  // fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection (G4double, const G4ThreeVector &, G4bool
130  // localCoordinates=true);
131  fastStep.SetPrimaryTrackFinalPolarization(pol, false);
132  fastStep.SetPrimaryTrackPathLength(len);
133 
134  if (FairTrajFilter::Instance()->IsAccepted(particle)) {
136  particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
137  }
138 
139  gGeoManager->SetCurrentPoint(particle->Vx(), particle->Vy(), particle->Vz());
140  }
141  ((FairGenericStack*)(gMC->GetStack()))->FastSimClearMovedIndex();
142 }
virtual G4bool IsApplicable([[gnu::unused]] const G4ParticleDefinition &)
virtual void Stepping()
virtual G4bool ModelTrigger([[gnu::unused]] const G4FastTrack &)
FairFastSimModel(G4String, G4Region *)
virtual void DoIt(const G4FastTrack &, G4FastStep &)
TGeoTrack * GetCurrentTrk()
static FairMCApplication * Instance()
static FairTrajFilter * Instance()