-
Notifications
You must be signed in to change notification settings - Fork 4.2k
/
GFlashEMShowerModel.cc
193 lines (154 loc) · 7.14 KB
/
GFlashEMShowerModel.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
//
// initial setup : E.Barberio & Joanna Weng
// big changes : Soon Jun & Dongwook Jang
// V.Ivanchenko rename the class, cleanup, and move
// to SimG4Core/Application - 2012/08/14
//
#include "SimG4Core/Application/interface/SteppingAction.h"
#include "SimG4Core/Application/interface/GFlashEMShowerModel.h"
#include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
#include "SimGeneral/GFlash/interface/GflashHit.h"
#include "G4Electron.hh"
#include "G4Positron.hh"
#include "G4VProcess.hh"
#include "G4VPhysicalVolume.hh"
#include "G4LogicalVolume.hh"
#include "G4TransportationManager.hh"
#include "G4EventManager.hh"
#include "G4FastSimulationManager.hh"
#include "G4TouchableHandle.hh"
#include "G4VSensitiveDetector.hh"
#include "G4SystemOfUnits.hh"
GFlashEMShowerModel::GFlashEMShowerModel(const G4String& modelName,
G4Envelope* envelope,
const edm::ParameterSet& parSet)
: G4VFastSimulationModel(modelName, envelope), theParSet(parSet)
{
theWatcherOn = parSet.getParameter<bool>("watcherOn");
theProfile = new GflashEMShowerProfile(parSet);
theRegion = const_cast<const G4Region*>(envelope);
theGflashStep = new G4Step();
theGflashTouchableHandle = new G4TouchableHistory();
theGflashNavigator = new G4Navigator();
}
// --------------------------------------------------------------------------
GFlashEMShowerModel::~GFlashEMShowerModel()
{
delete theProfile;
delete theGflashStep;
}
G4bool
GFlashEMShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
{
return ( &particleType == G4Electron::Electron() ||
&particleType == G4Positron::Positron() );
}
// --------------------------------------------------------------------------
G4bool GFlashEMShowerModel::ModelTrigger(const G4FastTrack & fastTrack )
{
// Mininum energy cutoff to parameterize
if(fastTrack.GetPrimaryTrack()->GetKineticEnergy() < Gflash::energyCutOff)
{ return false; }
if(excludeDetectorRegion(fastTrack)) { return false; }
// This will be changed accordingly when the way
// dealing with CaloRegion changes later.
G4TouchableHistory* touch =
(G4TouchableHistory*)(fastTrack.GetPrimaryTrack()->GetTouchable());
G4VPhysicalVolume* pCurrentVolume = touch->GetVolume();
if( pCurrentVolume == nullptr) { return false; }
G4LogicalVolume* lv = pCurrentVolume->GetLogicalVolume();
if(lv->GetRegion() != theRegion) { return false; }
//std::cout << "GFlashEMShowerModel::ModelTrigger: LV "
// << lv->GetRegion()->GetName() << std::endl;
// The parameterization starts inside crystals
//std::size_t pos1 = lv->GetName().find("EBRY");
//std::size_t pos2 = lv->GetName().find("EFRY");
//std::size_t pos3 = lv->GetName().find("HVQ");
//std::size_t pos4 = lv->GetName().find("HF");
//if(pos1 == std::string::npos && pos2 == std::string::npos &&
// pos3 == std::string::npos && pos4 == std::string::npos) return false;
return true;
}
// ---------------------------------------------------------------------------
void GFlashEMShowerModel::DoIt(const G4FastTrack& fastTrack,
G4FastStep& fastStep)
{
// Kill the parameterised particle:
fastStep.KillPrimaryTrack();
fastStep.ProposePrimaryTrackPathLength(0.0);
// Input variables for GFlashEMShowerProfile with showerType = 1,5
// Shower starts inside crystals
G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV;
G4double globalTime =
fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
G4double charge =
fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition()/cm;
G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum()/GeV;
G4int showerType = Gflash::findShowerType(position);
// Do actual parameterization
// The result of parameterization is gflashHitList
theProfile->initialize(showerType,energy,globalTime,charge,
position,momentum);
theProfile->parameterization();
// Make hits
makeHits(fastTrack);
}
// ---------------------------------------------------------------------------
void GFlashEMShowerModel::makeHits(const G4FastTrack& fastTrack)
{
std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
theGflashStep->SetTrack(const_cast<G4Track*>(fastTrack.GetPrimaryTrack()));
theGflashStep->GetPostStepPoint()
->SetProcessDefinedStep(const_cast<G4VProcess*>(fastTrack.GetPrimaryTrack()
->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
for(; spotIter != spotIterEnd; spotIter++){
// Put touchable for each hit so that touchable history
// keeps track of each step.
theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(spotIter->getPosition(),
G4ThreeVector(0,0,0),
theGflashTouchableHandle, false);
updateGflashStep(spotIter->getPosition(),spotIter->getTime());
// If there is a watcher defined in a job and the flag is turned on
if(theWatcherOn) {
SteppingAction* userSteppingAction =
(SteppingAction*) G4EventManager::GetEventManager()->GetUserSteppingAction();
userSteppingAction->m_g4StepSignal(theGflashStep);
}
// Send G4Step information to Hit/Digi if the volume is sensitive
// Copied from G4SteppingManager.cc
G4VPhysicalVolume* aCurrentVolume =
theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
if( aCurrentVolume == nullptr ) { continue; }
G4LogicalVolume* lv = aCurrentVolume->GetLogicalVolume();
if(lv->GetRegion() != theRegion) { continue; }
theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
G4VSensitiveDetector* aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
if( aSensitive == nullptr ) { continue; }
theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy());
aSensitive->Hit(theGflashStep);
}
}
// ---------------------------------------------------------------------------
void GFlashEMShowerModel::updateGflashStep(const G4ThreeVector& spotPosition,
G4double timeGlobal)
{
theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
}
// ---------------------------------------------------------------------------
G4bool GFlashEMShowerModel::excludeDetectorRegion(const G4FastTrack& fastTrack)
{
G4bool isExcluded=false;
//exclude regions where geometry are complicated
//+- one supermodule around the EB/EE boundary: 1.479 +- 0.0174*5
G4double eta = fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity();
if(std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566) { return true; }
return isExcluded;
}
// ---------------------------------------------------------------------------