7 #include <phparameter/PHParameters.h>
25 #include <calobase/TowerInfo.h>
26 #include <calobase/TowerInfoContainer.h>
27 #include <calobase/TowerInfoContainerv1.h>
28 #include <calobase/TowerInfoDefs.h>
36 #include <Geant4/G4ParticleDefinition.hh>
37 #include <Geant4/G4ReferenceCountedHandle.hh>
38 #include <Geant4/G4Step.hh>
39 #include <Geant4/G4StepPoint.hh>
40 #include <Geant4/G4StepStatus.hh>
41 #include <Geant4/G4String.hh>
42 #include <Geant4/G4SystemOfUnits.hh>
43 #include <Geant4/G4ThreeVector.hh>
44 #include <Geant4/G4TouchableHandle.hh>
45 #include <Geant4/G4Track.hh>
46 #include <Geant4/G4TrackStatus.hh>
47 #include <Geant4/G4TransportationManager.hh>
48 #include <Geant4/G4Types.hh>
49 #include <Geant4/G4VPhysicalVolume.hh>
50 #include <Geant4/G4VTouchable.hh>
51 #include <Geant4/G4VUserTrackInformation.hh>
65 , m_Detector(detector)
66 , m_Params(parameters)
67 , m_IsActive(m_Params->get_int_param(
"active"))
68 , m_IsBlackHole(m_Params->get_int_param(
"blackhole"))
69 , m_LightScintModelFlag(m_Params->get_int_param(
"light_scint_model"))
70 , m_doG4Hit(m_Params->get_int_param(
"saveg4hit"))
71 , m_tmin(m_Params->get_double_param(
"tmin"))
72 , m_tmax(m_Params->get_double_param(
"tmax"))
73 , m_dt(m_Params->get_double_param(
"dt"))
98 if (ihcalmapname.empty())
105 std::cout <<
"use empty filename to ignore mapfile" << std::endl;
108 TFile*
file = TFile::Open(ihcalmapname.c_str());
125 catch (std::exception&
e)
127 std::cout << e.what() << std::endl;
139 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
140 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
163 layer_id = std::get<1>(layer_tower);
164 tower_id = std::get<2>(layer_tower);
176 G4StepPoint* prePoint = aStep->GetPreStepPoint();
177 G4StepPoint* postPoint = aStep->GetPostStepPoint();
179 double pretime = prePoint->GetGlobalTime() / nanosecond;
180 double posttime = postPoint->GetGlobalTime() / nanosecond;
181 if (posttime < m_tmin || pretime >
m_tmax)
return false;
182 if ((posttime - pretime) >
m_dt)
return false;
183 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
184 const G4Track* aTrack = aStep->GetTrack();
186 double light_yield = eion;
194 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
195 const G4ThreeVector& worldPosition = postPoint->GetPosition();
196 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
197 float lx = localPosition.x() /
cm;
198 float ly = localPosition.y() /
cm;
201 int lcx = (int) (5.0 * lx) + 1;
202 int lcy = (int) (5.0 * (ly + 2.0)) + 1;
216 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
220 unsigned int ieta = tower_id;
221 unsigned int iphi = (
unsigned int) layer_id / 4;
227 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
246 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
247 G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
269 sector_id = std::get<0>(layer_tower);
270 layer_id = std::get<1>(layer_tower);
271 tower_id = std::get<2>(layer_tower);
280 G4double edep = aStep->GetTotalEnergyDeposit() /
GeV;
281 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
282 double 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:
324 std::cout <<
GetName() <<
": Bad step status combination for the same track " << std::endl;
330 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
331 std::cout <<
"phys pre vol: " << volume->GetName()
332 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
333 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
334 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
341 std::cout <<
GetName() <<
": New Hit for " << std::endl;
347 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
348 std::cout <<
"phys pre vol: " << volume->GetName()
349 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
350 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
351 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
367 m_Hit->
set_t(0, prePoint->GetGlobalTime() / nanosecond);
387 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
404 std::cout <<
GetName() <<
": hit was not created" << std::endl;
410 <<
", current trackid: " << aTrack->GetTrackID() << std::endl;
411 std::cout <<
"phys pre vol: " << volume->GetName()
412 <<
" post vol : " << touchpost->GetVolume()->GetName() << std::endl;
413 std::cout <<
" previous phys pre vol: " <<
m_SaveVolPre->GetName()
414 <<
" previous phys post vol: " <<
m_SaveVolPost->GetName() << std::endl;
420 std::cout <<
GetName() <<
": hits do not belong to the same track" << std::endl;
422 <<
", current trackid: " << aTrack->GetTrackID()
423 <<
", prestep status: " << prePoint->GetStepStatus()
440 m_Hit->
set_t(1, postPoint->GetGlobalTime() / nanosecond);
454 const G4TouchableHandle& theTouchable = prePoint->GetTouchableHandle();
455 const G4ThreeVector& worldPosition = postPoint->GetPosition();
456 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
457 float lx = localPosition.x() /
cm;
458 float ly = localPosition.y() /
cm;
461 int lcx = (int) (5.0 * lx) + 1;
462 int lcy = (int) (5.0 * (ly + 2.0)) + 1;
476 light_yield = light_yield *
GetLightCorrection(postPoint->GetPosition().x(), postPoint->GetPosition().y());
488 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
505 if (postPoint->GetStepStatus() == fGeomBoundary ||
506 postPoint->GetStepStatus() == fWorldBoundary ||
507 postPoint->GetStepStatus() == fAtRestDoItProc ||
508 aTrack->GetTrackStatus() == fStopAndKill)
549 std::cout <<
"PHG4IHCalSteppingAction::SetTopNode - unable to find " <<
m_HitNodeName << std::endl;
551 if (!m_AbsorberHitContainer)
555 std::cout <<
"PHG4IHcalSteppingAction::SetTopNode - unable to find " <<
m_AbsorberNodeName << std::endl;
567 else if (type ==
"G4HIT_ABSORBER")
572 std::cout <<
"Invalid output hit node type " << type << std::endl;
581 nodeItr.
findFirst(
"PHCompositeNode",
"DST"));
584 std::cout <<
"PHComposite node created: DST" << std::endl;