6 #include <phparameter/PHParameters.h>
19 #include <calobase/TowerInfo.h>
20 #include <calobase/TowerInfoContainer.h>
21 #include <calobase/TowerInfoContainerv1.h>
22 #include <calobase/TowerInfoDefs.h>
38 #include <Geant4/G4ParticleDefinition.hh>
39 #include <Geant4/G4ReferenceCountedHandle.hh>
40 #include <Geant4/G4Step.hh>
41 #include <Geant4/G4StepPoint.hh>
42 #include <Geant4/G4StepStatus.hh>
43 #include <Geant4/G4String.hh>
44 #include <Geant4/G4SystemOfUnits.hh>
45 #include <Geant4/G4ThreeVector.hh>
46 #include <Geant4/G4TouchableHandle.hh>
47 #include <Geant4/G4Track.hh>
48 #include <Geant4/G4TrackStatus.hh>
49 #include <Geant4/G4TransportationManager.hh>
50 #include <Geant4/G4Types.hh>
51 #include <Geant4/G4VPhysicalVolume.hh>
52 #include <Geant4/G4VTouchable.hh>
53 #include <Geant4/G4VUserTrackInformation.hh>
67 , m_Detector(detector)
68 , m_Params(parameters)
69 , m_IsActive(m_Params->get_int_param(
"active"))
70 , m_IsBlackHole(m_Params->get_int_param(
"blackhole"))
71 , m_LightScintModelFlag(m_Params->get_int_param(
"light_scint_model"))
72 , m_doG4Hit(m_Params->get_int_param(
"saveg4hit"))
73 , m_tmin(m_Params->get_double_param(
"tmin"))
74 , m_tmax(m_Params->get_double_param(
"tmax"))
75 , m_dt(m_Params->get_double_param(
"dt"))
100 if (ihcalmapname.empty())
107 std::cout <<
PHWHERE <<
" Could not locate " << url << std::endl;
108 std::cout <<
"use empty filename to ignore mapfile" << std::endl;
111 TFile*
file = TFile::Open(url.c_str());
129 catch (std::exception&
e)
131 std::cout << e.what() << std::endl;
143 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
144 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
165 layer_id = layer_tower.first;
166 tower_id = layer_tower.second;
178 G4StepPoint* prePoint = aStep->GetPreStepPoint();
179 G4StepPoint* postPoint = aStep->GetPostStepPoint();
181 double pretime = prePoint->GetGlobalTime() / nanosecond;
182 double posttime = postPoint->GetGlobalTime() / nanosecond;
183 if (posttime < m_tmin || pretime >
m_tmax)
return false;
184 if ((posttime - pretime) >
m_dt)
return false;
185 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
186 const G4Track* aTrack = aStep->GetTrack();
188 double light_yield = eion;
196 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
197 const G4ThreeVector& worldPosition = postPoint->GetPosition();
198 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
199 float lx = (localPosition.x() /
cm);
200 float lz = fabs(localPosition.z() /
cm);
202 int lcz = (int) (5.0 * lz) + 1;
203 int lcx = (int) (5.0 * (lx + 12.1)) + 1;
218 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
222 unsigned int ieta = tower_id;
223 unsigned int iphi = (
unsigned int) layer_id / 4;
229 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
248 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
249 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
270 layer_id = layer_tower.first;
271 tower_id = layer_tower.second;
277 layer_id = touch->GetCopyNumber();
280 G4double edep = aStep->GetTotalEnergyDeposit() /
GeV;
281 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
282 G4double light_yield = 0;
283 const G4Track* aTrack = aStep->GetTrack();
288 edep = aTrack->GetKineticEnergy() /
GeV;
289 G4Track* killtrack =
const_cast<G4Track*
>(aTrack);
290 killtrack->SetTrackStatus(fStopAndKill);
296 bool geantino =
false;
301 if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
302 aTrack->GetParticleDefinition()->GetParticleName().find(
"geantino") != std::string::npos)
306 G4StepPoint* prePoint = aStep->GetPreStepPoint();
307 G4StepPoint* postPoint = aStep->GetPostStepPoint();
311 switch (prePoint->GetStepStatus())
313 case fPostStepDoItProc:
320 std::cout <<
GetName() <<
": New Hit for " << std::endl;
326 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
327 std::cout <<
"phys pre vol: " << volume->GetName()
328 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
329 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
330 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
346 m_Hit->
set_t(0, prePoint->GetGlobalTime() / nanosecond);
365 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
382 std::cout <<
GetName() <<
": hit was not created" << std::endl;
388 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
389 std::cout <<
"phys pre vol: " << volume->GetName()
390 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
391 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
392 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
398 std::cout <<
GetName() <<
": hits do not belong to the same track" << std::endl;
400 <<
", current trackid: " << aTrack->GetTrackID()
401 <<
", prestep status: " << prePoint->GetStepStatus()
418 m_Hit->
set_t(1, postPoint->GetGlobalTime() / nanosecond);
432 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
433 const G4ThreeVector& worldPosition = postPoint->GetPosition();
434 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
435 float lx = (localPosition.x() /
cm);
436 float lz = fabs(localPosition.z() /
cm);
438 int lcz = (int) (5.0 * lz) + 1;
439 int lcx = (int) (5.0 * (lx + 12.1)) + 1;
454 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
466 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
483 if (postPoint->GetStepStatus() == fGeomBoundary ||
484 postPoint->GetStepStatus() == fWorldBoundary ||
485 postPoint->GetStepStatus() == fAtRestDoItProc ||
486 aTrack->GetTrackStatus() == fStopAndKill)
535 m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
536 m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename);
541 std::cout <<
"PHG4InnerHcalSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
547 std::cout <<
"PHG4HcalSteppingAction::SetTopNode - unable to find " << absorbernodename << std::endl;
556 nodeItr.
findFirst(
"PHCompositeNode",
"DST"));
559 std::cout <<
"PHComposite node created: DST" << std::endl;