18 #include "G4Electron.hh"
19 #include "G4Positron.hh"
20 #include "G4PionPlus.hh"
21 #include "G4PionMinus.hh"
27 #include "G4SDManager.hh"
36 const Gaudi::Property<std::string>& randomEngineName,
41 const Gaudi::Property<std::string>& CaloCellContainerSDName,
42 const Gaudi::Property<bool>& doG4Transport,
45 : G4VFastSimulationModel(
name),
46 m_rndmGenSvc(rndmGenSvc), m_randomEngineName(randomEngineName),
50 m_FastCaloSimSvc(FastCaloSimSvc),
51 m_CaloCellContainerSDName(CaloCellContainerSDName),
52 m_doG4Transport(doG4Transport),
78 bool isPositron = &
particleType == G4Positron::PositronDefinition();
85 const std::string pName =
particleType.GetParticleName();
86 G4cout<<
"[FastCaloSim::IsApplicable] Got " << pName <<G4endl;
87 if(
IsApplicable) G4cout<<
"[FastCaloSim::IsApplicable] APPLICABLE"<<G4endl;
88 else G4cout<<
"[FastCaloSim::IsApplicable] NOT APPLICABLE"<<G4endl;
99 G4cout<<
"[FastCaloSim::ModelTrigger] Got particle with " <<
"\n"
100 <<
" pdg=" <<fastTrack.GetPrimaryTrack() -> GetDefinition()->GetPDGEncoding() <<
"\n"
101 <<
" Ekin="<<fastTrack.GetPrimaryTrack() -> GetKineticEnergy() <<
"\n"
102 <<
" p=" <<fastTrack.GetPrimaryTrack() -> GetMomentum().mag() <<
"\n"
103 <<
" x=" <<fastTrack.GetPrimaryTrack() -> GetPosition().x() <<
"\n"
104 <<
" y=" <<fastTrack.GetPrimaryTrack() -> GetPosition().y() <<
"\n"
105 <<
" z=" <<fastTrack.GetPrimaryTrack() -> GetPosition().z() <<
"\n"
106 <<
" r=" <<fastTrack.GetPrimaryTrack() -> GetPosition().perp() <<
"\n"
107 <<
" eta=" <<fastTrack.GetPrimaryTrack() -> GetMomentum().eta() <<
"\n"
108 <<
" phi=" <<fastTrack.GetPrimaryTrack() -> GetMomentum().phi() <<
"\n"
114 if (fastTrack.GetPrimaryTrack() -> GetKineticEnergy() < 0.05) {
116 G4cout<<
"[FastCaloSim::ModelTrigger] Particle below 50 keV threshold. Passing to G4. "<<G4endl;
124 G4cout<<
"[FastCaloSim::ModelTrigger] Particle failed passedIDCaloBoundary z="<<fastTrack.GetPrimaryTrack() -> GetPosition().z()<<
" r="<<fastTrack.GetPrimaryTrack() -> GetPosition().perp()<<G4endl;
130 float minEkinPions = 200;
131 float minEkinOtherHadrons = 400;
134 const G4ParticleDefinition * G4Particle = fastTrack.GetPrimaryTrack() -> GetDefinition();
136 const float Ekin = fastTrack.GetPrimaryTrack() -> GetKineticEnergy();
139 bool isPhoton = G4Particle == G4Gamma::Definition();
140 bool isElectron = G4Particle == G4Electron::Definition();
141 bool isPositron = G4Particle == G4Positron::Definition();
142 bool isPionPlus = G4Particle == G4PionPlus::Definition();
143 bool isPionMinus = G4Particle == G4PionMinus::Definition();
148 G4cout<<
"[FastCaloSim::ModelTrigger] Model triggered"<<G4endl;
154 if (isPionPlus || isPionMinus){
155 bool passMinEkinPions = Ekin > minEkinPions;
158 if(passMinEkinPions) G4cout<<
"[FastCaloSim::ModelTrigger] Model triggered"<<G4endl;
159 else G4cout<<
"[FastCaloSim::ModelTrigger] Pion with Ekin="<<Ekin<<
" below the minimum "<<minEkinPions<<
" MeV threshold. Model not triggered."<<G4endl;
162 return passMinEkinPions;
167 bool passMinEkinOtherHadrons = Ekin > minEkinOtherHadrons;
170 if(passMinEkinOtherHadrons) G4cout<<
"[FastCaloSim::ModelTrigger] Model triggered"<<G4endl;
171 else G4cout<<
"[FastCaloSim::ModelTrigger] Other hadron with Ekin="<<Ekin<<
" below the minimum "<<minEkinOtherHadrons<<
" MeV threshold. Model not triggered."<<G4endl;
174 return passMinEkinOtherHadrons;
185 const G4Track * G4PrimaryTrack = fastTrack.GetPrimaryTrack();
187 const G4ParticleDefinition * G4Particle = G4PrimaryTrack -> GetDefinition();
189 signed int pdgID = G4Particle -> GetPDGEncoding();
192 if(G4PrimaryTrack -> GetKineticEnergy() < 10){
194 G4cout<<
"[FastCaloSim::DoIt] Skipping particle with Ekin: " << G4PrimaryTrack -> GetKineticEnergy() <<
" MeV. Below the 10 MeV threshold"<<G4endl;
196 fastStep.KillPrimaryTrack();
200 if(G4Particle == G4Electron::Definition() || G4Particle == G4Positron::Definition() || G4Particle == G4Gamma::Definition())
207 truthState.
set_pdgid(G4PionPlus::Definition() -> GetPDGEncoding());
211 truthState.SetPtEtaPhiM(G4PrimaryTrack -> GetMomentum().
perp(),
212 G4PrimaryTrack -> GetMomentum().
eta(),
213 G4PrimaryTrack -> GetMomentum().
phi(),
214 G4Particle -> GetPDGMass());
217 truthState.
set_vertex(G4PrimaryTrack -> GetPosition().
x(),
218 G4PrimaryTrack -> GetPosition().
y(),
219 G4PrimaryTrack -> GetPosition().
z());
225 if(pdgID == -2212 || pdgID == -2112) truthState.
set_Ekin_off(2 * G4Particle -> GetPDGMass());
240 G4cout<<
"[FastCaloSim::DoIt] Killing particle as extrapolation failed"<<G4endl;
242 fastStep.KillPrimaryTrack();
246 if(
m_FastCaloSimSvc->simulate(simState, &truthState, &extrapolState).isFailure()){
247 G4Exception(
"FastCaloSimSvc",
"FailedSimulationCall", FatalException,
"FastCaloSimSvc: Simulation call failed.");
252 G4cout<<
"[FastCaloSim::DoIt] Energy returned: " << simState.
E() << G4endl;
253 G4cout<<
"[FastCaloSim::DoIt] Energy fraction for layer: " << G4endl;
254 for (
int s = 0;
s < 24;
s++) G4cout<<
"[FastCaloSim::DoIt] Sampling " <<
s <<
" energy " << simState.
E(
s) << G4endl;
266 fastStep.KillPrimaryTrack();
267 fastStep.SetPrimaryTrackPathLength(0.0);
273 G4SDManager *sdm = G4SDManager::GetSDMpointer();
277 G4Exception(
"FastCaloSimSvc",
"FailedFindSensitiveDetector", FatalException,
"FastCaloSimSvc: Failed getting CaloCellContainer SD.");
283 if (!caloCellContainerSD){
284 G4Exception(
"FastCaloSimSvc",
"FailedCastSensitiveDetector", FatalException,
"FastCaloSimSvc: Failed casting G4VSensitiveDetector.");
287 return caloCellContainerSD;
304 const float barrelR = 1148;
305 const float barrelZ = 3549.5;
307 const float innerBeamPipeR = 120;
308 const float innerBeamPipeZ = 4587;
310 const float outerBeamPipeR = 41;
311 const float outerBeamPipeZ = 6783;
314 const G4ThreeVector particlePosition = fastTrack.GetPrimaryTrack() -> GetPosition();
316 const float r = particlePosition.perp();
317 const float z = particlePosition.z();
319 const G4ThreeVector particleDirection = fastTrack.GetPrimaryTrack() -> GetMomentum();
321 const G4ThreeVector helperLine = particlePosition + particleDirection;
325 const float triggerTolerance = 50;
328 if (std::abs(
z) <= barrelZ + triggerTolerance) {
329 if (
r >= barrelR &&
r < barrelR + triggerTolerance)
return helperLine.perp() >= barrelR;
332 if (std::abs(
z) >= barrelZ && std::abs(
z) <= innerBeamPipeZ + triggerTolerance){
333 if (
r >= innerBeamPipeR &&
r < innerBeamPipeR + triggerTolerance)
return helperLine.perp() >= innerBeamPipeR;
335 if (std::abs(
z) >= innerBeamPipeZ && std::abs(
z) <= outerBeamPipeZ){
336 if (
r >= outerBeamPipeR &&
r < outerBeamPipeR + triggerTolerance)
return helperLine.perp() >= outerBeamPipeR;
339 if (
r >= outerBeamPipeR &&
r<= innerBeamPipeR){
340 if (std::abs(
z) >= innerBeamPipeZ && std::abs(
z) < innerBeamPipeZ + triggerTolerance)
return std::abs(helperLine.z()) >= innerBeamPipeZ;
343 if (
r >= innerBeamPipeR &&
r <= barrelR){
344 if (std::abs(
z) >= barrelZ && std::abs(
z) < barrelZ + triggerTolerance)
return std::abs(helperLine.z()) >= barrelZ;