Monte Carlo LYSO
Geant4 simulation for the LYSO calorimeter prototype
stepping.cc
Go to the documentation of this file.
1 
5  #include "stepping.hh"
6 
7 MySteppingAction::MySteppingAction(MyEventAction *eventAction) : fEventAction(eventAction)
8 {}
9 
10 
11 
12 void MySteppingAction::UserSteppingAction(const G4Step *step)
13 {
14  // Get the logical volume containing the step
15  G4LogicalVolume *volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
16 
17  // Get the detector construction
18  const MyDetectorConstruction *detectorConstruction = static_cast<const MyDetectorConstruction*> (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
19 
20  // Settings depending on run mode type
21  G4int modeType = static_cast<const MyPrimaryGenerator*>(G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction())->GetModeType();
22  switch(modeType)
23  {
24  case 20: // 176Lu decay
25  case 21:
26  SetTimeOfDecay(step);
27  break;
28  case 22: // with trigger
29  SetTimeOfDecay(step);
30  SteppingForLuDecayBulkSignature(step, volume, detectorConstruction);
31  break;
32  case 30: // Cosmic rays
33  case 31:
34  SteppingForCosmicRaysDetectors(step, volume, detectorConstruction);
35  break;
36  }
37 
38  // Get the scoring volume (i.e. the scintillator logic volume)
39  G4LogicalVolume *fScoringVolume = detectorConstruction->GetScoringVolume();
40 
41 
42  // Here I'm selecting that I want data ONLY in the crystal
43  if(volume != fScoringVolume)
44  return;
45 
46  // Do not register optical photons (PDGEncoding == -22), they DON'T respect energy conservation
47  if(step->GetTrack()->GetParticleDefinition()->GetPDGEncoding()==-22)
48  return;
49 
50  // Store time and position of arrival of primary gamma (or daughters)
51  if(step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
52  {
53  G4double newtimein = step->GetPreStepPoint()->GetGlobalTime();
54  G4ThreeVector newposin = step->GetPreStepPoint()->GetPosition();
55 
56  G4double newtimeinter = step->GetPostStepPoint()->GetGlobalTime();
57  G4ThreeVector newposinter = step->GetPostStepPoint()->GetPosition();
58 
59  fEventAction->SetArrivalandFirstInteraction(newtimein, newposin, newtimeinter, newposinter);
60  }
61 
62  // Sum the energy deposited in the step
63  G4double edep = step->GetTotalEnergyDeposit();
64  fEventAction->AddEdep(edep);
65 
66  // Check if it is the maximum deposition of energy per unit length.
67  // If yes it will be stored with its position.
68  G4double dx = step->GetStepLength();
69  if(dx == 0) return;
70 
71  G4ThreeVector maxedeppos = (step->GetPreStepPoint()->GetPosition() + step->GetPostStepPoint()->GetPosition())/2;
72  fEventAction->SetMaxEdep(edep/dx, maxedeppos);
73 }
74 
75 
76 
77 void MySteppingAction::SteppingForCosmicRaysDetectors(const G4Step *step, G4LogicalVolume *volume, const MyDetectorConstruction *detectorConstruction)
78 {
79  // Get the scoring volume (i.e. the cosmic rays detector logic volume)
80  G4LogicalVolume *fCosmicVolume = detectorConstruction->GetCosmicTriggerVolume();
81 
82  // Here I'm selecting that I want data ONLY for cosmic rays detectors
83  if(!fCosmicVolume || volume != fCosmicVolume)
84  return;
85 
86  // Should see this condition. Now I'm selecting only primary muon
87  if(step->GetTrack()->GetTrackID() != 1)
88  return;
89 
90  // Up (0) or Bottom (1) detector
91  G4bool isUpOrBottom = step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber();
92 
93  if(!isUpOrBottom) // up
94  fEventAction->SetCosmicTriggerUp(true);
95  else // bottom
96  {
97  fEventAction->SetCosmicTriggerBottom(true);
98  step->GetTrack()->SetTrackStatus(fStopAndKill);
99  }
100 }
101 
102 
103 
104 void MySteppingAction::SetTimeOfDecay(const G4Step *step)
105 {
106  // Set time = 0 ns also for primary data. Instantaneous decay
107  if(step->GetTrack()->GetParentID() == 0 && step->IsFirstStepInVolume())
108  {
109  G4double newtimein = step->GetPreStepPoint()->GetGlobalTime();
110  G4ThreeVector newposin = step->GetPreStepPoint()->GetPosition();
111 
112  fEventAction->SetArrivalandFirstInteraction(newtimein, newposin, newtimein, newposin);
113  }
114 
115  // Find if the primary is killed by radioactive decay
116  if(step->GetTrack()->GetParentID() == 0 && step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessType() == 6)
117  {
118  const auto* fDecayProducts = step->GetSecondaryInCurrentStep();
119 
120  for(auto i = 0; i < fDecayProducts->size(); i++)
121  {
122  G4Track* modifiableTrack = const_cast<G4Track*>((*fDecayProducts)[i]);
123  modifiableTrack->SetGlobalTime(0.);
124  }
125  }
126 }
127 
128 
129 
130 void MySteppingAction::SteppingForLuDecayBulkSignature(const G4Step *step, G4LogicalVolume *volume, const MyDetectorConstruction *detectorConstruction)
131 {
132  // Get the scoring volume (i.e. the cosmic rays detector logic volume)
133  G4LogicalVolume *fSiVolume = detectorConstruction->GetDecayTriggerVolume();
134 
135  // Here I'm selecting that I want data ONLY for Si detectors
136  if(volume != fSiVolume)
137  return;
138 
139  // Save only electrons
140  if(step->GetTrack()->GetParticleDefinition()->GetPDGEncoding() != 11)
141  return;
142 
143  // Save only detector Front-57, default for 176Lu decay Si-Trigger runs
144  if(step->GetPreStepPoint()->GetTouchableHandle()->GetVolume(2)->GetTranslation().z() > GS::zScintillator || step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(2) != 57)
145  return;
146 
147  fEventAction->SetDecayTriggerSi(true);
148 }
Mandatory user initialization concrete class of G4VUserDetectorConstruction. It represents the whole ...
Definition: construction.hh:33
G4LogicalVolume * GetScoringVolume() const
Get the scoring volume. It will be used by MySteppingAction::UserSteppingAction().
Definition: construction.hh:46
User action concrete class of G4UserEventAction. In addition to defining procedures executed at the s...
Definition: event.hh:26
void SetArrivalandFirstInteraction(G4double newtimein, G4ThreeVector newposin, G4double newtimeinter, G4ThreeVector newposinter)
Stores the time and the position of arrival to the crystal of the primary gamma.
Definition: event.hh:55
void AddEdep(G4double edep)
For every G4Step inside the crystal it sums the energy deposit.
Definition: event.hh:75
void SetMaxEdep(G4double edepondx, G4ThreeVector maxedeppos)
Stores the maximum deposit of energy per unit length and its position inside the crystal.
Definition: event.hh:84
Mandatory user action concrete class of G4VUserPrimaryGeneratorAction. It defines the settings for th...
Definition: generator.hh:25
void UserSteppingAction(const G4Step *step) override
For every step updates and stores quantities about the physics inside the crystal by calling the meth...
Definition: stepping.cc:12
MyEventAction * fEventAction
Pointer to the MyEventAction object.
Definition: stepping.hh:43
MySteppingAction(MyEventAction *eventAction)
Constructor of the class.
Definition: stepping.cc:7
constexpr G4double zScintillator
z-position of the scintillator crystal.
Declaration of the class MySteppingAction.